Changeset 5306


Ignore:
Timestamp:
May 11, 2008, 8:40:01 AM (16 years ago)
Author:
steve
Message:

Working version of wave.py for Ole to work on.

Files:
1 added
8 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

    r5242 r5306  
    306306
    307307        if self.default_order == 1:
    308             self.set_timestepping_method('euler')
     308            pass
     309            #self.set_timestepping_method('euler')
    309310            #self.set_all_limiters(beta_euler)
    310311
    311312        if self.default_order == 2:
    312             self.set_timestepping_method('rk2')
     313            pass
     314            #self.set_timestepping_method('rk2')
    313315            #self.set_all_limiters(beta_rk2)
    314316       
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r5162 r5306  
    14091409        extrapolate_from_gradient(self)
    14101410       
    1411     def extrapolate_second_order_and_limit(self):
     1411    def extrapolate_second_order_and_limit_by_edge(self):
    14121412        #Call correct module function
    14131413        #(either from this module or C-extension)
    1414         extrapolate_second_order_and_limit(self)
     1414        extrapolate_second_order_and_limit_by_edge(self)
     1415
     1416
     1417    def extrapolate_second_order_and_limit_by_vertex(self):
     1418        #Call correct module function
     1419        #(either from this module or C-extension)
     1420        extrapolate_second_order_and_limit_by_vertex(self)
    14151421
    14161422    def bound_vertices_below_by_constant(self, bound):
     
    14681474         limit_gradient_by_neighbour,\
    14691475         extrapolate_from_gradient,\
    1470          extrapolate_second_order_and_limit,\
     1476         extrapolate_second_order_and_limit_by_edge,\
     1477         extrapolate_second_order_and_limit_by_vertex,\
    14711478         bound_vertices_below_by_constant,\
    14721479         bound_vertices_below_by_quantity,\
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r5162 r5306  
    266266                     double* vertex_values,
    267267                     double* edge_values,
    268                      long*   neighbours) {
     268                     long*   neighbours,
     269                     double* x_gradient,
     270                     double* y_gradient) {
     271
    269272
    270273        int i, k, k2, k3, k6;
     
    305308                    phi = min( min(r*beta, 1.0), phi);   
    306309                }
     310
     311                //Update gradient, vertex and edge values using phi limiter
     312                x_gradient[k] = x_gradient[k]*phi;
     313                y_gradient[k] = y_gradient[k]*phi;
    307314   
    308                 //Update vertex and edge values using phi limiter
    309315                vertex_values[k3+0] = qc + phi*dqa[0];
    310316                vertex_values[k3+1] = qc + phi*dqa[1];
     
    11261132
    11271133
    1128 PyObject *extrapolate_second_order_and_limit(PyObject *self, PyObject *args) {
     1134PyObject *extrapolate_second_order_and_limit_by_edge(PyObject *self, PyObject *args) {
    11291135    /* Compute edge values using second order approximation and limit values
    11301136       so that edge values are limited by the two corresponding centroid values
     
    12521258
    12531259
     1260PyObject *extrapolate_second_order_and_limit_by_vertex(PyObject *self, PyObject *args) {
     1261    /* Compute edge values using second order approximation and limit values
     1262       so that edge values are limited by the two corresponding centroid values
     1263       
     1264       Python Call:
     1265       extrapolate_second_order_and_limit(domain,quantity,beta)
     1266    */
     1267
     1268    PyObject *quantity, *domain;
     1269       
     1270    PyArrayObject
     1271        *domain_centroids,              //Coordinates at centroids
     1272        *domain_vertex_coordinates,     //Coordinates at vertices
     1273        *domain_number_of_boundaries,   //Number of boundaries for each triangle
     1274        *domain_surrogate_neighbours,   //True neighbours or - if one missing - self
     1275        *domain_neighbours,             //True neighbours, or if negative a link to boundary
     1276
     1277        *quantity_centroid_values,      //Values at centroids   
     1278        *quantity_vertex_values,        //Values at vertices
     1279        *quantity_edge_values,          //Values at edges
     1280        *quantity_phi,                  //limiter phi values
     1281        *quantity_x_gradient,           //x gradient
     1282        *quantity_y_gradient;           //y gradient
     1283   
     1284
     1285  // Local variables
     1286  int ntri;
     1287  double beta;
     1288  int err;
     1289   
     1290  // Convert Python arguments to C
     1291  if (!PyArg_ParseTuple(args, "O",&quantity)) {
     1292      PyErr_SetString(PyExc_RuntimeError,
     1293                      "quantity_ext.c: extrapolate_second_order_and_limit could not parse input");     
     1294      return NULL;
     1295  }
     1296
     1297  domain = PyObject_GetAttrString(quantity, "domain");
     1298  if (!domain) {
     1299      PyErr_SetString(PyExc_RuntimeError,
     1300                      "quantity_ext.c: extrapolate_second_order_and_limit could not obtain domain object from quantity");       
     1301      return NULL;
     1302  }
     1303
     1304       
     1305  // Get pertinent variables
     1306  domain_centroids              = get_consecutive_array(domain,   "centroid_coordinates");
     1307  domain_surrogate_neighbours   = get_consecutive_array(domain,   "surrogate_neighbours");
     1308  domain_number_of_boundaries   = get_consecutive_array(domain,   "number_of_boundaries");
     1309  domain_vertex_coordinates     = get_consecutive_array(domain,   "vertex_coordinates");
     1310  domain_neighbours             = get_consecutive_array(domain,   "neighbours");
     1311
     1312  quantity_centroid_values      = get_consecutive_array(quantity, "centroid_values");
     1313  quantity_vertex_values        = get_consecutive_array(quantity, "vertex_values");
     1314  quantity_edge_values          = get_consecutive_array(quantity, "edge_values");
     1315  quantity_phi                  = get_consecutive_array(quantity, "phi");
     1316  quantity_x_gradient           = get_consecutive_array(quantity, "x_gradient");
     1317  quantity_y_gradient           = get_consecutive_array(quantity, "y_gradient");
     1318
     1319  beta = get_python_double(quantity,"beta");
     1320
     1321  ntri = quantity_centroid_values -> dimensions[0];
     1322
     1323  err = _compute_gradients(ntri,
     1324                           (double*) domain_centroids -> data,
     1325                           (double*) quantity_centroid_values -> data,
     1326                           (long*)   domain_number_of_boundaries -> data,
     1327                           (long*)   domain_surrogate_neighbours -> data,
     1328                           (double*) quantity_x_gradient -> data,
     1329                           (double*) quantity_y_gradient -> data);
     1330
     1331  if (err != 0) {
     1332      PyErr_SetString(PyExc_RuntimeError,
     1333                          "quantity_ext.c: Internal function _compute_gradient failed");
     1334      return NULL;
     1335  }
     1336
     1337
     1338  err = _extrapolate_from_gradient(ntri,
     1339                                   (double*) domain_centroids -> data,
     1340                                   (double*) quantity_centroid_values -> data,
     1341                                   (double*) domain_vertex_coordinates -> data,
     1342                                   (double*) quantity_vertex_values -> data,
     1343                                   (double*) quantity_edge_values -> data,
     1344                                   (double*) quantity_x_gradient -> data,
     1345                                   (double*) quantity_y_gradient -> data);
     1346
     1347  if (err != 0) {
     1348      PyErr_SetString(PyExc_RuntimeError,
     1349                          "quantity_ext.c: Internal function _extrapolate_from_gradient failed");
     1350      return NULL;
     1351  }
     1352
     1353
     1354  err = _limit_vertices_by_all_neighbours(ntri, beta,
     1355                                       (double*) quantity_centroid_values -> data,
     1356                                       (double*) quantity_vertex_values -> data,
     1357                                       (double*) quantity_edge_values -> data,
     1358                                       (long*)   domain_neighbours -> data,
     1359                                       (double*) quantity_x_gradient -> data,
     1360                                       (double*) quantity_y_gradient -> data);
     1361
     1362  if (err != 0) {
     1363      PyErr_SetString(PyExc_RuntimeError,
     1364                          "quantity_ext.c: Internal function _limit_vertices_by_all_neighbours failed");
     1365      return NULL;
     1366  }
     1367
     1368
     1369  // Release
     1370  Py_DECREF(domain_centroids);
     1371  Py_DECREF(domain_surrogate_neighbours);
     1372  Py_DECREF(domain_number_of_boundaries);
     1373  Py_DECREF(domain_vertex_coordinates);
     1374
     1375  Py_DECREF(quantity_centroid_values);
     1376  Py_DECREF(quantity_vertex_values);
     1377  Py_DECREF(quantity_edge_values);
     1378  Py_DECREF(quantity_phi);
     1379  Py_DECREF(quantity_x_gradient);
     1380  Py_DECREF(quantity_y_gradient);
     1381
     1382  return Py_BuildValue("");
     1383}
     1384
     1385
    12541386
    12551387PyObject *compute_gradients(PyObject *self, PyObject *args) {
     
    14451577            *centroid_values, //Conserved quantities at centroids
    14461578            *edge_values,     //Conserved quantities at edges
    1447             *neighbours;
     1579            *neighbours,
     1580            *x_gradient,
     1581            *y_gradient;
    14481582
    14491583        double beta_w; //Safety factor
     
    14811615        vertex_values    = get_consecutive_array(quantity, "vertex_values");
    14821616        edge_values      = get_consecutive_array(quantity, "edge_values");
    1483         beta_w           = PyFloat_AsDouble(Tmp);
     1617        x_gradient       = get_consecutive_array(quantity, "x_gradient");
     1618        y_gradient       = get_consecutive_array(quantity, "y_gradient");
     1619        beta_w           = get_python_double(domain,"beta_w");
     1620
    14841621
    14851622
     
    14901627                                                (double*) vertex_values -> data,
    14911628                                                (double*) edge_values -> data,
    1492                                                 (long*)   neighbours -> data);
     1629                                                (long*)   neighbours -> data,
     1630                                                (double*) x_gradient -> data,
     1631                                                (double*) y_gradient -> data);
     1632
     1633
    14931634       
    14941635        if (err != 0) {
     
    15041645        Py_DECREF(vertex_values);
    15051646        Py_DECREF(edge_values);
     1647        Py_DECREF(x_gradient);
     1648        Py_DECREF(y_gradient);
    15061649        Py_DECREF(Tmp);
    15071650
     
    15571700        beta_w           = get_python_double(domain,"beta_w");
    15581701
    1559         Py_DECREF(domain);
     1702
    15601703
    15611704        N = centroid_values -> dimensions[0];
     
    15811724        Py_DECREF(vertex_values);
    15821725        Py_DECREF(edge_values);
     1726        Py_DECREF(x_gradient);
     1727        Py_DECREF(y_gradient);
    15831728
    15841729
     
    16251770
    16261771
    1627         Py_DECREF(domain);
    16281772
    16291773        N = centroid_values -> dimensions[0];
     
    19212065        {"extrapolate_from_gradient", extrapolate_from_gradient,
    19222066                METH_VARARGS, "Print out"},
    1923         {"extrapolate_second_order_and_limit", extrapolate_second_order_and_limit,
     2067        {"extrapolate_second_order_and_limit_by_edge", extrapolate_second_order_and_limit_by_edge,
     2068                METH_VARARGS, "Print out"},
     2069        {"extrapolate_second_order_and_limit_by_vertex", extrapolate_second_order_and_limit_by_vertex,
    19242070                METH_VARARGS, "Print out"},
    19252071        {"interpolate_from_vertices_to_edges",
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r5303 r5306  
    401401    def distribute_to_vertices_and_edges(self):
    402402        # Call correct module function
    403         # (either from this module or C-extension)
    404403        if self.use_edge_limiter:
    405             protect_against_infinitesimal_and_negative_heights(self)
    406             for name in self.conserved_quantities:
    407                 Q = self.quantities[name]
    408                 if self._order_ == 1:
    409                     Q.extrapolate_first_order()
    410                 elif self._order_ == 2:
    411                     Q.extrapolate_second_order_and_limit()
    412                 else:
    413                     raise 'Unknown order'
    414             balance_deep_and_shallow(self)
    415 
    416             #Compute edge values by interpolation
    417             for name in self.conserved_quantities:
    418                 Q = self.quantities[name]
    419                 Q.interpolate_from_vertices_to_edges()
     404            distribute_using_edge_limiter(self)           
    420405        else:
    421             distribute_to_vertices_and_edges(self)
     406            distribute_using_vertex_limiter(self)
    422407
    423408
     
    751736
    752737
    753 def distribute_to_vertices_and_edges(domain):
     738def distribute_using_vertex_limiter(domain):
    754739    """Distribution from centroids to vertices specific to the
    755740    shallow water wave
     
    801786                Q.extrapolate_first_order()
    802787            elif domain._order_ == 2:
    803                 Q.extrapolate_second_order()
    804                 Q.limit()
     788                Q.extrapolate_second_order_and_limit_by_vertex()
    805789            else:
    806790                raise 'Unknown order'
     
    808792
    809793    # Take bed elevation into account when water heights are small
     794    balance_deep_and_shallow(domain)
     795
     796    # Compute edge values by interpolation
     797    for name in domain.conserved_quantities:
     798        Q = domain.quantities[name]
     799        Q.interpolate_from_vertices_to_edges()
     800
     801
     802
     803def distribute_using_edge_limiter(domain):
     804    """Distribution from centroids to edges specific to the
     805    shallow water wave
     806    equation.
     807
     808    It will ensure that h (w-z) is always non-negative even in the
     809    presence of steep bed-slopes by taking a weighted average between shallow
     810    and deep cases.
     811
     812    In addition, all conserved quantities get distributed as per either a
     813    constant (order==1) or a piecewise linear function (order==2).
     814
     815
     816    Precondition:
     817      All quantities defined at centroids and bed elevation defined at
     818      vertices.
     819
     820    Postcondition
     821      Conserved quantities defined at vertices
     822
     823    """
     824
     825    # Remove very thin layers of water
     826    protect_against_infinitesimal_and_negative_heights(domain)
     827
     828
     829    for name in domain.conserved_quantities:
     830        Q = domain.quantities[name]
     831        if domain._order_ == 1:
     832            Q.extrapolate_first_order()
     833        elif domain._order_ == 2:
     834            Q.extrapolate_second_order_and_limit_by_edge()
     835        else:
     836            raise 'Unknown order'
     837
    810838    balance_deep_and_shallow(domain)
    811839
  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r5297 r5306  
    18231823}
    18241824
    1825 
    1826 PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
     1825//=========================================================================
     1826// Python Glue
     1827//=========================================================================
     1828
     1829PyObject *compute_fluxes_ext_central_new(PyObject *self, PyObject *args) {
    18271830  /*Compute all fluxes and the timestep suitable for all volumes
    18281831    in domain.
     
    18401843
    18411844    Python call:
    1842     domain.timestep = compute_fluxes(timestep,
    1843                                      domain.epsilon,
    1844                                      domain.H0,
    1845                                      domain.g,
    1846                                      domain.neighbours,
    1847                                      domain.neighbour_edges,
    1848                                      domain.normals,
    1849                                      domain.edgelengths,
    1850                                      domain.radii,
    1851                                      domain.areas,
    1852                                      tri_full_flag,
    1853                                      Stage.edge_values,
    1854                                      Xmom.edge_values,
    1855                                      Ymom.edge_values,
    1856                                      Bed.edge_values,
    1857                                      Stage.boundary_values,
    1858                                      Xmom.boundary_values,
    1859                                      Ymom.boundary_values,
    1860                                      Stage.explicit_update,
    1861                                      Xmom.explicit_update,
    1862                                      Ymom.explicit_update,
    1863                                      already_computed_flux,
    1864                                      optimise_dry_cells)                                     
     1845    timestep = compute_fluxes(timestep, domain, stage, xmom, ymom, bed)
    18651846
    18661847
    18671848    Post conditions:
    18681849      domain.explicit_update is reset to computed flux values
    1869       domain.timestep is set to the largest step satisfying all volumes.
     1850      returns timestep which is the largest step satisfying all volumes.
    18701851
    18711852
    18721853  */
    18731854
    1874 
    1875   PyArrayObject *neighbours, *neighbour_edges,
    1876     *normals, *edgelengths, *radii, *areas,
    1877     *tri_full_flag,
    1878     *stage_edge_values,
    1879     *xmom_edge_values,
    1880     *ymom_edge_values,
    1881     *bed_edge_values,
    1882     *stage_boundary_values,
    1883     *xmom_boundary_values,
    1884     *ymom_boundary_values,
    1885     *stage_explicit_update,
    1886     *xmom_explicit_update,
    1887     *ymom_explicit_update,
    1888     *already_computed_flux, //Tracks whether the flux across an edge has already been computed
    1889     *max_speed_array; //Keeps track of max speeds for each triangle
    1890 
    1891    
    1892   double timestep, epsilon, H0, g;
    1893   int optimise_dry_cells;
    1894    
    1895   // Convert Python arguments to C
    1896   if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
    1897                         &timestep,
    1898                         &epsilon,
    1899                         &H0,
    1900                         &g,
    1901                         &neighbours,
    1902                         &neighbour_edges,
    1903                         &normals,
    1904                         &edgelengths, &radii, &areas,
    1905                         &tri_full_flag,
    1906                         &stage_edge_values,
    1907                         &xmom_edge_values,
    1908                         &ymom_edge_values,
    1909                         &bed_edge_values,
    1910                         &stage_boundary_values,
    1911                         &xmom_boundary_values,
    1912                         &ymom_boundary_values,
    1913                         &stage_explicit_update,
    1914                         &xmom_explicit_update,
    1915                         &ymom_explicit_update,
    1916                         &already_computed_flux,
    1917                         &max_speed_array,
    1918                         &optimise_dry_cells)) {
    1919     PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    1920     return NULL;
    1921   }
    1922 
    1923  
     1855    PyObject
     1856        *domain,
     1857        *stage,
     1858        *xmom,
     1859        *ymom,
     1860        *bed;
     1861
     1862    PyArrayObject
     1863        *neighbours,
     1864        *neighbour_edges,
     1865        *normals,
     1866        *edgelengths,
     1867        *radii,
     1868        *areas,
     1869        *tri_full_flag,
     1870        *stage_edge_values,
     1871        *xmom_edge_values,
     1872        *ymom_edge_values,
     1873        *bed_edge_values,
     1874        *stage_boundary_values,
     1875        *xmom_boundary_values,
     1876        *ymom_boundary_values,
     1877        *stage_explicit_update,
     1878        *xmom_explicit_update,
     1879        *ymom_explicit_update,
     1880        *already_computed_flux, //Tracks whether the flux across an edge has already been computed
     1881        *max_speed_array; //Keeps track of max speeds for each triangle
     1882
     1883   
     1884    double timestep, epsilon, H0, g;
     1885    int optimise_dry_cells;
     1886   
     1887    // Convert Python arguments to C
     1888    if (!PyArg_ParseTuple(args, "dOOOO", &timestep, &domain, &stage, &xmom, &ymom, &bed )) {
     1889        PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     1890        return NULL;
     1891    }
     1892
     1893    epsilon           = get_python_double(domain,"epsilon");
     1894    H0                = get_python_double(domain,"H0");
     1895    g                 = get_python_double(domain,"g");
     1896    optimise_dry_cells = get_python_integer(domain,"optimse_dry_cells");
     1897   
     1898    neighbours        = get_consecutive_array(domain, "neighbours");
     1899    neighbour_edges   = get_consecutive_array(domain, "neighbour_edges");
     1900    normals           = get_consecutive_array(domain, "normals");
     1901    edgelengths       = get_consecutive_array(domain, "edge_lengths");   
     1902    radii             = get_consecutive_array(domain, "radii");   
     1903    areas             = get_consecutive_array(domain, "areas");   
     1904    tri_full_flag     = get_consecutive_array(domain, "normals");
     1905    already_computed_flux  = get_consecutive_array(domain, "already_computed_flux");
     1906    max_speed_array   = get_consecutive_array(domain, "max_speed");
     1907   
     1908    stage_edge_values = get_consecutive_array(stage, "edge_values");   
     1909    xmom_edge_values  = get_consecutive_array(xmom, "edge_values");   
     1910    ymom_edge_values  = get_consecutive_array(ymom, "edge_values");
     1911    bed_edge_values   = get_consecutive_array(bed, "edge_values");   
     1912
     1913    stage_boundary_values = get_consecutive_array(stage, "boundary_values");   
     1914    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");   
     1915    ymom_boundary_values  = get_consecutive_array(ymom, "boundary_values");
     1916
     1917    stage_explicit_update = get_consecutive_array(stage, "explicit_update");   
     1918    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");   
     1919    ymom_explicit_update  = get_consecutive_array(ymom, "explicit_update");
     1920
     1921
    19241922  int number_of_elements = stage_edge_values -> dimensions[0];
    19251923
     
    19511949                                     (double*) max_speed_array -> data,
    19521950                                     optimise_dry_cells);
     1951
     1952  Py_DECREF(neighbours);
     1953  Py_DECREF(neighbour_edges);
     1954  Py_DECREF(normals);
     1955  Py_DECREF(edgelengths);
     1956  Py_DECREF(radii);
     1957  Py_DECREF(areas);
     1958  Py_DECREF(tri_full_flag);
     1959  Py_DECREF(already_computed_flux);
     1960  Py_DECREF(max_speed_array);
     1961  Py_DECREF(stage_edge_values);
     1962  Py_DECREF(xmom_edge_values);
     1963  Py_DECREF(ymom_edge_values);
     1964  Py_DECREF(bed_edge_values);
     1965  Py_DECREF(stage_boundary_values);
     1966  Py_DECREF(xmom_boundary_values);
     1967  Py_DECREF(ymom_boundary_values);
     1968  Py_DECREF(stage_explicit_update);
     1969  Py_DECREF(xmom_explicit_update);
     1970  Py_DECREF(ymom_explicit_update);
     1971
    19531972 
    19541973  // Return updated flux timestep
    19551974  return Py_BuildValue("d", timestep);
    19561975}
     1976
     1977
     1978
     1979
     1980
     1981
     1982PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
     1983  /*Compute all fluxes and the timestep suitable for all volumes
     1984    in domain.
     1985
     1986    Compute total flux for each conserved quantity using "flux_function_central"
     1987
     1988    Fluxes across each edge are scaled by edgelengths and summed up
     1989    Resulting flux is then scaled by area and stored in
     1990    explicit_update for each of the three conserved quantities
     1991    stage, xmomentum and ymomentum
     1992
     1993    The maximal allowable speed computed by the flux_function for each volume
     1994    is converted to a timestep that must not be exceeded. The minimum of
     1995    those is computed as the next overall timestep.
     1996
     1997    Python call:
     1998    domain.timestep = compute_fluxes(timestep,
     1999                                     domain.epsilon,
     2000                                     domain.H0,
     2001                                     domain.g,
     2002                                     domain.neighbours,
     2003                                     domain.neighbour_edges,
     2004                                     domain.normals,
     2005                                     domain.edgelengths,
     2006                                     domain.radii,
     2007                                     domain.areas,
     2008                                     tri_full_flag,
     2009                                     Stage.edge_values,
     2010                                     Xmom.edge_values,
     2011                                     Ymom.edge_values,
     2012                                     Bed.edge_values,
     2013                                     Stage.boundary_values,
     2014                                     Xmom.boundary_values,
     2015                                     Ymom.boundary_values,
     2016                                     Stage.explicit_update,
     2017                                     Xmom.explicit_update,
     2018                                     Ymom.explicit_update,
     2019                                     already_computed_flux,
     2020                                     optimise_dry_cells)                                     
     2021
     2022
     2023    Post conditions:
     2024      domain.explicit_update is reset to computed flux values
     2025      domain.timestep is set to the largest step satisfying all volumes.
     2026
     2027
     2028  */
     2029
     2030
     2031  PyArrayObject *neighbours, *neighbour_edges,
     2032    *normals, *edgelengths, *radii, *areas,
     2033    *tri_full_flag,
     2034    *stage_edge_values,
     2035    *xmom_edge_values,
     2036    *ymom_edge_values,
     2037    *bed_edge_values,
     2038    *stage_boundary_values,
     2039    *xmom_boundary_values,
     2040    *ymom_boundary_values,
     2041    *stage_explicit_update,
     2042    *xmom_explicit_update,
     2043    *ymom_explicit_update,
     2044    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
     2045    *max_speed_array; //Keeps track of max speeds for each triangle
     2046
     2047   
     2048  double timestep, epsilon, H0, g;
     2049  int optimise_dry_cells;
     2050   
     2051  // Convert Python arguments to C
     2052  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
     2053                        &timestep,
     2054                        &epsilon,
     2055                        &H0,
     2056                        &g,
     2057                        &neighbours,
     2058                        &neighbour_edges,
     2059                        &normals,
     2060                        &edgelengths, &radii, &areas,
     2061                        &tri_full_flag,
     2062                        &stage_edge_values,
     2063                        &xmom_edge_values,
     2064                        &ymom_edge_values,
     2065                        &bed_edge_values,
     2066                        &stage_boundary_values,
     2067                        &xmom_boundary_values,
     2068                        &ymom_boundary_values,
     2069                        &stage_explicit_update,
     2070                        &xmom_explicit_update,
     2071                        &ymom_explicit_update,
     2072                        &already_computed_flux,
     2073                        &max_speed_array,
     2074                        &optimise_dry_cells)) {
     2075    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     2076    return NULL;
     2077  }
     2078
     2079 
     2080  int number_of_elements = stage_edge_values -> dimensions[0];
     2081
     2082  // Call underlying flux computation routine and update
     2083  // the explicit update arrays
     2084  timestep = _compute_fluxes_central(number_of_elements,
     2085                                     timestep,
     2086                                     epsilon,
     2087                                     H0,
     2088                                     g,
     2089                                     (long*) neighbours -> data,
     2090                                     (long*) neighbour_edges -> data,
     2091                                     (double*) normals -> data,
     2092                                     (double*) edgelengths -> data,
     2093                                     (double*) radii -> data,
     2094                                     (double*) areas -> data,
     2095                                     (long*) tri_full_flag -> data,
     2096                                     (double*) stage_edge_values -> data,
     2097                                     (double*) xmom_edge_values -> data,
     2098                                     (double*) ymom_edge_values -> data,
     2099                                     (double*) bed_edge_values -> data,
     2100                                     (double*) stage_boundary_values -> data,
     2101                                     (double*) xmom_boundary_values -> data,
     2102                                     (double*) ymom_boundary_values -> data,
     2103                                     (double*) stage_explicit_update -> data,
     2104                                     (double*) xmom_explicit_update -> data,
     2105                                     (double*) ymom_explicit_update -> data,
     2106                                     (long*) already_computed_flux -> data,
     2107                                     (double*) max_speed_array -> data,
     2108                                     optimise_dry_cells);
     2109 
     2110  // Return updated flux timestep
     2111  return Py_BuildValue("d", timestep);
     2112}
     2113
     2114
    19572115
    19582116
     
    24902648  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
    24912649  {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
     2650  {"compute_fluxes_ext_central_new", compute_fluxes_ext_central_new, METH_VARARGS, "Print out"},
    24922651  {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
    24932652  {"gravity", gravity, METH_VARARGS, "Print out"},
  • anuga_core/source/anuga/utilities/util_ext.h

    r5162 r5306  
    298298}
    299299
     300int get_python_integer(PyObject *O, char *name) {
     301  PyObject *TObject;
     302  int tmp;
     303 
     304
     305  //Get double from attribute
     306  TObject = PyObject_GetAttrString(O, name);
     307  if (!TObject) {
     308    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_integer could not obtain double from object");
     309    return 0;
     310  } 
     311 
     312  tmp = PyFloat_AsDouble(TObject);
     313 
     314  Py_DECREF(TObject);
     315 
     316  return tmp;
     317}
     318
    300319
    301320PyObject *get_python_object(PyObject *O, char *name) {
  • anuga_core/source/anuga_parallel/run_parallel_advection.py

    r5243 r5306  
    102102    t0 = time.time()
    103103
    104 # Start the evolve computions
    105 import hotshot
    106 profiler = hotshot.Profile("hotshot." + str(numprocs) + "." + str(myid) + ".prof")
    107 s = '''for t in domain.evolve(yieldstep = 0.1, finaltime = 2.0):
    108   if myid == 0:
    109     domain.write_time()
    110 '''
    111    
    112 
    113 result = profiler.runctx(s, globals(), locals())
    114 profiler.close()
    115 
    116 # for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
    117 #     if myid == 0:
    118 #         domain.write_time()
     104for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
     105    if myid == 0:
     106        domain.write_time()
    119107
    120108# Output some computation statistics
  • anuga_validation/convergence_study/wave.py

    r5300 r5306  
    3939#start_screen_catcher(output_dir+sep)
    4040
     41interactive_visualisation = False
     42
    4143#------------------------------------------------------------------------------
    4244# Setup domain
     
    7779domain.beta_h    = 0.0
    7880
    79 interactive_visualisation = True
    8081
    8182
Note: See TracChangeset for help on using the changeset viewer.