Changeset 5306

Timestamp:

May 11, 2008, 8:40:01 AM (17 years ago)

Author:

steve

Message:

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

Files:

• anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py (modified) (1 diff)
• anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py (modified) (2 diffs)
• anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c (modified) (12 diffs)
• anuga_core/source/anuga/shallow_water/shallow_water_domain.py (modified) (4 diffs)
• anuga_core/source/anuga/shallow_water/shallow_water_ext.c (modified) (4 diffs)
• anuga_core/source/anuga/utilities/util_ext.h (modified) (1 diff)
• anuga_core/source/anuga_parallel/run_parallel_advection.py (modified) (1 diff)
• anuga_validation/convergence_study/plot_timeseries.py (added)
• anuga_validation/convergence_study/wave.py (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -306 +306 @@
 
         if self.default_order == 1:
-            self.set_timestepping_method('euler')
+            pass
+            #self.set_timestepping_method('euler')
             #self.set_all_limiters(beta_euler)
 
         if self.default_order == 2:
-            self.set_timestepping_method('rk2')
+            pass
+            #self.set_timestepping_method('rk2')
             #self.set_all_limiters(beta_rk2)
         

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -1409 +1409 @@
         extrapolate_from_gradient(self)
         
-    def extrapolate_second_order_and_limit(self):
+    def extrapolate_second_order_and_limit_by_edge(self):
         #Call correct module function
         #(either from this module or C-extension)
-        extrapolate_second_order_and_limit(self)
+        extrapolate_second_order_and_limit_by_edge(self)
+
+
+    def extrapolate_second_order_and_limit_by_vertex(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        extrapolate_second_order_and_limit_by_vertex(self)
 
     def bound_vertices_below_by_constant(self, bound):
@@ -1468 +1474 @@
          limit_gradient_by_neighbour,\
          extrapolate_from_gradient,\
-         extrapolate_second_order_and_limit,\
+         extrapolate_second_order_and_limit_by_edge,\
+         extrapolate_second_order_and_limit_by_vertex,\
          bound_vertices_below_by_constant,\
          bound_vertices_below_by_quantity,\

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

@@ -266 +266 @@
                      double* vertex_values,
                      double* edge_values,
-                     long*   neighbours) {
+                     long*   neighbours,
+                     double* x_gradient,
+                     double* y_gradient) {
+
 
         int i, k, k2, k3, k6;
@@ -305 +308 @@
                     phi = min( min(r*beta, 1.0), phi);    
                 }
+
+                //Update gradient, vertex and edge values using phi limiter
+                x_gradient[k] = x_gradient[k]*phi;
+                y_gradient[k] = y_gradient[k]*phi;
     
-                //Update vertex and edge values using phi limiter
                 vertex_values[k3+0] = qc + phi*dqa[0];
                 vertex_values[k3+1] = qc + phi*dqa[1];
@@ -1126 +1132 @@
 
 
-PyObject *extrapolate_second_order_and_limit(PyObject *self, PyObject *args) {
+PyObject *extrapolate_second_order_and_limit_by_edge(PyObject *self, PyObject *args) {
     /* Compute edge values using second order approximation and limit values 
        so that edge values are limited by the two corresponding centroid values
@@ -1252 +1258 @@
 
 
+PyObject *extrapolate_second_order_and_limit_by_vertex(PyObject *self, PyObject *args) {
+    /* Compute edge values using second order approximation and limit values 
+       so that edge values are limited by the two corresponding centroid values
+       
+       Python Call:
+       extrapolate_second_order_and_limit(domain,quantity,beta)
+    */
+
+    PyObject *quantity, *domain;
+        
+    PyArrayObject
+        *domain_centroids,              //Coordinates at centroids
+        *domain_vertex_coordinates,     //Coordinates at vertices
+        *domain_number_of_boundaries,   //Number of boundaries for each triangle
+        *domain_surrogate_neighbours,   //True neighbours or - if one missing - self
+        *domain_neighbours,             //True neighbours, or if negative a link to boundary
+
+        *quantity_centroid_values,      //Values at centroids   
+        *quantity_vertex_values,        //Values at vertices
+        *quantity_edge_values,          //Values at edges
+        *quantity_phi,                  //limiter phi values
+        *quantity_x_gradient,           //x gradient
+        *quantity_y_gradient;           //y gradient
+    
+
+  // Local variables
+  int ntri;
+  double beta;
+  int err;
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "O",&quantity)) {
+      PyErr_SetString(PyExc_RuntimeError, 
+                      "quantity_ext.c: extrapolate_second_order_and_limit could not parse input");      
+      return NULL;
+  }
+
+  domain = PyObject_GetAttrString(quantity, "domain");
+  if (!domain) {
+      PyErr_SetString(PyExc_RuntimeError, 
+                      "quantity_ext.c: extrapolate_second_order_and_limit could not obtain domain object from quantity");       
+      return NULL;
+  }
+
+        
+  // Get pertinent variables
+  domain_centroids              = get_consecutive_array(domain,   "centroid_coordinates");
+  domain_surrogate_neighbours   = get_consecutive_array(domain,   "surrogate_neighbours");
+  domain_number_of_boundaries   = get_consecutive_array(domain,   "number_of_boundaries");
+  domain_vertex_coordinates     = get_consecutive_array(domain,   "vertex_coordinates");
+  domain_neighbours             = get_consecutive_array(domain,   "neighbours");
+
+  quantity_centroid_values      = get_consecutive_array(quantity, "centroid_values");
+  quantity_vertex_values        = get_consecutive_array(quantity, "vertex_values");
+  quantity_edge_values          = get_consecutive_array(quantity, "edge_values");
+  quantity_phi                  = get_consecutive_array(quantity, "phi");
+  quantity_x_gradient           = get_consecutive_array(quantity, "x_gradient");
+  quantity_y_gradient           = get_consecutive_array(quantity, "y_gradient");
+
+  beta = get_python_double(quantity,"beta");
+
+  ntri = quantity_centroid_values -> dimensions[0];
+
+  err = _compute_gradients(ntri,
+                           (double*) domain_centroids -> data,
+                           (double*) quantity_centroid_values -> data,
+                           (long*)   domain_number_of_boundaries -> data,
+                           (long*)   domain_surrogate_neighbours -> data,
+                           (double*) quantity_x_gradient -> data,
+                           (double*) quantity_y_gradient -> data);
+
+  if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+                          "quantity_ext.c: Internal function _compute_gradient failed");
+      return NULL;
+  }
+
+
+  err = _extrapolate_from_gradient(ntri, 
+                                   (double*) domain_centroids -> data,
+                                   (double*) quantity_centroid_values -> data,
+                                   (double*) domain_vertex_coordinates -> data,
+                                   (double*) quantity_vertex_values -> data,
+                                   (double*) quantity_edge_values -> data,
+                                   (double*) quantity_x_gradient -> data,
+                                   (double*) quantity_y_gradient -> data);
+
+  if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+                          "quantity_ext.c: Internal function _extrapolate_from_gradient failed");
+      return NULL;
+  }
+
+
+  err = _limit_vertices_by_all_neighbours(ntri, beta,
+                                       (double*) quantity_centroid_values -> data,
+                                       (double*) quantity_vertex_values -> data,
+                                       (double*) quantity_edge_values -> data,
+                                       (long*)   domain_neighbours -> data,
+                                       (double*) quantity_x_gradient -> data,
+                                       (double*) quantity_y_gradient -> data);
+
+  if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+                          "quantity_ext.c: Internal function _limit_vertices_by_all_neighbours failed");
+      return NULL;
+  }
+
+
+  // Release
+  Py_DECREF(domain_centroids);
+  Py_DECREF(domain_surrogate_neighbours);
+  Py_DECREF(domain_number_of_boundaries);
+  Py_DECREF(domain_vertex_coordinates);
+
+  Py_DECREF(quantity_centroid_values);
+  Py_DECREF(quantity_vertex_values);
+  Py_DECREF(quantity_edge_values);
+  Py_DECREF(quantity_phi);
+  Py_DECREF(quantity_x_gradient);
+  Py_DECREF(quantity_y_gradient);
+
+  return Py_BuildValue("");
+}
+
+
 
 PyObject *compute_gradients(PyObject *self, PyObject *args) {
@@ -1445 +1577 @@
             *centroid_values, //Conserved quantities at centroids
             *edge_values,     //Conserved quantities at edges
-            *neighbours;
+            *neighbours,
+            *x_gradient,
+            *y_gradient;
 
         double beta_w; //Safety factor
@@ -1481 +1615 @@
         vertex_values    = get_consecutive_array(quantity, "vertex_values");
         edge_values      = get_consecutive_array(quantity, "edge_values");
-        beta_w           = PyFloat_AsDouble(Tmp);
+        x_gradient       = get_consecutive_array(quantity, "x_gradient");
+        y_gradient       = get_consecutive_array(quantity, "y_gradient");
+        beta_w           = get_python_double(domain,"beta_w");
+
 
 
@@ -1490 +1627 @@
                                                 (double*) vertex_values -> data,
                                                 (double*) edge_values -> data,
-                                                (long*)   neighbours -> data);
+                                                (long*)   neighbours -> data,
+                                                (double*) x_gradient -> data,
+                                                (double*) y_gradient -> data);
+
+
         
         if (err != 0) {
@@ -1504 +1645 @@
         Py_DECREF(vertex_values);
         Py_DECREF(edge_values);
+        Py_DECREF(x_gradient);
+        Py_DECREF(y_gradient);
         Py_DECREF(Tmp);
 
@@ -1557 +1700 @@
         beta_w           = get_python_double(domain,"beta_w");
 
-        Py_DECREF(domain);
+
 
         N = centroid_values -> dimensions[0];
@@ -1581 +1724 @@
         Py_DECREF(vertex_values);
         Py_DECREF(edge_values);
+        Py_DECREF(x_gradient);
+        Py_DECREF(y_gradient);
 
 
@@ -1625 +1770 @@
 
 
-        Py_DECREF(domain);
 
         N = centroid_values -> dimensions[0];
@@ -1921 +2065 @@
         {"extrapolate_from_gradient", extrapolate_from_gradient,
                 METH_VARARGS, "Print out"},
-        {"extrapolate_second_order_and_limit", extrapolate_second_order_and_limit,
+        {"extrapolate_second_order_and_limit_by_edge", extrapolate_second_order_and_limit_by_edge,
+                METH_VARARGS, "Print out"},
+        {"extrapolate_second_order_and_limit_by_vertex", extrapolate_second_order_and_limit_by_vertex,
                 METH_VARARGS, "Print out"},
         {"interpolate_from_vertices_to_edges",

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -401 +401 @@
     def distribute_to_vertices_and_edges(self):
         # Call correct module function
-        # (either from this module or C-extension)
         if self.use_edge_limiter:
-            protect_against_infinitesimal_and_negative_heights(self)
-            for name in self.conserved_quantities:
-                Q = self.quantities[name]
-                if self._order_ == 1:
-                    Q.extrapolate_first_order()
-                elif self._order_ == 2:
-                    Q.extrapolate_second_order_and_limit()
-                else:
-                    raise 'Unknown order'
-            balance_deep_and_shallow(self)
-
-            #Compute edge values by interpolation
-            for name in self.conserved_quantities:
-                Q = self.quantities[name]
-                Q.interpolate_from_vertices_to_edges()
+            distribute_using_edge_limiter(self)            
         else:
-            distribute_to_vertices_and_edges(self)
+            distribute_using_vertex_limiter(self)
 
 
@@ -751 +736 @@
 
 
-def distribute_to_vertices_and_edges(domain):
+def distribute_using_vertex_limiter(domain):
     """Distribution from centroids to vertices specific to the
     shallow water wave
@@ -801 +786 @@
                 Q.extrapolate_first_order()
             elif domain._order_ == 2:
-                Q.extrapolate_second_order()
-                Q.limit()
+                Q.extrapolate_second_order_and_limit_by_vertex()
             else:
                 raise 'Unknown order'
@@ -808 +792 @@
 
     # Take bed elevation into account when water heights are small
+    balance_deep_and_shallow(domain)
+
+    # Compute edge values by interpolation
+    for name in domain.conserved_quantities:
+        Q = domain.quantities[name]
+        Q.interpolate_from_vertices_to_edges()
+
+
+
+def distribute_using_edge_limiter(domain):
+    """Distribution from centroids to edges specific to the
+    shallow water wave
+    equation.
+
+    It will ensure that h (w-z) is always non-negative even in the
+    presence of steep bed-slopes by taking a weighted average between shallow
+    and deep cases.
+
+    In addition, all conserved quantities get distributed as per either a
+    constant (order==1) or a piecewise linear function (order==2).
+
+
+    Precondition:
+      All quantities defined at centroids and bed elevation defined at
+      vertices.
+
+    Postcondition
+      Conserved quantities defined at vertices
+
+    """
+
+    # Remove very thin layers of water
+    protect_against_infinitesimal_and_negative_heights(domain)
+
+
+    for name in domain.conserved_quantities:
+        Q = domain.quantities[name]
+        if domain._order_ == 1:
+            Q.extrapolate_first_order()
+        elif domain._order_ == 2:
+            Q.extrapolate_second_order_and_limit_by_edge()
+        else:
+            raise 'Unknown order'
+
     balance_deep_and_shallow(domain)
 

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -1823 +1823 @@
 }
 
-
-PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
+//=========================================================================
+// Python Glue
+//=========================================================================
+
+PyObject *compute_fluxes_ext_central_new(PyObject *self, PyObject *args) {
   /*Compute all fluxes and the timestep suitable for all volumes
     in domain.
@@ -1840 +1843 @@
 
     Python call:
-    domain.timestep = compute_fluxes(timestep,
-                                     domain.epsilon,
-                                     domain.H0,
-                                     domain.g,
-                                     domain.neighbours,
-                                     domain.neighbour_edges,
-                                     domain.normals,
-                                     domain.edgelengths,
-                                     domain.radii,
-                                     domain.areas,
-                                     tri_full_flag,
-                                     Stage.edge_values,
-                                     Xmom.edge_values,
-                                     Ymom.edge_values,
-                                     Bed.edge_values,
-                                     Stage.boundary_values,
-                                     Xmom.boundary_values,
-                                     Ymom.boundary_values,
-                                     Stage.explicit_update,
-                                     Xmom.explicit_update,
-                                     Ymom.explicit_update,
-                                     already_computed_flux,
-                                     optimise_dry_cells)                                     
+    timestep = compute_fluxes(timestep, domain, stage, xmom, ymom, bed)
 
 
     Post conditions:
       domain.explicit_update is reset to computed flux values
-      domain.timestep is set to the largest step satisfying all volumes.
+      returns timestep which is the largest step satisfying all volumes.
 
 
   */
 
-
-  PyArrayObject *neighbours, *neighbour_edges,
-    *normals, *edgelengths, *radii, *areas,
-    *tri_full_flag,
-    *stage_edge_values,
-    *xmom_edge_values,
-    *ymom_edge_values,
-    *bed_edge_values,
-    *stage_boundary_values,
-    *xmom_boundary_values,
-    *ymom_boundary_values,
-    *stage_explicit_update,
-    *xmom_explicit_update,
-    *ymom_explicit_update,
-    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
-    *max_speed_array; //Keeps track of max speeds for each triangle
-
-    
-  double timestep, epsilon, H0, g;
-  int optimise_dry_cells;
-    
-  // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
-                        &timestep,
-                        &epsilon,
-                        &H0,
-                        &g,
-                        &neighbours,
-                        &neighbour_edges,
-                        &normals,
-                        &edgelengths, &radii, &areas,
-                        &tri_full_flag,
-                        &stage_edge_values,
-                        &xmom_edge_values,
-                        &ymom_edge_values,
-                        &bed_edge_values,
-                        &stage_boundary_values,
-                        &xmom_boundary_values,
-                        &ymom_boundary_values,
-                        &stage_explicit_update,
-                        &xmom_explicit_update,
-                        &ymom_explicit_update,
-                        &already_computed_flux,
-                        &max_speed_array,
-                        &optimise_dry_cells)) {
-    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
-    return NULL;
-  }
-
- 
+    PyObject 
+        *domain,
+        *stage, 
+        *xmom, 
+        *ymom, 
+        *bed;
+
+    PyArrayObject 
+        *neighbours, 
+        *neighbour_edges,
+        *normals, 
+        *edgelengths, 
+        *radii, 
+        *areas,
+        *tri_full_flag,
+        *stage_edge_values,
+        *xmom_edge_values,
+        *ymom_edge_values,
+        *bed_edge_values,
+        *stage_boundary_values,
+        *xmom_boundary_values,
+        *ymom_boundary_values,
+        *stage_explicit_update,
+        *xmom_explicit_update,
+        *ymom_explicit_update,
+        *already_computed_flux, //Tracks whether the flux across an edge has already been computed
+        *max_speed_array; //Keeps track of max speeds for each triangle
+
+    
+    double timestep, epsilon, H0, g;
+    int optimise_dry_cells;
+    
+    // Convert Python arguments to C
+    if (!PyArg_ParseTuple(args, "dOOOO", &timestep, &domain, &stage, &xmom, &ymom, &bed )) {
+        PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
+        return NULL;
+    }
+
+    epsilon           = get_python_double(domain,"epsilon");
+    H0                = get_python_double(domain,"H0");
+    g                 = get_python_double(domain,"g");
+    optimise_dry_cells = get_python_integer(domain,"optimse_dry_cells");
+    
+    neighbours        = get_consecutive_array(domain, "neighbours");
+    neighbour_edges   = get_consecutive_array(domain, "neighbour_edges"); 
+    normals           = get_consecutive_array(domain, "normals");
+    edgelengths       = get_consecutive_array(domain, "edge_lengths");    
+    radii             = get_consecutive_array(domain, "radii");    
+    areas             = get_consecutive_array(domain, "areas");    
+    tri_full_flag     = get_consecutive_array(domain, "normals");
+    already_computed_flux  = get_consecutive_array(domain, "already_computed_flux");
+    max_speed_array   = get_consecutive_array(domain, "max_speed");
+    
+    stage_edge_values = get_consecutive_array(stage, "edge_values");    
+    xmom_edge_values  = get_consecutive_array(xmom, "edge_values");    
+    ymom_edge_values  = get_consecutive_array(ymom, "edge_values"); 
+    bed_edge_values   = get_consecutive_array(bed, "edge_values");    
+
+    stage_boundary_values = get_consecutive_array(stage, "boundary_values");    
+    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");    
+    ymom_boundary_values  = get_consecutive_array(ymom, "boundary_values"); 
+
+    stage_explicit_update = get_consecutive_array(stage, "explicit_update");    
+    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");    
+    ymom_explicit_update  = get_consecutive_array(ymom, "explicit_update"); 
+
+
   int number_of_elements = stage_edge_values -> dimensions[0];
 
@@ -1951 +1949 @@
                                      (double*) max_speed_array -> data,
                                      optimise_dry_cells);
+
+  Py_DECREF(neighbours);
+  Py_DECREF(neighbour_edges);
+  Py_DECREF(normals);
+  Py_DECREF(edgelengths);
+  Py_DECREF(radii);
+  Py_DECREF(areas);
+  Py_DECREF(tri_full_flag);
+  Py_DECREF(already_computed_flux);
+  Py_DECREF(max_speed_array);
+  Py_DECREF(stage_edge_values);
+  Py_DECREF(xmom_edge_values);
+  Py_DECREF(ymom_edge_values);
+  Py_DECREF(bed_edge_values);
+  Py_DECREF(stage_boundary_values);
+  Py_DECREF(xmom_boundary_values);
+  Py_DECREF(ymom_boundary_values);
+  Py_DECREF(stage_explicit_update);
+  Py_DECREF(xmom_explicit_update);
+  Py_DECREF(ymom_explicit_update);
+
   
   // Return updated flux timestep
   return Py_BuildValue("d", timestep);
 }
+
+
+
+
+
+
+PyObject *compute_fluxes_ext_central(PyObject *self, PyObject *args) {
+  /*Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Compute total flux for each conserved quantity using "flux_function_central"
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for each of the three conserved quantities
+    stage, xmomentum and ymomentum
+
+    The maximal allowable speed computed by the flux_function for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Python call:
+    domain.timestep = compute_fluxes(timestep,
+                                     domain.epsilon,
+                                     domain.H0,
+                                     domain.g,
+                                     domain.neighbours,
+                                     domain.neighbour_edges,
+                                     domain.normals,
+                                     domain.edgelengths,
+                                     domain.radii,
+                                     domain.areas,
+                                     tri_full_flag,
+                                     Stage.edge_values,
+                                     Xmom.edge_values,
+                                     Ymom.edge_values,
+                                     Bed.edge_values,
+                                     Stage.boundary_values,
+                                     Xmom.boundary_values,
+                                     Ymom.boundary_values,
+                                     Stage.explicit_update,
+                                     Xmom.explicit_update,
+                                     Ymom.explicit_update,
+                                     already_computed_flux,
+                                     optimise_dry_cells)                                     
+
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+
+
+  */
+
+
+  PyArrayObject *neighbours, *neighbour_edges,
+    *normals, *edgelengths, *radii, *areas,
+    *tri_full_flag,
+    *stage_edge_values,
+    *xmom_edge_values,
+    *ymom_edge_values,
+    *bed_edge_values,
+    *stage_boundary_values,
+    *xmom_boundary_values,
+    *ymom_boundary_values,
+    *stage_explicit_update,
+    *xmom_explicit_update,
+    *ymom_explicit_update,
+    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
+    *max_speed_array; //Keeps track of max speeds for each triangle
+
+    
+  double timestep, epsilon, H0, g;
+  int optimise_dry_cells;
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
+                        &timestep,
+                        &epsilon,
+                        &H0,
+                        &g,
+                        &neighbours,
+                        &neighbour_edges,
+                        &normals,
+                        &edgelengths, &radii, &areas,
+                        &tri_full_flag,
+                        &stage_edge_values,
+                        &xmom_edge_values,
+                        &ymom_edge_values,
+                        &bed_edge_values,
+                        &stage_boundary_values,
+                        &xmom_boundary_values,
+                        &ymom_boundary_values,
+                        &stage_explicit_update,
+                        &xmom_explicit_update,
+                        &ymom_explicit_update,
+                        &already_computed_flux,
+                        &max_speed_array,
+                        &optimise_dry_cells)) {
+    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
+    return NULL;
+  }
+
+ 
+  int number_of_elements = stage_edge_values -> dimensions[0];
+
+  // Call underlying flux computation routine and update 
+  // the explicit update arrays 
+  timestep = _compute_fluxes_central(number_of_elements,
+                                     timestep,
+                                     epsilon,
+                                     H0,
+                                     g,
+                                     (long*) neighbours -> data,
+                                     (long*) neighbour_edges -> data,
+                                     (double*) normals -> data,
+                                     (double*) edgelengths -> data, 
+                                     (double*) radii -> data, 
+                                     (double*) areas -> data,
+                                     (long*) tri_full_flag -> data,
+                                     (double*) stage_edge_values -> data,
+                                     (double*) xmom_edge_values -> data,
+                                     (double*) ymom_edge_values -> data,
+                                     (double*) bed_edge_values -> data,
+                                     (double*) stage_boundary_values -> data,
+                                     (double*) xmom_boundary_values -> data,
+                                     (double*) ymom_boundary_values -> data,
+                                     (double*) stage_explicit_update -> data,
+                                     (double*) xmom_explicit_update -> data,
+                                     (double*) ymom_explicit_update -> data,
+                                     (long*) already_computed_flux -> data,
+                                     (double*) max_speed_array -> data,
+                                     optimise_dry_cells);
+  
+  // Return updated flux timestep
+  return Py_BuildValue("d", timestep);
+}
+
+
 
 
@@ -2490 +2648 @@
   {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
   {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
+  {"compute_fluxes_ext_central_new", compute_fluxes_ext_central_new, METH_VARARGS, "Print out"},
   {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
   {"gravity", gravity, METH_VARARGS, "Print out"},

anuga_core/source/anuga/utilities/util_ext.h

@@ -298 +298 @@
 }
 
+int get_python_integer(PyObject *O, char *name) {
+  PyObject *TObject;
+  int tmp;
+  
+
+  //Get double from attribute
+  TObject = PyObject_GetAttrString(O, name);
+  if (!TObject) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_integer could not obtain double from object");
+    return 0;
+  }  
+  
+  tmp = PyFloat_AsDouble(TObject);
+  
+  Py_DECREF(TObject);
+  
+  return tmp;
+}
+
 
 PyObject *get_python_object(PyObject *O, char *name) {

anuga_core/source/anuga_parallel/run_parallel_advection.py

@@ -102 +102 @@
     t0 = time.time()
 
-# Start the evolve computions
-import hotshot
-profiler = hotshot.Profile("hotshot." + str(numprocs) + "." + str(myid) + ".prof")
-s = '''for t in domain.evolve(yieldstep = 0.1, finaltime = 2.0):
-  if myid == 0:
-    domain.write_time()
-'''
-    
-
-result = profiler.runctx(s, globals(), locals())
-profiler.close()
-
-# for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
-#     if myid == 0:
-#         domain.write_time()
+for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
+    if myid == 0:
+        domain.write_time()
 
 # Output some computation statistics

anuga_validation/convergence_study/wave.py

@@ -39 +39 @@
 #start_screen_catcher(output_dir+sep)
 
+interactive_visualisation = False
+
 #------------------------------------------------------------------------------
 # Setup domain
@@ -77 +79 @@
 domain.beta_h    = 0.0
 
-interactive_visualisation = True