Ignore:
Timestamp:
Sep 12, 2007, 4:13:11 PM (17 years ago)
Author:
ole
Message:

Further optimisation when beta_h == 0 (gained about 5s bring the
total time of 30s down to 25s in the okushiri performance test).

Also cleared out some obsolete code

File:
1 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r4731 r4733  
    1818#include <stdio.h>
    1919
    20 //Shared code snippets
     20// Shared code snippets
    2121#include "util_ext.h"
    2222
    … …  
    316316
    317317// Computational function for flux computation (using stage w=z+h)
    318 // FIXME (Ole): Is this used anywhere??
    319318int flux_function_kinetic(double *q_left, double *q_right,
    320319                  double z_left, double z_right,
    … …  
    471470
    472471int _balance_deep_and_shallow(int N,
     472                              double beta_h,
    473473                              double* wc,
    474474                              double* zc,
    475                               double* hc,
    476475                              double* wv,
    477476                              double* zv,
    478                               double* hv,
    479                               double* hvbar,
     477                              double* hvbar, // Retire this
    480478                              double* xmomc,
    481479                              double* ymomc,
    … …  
    487485
    488486  int k, k3, i;
    489   double dz, hmin, alpha, h_diff;
     487  double dz, hmin, alpha, h_diff, hc_k;
     488  double epsilon = 1.0e-6; // Temporary measure
     489  double hv[3]; // Depths at vertices
    490490
    491491  // Compute linear combination between w-limited stages and
    … …  
    501501
    502502    k3 = 3*k;
    503 
     503    hc_k = wc[k] - zc[k]; // Centroid value at triangle k
    504504   
    505505    dz = 0.0;
    … …  
    511511    }
    512512
     513    // Calculate depth at vertices
     514    hv[0] = wv[k3] -   zv[k3];
     515    hv[1] = wv[k3+1] - zv[k3+1];
     516    hv[2] = wv[k3+2] - zv[k3+2];       
     517   
    513518    // Calculate minimal depth across all three vertices
    514     hmin = min(hv[k3], min(hv[k3+1], hv[k3+2]));
    515    
     519    hmin = min(hv[0], min(hv[1], hv[2]));
    516520   
    517521
    … …  
    544548        for (i=0; i<3; i++) {
    545549
    546           h_diff = hvbar[k3+i] - hv[k3+i];
     550          // FIXME (Ole): Simplify when (if) hvbar gets retired
     551          if (beta_h > epsilon) {
     552            h_diff = hvbar[k3+i] - hv[i];
     553          } else {
     554            h_diff = hc_k - hv[i];       
     555          }
    547556       
    548557          if (h_diff <= 0) {
    … …  
    554563            // h-limited stage.
    555564           
    556             alpha = min(alpha, (hvbar[k3+i] - H0)/h_diff);
     565            // FIXME (Ole): Simplify when (if) hvbar gets retired           
     566            if (beta_h > epsilon) {       
     567              alpha = min(alpha, (hvbar[k3+i] - H0)/h_diff);
     568            } else {
     569              alpha = min(alpha, (hc_k - H0)/h_diff);       
     570            }
    557571          }
    558572        }
    … …  
    594608    if (alpha < 1) {
    595609      for (i=0; i<3; i++) {
    596         wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i];
     610     
     611        // FIXME (Ole): Simplify when (if) hvbar gets retired       
     612        if (beta_h > epsilon) {   
     613          wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[i];
     614        } else {
     615          wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];   
     616        }
    597617
    598618        // Update momentum as a linear combination of
    … …  
    647667        //New code: Adjust momentum to guarantee speeds are physical
    648668        //          ensure h is non negative
    649         //FIXME (Ole): This is only implemented in this C extension and
    650         //             has no Python equivalent
    651669             
    652670        if (hc <= 0.0) {
    … …  
    14511469
    14521470
    1453 
    1454 // FIXME (Ole): This function is obsolete as of 12 Sep 2007
    1455 PyObject *extrapolate_second_order_sw_original(PyObject *self, PyObject *args) {
    1456   /*Compute the vertex values based on a linear reconstruction on each triangle
    1457     These values are calculated as follows:
    1458     1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle
    1459     formed by the centroids of its three neighbours.
    1460     2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product
    1461     of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average
    1462     of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the
    1463     original triangle has gradient (a,b).
    1464     3) Provisional vertex jumps dqv[0,1,2] are computed and these are then limited by calling the functions
    1465     find_qmin_and_qmax and limit_gradient
    1466 
    1467     Python call:
    1468     extrapolate_second_order_sw(domain.surrogate_neighbours,
    1469                                 domain.number_of_boundaries
    1470                                 domain.centroid_coordinates,
    1471                                 Stage.centroid_values
    1472                                 Xmom.centroid_values
    1473                                 Ymom.centroid_values
    1474                                 domain.vertex_coordinates,
    1475                                 Stage.vertex_values,
    1476                                 Xmom.vertex_values,
    1477                                 Ymom.vertex_values)
    1478 
    1479     Post conditions:
    1480             The vertices of each triangle have values from a limited linear reconstruction
    1481             based on centroid values
    1482 
    1483   */
    1484   PyArrayObject *surrogate_neighbours,
    1485     *number_of_boundaries,
    1486     *centroid_coordinates,
    1487     *stage_centroid_values,
    1488     *xmom_centroid_values,
    1489     *ymom_centroid_values,
    1490         *elevation_centroid_values,
    1491     *vertex_coordinates,
    1492     *stage_vertex_values,
    1493     *xmom_vertex_values,
    1494     *ymom_vertex_values,
    1495         *elevation_vertex_values;
    1496   PyObject *domain, *Tmp;
    1497   //Local variables
    1498   double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
    1499   int number_of_elements,k,k0,k1,k2,k3,k6,coord_index,i;
    1500   double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2;//vertices of the auxiliary triangle
    1501   double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
    1502   double dqv[3], qmin, qmax, hmin, hmax;
    1503   double hc, h0, h1, h2;
    1504   double epsilon=1.0e-12;
    1505   int optimise_dry_cells=1; // Optimisation flag   
    1506   double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, beta_tmp;
    1507   double minimum_allowed_height;
    1508  
    1509   // Provisional jumps from centroids to v'tices and safety factor re limiting
    1510   // by which these jumps are limited
    1511   // Convert Python arguments to C
    1512   if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
    1513                         &domain,
    1514                         &surrogate_neighbours,
    1515                         &number_of_boundaries,
    1516                         &centroid_coordinates,
    1517                         &stage_centroid_values,
    1518                         &xmom_centroid_values,
    1519                         &ymom_centroid_values,
    1520                         &elevation_centroid_values,
    1521                         &vertex_coordinates,
    1522                         &stage_vertex_values,
    1523                         &xmom_vertex_values,
    1524                         &ymom_vertex_values,
    1525                         &elevation_vertex_values,
    1526                         &optimise_dry_cells)) {                 
    1527                        
    1528     PyErr_SetString(PyExc_RuntimeError, "Input arguments to extrapolate_second_order_sw failed");
    1529     return NULL;
    1530   }
    1531 
    1532  
    1533   // Get the safety factor beta_w, set in the config.py file. This is used in the limiting process
    1534   Tmp = PyObject_GetAttrString(domain, "beta_w");
    1535   if (!Tmp) {
    1536     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w from domain");
    1537     return NULL;
    1538   } 
    1539   beta_w = PyFloat_AsDouble(Tmp);
    1540   Py_DECREF(Tmp);
    1541  
    1542   Tmp = PyObject_GetAttrString(domain, "beta_w_dry");
    1543   if (!Tmp) {
    1544     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w_dry from domain");
    1545     return NULL;
    1546   } 
    1547   beta_w_dry = PyFloat_AsDouble(Tmp);
    1548   Py_DECREF(Tmp);
    1549  
    1550   Tmp = PyObject_GetAttrString(domain, "beta_uh");
    1551   if (!Tmp) {
    1552     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh from domain");
    1553     return NULL;
    1554   } 
    1555   beta_uh = PyFloat_AsDouble(Tmp);
    1556   Py_DECREF(Tmp);
    1557  
    1558   Tmp = PyObject_GetAttrString(domain, "beta_uh_dry");
    1559   if (!Tmp) {
    1560     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh_dry from domain");
    1561     return NULL;
    1562   } 
    1563   beta_uh_dry = PyFloat_AsDouble(Tmp);
    1564   Py_DECREF(Tmp);
    1565 
    1566   Tmp = PyObject_GetAttrString(domain, "beta_vh");
    1567   if (!Tmp) {
    1568     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh from domain");
    1569     return NULL;
    1570   } 
    1571   beta_vh = PyFloat_AsDouble(Tmp);
    1572   Py_DECREF(Tmp);
    1573  
    1574   Tmp = PyObject_GetAttrString(domain, "beta_vh_dry");
    1575   if (!Tmp) {
    1576     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh_dry from domain");
    1577     return NULL;
    1578   } 
    1579   beta_vh_dry = PyFloat_AsDouble(Tmp);
    1580   Py_DECREF(Tmp);
    1581  
    1582   Tmp = PyObject_GetAttrString(domain, "minimum_allowed_height");
    1583   if (!Tmp) {
    1584     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object minimum_allowed_height");
    1585     return NULL;
    1586   } 
    1587   minimum_allowed_height = PyFloat_AsDouble(Tmp);
    1588   Py_DECREF(Tmp); 
    1589 
    1590   Tmp = PyObject_GetAttrString(domain, "epsilon");
    1591   if (!Tmp) {
    1592     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object epsilon");
    1593     return NULL;
    1594   } 
    1595   epsilon = PyFloat_AsDouble(Tmp);
    1596   Py_DECREF(Tmp); 
    1597    
    1598   number_of_elements = stage_centroid_values -> dimensions[0];
    1599   for (k=0; k<number_of_elements; k++) {
    1600     k3=k*3;
    1601     k6=k*6;
    1602 
    1603    
    1604     if (((long *) number_of_boundaries->data)[k]==3){
    1605       // No neighbours, set gradient on the triangle to zero*/
    1606       ((double *) stage_vertex_values->data)[k3]=((double *)stage_centroid_values->data)[k];
    1607       ((double *) stage_vertex_values->data)[k3+1]=((double *)stage_centroid_values->data)[k];
    1608       ((double *) stage_vertex_values->data)[k3+2]=((double *)stage_centroid_values->data)[k];
    1609       ((double *) xmom_vertex_values->data)[k3]=((double *)xmom_centroid_values->data)[k];
    1610       ((double *) xmom_vertex_values->data)[k3+1]=((double *)xmom_centroid_values->data)[k];
    1611       ((double *) xmom_vertex_values->data)[k3+2]=((double *)xmom_centroid_values->data)[k];
    1612       ((double *) ymom_vertex_values->data)[k3]=((double *)ymom_centroid_values->data)[k];
    1613       ((double *) ymom_vertex_values->data)[k3+1]=((double *)ymom_centroid_values->data)[k];
    1614       ((double *) ymom_vertex_values->data)[k3+2]=((double *)ymom_centroid_values->data)[k];
    1615       continue;
    1616     }
    1617     else {
    1618       // Triangle k has one or more neighbours.
    1619       // Get centroid and vertex coordinates of the triangle
    1620      
    1621       // Get the vertex coordinates
    1622       xv0=((double *)vertex_coordinates->data)[k6]; yv0=((double *)vertex_coordinates->data)[k6+1];
    1623       xv1=((double *)vertex_coordinates->data)[k6+2]; yv1=((double *)vertex_coordinates->data)[k6+3];
    1624       xv2=((double *)vertex_coordinates->data)[k6+4]; yv2=((double *)vertex_coordinates->data)[k6+5];
    1625      
    1626       // Get the centroid coordinates
    1627       coord_index=2*k;
    1628       x=((double *)centroid_coordinates->data)[coord_index];
    1629       y=((double *)centroid_coordinates->data)[coord_index+1];
    1630      
    1631       // Store x- and y- differentials for the vertices of the FV triangle relative to the centroid
    1632       dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
    1633       dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
    1634     }
    1635 
    1636 
    1637    
    1638    
    1639    
    1640            
    1641     if (((long *)number_of_boundaries->data)[k]<=1){
    1642    
    1643       //==============================================
    1644       // Number of boundaries <= 1
    1645       //==============================================   
    1646    
    1647    
    1648       // If no boundaries, auxiliary triangle is formed from the centroids of the three neighbours
    1649       // If one boundary, auxiliary triangle is formed from this centroid and its two neighbours
    1650       k0=((long *)surrogate_neighbours->data)[k3];
    1651       k1=((long *)surrogate_neighbours->data)[k3+1];
    1652       k2=((long *)surrogate_neighbours->data)[k3+2];
    1653      
    1654       // Get the auxiliary triangle's vertex coordinates (really the centroids of neighbouring triangles)
    1655       coord_index=2*k0;
    1656       x0=((double *)centroid_coordinates->data)[coord_index];
    1657       y0=((double *)centroid_coordinates->data)[coord_index+1];
    1658       coord_index=2*k1;
    1659       x1=((double *)centroid_coordinates->data)[coord_index];
    1660       y1=((double *)centroid_coordinates->data)[coord_index+1];
    1661       coord_index=2*k2;
    1662       x2=((double *)centroid_coordinates->data)[coord_index];
    1663       y2=((double *)centroid_coordinates->data)[coord_index+1];
    1664      
    1665       // Store x- and y- differentials for the vertices of the auxiliary triangle
    1666       dx1=x1-x0; dx2=x2-x0;
    1667       dy1=y1-y0; dy2=y2-y0;
    1668      
    1669       // Calculate 2*area of the auxiliary triangle
    1670       area2 = dy2*dx1 - dy1*dx2;//the triangle is guaranteed to be counter-clockwise
    1671      
    1672       // If the mesh is 'weird' near the boundary, the triangle might be flat or clockwise:
    1673       if (area2<=0) {
    1674         PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: negative triangle area encountered");
    1675         return NULL;
    1676       } 
    1677      
    1678       // Calculate heights of neighbouring cells
    1679       hc = ((double *)stage_centroid_values->data)[k]  - ((double *)elevation_centroid_values->data)[k];
    1680       h0 = ((double *)stage_centroid_values->data)[k0] - ((double *)elevation_centroid_values->data)[k0];
    1681       h1 = ((double *)stage_centroid_values->data)[k1] - ((double *)elevation_centroid_values->data)[k1];
    1682       h2 = ((double *)stage_centroid_values->data)[k2] - ((double *)elevation_centroid_values->data)[k2];
    1683       hmin = min(hc,min(h0,min(h1,h2)));  // FIXME Don't need to include hc
    1684      
    1685      
    1686       if (optimise_dry_cells) {     
    1687         // Check if linear reconstruction is necessary for triangle k
    1688         // This check will exclude dry cells.
    1689 
    1690         hmax = max(h0,max(h1,h2));     
    1691         if (hmax < epsilon) {
    1692           continue;
    1693         }
    1694       }
    1695 
    1696            
    1697       //-----------------------------------
    1698       // stage
    1699       //-----------------------------------     
    1700      
    1701       // Calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
    1702       dq0=((double *)stage_centroid_values->data)[k0]-((double *)stage_centroid_values->data)[k];
    1703      
    1704       // Calculate differentials between the vertices of the auxiliary triangle
    1705       dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k0];
    1706       dq2=((double *)stage_centroid_values->data)[k2]-((double *)stage_centroid_values->data)[k0];
    1707      
    1708       // Calculate the gradient of stage on the auxiliary triangle
    1709       a = dy2*dq1 - dy1*dq2;
    1710       a /= area2;
    1711       b = dx1*dq2 - dx2*dq1;
    1712       b /= area2;
    1713      
    1714       // Calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
    1715       dqv[0]=a*dxv0+b*dyv0;
    1716       dqv[1]=a*dxv1+b*dyv1;
    1717       dqv[2]=a*dxv2+b*dyv2;
    1718      
    1719       // Now we want to find min and max of the centroid and the vertices of the auxiliary triangle
    1720       // and compute jumps from the centroid to the min and max
    1721       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
    1722      
    1723       // Playing with dry wet interface
    1724       hmin = qmin;
    1725       beta_tmp = beta_w;
    1726       if (hmin<minimum_allowed_height)
    1727         beta_tmp = beta_w_dry;
    1728       limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
    1729       for (i=0;i<3;i++)
    1730         ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
    1731      
    1732      
    1733       //-----------------------------------
    1734       // xmomentum
    1735       //-----------------------------------           
    1736 
    1737       // Calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
    1738       dq0=((double *)xmom_centroid_values->data)[k0]-((double *)xmom_centroid_values->data)[k];
    1739      
    1740       // Calculate differentials between the vertices of the auxiliary triangle
    1741       dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k0];
    1742       dq2=((double *)xmom_centroid_values->data)[k2]-((double *)xmom_centroid_values->data)[k0];
    1743      
    1744       // Calculate the gradient of xmom on the auxiliary triangle
    1745       a = dy2*dq1 - dy1*dq2;
    1746       a /= area2;
    1747       b = dx1*dq2 - dx2*dq1;
    1748       b /= area2;
    1749      
    1750       // Calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
    1751       dqv[0]=a*dxv0+b*dyv0;
    1752       dqv[1]=a*dxv1+b*dyv1;
    1753       dqv[2]=a*dxv2+b*dyv2;
    1754      
    1755       // Now we want to find min and max of the centroid and the vertices of the auxiliary triangle
    1756       // and compute jumps from the centroid to the min and max
    1757       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
    1758       beta_tmp = beta_uh;
    1759       if (hmin<minimum_allowed_height)
    1760         beta_tmp = beta_uh_dry;
    1761       limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
    1762       for (i=0;i<3;i++)
    1763         ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
    1764      
    1765      
    1766       //-----------------------------------
    1767       // ymomentum
    1768       //-----------------------------------                 
    1769 
    1770       // Calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
    1771       dq0=((double *)ymom_centroid_values->data)[k0]-((double *)ymom_centroid_values->data)[k];
    1772      
    1773       // Calculate differentials between the vertices of the auxiliary triangle
    1774       dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k0];
    1775       dq2=((double *)ymom_centroid_values->data)[k2]-((double *)ymom_centroid_values->data)[k0];
    1776      
    1777       // Calculate the gradient of xmom on the auxiliary triangle
    1778       a = dy2*dq1 - dy1*dq2;
    1779       a /= area2;
    1780       b = dx1*dq2 - dx2*dq1;
    1781       b /= area2;
    1782      
    1783       // Calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
    1784       dqv[0]=a*dxv0+b*dyv0;
    1785       dqv[1]=a*dxv1+b*dyv1;
    1786       dqv[2]=a*dxv2+b*dyv2;
    1787      
    1788       // Now we want to find min and max of the centroid and the vertices of the auxiliary triangle
    1789       // and compute jumps from the centroid to the min and max
    1790       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
    1791       beta_tmp = beta_vh;
    1792       if (hmin<minimum_allowed_height)
    1793         beta_tmp = beta_vh_dry;
    1794       limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
    1795       for (i=0;i<3;i++)
    1796         ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
    1797     } // End number_of_boundaries <=1
    1798     else{
    1799 
    1800       //==============================================
    1801       // Number of boundaries == 2
    1802       //==============================================       
    1803        
    1804       // One internal neighbour and gradient is in direction of the neighbour's centroid
    1805      
    1806       // Find the only internal neighbour
    1807       for (k2=k3;k2<k3+3;k2++){//k2 just indexes the edges of triangle k
    1808         if (((long *)surrogate_neighbours->data)[k2]!=k)//find internal neighbour of triabngle k
    1809           break;
    1810       }
    1811       if ((k2==k3+3)) {//if we didn't find an internal neighbour
    1812         PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: Internal neighbour not found");     
    1813         return NULL;//error
    1814       }
    1815      
    1816       k1=((long *)surrogate_neighbours->data)[k2];
    1817      
    1818       // The coordinates of the triangle are already (x,y). Get centroid of the neighbour (x1,y1)
    1819       coord_index=2*k1;
    1820       x1=((double *)centroid_coordinates->data)[coord_index];
    1821       y1=((double *)centroid_coordinates->data)[coord_index+1];
    1822      
    1823       // Compute x- and y- distances between the centroid of the FV triangle and that of its neighbour
    1824       dx1=x1-x; dy1=y1-y;
    1825      
    1826       // Set area2 as the square of the distance
    1827       area2=dx1*dx1+dy1*dy1;
    1828      
    1829       // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
    1830       // respectively correspond to the x- and y- gradients of the conserved quantities
    1831       dx2=1.0/area2;
    1832       dy2=dx2*dy1;
    1833       dx2*=dx1;
    1834      
    1835      
    1836      
    1837       //-----------------------------------
    1838       // stage
    1839       //-----------------------------------           
    1840 
    1841       // Compute differentials
    1842       dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k];
    1843      
    1844       // Calculate the gradient between the centroid of the FV triangle and that of its neighbour
    1845       a=dq1*dx2;
    1846       b=dq1*dy2;
    1847      
    1848       // Calculate provisional vertex jumps, to be limited
    1849       dqv[0]=a*dxv0+b*dyv0;
    1850       dqv[1]=a*dxv1+b*dyv1;
    1851       dqv[2]=a*dxv2+b*dyv2;
    1852      
    1853       // Now limit the jumps
    1854       if (dq1>=0.0){
    1855         qmin=0.0;
    1856         qmax=dq1;
    1857       }
    1858       else{
    1859         qmin=dq1;
    1860         qmax=0.0;
    1861       }
    1862      
    1863      
    1864       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
    1865       for (i=0;i<3;i++)
    1866         ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
    1867      
    1868       //-----------------------------------
    1869       // xmomentum
    1870       //-----------------------------------                       
    1871      
    1872       // Compute differentials
    1873       dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k];
    1874      
    1875       // Calculate the gradient between the centroid of the FV triangle and that of its neighbour
    1876       a=dq1*dx2;
    1877       b=dq1*dy2;
    1878      
    1879       // Calculate provisional vertex jumps, to be limited
    1880       dqv[0]=a*dxv0+b*dyv0;
    1881       dqv[1]=a*dxv1+b*dyv1;
    1882       dqv[2]=a*dxv2+b*dyv2;
    1883      
    1884       // Now limit the jumps
    1885       if (dq1>=0.0){
    1886         qmin=0.0;
    1887         qmax=dq1;
    1888       }
    1889       else{
    1890         qmin=dq1;
    1891         qmax=0.0;
    1892       }
    1893       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
    1894       for (i=0;i<3;i++)
    1895         ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
    1896      
    1897       //-----------------------------------
    1898       // ymomentum
    1899       //-----------------------------------                       
    1900 
    1901       // Compute differentials
    1902       dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k];
    1903      
    1904       // Calculate the gradient between the centroid of the FV triangle and that of its neighbour
    1905       a=dq1*dx2;
    1906       b=dq1*dy2;
    1907      
    1908       // Calculate provisional vertex jumps, to be limited
    1909       dqv[0]=a*dxv0+b*dyv0;
    1910       dqv[1]=a*dxv1+b*dyv1;
    1911       dqv[2]=a*dxv2+b*dyv2;
    1912      
    1913       // Now limit the jumps
    1914       if (dq1>=0.0){
    1915         qmin=0.0;
    1916         qmax=dq1;
    1917       }
    1918       else{
    1919         qmin=dq1;
    1920         qmax=0.0;
    1921       }
    1922       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
    1923       for (i=0;i<3;i++)
    1924         ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
    1925     }//else [number_of_boundaries==2]
    1926   }//for k=0 to number_of_elements-1
    1927  
    1928   return Py_BuildValue("");
    1929 }//extrapolate_second-order_sw
    19301471
    19311472
    … …  
    22951836}
    22961837
    2297 
    2298 
    2299 
    2300 // THIS FUNCTION IS NOW OBSOLETE
    2301 PyObject *compute_fluxes_ext_central_original(PyObject *self, PyObject *args) {
    2302   /*Compute all fluxes and the timestep suitable for all volumes
    2303     in domain.
    2304 
    2305     Compute total flux for each conserved quantity using "flux_function_central"
    2306 
    2307     Fluxes across each edge are scaled by edgelengths and summed up
    2308     Resulting flux is then scaled by area and stored in
    2309     explicit_update for each of the three conserved quantities
    2310     stage, xmomentum and ymomentum
    2311 
    2312     The maximal allowable speed computed by the flux_function for each volume
    2313     is converted to a timestep that must not be exceeded. The minimum of
    2314     those is computed as the next overall timestep.
    2315 
    2316     Python call:
    2317     domain.timestep = compute_fluxes(timestep,
    2318                                      domain.epsilon,
    2319                                      domain.H0,
    2320                                      domain.g,
    2321                                      domain.neighbours,
    2322                                      domain.neighbour_edges,
    2323                                      domain.normals,
    2324                                      domain.edgelengths,
    2325                                      domain.radii,
    2326                                      domain.areas,
    2327                                      tri_full_flag,
    2328                                      Stage.edge_values,
    2329                                      Xmom.edge_values,
    2330                                      Ymom.edge_values,
    2331                                      Bed.edge_values,
    2332                                      Stage.boundary_values,
    2333                                      Xmom.boundary_values,
    2334                                      Ymom.boundary_values,
    2335                                      Stage.explicit_update,
    2336                                      Xmom.explicit_update,
    2337                                      Ymom.explicit_update,
    2338                                      already_computed_flux,
    2339                                      optimise_dry_cells)                                     
    2340 
    2341 
    2342     Post conditions:
    2343       domain.explicit_update is reset to computed flux values
    2344       domain.timestep is set to the largest step satisfying all volumes.
    2345 
    2346 
    2347   */
    2348 
    2349 
    2350   PyArrayObject *neighbours, *neighbour_edges,
    2351     *normals, *edgelengths, *radii, *areas,
    2352     *tri_full_flag,
    2353     *stage_edge_values,
    2354     *xmom_edge_values,
    2355     *ymom_edge_values,
    2356     *bed_edge_values,
    2357     *stage_boundary_values,
    2358     *xmom_boundary_values,
    2359     *ymom_boundary_values,
    2360     *stage_explicit_update,
    2361     *xmom_explicit_update,
    2362     *ymom_explicit_update,
    2363     *already_computed_flux, //Tracks whether the flux across an edge has already been computed
    2364     *max_speed_array; //Keeps track of max speeds for each triangle
    2365 
    2366 
    2367   // Local variables
    2368   double timestep, max_speed, epsilon, g, H0, length, area;
    2369   int optimise_dry_cells=0; // Optimisation flag 
    2370   double normal[2], ql[3], qr[3], zl, zr;
    2371   double edgeflux[3]; // Work array for summing up fluxes
    2372 
    2373   int number_of_elements, k, i, m, n;
    2374 
    2375   int ki, nm=0, ki2; // Index shorthands
    2376   static long call=1; // Static local variable flagging already computed flux
    2377 
    2378 
    2379   // Convert Python arguments to C
    2380   if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
    2381                         &timestep,
    2382                         &epsilon,
    2383                         &H0,
    2384                         &g,
    2385                         &neighbours,
    2386                         &neighbour_edges,
    2387                         &normals,
    2388                         &edgelengths, &radii, &areas,
    2389                         &tri_full_flag,
    2390                         &stage_edge_values,
    2391                         &xmom_edge_values,
    2392                         &ymom_edge_values,
    2393                         &bed_edge_values,
    2394                         &stage_boundary_values,
    2395                         &xmom_boundary_values,
    2396                         &ymom_boundary_values,
    2397                         &stage_explicit_update,
    2398                         &xmom_explicit_update,
    2399                         &ymom_explicit_update,
    2400                         &already_computed_flux,
    2401                         &max_speed_array,
    2402                         &optimise_dry_cells)) {
    2403     PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    2404     return NULL;
    2405   }
    2406  
    2407  
    2408   number_of_elements = stage_edge_values -> dimensions[0];
    2409 
    2410   call++; // Flag 'id' of flux calculation for this timestep
    2411  
    2412   // Set explicit_update to zero for all conserved_quantities.
    2413   // This assumes compute_fluxes called before forcing terms
    2414   for (k=0; k<number_of_elements; k++) {
    2415     ((double *) stage_explicit_update -> data)[k]=0.0;
    2416     ((double *) xmom_explicit_update -> data)[k]=0.0;
    2417     ((double *) ymom_explicit_update -> data)[k]=0.0; 
    2418   }
    2419  
    2420   // For all triangles
    2421   for (k=0; k<number_of_elements; k++) {
    2422  
    2423     // Loop through neighbours and compute edge flux for each 
    2424     for (i=0; i<3; i++) {
    2425       ki = k*3+i; // Linear index (triangle k, edge i)
    2426      
    2427       if (((long *) already_computed_flux->data)[ki] == call)
    2428         // We've already computed the flux across this edge
    2429         continue;
    2430        
    2431        
    2432       ql[0] = ((double *) stage_edge_values -> data)[ki];
    2433       ql[1] = ((double *) xmom_edge_values -> data)[ki];
    2434       ql[2] = ((double *) ymom_edge_values -> data)[ki];
    2435       zl =    ((double *) bed_edge_values -> data)[ki];
    2436 
    2437       // Quantities at neighbour on nearest face
    2438       n = ((long *) neighbours -> data)[ki];
    2439       if (n < 0) {
    2440         m = -n-1; // Convert negative flag to index
    2441        
    2442         qr[0] = ((double *) stage_boundary_values -> data)[m];
    2443         qr[1] = ((double *) xmom_boundary_values -> data)[m];
    2444         qr[2] = ((double *) ymom_boundary_values -> data)[m];
    2445         zr = zl; //Extend bed elevation to boundary
    2446       } else {
    2447         m = ((long *) neighbour_edges -> data)[ki];
    2448         nm = n*3+m; // Linear index (triangle n, edge m)
    2449        
    2450         qr[0] = ((double *) stage_edge_values -> data)[nm];
    2451         qr[1] = ((double *) xmom_edge_values -> data)[nm];
    2452         qr[2] = ((double *) ymom_edge_values -> data)[nm];
    2453         zr =    ((double *) bed_edge_values -> data)[nm];
    2454       }
    2455      
    2456      
    2457       if (optimise_dry_cells) {     
    2458         // Check if flux calculation is necessary across this edge
    2459         // This check will exclude dry cells.
    2460         // This will also optimise cases where zl != zr as
    2461         // long as both are dry
    2462 
    2463         if ( fabs(ql[0] - zl) < epsilon &&
    2464              fabs(qr[0] - zr) < epsilon ) {
    2465           // Cell boundary is dry
    2466          
    2467           ((long *) already_computed_flux -> data)[ki] = call; // #k Done       
    2468           if (n>=0)
    2469             ((long *) already_computed_flux -> data)[nm] = call; // #n Done
    2470        
    2471           max_speed = 0.0;
    2472           continue;
    2473         }
    2474       }
    2475      
    2476            
    2477       // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
    2478       ki2 = 2*ki; //k*6 + i*2
    2479       normal[0] = ((double *) normals -> data)[ki2];
    2480       normal[1] = ((double *) normals -> data)[ki2+1];
    2481      
    2482 
    2483       // Edge flux computation (triangle k, edge i)
    2484       flux_function_central(ql, qr, zl, zr,
    2485                             normal[0], normal[1],
    2486                             epsilon, H0, g,
    2487                             edgeflux, &max_speed);
    2488      
    2489      
    2490       // Multiply edgeflux by edgelength
    2491       length = ((double *) edgelengths -> data)[ki];
    2492       edgeflux[0] *= length;           
    2493       edgeflux[1] *= length;           
    2494       edgeflux[2] *= length;                       
    2495      
    2496      
    2497       // Update triangle k with flux from edge i
    2498       ((double *) stage_explicit_update -> data)[k] -= edgeflux[0];
    2499       ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1];
    2500       ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2];
    2501      
    2502       ((long *) already_computed_flux -> data)[ki] = call; // #k Done
    2503      
    2504      
    2505       // Update neighbour n with same flux but reversed sign
    2506       if (n>=0){
    2507         ((double *) stage_explicit_update -> data)[n] += edgeflux[0];
    2508         ((double *) xmom_explicit_update -> data)[n] += edgeflux[1];
    2509         ((double *) ymom_explicit_update -> data)[n] += edgeflux[2];
    2510        
    2511         ((long *) already_computed_flux -> data)[nm] = call; // #n Done
    2512       }
    2513      
    2514      
    2515       // Update timestep based on edge i and possibly neighbour n
    2516       if ( ((long *) tri_full_flag -> data)[k] == 1) {
    2517         if (max_speed > epsilon) {
    2518           timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
    2519           if (n>=0)
    2520             timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
    2521         }
    2522       }
    2523      
    2524     } // End edge i
    2525    
    2526    
    2527     // Normalise triangle k by area and store for when all conserved
    2528     // quantities get updated
    2529     area = ((double *) areas -> data)[k];
    2530     ((double *) stage_explicit_update -> data)[k] /= area;
    2531     ((double *) xmom_explicit_update -> data)[k] /= area;
    2532     ((double *) ymom_explicit_update -> data)[k] /= area;
    2533    
    2534    
    2535     // Keep track of maximal speeds
    2536     ((double *) max_speed_array -> data)[k] = max_speed;   
    2537    
    2538   } // End triangle k
    2539  
    2540   return Py_BuildValue("d", timestep);
    2541 }
    25421838
    25431839
    … …  
    27652061PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
    27662062  //
    2767   //    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
     2063  //    balance_deep_and_shallow(beta_h, wc, zc, wv, zv,
    27682064  //                             xmomc, ymomc, xmomv, ymomv)
    27692065
    … …  
    27722068    *wc,            //Stage at centroids
    27732069    *zc,            //Elevation at centroids
    2774     *hc,            //Height at centroids
    27752070    *wv,            //Stage at vertices
    27762071    *zv,            //Elevation at vertices
    2777     *hv,            //Depths at vertices
    27782072    *hvbar,         //h-Limited depths at vertices
    27792073    *xmomc,         //Momentums at centroids and vertices
    … …  
    27852079   
    27862080  double alpha_balance = 2.0;
    2787   double H0;
     2081  double H0, beta_h;
    27882082
    27892083  int N, tight_slope_limiters; //, err;
    27902084
    27912085  // Convert Python arguments to C
    2792   if (!PyArg_ParseTuple(args, "OOOOOOOOOOOO",
     2086  if (!PyArg_ParseTuple(args, "OdOOOOOOOOO",
    27932087                        &domain,
    2794                         &wc, &zc, &hc,
    2795                         &wv, &zv, &hv, &hvbar,
     2088                        &beta_h,
     2089                        &wc, &zc,
     2090                        &wv, &zv, &hvbar,
    27962091                        &xmomc, &ymomc, &xmomv, &ymomv)) {
    2797     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
     2092    PyErr_SetString(PyExc_RuntimeError,
     2093                    "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
    27982094    return NULL;
    27992095  } 
    … …  
    28322128 
    28332129 
    2834      
    2835   //alpha_balance = 2.0;
    2836   //H0 = 0.001;
    2837   //tight_slope_limiters = 1;
    2838  
    28392130  N = wc -> dimensions[0];
    2840 
    28412131  _balance_deep_and_shallow(N,
     2132                            beta_h,
    28422133                            (double*) wc -> data,
    28432134                            (double*) zc -> data,
    2844                             (double*) hc -> data,
    28452135                            (double*) wv -> data,
    28462136                            (double*) zv -> data,
    2847                             (double*) hv -> data,
    28482137                            (double*) hvbar -> data,
    28492138                            (double*) xmomc -> data,
    … …  
    30842373  Py_InitModule("shallow_water_ext", MethodTable);
    30852374
    3086   import_array();     //Necessary for handling of NumPY structures
    3087 }
     2375  import_array(); // Necessary for handling of NumPY structures
     2376}
Note: See TracChangeset for help on using the changeset viewer.