Changeset 8376

Timestamp:

Apr 1, 2012, 11:46:38 PM (13 years ago)

Author:

steve

Message:

Added well balanced flux calc

Location:

trunk/anuga_core/source/anuga/shallow_water

Files:

• shallow_water_domain.py (modified) (8 diffs)
• shallow_water_ext.c (modified) (11 diffs)
• sw_domain.h (modified) (3 diffs)
• test_shallow_water_domain.py (modified) (14 diffs)

trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -64 +64 @@
 SeeAlso:
     TRAC administration of ANUGA (User Manuals etc) at
-    https://datamining.anu.edu.au/anuga and Subversion repository at
-    $HeadURL: https://datamining.anu.edu.au/svn/anuga/trunk/anuga_core/source/
+    https://anuga.anu.edu.au and Subversion repository at
+    $HeadURL: https://anuga.anu.edu.au/svn/anuga/trunk/anuga_core/source/
                 anuga/shallow_water/shallow_water_domain.py $
 
@@ -165 +165 @@
         self.set_defaults()
 
+
         #-------------------------------
         # Operators and Forcing Terms
         #-------------------------------
         self.forcing_terms.append(manning_friction_implicit)
-        #self.forcing_terms.append(gravity_)
+        #self.forcing_terms.append(gravity_wb)
         self.forcing_terms.append(gravity)
 
@@ -219 +220 @@
         from anuga.config import use_edge_limiter
         from anuga.config import use_centroid_velocities
+        from anuga.config import compute_fluxes_method
         
 
@@ -246 +248 @@
         self.use_edge_limiter = use_edge_limiter
         self.optimised_gradient_limiter = optimised_gradient_limiter
-        self.use_centroid_velocities = use_centroid_velocities       
+        self.use_centroid_velocities = use_centroid_velocities
+
+        self.set_compute_fluxes_method(compute_fluxes_method)
 
 
@@ -274 +278 @@
 
 
-    def set_use_edge_limiter(self, flag=True):
-        """Cludge to allow unit test to pass, but to
-        also introduce new edge limiting. The flag is
-        tested in distribute_to_vertices_and_edges.
-        """
-        if flag:
-            self.use_edge_limiter = True
+    def set_compute_fluxes_method(self, flag='original'):
+        """Set method for computing fluxes.
+
+        Currently
+           'original'
+           'well_balanced_1
+        """
+        compute_fluxes_methods = ['original', 'well_balanced_1']
+
+        if flag in compute_fluxes_methods:
+            self.compute_fluxes_method = flag
         else:
-            self.use_edge_limiter = False
+            raise Exception('Unknown compute_fluxes_method')
+
+
+    def get_compute_fluxes_method(self):
+        """Get method for computing fluxes.
+
+        Currently
+           'original'
+           'well_balanced_1
+        """
+
+        return self.compute_fluxes_method
+
+
 
 
@@ -580 +601 @@
 
     def compute_fluxes(self):
-        """Call correct module function
-            (either from this module or C-extension)"""
-        compute_fluxes(self)
+        """Compute fluxes and timestep suitable for all volumes in domain.
+
+        Compute total flux for each conserved quantity using "flux_function"
+
+        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.
+
+        Post conditions:
+          domain.explicit_update is reset to computed flux values
+          domain.flux_timestep is set to the largest step satisfying all volumes.
+
+        This wrapper calls the underlying C version of compute fluxes
+        """
+
+        if self.compute_fluxes_method == 'original':
+            from shallow_water_ext import compute_fluxes_ext_central_structure
+            self.flux_timestep = compute_fluxes_ext_central_structure(self)
+        elif self.compute_fluxes_method == 'well_balanced_1':
+            from shallow_water_ext import compute_fluxes_ext_wb
+            self.flux_timestep = compute_fluxes_ext_wb(self)
+        else:
+            raise Exception('unknown compute_fluxes_method')
 
 
@@ -1137 +1183 @@
 
 
-    compute_fluxes_ext_central_structure(domain)
+    domain.flux_timestep = compute_fluxes_ext_central_structure(domain)
 
 
@@ -1358 +1404 @@
     """
 
-    from shallow_water_ext import gravity_wb_c
+    from shallow_water_ext import gravity_wb as gravity_wb_c
 
     gravity_wb_c(domain)

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

@@ -400 +400 @@
   }
   
+  return 0;
+}
+
+// Innermost flux function (using stage w=z+h)
+int _flux_function_central_wb(double *q_left, double *q_right,
+                    double z_left, double h_left, double h1_left, double h2_left,
+                    double z_right, double h_right, double h1_right, double h2_right,
+                    double n1, double n2,
+                    double epsilon,
+                    double h0,
+                    double limiting_threshold,
+                    double g,
+                    double *edgeflux,
+                    double *max_speed)
+{
+
+  /*Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of the 'stage', w = h+z using
+    the 'central scheme' as described in
+
+    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
+    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
+    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
+
+    The implemented formula is given in equation (3.15) on page 714
+  */
+
+
+  int i;
+  double hl, hr;
+  double w_left, uh_left, vh_left, u_left;
+  double w_right, uh_right, vh_right, u_right;
+  double s_min, s_max, soundspeed_left, soundspeed_right;
+  double denom, inverse_denominator, z;
+  double p_left, p_right;
+
+  // Workspace (allocate once, use many)
+  static double q_left_rotated[3], q_right_rotated[3], flux_right[3], flux_left[3];
+
+  // Copy conserved quantities to protect from modification
+  q_left_rotated[0] = q_left[0];
+  q_right_rotated[0] = q_right[0];
+  q_left_rotated[1] = q_left[1];
+  q_right_rotated[1] = q_right[1];
+  q_left_rotated[2] = q_left[2];
+  q_right_rotated[2] = q_right[2];
+
+  // Align x- and y-momentum with x-axis
+  _rotate(q_left_rotated, n1, n2);
+  _rotate(q_right_rotated, n1, n2);
+
+  z = 0.5*(z_left + z_right); // Average elevation values.
+                              // Even though this will nominally allow
+                  // for discontinuities in the elevation data,
+                  // there is currently no numerical support for
+                  // this so results may be strange near
+                  // jumps in the bed.
+
+/*
+  printf("ql[0] %g\n",q_left[0]);
+  printf("ql[1] %g\n",q_left[1]);
+  printf("ql[2] %g\n",q_left[2]);
+  printf("h_left %g\n",h_left);
+  printf("z_left %g\n",z_left);
+  printf("z %g\n",z);
+  printf("qlr[0] %g\n",q_left_rotated[0]);
+  printf("qlr[1] %g\n",q_left_rotated[1]);
+  printf("qlr[2] %g\n",q_left_rotated[2]);
+*/
+  
+  // Compute speeds in x-direction
+  w_left = q_left_rotated[0];
+  hl = w_left - z;
+/*
+  printf("w_left %g\n",w_left);
+  printf("hl %g\n",hl);
+  printf("hl-h_left %g\n",hl-h_left);
+*/
+  uh_left = q_left_rotated[1];
+  u_left = _compute_speed(&uh_left, &h_left,
+              epsilon, h0, limiting_threshold);
+
+  w_right = q_right_rotated[0];
+  h_right = w_right - z;
+  uh_right = q_right_rotated[1];
+  u_right = _compute_speed(&uh_right, &h_right,
+               epsilon, h0, limiting_threshold);
+
+  // Momentum in y-direction
+  vh_left  = q_left_rotated[2];
+  vh_right = q_right_rotated[2];
+
+  // Limit y-momentum if necessary
+  // Leaving this out, improves speed significantly (Ole 27/5/2009)
+  // All validation tests pass, so do we really need it anymore?
+  _compute_speed(&vh_left, &h_left,
+         epsilon, h0, limiting_threshold);
+  _compute_speed(&vh_right, &h_right,
+         epsilon, h0, limiting_threshold);
+
+  // Maximal and minimal wave speeds
+  soundspeed_left  = sqrt(g*h_left);
+  soundspeed_right = sqrt(g*h_right);
+
+  // Code to use fast square root optimisation if desired.
+  // Timings on AMD 64 for the Okushiri profile gave the following timings
+  //
+  // SQRT           Total    Flux
+  //=============================
+  //
+  // Ref            405s     152s
+  // Fast (dbl)     453s     173s
+  // Fast (sng)     437s     171s
+  //
+  // Consequently, there is currently (14/5/2009) no reason to use this
+  // approximation.
+
+  //soundspeed_left  = fast_squareroot_approximation(g*h_left);
+  //soundspeed_right = fast_squareroot_approximation(g*h_right);
+
+  s_max = max(u_left + soundspeed_left, u_right + soundspeed_right);
+  if (s_max < 0.0)
+  {
+    s_max = 0.0;
+  }
+
+  s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
+  if (s_min > 0.0)
+  {
+    s_min = 0.0;
+  }
+
+
+
+  //p_left = 0.5*g*h_left*h_left;
+  //p_right = 0.5*g*h_right*h_right;
+
+  p_left  = 0.5*g/6.0*(h1_left*h1_left   + 4.0*h_left*h_left   + h2_left*h2_left);
+  p_right = 0.5*g/6.0*(h1_right*h1_right + 4.0*h_right*h_right + h2_right*h2_right);
+
+  // Flux formulas
+  flux_left[0] = u_left*h_left;
+  flux_left[1] = u_left*uh_left + p_left;
+  flux_left[2] = u_left*vh_left;
+
+  flux_right[0] = u_right*h_right;
+  flux_right[1] = u_right*uh_right + p_right;
+  flux_right[2] = u_right*vh_right;
+
+  // Flux computation
+  denom = s_max - s_min;
+  if (denom < epsilon)
+  { // FIXME (Ole): Try using h0 here
+    memset(edgeflux, 0, 3*sizeof(double));
+    *max_speed = 0.0;
+  }
+  else
+  {
+    inverse_denominator = 1.0/denom;
+    for (i = 0; i < 3; i++)
+    {
+      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
+      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i]);
+      edgeflux[i] *= inverse_denominator;
+    }
+
+    // Maximal wavespeed
+    *max_speed = max(fabs(s_max), fabs(s_min));
+
+    // Rotate back
+    _rotate(edgeflux, n1, -n2);
+  }
+
+  return 0;
+}
+
+
+// Innermost flux function (using stage w=z+h)
+int _flux_function_central_wb_old(double *q_left, double *q_right,
+                    double z_left, double h_left, double h1_left, double h2_left,
+                    double z_right, double h_right, double h1_right, double h2_right,
+                    double n1, double n2,
+                    double epsilon,
+                    double h0,
+                    double limiting_threshold,
+                    double g,
+                    double *edgeflux,
+                    double *max_speed)
+{
+
+  /*Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of the 'stage', w = h+z using
+    the 'central scheme' as described in
+
+    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
+    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
+    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
+
+    The implemented formula is given in equation (3.15) on page 714
+  */
+
+  int i;
+
+  double uh_left, vh_left, u_left;
+  double uh_right, vh_right, u_right;
+  double p_left, p_right;
+  double s_min, s_max, soundspeed_left, soundspeed_right;
+  double denom, inverse_denominator;
+
+  // Workspace (allocate once, use many)
+  static double q_left_rotated[3], q_right_rotated[3], flux_right[3], flux_left[3];
+
+  // Copy conserved quantities to protect from modification
+  q_left_rotated[0] = q_left[0];
+  q_right_rotated[0] = q_right[0];
+  q_left_rotated[1] = q_left[1];
+  q_right_rotated[1] = q_right[1];
+  q_left_rotated[2] = q_left[2];
+  q_right_rotated[2] = q_right[2];
+
+  // Align x- and y-momentum with x-axis
+  _rotate(q_left_rotated, n1, n2);
+  _rotate(q_right_rotated, n1, n2);
+
+  //printf("zr-zl %g\n",z_right-z_left);
+
+  // Compute left and right speeds  in x-direction
+  uh_left = q_left_rotated[1];
+  u_left = _compute_speed(&uh_left, &h_left,
+              epsilon, h0, limiting_threshold);
+
+  uh_right = q_right_rotated[1];
+  u_right = _compute_speed(&uh_right, &h_right,
+               epsilon, h0, limiting_threshold);
+
+  // Momentum in y-direction
+  vh_left  = q_left_rotated[2];
+  vh_right = q_right_rotated[2];
+
+
+  // Maximal and minimal wave speeds
+  soundspeed_left  = sqrt(g*h_left);
+  soundspeed_right = sqrt(g*h_right);
+
+
+  s_max = max(u_left + soundspeed_left, u_right + soundspeed_right);
+  if (s_max < 0.0)
+  {
+    s_max = 0.0;
+  }
+
+  s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
+  if (s_min > 0.0)
+  {
+    s_min = 0.0;
+  }
+
+  // Calculate pressure terms using Simpson's rule (exact for pw linear)
+  p_left  = 0.5*g/6.0*(h1_left*h1_left   + 4.0*h_left*h_left   + h2_left*h2_left);
+  p_right = 0.5*g/6.0*(h1_right*h1_right + 4.0*h_right*h_right + h2_right*h2_right);
+
+  p_left = 0.5*g*h_left*h_left;
+  p_right = 0.5*g*h_right*h_right;
+
+  // Flux formulas
+  flux_left[0] = u_left*h_left;
+  flux_left[1] = u_left*uh_left + p_left;
+  flux_left[2] = u_left*vh_left;
+
+  flux_right[0] = u_right*h_right;
+  flux_right[1] = u_right*uh_right + p_right;
+  flux_right[2] = u_right*vh_right;
+
+  // Flux computation
+  denom = s_max - s_min;
+  if (denom < epsilon)
+  { // FIXME (Ole): Try using h0 here
+    memset(edgeflux, 0, 3*sizeof(double));
+    *max_speed = 0.0;
+  }
+  else
+  {
+    inverse_denominator = 1.0/denom;
+    for (i = 0; i < 3; i++)
+    {
+      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
+      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i]);
+      edgeflux[i] *= inverse_denominator;
+    }
+
+    // Maximal wavespeed
+    *max_speed = max(fabs(s_max), fabs(s_min));
+
+    // Rotate back
+    _rotate(edgeflux, n1, -n2);
+  }
+
   return 0;
 }
@@ -1243 +1540 @@
     double g, avg_h, wx, wy, fact;
     double x0, y0, x1, y1, x2, y2;
-    double h[3];
-    double w[3];
+    double h0,h1,h2;
+    double hh[3];
+    double w0,w1,w2;
     double sidex, sidey;
     double n0, n1;
@@ -1269 +1567 @@
         //------------------------------------
 
-        // Get vertex bed values for gradient calculation
-        w[0] = D.stage_vertex_values[k3 + 0];
-        w[1] = D.stage_vertex_values[k3 + 1];
-        w[2] = D.stage_vertex_values[k3 + 2];
+        // Get vertex stage values for gradient calculation
+        w0 = D.stage_vertex_values[k3 + 0];
+        w1 = D.stage_vertex_values[k3 + 1];
+        w2 = D.stage_vertex_values[k3 + 2];
 
         // Compute stage slope
@@ -1285 +1583 @@
 
         //printf("x0 %g, y0 %g, x1 %g, y1 %g, x2 %g, y2 %g \n",x0,y0,x1,y1,x2,y2);
-        _gradient(x0, y0, x1, y1, x2, y2, w[0], w[1], w[2], &wx, &wy);
+        _gradient(x0, y0, x1, y1, x2, y2, w0, w1, w2, &wx, &wy);
 
         avg_h = D.stage_centroid_values[k] - D.bed_centroid_values[k];
@@ -1298 +1596 @@
         
         // Get edge depths
-        h[0] = D.stage_edge_values[k3 + 0] - D.bed_edge_values[k3 + 0];
-        h[1] = D.stage_edge_values[k3 + 1] - D.bed_edge_values[k3 + 1];
-        h[2] = D.stage_edge_values[k3 + 2] - D.bed_edge_values[k3 + 2];
-
+        hh[0] = D.stage_edge_values[k3 + 0] - D.bed_edge_values[k3 + 0];
+        hh[1] = D.stage_edge_values[k3 + 1] - D.bed_edge_values[k3 + 1];
+        hh[2] = D.stage_edge_values[k3 + 2] - D.bed_edge_values[k3 + 2];
+
+
+        printf("h0,1,2 %f %f %f\n",hh[0],hh[1],hh[2]);
+        
         // Calculate the side correction term
         sidex = 0.0;
@@ -1308 +1609 @@
             n0 = D.normals[k6+2*i];
             n1 = D.normals[k6+2*i+1];
-            fact = 0.5*g*D.edgelengths[k3+i]*h[i]*h[i];
-            sidex += fact*n0;
-            sidey += fact*n1;
+
+            printf("n0, n1 %i %g %g\n",i,n0,n1);
+            fact = 0.5*g*hh[i]*hh[i]*D.edgelengths[k3+i];
+            sidex = sidex + fact*n0;
+            sidey = sidey + fact*n1;
         }
 
         // Update momentum with side terms
-        D.xmom_explicit_update[k] += sidex;
-        D.ymom_explicit_update[k] += sidey;
+        D.xmom_explicit_update[k] =- sidex;
+        D.ymom_explicit_update[k] =- sidey;
 
     }
@@ -2431 +2734 @@
 
 
-int _compute_fluxes_central_structure(struct domain *D) {
+double _compute_fluxes_central_structure(struct domain *D) {
 
     // Local variables
@@ -2603 +2906 @@
 
 
-    D->flux_timestep = timestep;
-
-    return 0;
+    return timestep;
 }
+
+double _compute_fluxes_central_wb(struct domain *D) {
+
+    // Local variables
+    double max_speed, length, inv_area, zl, zr;
+    double timestep = 1.0e30;
+    double h0 = D->H0*D->H0; // This ensures a good balance when h approaches H0.
+
+    double limiting_threshold = 10 * D->H0; // Avoid applying limiter below this
+    // threshold for performance reasons.
+    // See ANUGA manual under flux limiting
+    int k, i, m, n;
+    int k3, k3i, k3i1, k3i2, k2i; // Index short hands
+
+    int n3m = 0, n3m1, n3m2;  // Index short hand for neightbours
+
+    double hl, hl1, hl2;
+    double hr, hr1, hr2;
+
+    // Workspace (making them static actually made function slightly slower (Ole))
+    double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
+
+    static long call = 1; // Static local variable flagging already computed flux
+
+    // Start computation
+    call++; // Flag 'id' of flux calculation for this timestep
+
+
+    timestep = D->evolve_max_timestep;
+
+    // Set explicit_update to zero for all conserved_quantities.
+    // This assumes compute_fluxes called before forcing terms
+    memset((char*) D->stage_explicit_update, 0, D->number_of_elements * sizeof (double));
+    memset((char*) D->xmom_explicit_update, 0, D->number_of_elements * sizeof (double));
+    memset((char*) D->ymom_explicit_update, 0, D->number_of_elements * sizeof (double));
+
+    // For all triangles
+    for (k = 0; k < D->number_of_elements; k++) {
+        // Loop through neighbours and compute edge flux for each
+        for (i = 0; i < 3; i++) {
+            k3 = 3*k;
+            k3i = k3 + i; // Linear index to edge i of triangle k
+            k3i1 = k3 + (i+1)%3;
+            k3i2 = k3 + (i+2)%3;
+
+            if (D->already_computed_flux[k3i] == call) {
+                // We've already computed the flux across this edge
+                continue;
+            }
+
+
+            // Get the inside values at the vertices from the triangle k, edge i
+            ql[0] = D->stage_edge_values[k3i];
+            ql[1] = D->xmom_edge_values[k3i];
+            ql[2] = D->ymom_edge_values[k3i];
+
+            zl = D->bed_edge_values[k3i];
+            hl = D->stage_edge_values[k3i] - zl;
+
+            hl1 = D->stage_vertex_values[k3i1] - D->bed_vertex_values[k3i1];
+            hl2 = D->stage_vertex_values[k3i2] - D->bed_vertex_values[k3i2];
+
+            // Get right hand side values either from neighbouring triangle
+            // or from boundary array (Quantities at neighbour on nearest face).
+            n = D->neighbours[k3i];
+            if (n < 0) {
+                // Neighbour is a boundary condition
+                m = -n - 1; // Convert negative flag to boundary index
+
+                qr[0] = D->stage_boundary_values[m];
+                qr[1] = D->xmom_boundary_values[m];
+                qr[2] = D->ymom_boundary_values[m];
+
+                zr = zl; // Extend bed elevation to boundary and assume flat
+                hr = D->stage_boundary_values[m] - zr;
+                hr1 = hr;
+                hr2 = hr;
+
+            }
+            else {
+                // Neighbour is a real triangle
+                m = D->neighbour_edges[k3i];
+                n3m  = 3*n + m; // Linear index (triangle n, edge m)
+                n3m1 = 3*n + (m+1)%3;
+                n3m2 = 3*n + (m+2)%3;
+
+                qr[0] = D->stage_edge_values[n3m];
+                qr[1] = D->xmom_edge_values[n3m];
+                qr[2] = D->ymom_edge_values[n3m];
+
+                zr = D->bed_edge_values[n3m];
+                hr = D->stage_edge_values[n3m] - zr;
+
+                hr1 = D->stage_vertex_values[n3m2] - D->bed_vertex_values[n3m2];
+                hr2 = D->stage_vertex_values[n3m1] - D->bed_vertex_values[n3m1];
+
+            }
+
+            // Now we have values for this edge - both from left and right side.
+
+            if (D->optimise_dry_cells) {
+                // Check if flux calculation is necessary across this edge
+                // This check will exclude dry cells.
+                // This will also optimise cases where zl != zr as
+                // long as both are dry
+
+                if (fabs(ql[0] - zl) < D->epsilon &&
+                        fabs(qr[0] - zr) < D->epsilon) {
+                    // Cell boundary is dry
+
+                    D->already_computed_flux[k3i] = call; // #k Done
+                    if (n >= 0) {
+                        D->already_computed_flux[n3m] = call; // #n Done
+                    }
+
+                    max_speed = 0.0;
+                    continue;
+                }
+            }
+
+
+            if (fabs(zl-zr)>1.0e-10) {
+                report_python_error(AT,"Discontinuous Elevation");
+                return 0.0;
+            }
+
+            // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
+            k2i = 2 * k3i; //k*6 + i*2
+
+
+
+/*
+            _flux_function_central(ql, qr, zl, zr,
+                    D->normals[k2i], D->normals[k2i + 1],
+                    D->epsilon, h0, limiting_threshold, D->g,
+                    edgeflux, &max_speed);
+*/
+
+            // Edge flux computation (triangle k, edge i)
+
+            _flux_function_central_wb(ql, qr,
+                    zl, hl, hl1, hl2,
+                    zr, hr, hr1, hr2,
+                    D->normals[k2i], D->normals[k2i + 1],
+                    D->epsilon, h0, limiting_threshold, D->g,
+                    edgeflux, &max_speed);
+
+
+
+            // Multiply edgeflux by edgelength
+            length = D->edgelengths[k3i];
+            edgeflux[0] *= length;
+            edgeflux[1] *= length;
+            edgeflux[2] *= length;
+
+
+            // Update triangle k with flux from edge i
+            D->stage_explicit_update[k] -= edgeflux[0];
+            D->xmom_explicit_update[k]  -= edgeflux[1];
+            D->ymom_explicit_update[k]  -= edgeflux[2];
+
+            D->already_computed_flux[k3i] = call; // #k Done
+
+
+            // Update neighbour n with same flux but reversed sign
+            if (n >= 0) {
+                D->stage_explicit_update[n] += edgeflux[0];
+                D->xmom_explicit_update[n]  += edgeflux[1];
+                D->ymom_explicit_update[n]  += edgeflux[2];
+
+                D->already_computed_flux[n3m] = call; // #n Done
+            }
+
+            // Update timestep based on edge i and possibly neighbour n
+            if (D->tri_full_flag[k] == 1) {
+                if (max_speed > D->epsilon) {
+                    // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
+
+                    // CFL for triangle k
+                    timestep = min(timestep, D->radii[k] / max_speed);
+
+                    if (n >= 0) {
+                        // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
+                        timestep = min(timestep, D->radii[n] / max_speed);
+                    }
+
+                    // Ted Rigby's suggested less conservative version
+                    //if (n>=0) {
+                    //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
+                    //} else {
+                    //  timestep = min(timestep, radii[k]/max_speed);
+                    // }
+                }
+            }
+
+        } // End edge i (and neighbour n)
+
+
+        // Normalise triangle k by area and store for when all conserved
+        // quantities get updated
+        inv_area = 1.0 / D->areas[k];
+        D->stage_explicit_update[k] *= inv_area;
+        D->xmom_explicit_update[k]  *= inv_area;
+        D->ymom_explicit_update[k]  *= inv_area;
+
+
+        // Keep track of maximal speeds
+        D->max_speed[k] = max_speed;
+
+    } // End triangle k
+
+
+    return timestep;
+}
+
 
 
@@ -2949 +3465 @@
     struct domain D;
     PyObject *domain;
+    double timestep;
     
     if (!PyArg_ParseTuple(args, "O", &domain)) {
@@ -2961 +3478 @@
     // Call underlying flux computation routine and update
     // the explicit update arrays
-    _compute_fluxes_central_structure(&D);
-
-
-    return Py_BuildValue("");
+    timestep = _compute_fluxes_central_structure(&D);
+
+
+    return Py_BuildValue("d", timestep);
 }
+
+
+
+PyObject *compute_fluxes_ext_wb(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
+
+      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_ext_wb(domain, timestep)
+
+
+      Post conditions:
+        domain.explicit_update is reset to computed flux values
+
+      Returns:
+        timestep which is the largest step satisfying all volumes.
+
+
+     */
+
+    struct domain D;
+    PyObject *domain;
+    double timestep;
+
+    if (!PyArg_ParseTuple(args, "O", &domain)) {
+        report_python_error(AT, "could not parse input arguments");
+        return NULL;
+    }
+
+    // populate the C domain structure with pointers
+    // to the python domain data
+    get_python_domain(&D,domain);
+
+    // Call underlying flux computation routine and update
+    // the explicit update arrays
+    timestep = _compute_fluxes_central_wb(&D);
+
+
+    return Py_BuildValue("d", timestep);
+}
+
 
 
@@ -3355 +3926 @@
   {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
   {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
+  {"compute_fluxes_ext_wb", compute_fluxes_ext_wb, METH_VARARGS, "Print out"},
   {"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
   {"compute_fluxes_ext_central_structure", compute_fluxes_ext_central_structure, METH_VARARGS, "Print out"},

trunk/anuga_core/source/anuga/shallow_water/sw_domain.h

@@ -11 +11 @@
 // structure
 struct domain {
+    // Changing these don't change the data in python object
     long    number_of_elements;
     double  epsilon;
@@ -17 +18 @@
     long    optimise_dry_cells;
     double  evolve_max_timestep;
-    double  flux_timestep;
 
+    // The values in the python object will be changed
     long*   neighbours;
     long*   neighbour_edges;
@@ -82 +83 @@
     D->evolve_max_timestep = get_python_double(domain, "evolve_max_timestep");
 
-    D->flux_timestep = get_python_double(domain, "flux_timestep");
-
     neighbours = get_consecutive_array(domain, "neighbours");
     D->neighbours = (long *) neighbours->data;

trunk/anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -692 +692 @@
             assert num.allclose(t, q[0]), msg
 
-    def test_compute_fluxes0(self):
+    def test_compute_fluxes_structure_0(self):
         # Do a full triangle and check that fluxes cancel out for
         # the constant stage case
@@ -775 +775 @@
         assert num.allclose(e0*edgeflux0 + e1*edgeflux1 + e2*edgeflux2, 0.)
 
+
+        from shallow_water_ext import compute_fluxes_ext_central_structure as compute_fluxes
         # Now check that compute_flux yields zeros as well
-        domain.compute_fluxes()
+        compute_fluxes(domain)
 
         for name in ['stage', 'xmomentum', 'ymomentum']:
             assert num.allclose(domain.quantities[name].explicit_update[1], 0)
 
-    def test_compute_fluxes1(self):
+    def test_compute_fluxes_structure_1(self):
         #Use values from previous version
         a = [0.0, 0.0]
@@ -868 +870 @@
         assert num.allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
 
-        domain.compute_fluxes()
+
+
+        from shallow_water_ext import compute_fluxes_ext_central_structure as compute_fluxes
+        # Now check that compute_flux yields zeros as well
+        compute_fluxes(domain)
+
 
         for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
@@ -881 +888 @@
                             [-69.68888889, -35.93333333, 0., 69.68888889])
 
-    def test_compute_fluxes2(self):
+    def test_compute_fluxes_structure_2(self):
         #Random values, incl momentum
         a = [0.0, 0.0]
@@ -952 +959 @@
                        e2*edgeflux2) / domain.areas[1]
 
-        domain.compute_fluxes()
+
+        from shallow_water_ext import compute_fluxes_ext_central_structure as compute_fluxes
+        # Now check that compute_flux yields zeros as well
+        compute_fluxes(domain)
 
         for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
@@ -959 +969 @@
 
     # FIXME (Ole): Need test like this for fluxes in very shallow water.
-    def test_compute_fluxes3(self):
+    def test_compute_fluxes_structure_3(self):
         #Random values, incl momentum
         a = [0.0, 0.0]
@@ -980 +990 @@
 
         zl = zr = -3.75    # Assume constant bed (must be less than stage)
-        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
+        domain.set_quantity('elevation', zl*num.ones((4, 3), num.float)) #array default#
 
         edgeflux = num.zeros(3, num.float)
@@ -1031 +1041 @@
                        e2*edgeflux2) / domain.areas[1]
 
+
+        from shallow_water_ext import compute_fluxes_ext_central_structure as compute_fluxes
+        # Now check that compute_flux yields zeros as well
+        flux_timestep = compute_fluxes(domain)
+
+        #domain.compute_fluxes()
+
+        #print domain.flux_timestep
+        assert num.allclose(flux_timestep, 0.0426244319785)
+
+
+
+        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
+            #print total_flux[i]
+            assert num.allclose(total_flux[i],
+                                domain.quantities[name].explicit_update[1])
+
+
+    def test_compute_fluxes_old_0(self):
+        # Do a full triangle and check that fluxes cancel out for
+        # the constant stage case
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+        domain.set_quantity('stage', [[2,2,2], [2,2,2], [2,2,2], [2,2,2]])
+        domain.check_integrity()
+
+        assert num.allclose(domain.neighbours,
+                            [[-1,1,-2], [2,3,0], [-3,-4,1],[1,-5,-6]])
+        assert num.allclose(domain.neighbour_edges,
+                            [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
+
+        zl = zr = 0.     # Assume flat bed
+
+        edgeflux = num.zeros(3, num.float)
+        edgeflux0 = num.zeros(3, num.float)
+        edgeflux1 = num.zeros(3, num.float)
+        edgeflux2 = num.zeros(3, num.float)
+        H0 = 0.0
+
+        # Flux across right edge of volume 1
+        normal = domain.get_normal(1, 0)
+        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
+        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
+                                  epsilon, g, H0)
+
+        # Check that flux seen from other triangles is inverse
+        (ql, qr) = (qr, ql)
+        normal = domain.get_normal(2, 2)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
+                                  epsilon, g, H0)
+
+        assert num.allclose(edgeflux0 + edgeflux, 0.)
+
+        # Flux across upper edge of volume 1
+        normal = domain.get_normal(1, 1)
+        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
+        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
+                                  epsilon, g, H0)
+
+        # Check that flux seen from other triangles is inverse
+        (ql, qr) = (qr, ql)
+        normal = domain.get_normal(3, 0)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
+                                  epsilon, g, H0)
+
+        assert num.allclose(edgeflux1 + edgeflux, 0.)
+
+        # Flux across lower left hypotenuse of volume 1
+        normal = domain.get_normal(1, 2)
+        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
+        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
+                                  epsilon, g, H0)
+
+        # Check that flux seen from other triangles is inverse
+        (ql, qr) = (qr, ql)
+        normal = domain.get_normal(0, 1)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
+                                  epsilon, g, H0)
+        assert num.allclose(edgeflux2 + edgeflux, 0.)
+
+        # Scale by edgelengths, add up anc check that total flux is zero
+        e0 = domain.edgelengths[1, 0]
+        e1 = domain.edgelengths[1, 1]
+        e2 = domain.edgelengths[1, 2]
+
+        assert num.allclose(e0*edgeflux0 + e1*edgeflux1 + e2*edgeflux2, 0.)
+
+        # Now check that compute_flux yields zeros as well
+        domain.compute_fluxes()
+
+        for name in ['stage', 'xmomentum', 'ymomentum']:
+            assert num.allclose(domain.quantities[name].explicit_update[1], 0)
+
+    def test_compute_fluxes_old_1(self):
+        #Use values from previous version
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+        val0 = 2. + 2.0/3
+        val1 = 4. + 4.0/3
+        val2 = 8. + 2.0/3
+        val3 = 2. + 8.0/3
+
+        domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1],
+                                      [val2, val2, val2], [val3, val3, val3]])
+        domain.check_integrity()
+
+        zl = zr = 0.    # Assume flat bed
+
+        edgeflux = num.zeros(3, num.float)
+        edgeflux0 = num.zeros(3, num.float)
+        edgeflux1 = num.zeros(3, num.float)
+        edgeflux2 = num.zeros(3, num.float)
+        H0 = 0.0
+
+        # Flux across right edge of volume 1
+        normal = domain.get_normal(1, 0)    # Get normal 0 of triangle 1
+        assert num.allclose(normal, [1, 0])
+
+        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
+        assert num.allclose(ql, [val1, 0, 0])
+
+        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
+        assert num.allclose(qr, [val2, 0, 0])
+
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
+                                  epsilon, g, H0)
+
+        # Flux across edge in the east direction (as per normal vector)
+        assert num.allclose(edgeflux0, [-15.3598804, 253.71111111, 0.])
+        assert num.allclose(max_speed, 9.21592824046)
+
+        #Flux across edge in the west direction (opposite sign for xmomentum)
+        normal_opposite = domain.get_normal(2, 2)   # Get normal 2 of triangle 2
+        assert num.allclose(normal_opposite, [-1, 0])
+
+        max_speed = flux_function(normal_opposite, ql, qr, zl, zr, edgeflux,
+                                  epsilon, g, H0)
+        assert num.allclose(edgeflux, [-15.3598804, -253.71111111, 0.])
+
+        #Flux across upper edge of volume 1
+        normal = domain.get_normal(1, 1)
+        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
+        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
+                                  epsilon, g, H0)
+
+        assert num.allclose(edgeflux1, [2.4098563, 0., 123.04444444])
+        assert num.allclose(max_speed, 7.22956891292)
+
+        #Flux across lower left hypotenuse of volume 1
+        normal = domain.get_normal(1, 2)
+        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
+        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
+                                  epsilon, g, H0)
+
+        assert num.allclose(edgeflux2, [9.63942522, -61.59685738, -61.59685738])
+        assert num.allclose(max_speed, 7.22956891292)
+
+        #Scale, add up and check that compute_fluxes is correct for vol 1
+        e0 = domain.edgelengths[1, 0]
+        e1 = domain.edgelengths[1, 1]
+        e2 = domain.edgelengths[1, 2]
+
+        total_flux = -(e0*edgeflux0 +
+                       e1*edgeflux1 +
+                       e2*edgeflux2) / domain.areas[1]
+
+        assert num.allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
+
         domain.compute_fluxes()
 
@@ -1036 +1240 @@
             assert num.allclose(total_flux[i],
                                 domain.quantities[name].explicit_update[1])
+
+        assert num.allclose(domain.quantities['stage'].explicit_update,
+                            [0., -0.68218178, -111.77316251, -35.68522449])
+        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
+                            [-69.68888889, -166.6, 69.68888889, 0])
+        assert num.allclose(domain.quantities['ymomentum'].explicit_update,
+                            [-69.68888889, -35.93333333, 0., 69.68888889])
+
+    def test_compute_fluxes_old_2(self):
+        #Random values, incl momentum
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+        val0 = 2. + 2.0/3
+        val1 = 4. + 4.0/3
+        val2 = 8. + 2.0/3
+        val3 = 2. + 8.0/3
+
+        zl = zr = 0    # Assume flat zero bed
+        edgeflux = num.zeros(3, num.float)
+        edgeflux0 = num.zeros(3, num.float)
+        edgeflux1 = num.zeros(3, num.float)
+        edgeflux2 = num.zeros(3, num.float)
+        H0 = 0.0
+
+        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
+
+        domain.set_quantity('stage', [[val0, val0-1, val0-2],
+                                      [val1, val1+1, val1],
+                                      [val2, val2-2, val2],
+                                      [val3-0.5, val3, val3]])
+
+        domain.set_quantity('xmomentum',
+                            [[1,2,3], [3,4,5], [1,-1,0], [0,-2,2]])
+
+        domain.set_quantity('ymomentum',
+                            [[1,-1,0], [0,-3,2], [0,1,0], [-1,2,2]])
+
+        domain.check_integrity()
+
+        # Flux across right edge of volume 1
+        normal = domain.get_normal(1, 0)
+        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
+        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
+                                  epsilon, g, H0)
+
+        # Flux across upper edge of volume 1
+        normal = domain.get_normal(1, 1)
+        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
+        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
+                                  epsilon, g, H0)
+
+        # Flux across lower left hypotenuse of volume 1
+        normal = domain.get_normal(1, 2)
+        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
+        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
+                                  epsilon, g, H0)
+
+        # Scale, add up and check that compute_fluxes is correct for vol 1
+        e0 = domain.edgelengths[1, 0]
+        e1 = domain.edgelengths[1, 1]
+        e2 = domain.edgelengths[1, 2]
+
+        total_flux = -(e0*edgeflux0 +
+                       e1*edgeflux1 +
+                       e2*edgeflux2) / domain.areas[1]
+
+        domain.compute_fluxes()
+
+        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
+            assert num.allclose(total_flux[i],
+                                domain.quantities[name].explicit_update[1])
+
+    # FIXME (Ole): Need test like this for fluxes in very shallow water.
+    def test_compute_fluxes_old_3(self):
+        #Random values, incl momentum
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        val0 = 2.+2.0/3
+        val1 = 4.+4.0/3
+        val2 = 8.+2.0/3
+        val3 = 2.+8.0/3
+
+        zl = zr = -3.75    # Assume constant bed (must be less than stage)
+        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
+
+        edgeflux = num.zeros(3, num.float)
+        edgeflux0 = num.zeros(3, num.float)
+        edgeflux1 = num.zeros(3, num.float)
+        edgeflux2 = num.zeros(3, num.float)
+        H0 = 0.0
+
+        domain.set_quantity('stage', [[val0, val0-1, val0-2],
+                                      [val1, val1+1, val1],
+                                      [val2, val2-2, val2],
+                                      [val3-0.5, val3, val3]])
+
+        domain.set_quantity('xmomentum',
+                            [[1,2,3], [3,4,5], [1,-1,0], [0,-2,2]])
+
+        domain.set_quantity('ymomentum',
+                            [[1,-1,0], [0,-3,2], [0,1,0], [-1,2,2]])
+
+        domain.check_integrity()
+
+        # Flux across right edge of volume 1
+        normal = domain.get_normal(1, 0)
+        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
+        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
+                                  epsilon, g, H0)
+
+        # Flux across upper edge of volume 1
+        normal = domain.get_normal(1, 1)
+        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
+        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
+                                  epsilon, g, H0)
+
+        # Flux across lower left hypotenuse of volume 1
+        normal = domain.get_normal(1, 2)
+        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
+        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
+                                  epsilon, g, H0)
+
+        # Scale, add up and check that compute_fluxes is correct for vol 1
+        e0 = domain.edgelengths[1, 0]
+        e1 = domain.edgelengths[1, 1]
+        e2 = domain.edgelengths[1, 2]
+
+        total_flux = -(e0*edgeflux0 +
+                       e1*edgeflux1 +
+                       e2*edgeflux2) / domain.areas[1]
+
+        domain.compute_fluxes()
+
+        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
+            assert num.allclose(total_flux[i],
+                                domain.quantities[name].explicit_update[1])
+
+
 
     def xtest_catching_negative_heights(self):
@@ -1789 +2158 @@
         h = 0.1
         def stage(x, y):
-            return slope(x, y) + h
+            return slope(x,y) + h
 
         domain.set_quantity('elevation', slope)
@@ -1799 +2168 @@
 
         domain.compute_forcing_terms()
+        #from shallow_water_ext import gravity
+        #gravity(domain)
+
+        #print domain.quantities['xmomentum'].explicit_update
+        #print domain.quantities['ymomentum'].explicit_update
+
 
         assert num.allclose(domain.quantities['stage'].explicit_update, 0)
@@ -1804 +2179 @@
                             -g*h*3)
         assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
+
+
+
+    def Xtest_gravity_wb(self):
+        #Assuming no friction
+
+        from anuga.config import g
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        #Set up for a gradient of (3,0) at mid triangle (bce)
+        def slope(x, y):
+            return 3*x
+
+        h = 0.1
+        def stage(x, y):
+            return slope(x,y)+h
+
+        domain.set_quantity('elevation', slope)
+        domain.set_quantity('stage', stage)
+
+        for name in domain.conserved_quantities:
+            assert num.allclose(domain.quantities[name].explicit_update, 0)
+            assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
+
+        from shallow_water_ext import gravity_wb
+        gravity_wb(domain)
+
+
+        print domain.quantities['xmomentum'].explicit_update
+        print domain.quantities['ymomentum'].explicit_update
+        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
+        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
+                            -g*h*3)
+        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
+
+
+    def Xtest_gravity_wb_2(self):
+        #Assuming no friction
+
+        from anuga.config import g
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        #Set up for a gradient of (3,0) at mid triangle (bce)
+        def slope(x, y):
+            return 3*x
+
+        h = 15.0
+        def stage(x, y):
+            return h
+
+        domain.set_quantity('elevation', slope)
+        domain.set_quantity('stage', stage)
+
+        for name in domain.conserved_quantities:
+            assert num.allclose(domain.quantities[name].explicit_update, 0)
+            assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
+
+        from shallow_water_ext import gravity_wb
+        gravity_wb(domain)
+
+
+        print domain.quantities['xmomentum'].explicit_update
+        print domain.quantities['ymomentum'].explicit_update
+        
+        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
+        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
+                            -g*h*3)
+        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
+
 
     def test_manning_friction(self):
@@ -6200 +6668 @@
 
         
-        verbose = False
+        verbose = True
         
         #----------------------------------------------------------------------
@@ -6294 +6762 @@
                 #--------------------------------------------------------------
 
-                for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
-                    pass
-                    #if verbose :
-                    #    print domain.timestepping_statistics()
-
-                    #    print domain.volumetric_balance_statistics()                                    
+                for t in domain.evolve(yieldstep=10.0, finaltime=finaltime):
+                    #pass
+                    if verbose :
+                        print domain.timestepping_statistics()
+                        print domain.volumetric_balance_statistics()