Changeset 8456

Timestamp:

Jul 6, 2012, 5:45:58 PM (13 years ago)

Author:

steve

Message:

Adding in some of gareth's methods

Location:

trunk/anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/generic_domain.py (modified) (1 diff)
• operators/kinematic_viscosity_operator.py (modified) (2 diffs)
• shallow_water/boundaries.py (modified) (1 diff)
• shallow_water/shallow_water_domain.py (modified) (4 diffs)
• shallow_water/shallow_water_ext.c (modified) (9 diffs)
• shallow_water/sw_domain.h (modified) (6 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (2 diffs)

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py

@@ -100 +100 @@
         self.centroid_coordinates = self.mesh.centroid_coordinates
         self.vertex_coordinates = self.mesh.vertex_coordinates
+        self.edge_coordinates = self.mesh.edge_midpoint_coordinates
         self.boundary = self.mesh.boundary
         self.boundary_enumeration = self.mesh.boundary_enumeration

trunk/anuga_core/source/anuga/operators/kinematic_viscosity_operator.py

@@ -27 +27 @@
     """
 
-    def __init__(self, domain, use_triangle_areas=True, verbose=False):
+    def __init__(self, domain, diffusivity=None, use_triangle_areas=True, verbose=False):
+
         if verbose: log.critical('Kinematic Viscosity: Beginning Initialisation')
         
@@ -37 +38 @@
         self.boundary = domain.boundary
         self.boundary_enumeration = domain.boundary_enumeration
-        
+
+        self.diffusivity = diffusivity
         # Setup a quantity as diffusivity
-        # FIXME SR: Could/Should pass a quantity which already exists
-        self.diffusivity = Quantity(self.domain)
-        self.diffusivity.set_values(1.0)
-        self.diffusivity.set_boundary_values(1.0)
+        if diffusivity is None:
+            self.diffusivity = Quantity(self.domain)
+            self.diffusivity.set_values(1.0)
+            self.diffusivity.set_boundary_values(1.0)
+
+        if isinstance(diffusivity, (int, float)):
+            self.diffusivity = Quantity(self.domain)
+            self.diffusivity.set_values(diffusivity)
+            self.diffusivity.set_boundary_values(diffusivity)
+
+    
+        assert isinstance(self.diffusivity, Quantity)
         
 

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

@@ -117 +117 @@
 
         self.domain = domain
+
+        if isinstance(function, (int, float)):
+            tmp = function
+            function = lambda t: tmp
+            
         self.function = function
 

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

@@ -798 +798 @@
         """Call correct module function
             (either from this module or C-extension)"""
-        extrapolate_second_order_sw(self)
+        extrapolate_second_order_sw_old(self)
 
     def compute_fluxes(self):
@@ -1436 +1436 @@
 ################################################################################
 
-def extrapolate_second_order_sw(domain):
+def extrapolate_second_order_sw_old(domain):
     """Wrapper calling C version of extrapolate_second_order_sw.
 
@@ -1445 +1445 @@
     """
 
-    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
+    from shallow_water_ext import extrapolate_second_order_sw_old as extrapol2
 
     # Shortcuts
@@ -1468 +1468 @@
               int(domain.optimise_dry_cells),
               int(domain.extrapolate_velocity_second_order))
+
+
+def extrapolate_second_order_sw(domain):
+    """Wrapper calling C version of extrapolate_second_order_sw.
+
+    domain  the domain to operate on
+
+    Note MH090605: The following method belongs to the shallow_water domain
+    class, see comments in the corresponding method in shallow_water_ext.c
+    """
+
+    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
+    extrapol2(domain)
+
 
 def distribute_using_vertex_limiter(domain):

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

@@ -2060 +2060 @@
 // Computational routine
 
+int _extrapolate_second_order_sw_old(int number_of_elements,
+        double epsilon,
+        double minimum_allowed_height,
+        double beta_w,
+        double beta_w_dry,
+        double beta_uh,
+        double beta_uh_dry,
+        double beta_vh,
+        double beta_vh_dry,
+        long* surrogate_neighbours,
+        long* number_of_boundaries,
+        double* centroid_coordinates,
+        double* stage_centroid_values,
+        double* xmom_centroid_values,
+        double* ymom_centroid_values,
+        double* elevation_centroid_values,
+        double* vertex_coordinates,
+        double* stage_vertex_values,
+        double* xmom_vertex_values,
+        double* ymom_vertex_values,
+        double* elevation_vertex_values,
+        int optimise_dry_cells,
+        int extrapolate_velocity_second_order) {
+
+
+
+    // Local variables
+    double a, b; // Gradient vector used to calculate vertex values from centroids
+    int k, k0, k1, k2, k3, k6, coord_index, i;
+    double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
+    double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
+    double dqv[3], qmin, qmax, hmin, hmax;
+    double hc, h0, h1, h2, beta_tmp, hfactor;
+    double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements], dk, dv0, dv1, dv2;
+
+    if (extrapolate_velocity_second_order == 1) {
+        // Replace momentum centroid with velocity centroid to allow velocity
+        // extrapolation This will be changed back at the end of the routine
+        for (k = 0; k < number_of_elements; k++) {
+
+            dk = max(stage_centroid_values[k] - elevation_centroid_values[k], minimum_allowed_height);
+            xmom_centroid_store[k] = xmom_centroid_values[k];
+            xmom_centroid_values[k] = xmom_centroid_values[k] / dk;
+
+            ymom_centroid_store[k] = ymom_centroid_values[k];
+            ymom_centroid_values[k] = ymom_centroid_values[k] / dk;
+        }
+    }
+
+    // Begin extrapolation routine
+    for (k = 0; k < number_of_elements; k++) {
+        k3 = k * 3;
+        k6 = k * 6;
+
+        if (number_of_boundaries[k] == 3) {
+            // No neighbours, set gradient on the triangle to zero
+
+            stage_vertex_values[k3] = stage_centroid_values[k];
+            stage_vertex_values[k3 + 1] = stage_centroid_values[k];
+            stage_vertex_values[k3 + 2] = stage_centroid_values[k];
+            xmom_vertex_values[k3] = xmom_centroid_values[k];
+            xmom_vertex_values[k3 + 1] = xmom_centroid_values[k];
+            xmom_vertex_values[k3 + 2] = xmom_centroid_values[k];
+            ymom_vertex_values[k3] = ymom_centroid_values[k];
+            ymom_vertex_values[k3 + 1] = ymom_centroid_values[k];
+            ymom_vertex_values[k3 + 2] = ymom_centroid_values[k];
+
+            continue;
+        } else {
+            // Triangle k has one or more neighbours.
+            // Get centroid and vertex coordinates of the triangle
+
+            // Get the vertex coordinates
+            xv0 = vertex_coordinates[k6];
+            yv0 = vertex_coordinates[k6 + 1];
+            xv1 = vertex_coordinates[k6 + 2];
+            yv1 = vertex_coordinates[k6 + 3];
+            xv2 = vertex_coordinates[k6 + 4];
+            yv2 = vertex_coordinates[k6 + 5];
+
+            // Get the centroid coordinates
+            coord_index = 2 * k;
+            x = centroid_coordinates[coord_index];
+            y = centroid_coordinates[coord_index + 1];
+
+            // Store x- and y- differentials for the vertices of
+            // triangle k relative to the centroid
+            dxv0 = xv0 - x;
+            dxv1 = xv1 - x;
+            dxv2 = xv2 - x;
+            dyv0 = yv0 - y;
+            dyv1 = yv1 - y;
+            dyv2 = yv2 - y;
+        }
+
+
+
+        if (number_of_boundaries[k] <= 1) {
+            //==============================================
+            // Number of boundaries <= 1
+            //==============================================
+
+
+            // If no boundaries, auxiliary triangle is formed
+            // from the centroids of the three neighbours
+            // If one boundary, auxiliary triangle is formed
+            // from this centroid and its two neighbours
+
+            k0 = surrogate_neighbours[k3];
+            k1 = surrogate_neighbours[k3 + 1];
+            k2 = surrogate_neighbours[k3 + 2];
+
+            // Get the auxiliary triangle's vertex coordinates
+            // (really the centroids of neighbouring triangles)
+            coord_index = 2 * k0;
+            x0 = centroid_coordinates[coord_index];
+            y0 = centroid_coordinates[coord_index + 1];
+
+            coord_index = 2 * k1;
+            x1 = centroid_coordinates[coord_index];
+            y1 = centroid_coordinates[coord_index + 1];
+
+            coord_index = 2 * k2;
+            x2 = centroid_coordinates[coord_index];
+            y2 = centroid_coordinates[coord_index + 1];
+
+            // Store x- and y- differentials for the vertices
+            // of the auxiliary triangle
+            dx1 = x1 - x0;
+            dx2 = x2 - x0;
+            dy1 = y1 - y0;
+            dy2 = y2 - y0;
+
+            // Calculate 2*area of the auxiliary triangle
+            // The triangle is guaranteed to be counter-clockwise
+            area2 = dy2 * dx1 - dy1*dx2;
+
+            // If the mesh is 'weird' near the boundary,
+            // the triangle might be flat or clockwise
+            // Default to zero gradient
+            if (area2 <= 0) {
+                //printf("Error negative triangle area \n");
+                //report_python_error(AT, "Negative triangle area");
+                //return -1;
+
+                stage_vertex_values[k3] = stage_centroid_values[k];
+                stage_vertex_values[k3 + 1] = stage_centroid_values[k];
+                stage_vertex_values[k3 + 2] = stage_centroid_values[k];
+                xmom_vertex_values[k3] = xmom_centroid_values[k];
+                xmom_vertex_values[k3 + 1] = xmom_centroid_values[k];
+                xmom_vertex_values[k3 + 2] = xmom_centroid_values[k];
+                ymom_vertex_values[k3] = ymom_centroid_values[k];
+                ymom_vertex_values[k3 + 1] = ymom_centroid_values[k];
+                ymom_vertex_values[k3 + 2] = ymom_centroid_values[k];
+
+                continue;
+            }
+
+            // Calculate heights of neighbouring cells
+            hc = stage_centroid_values[k] - elevation_centroid_values[k];
+            h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
+            h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
+            h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
+            hmin = min(min(h0, min(h1, h2)), hc);
+            //hfactor = hc/(hc + 1.0);
+
+            hfactor = 0.0;
+            if (hmin > 0.001) {
+                hfactor = (hmin - 0.001) / (hmin + 0.004);
+            }
+
+            if (optimise_dry_cells) {
+                // Check if linear reconstruction is necessary for triangle k
+                // This check will exclude dry cells.
+
+                hmax = max(h0, max(h1, h2));
+                if (hmax < epsilon) {
+                    continue;
+                }
+            }
+
+            //-----------------------------------
+            // stage
+            //-----------------------------------
+
+            // Calculate the difference between vertex 0 of the auxiliary
+            // triangle and the centroid of triangle k
+            dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
+
+            // Calculate differentials between the vertices
+            // of the auxiliary triangle (centroids of neighbouring triangles)
+            dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
+            dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
+
+            inv_area2 = 1.0 / area2;
+            // Calculate the gradient of stage on the auxiliary triangle
+            a = dy2 * dq1 - dy1*dq2;
+            a *= inv_area2;
+            b = dx1 * dq2 - dx2*dq1;
+            b *= inv_area2;
+
+            // Calculate provisional jumps in stage from the centroid
+            // of triangle k to its vertices, to be limited
+            dqv[0] = a * dxv0 + b*dyv0;
+            dqv[1] = a * dxv1 + b*dyv1;
+            dqv[2] = a * dxv2 + b*dyv2;
+
+            // Now we want to find min and max of the centroid and the
+            // vertices of the auxiliary triangle and compute jumps
+            // from the centroid to the min and max
+            find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+
+            // Playing with dry wet interface
+            //hmin = qmin;
+            //beta_tmp = beta_w_dry;
+            //if (hmin>minimum_allowed_height)
+            beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
+
+            //printf("min_alled_height = %f\n",minimum_allowed_height);
+            //printf("hmin = %f\n",hmin);
+            //printf("beta_w = %f\n",beta_w);
+            //printf("beta_tmp = %f\n",beta_tmp);
+            // Limit the gradient
+            limit_gradient(dqv, qmin, qmax, beta_tmp);
+
+            //for (i=0;i<3;i++)
+            stage_vertex_values[k3 + 0] = stage_centroid_values[k] + dqv[0];
+            stage_vertex_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
+            stage_vertex_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
+
+
+            //-----------------------------------
+            // xmomentum
+            //-----------------------------------
+
+            // Calculate the difference between vertex 0 of the auxiliary
+            // triangle and the centroid of triangle k
+            dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
+
+            // Calculate differentials between the vertices
+            // of the auxiliary triangle
+            dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
+            dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
+
+            // Calculate the gradient of xmom on the auxiliary triangle
+            a = dy2 * dq1 - dy1*dq2;
+            a *= inv_area2;
+            b = dx1 * dq2 - dx2*dq1;
+            b *= inv_area2;
+
+            // Calculate provisional jumps in stage from the centroid
+            // of triangle k to its vertices, to be limited
+            dqv[0] = a * dxv0 + b*dyv0;
+            dqv[1] = a * dxv1 + b*dyv1;
+            dqv[2] = a * dxv2 + b*dyv2;
+
+            // Now we want to find min and max of the centroid and the
+            // vertices of the auxiliary triangle and compute jumps
+            // from the centroid to the min and max
+            find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+            //beta_tmp = beta_uh;
+            //if (hmin<minimum_allowed_height)
+            //beta_tmp = beta_uh_dry;
+            beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
+
+            // Limit the gradient
+            limit_gradient(dqv, qmin, qmax, beta_tmp);
+
+            for (i = 0; i < 3; i++) {
+                xmom_vertex_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
+            }
+
+            //-----------------------------------
+            // ymomentum
+            //-----------------------------------
+
+            // Calculate the difference between vertex 0 of the auxiliary
+            // triangle and the centroid of triangle k
+            dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
+
+            // Calculate differentials between the vertices
+            // of the auxiliary triangle
+            dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
+            dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
+
+            // Calculate the gradient of xmom on the auxiliary triangle
+            a = dy2 * dq1 - dy1*dq2;
+            a *= inv_area2;
+            b = dx1 * dq2 - dx2*dq1;
+            b *= inv_area2;
+
+            // Calculate provisional jumps in stage from the centroid
+            // of triangle k to its vertices, to be limited
+            dqv[0] = a * dxv0 + b*dyv0;
+            dqv[1] = a * dxv1 + b*dyv1;
+            dqv[2] = a * dxv2 + b*dyv2;
+
+            // Now we want to find min and max of the centroid and the
+            // vertices of the auxiliary triangle and compute jumps
+            // from the centroid to the min and max
+            find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+
+            //beta_tmp = beta_vh;
+            //
+            //if (hmin<minimum_allowed_height)
+            //beta_tmp = beta_vh_dry;
+            beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;
+
+            // Limit the gradient
+            limit_gradient(dqv, qmin, qmax, beta_tmp);
+
+            for (i = 0; i < 3; i++) {
+                ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
+            }
+        }// End number_of_boundaries <=1
+        else {
+
+            //==============================================
+            // Number of boundaries == 2
+            //==============================================
+
+            // One internal neighbour and gradient is in direction of the neighbour's centroid
+
+            // Find the only internal neighbour (k1?)
+            for (k2 = k3; k2 < k3 + 3; k2++) {
+                // Find internal neighbour of triangle k
+                // k2 indexes the edges of triangle k
+
+                if (surrogate_neighbours[k2] != k) {
+                    break;
+                }
+            }
+
+            if ((k2 == k3 + 3)) {
+                // If we didn't find an internal neighbour
+                report_python_error(AT, "Internal neighbour not found");
+                return -1;
+            }
+
+            k1 = surrogate_neighbours[k2];
+
+            // The coordinates of the triangle are already (x,y).
+            // Get centroid of the neighbour (x1,y1)
+            coord_index = 2 * k1;
+            x1 = centroid_coordinates[coord_index];
+            y1 = centroid_coordinates[coord_index + 1];
+
+            // Compute x- and y- distances between the centroid of
+            // triangle k and that of its neighbour
+            dx1 = x1 - x;
+            dy1 = y1 - y;
+
+            // Set area2 as the square of the distance
+            area2 = dx1 * dx1 + dy1*dy1;
+
+            // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2)
+            // and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
+            // respectively correspond to the x- and y- gradients
+            // of the conserved quantities
+            dx2 = 1.0 / area2;
+            dy2 = dx2*dy1;
+            dx2 *= dx1;
+
+
+            //-----------------------------------
+            // stage
+            //-----------------------------------
+
+            // Compute differentials
+            dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
+
+            // Calculate the gradient between the centroid of triangle k
+            // and that of its neighbour
+            a = dq1*dx2;
+            b = dq1*dy2;
+
+            // Calculate provisional vertex jumps, to be limited
+            dqv[0] = a * dxv0 + b*dyv0;
+            dqv[1] = a * dxv1 + b*dyv1;
+            dqv[2] = a * dxv2 + b*dyv2;
+
+            // Now limit the jumps
+            if (dq1 >= 0.0) {
+                qmin = 0.0;
+                qmax = dq1;
+            } else {
+                qmin = dq1;
+                qmax = 0.0;
+            }
+
+            // Limit the gradient
+            limit_gradient(dqv, qmin, qmax, beta_w);
+
+            //for (i=0; i < 3; i++)
+            //{
+            stage_vertex_values[k3] = stage_centroid_values[k] + dqv[0];
+            stage_vertex_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
+            stage_vertex_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
+            //}
+
+            //-----------------------------------
+            // xmomentum
+            //-----------------------------------
+
+            // Compute differentials
+            dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
+
+            // Calculate the gradient between the centroid of triangle k
+            // and that of its neighbour
+            a = dq1*dx2;
+            b = dq1*dy2;
+
+            // Calculate provisional vertex jumps, to be limited
+            dqv[0] = a * dxv0 + b*dyv0;
+            dqv[1] = a * dxv1 + b*dyv1;
+            dqv[2] = a * dxv2 + b*dyv2;
+
+            // Now limit the jumps
+            if (dq1 >= 0.0) {
+                qmin = 0.0;
+                qmax = dq1;
+            } else {
+                qmin = dq1;
+                qmax = 0.0;
+            }
+
+            // Limit the gradient
+            limit_gradient(dqv, qmin, qmax, beta_w);
+
+            //for (i=0;i<3;i++)
+            //xmom_vertex_values[k3] = xmom_centroid_values[k] + dqv[0];
+            //xmom_vertex_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
+            //xmom_vertex_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
+
+            for (i = 0; i < 3; i++) {
+                xmom_vertex_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
+            }
+
+            //-----------------------------------
+            // ymomentum
+            //-----------------------------------
+
+            // Compute differentials
+            dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
+
+            // Calculate the gradient between the centroid of triangle k
+            // and that of its neighbour
+            a = dq1*dx2;
+            b = dq1*dy2;
+
+            // Calculate provisional vertex jumps, to be limited
+            dqv[0] = a * dxv0 + b*dyv0;
+            dqv[1] = a * dxv1 + b*dyv1;
+            dqv[2] = a * dxv2 + b*dyv2;
+
+            // Now limit the jumps
+            if (dq1 >= 0.0) {
+                qmin = 0.0;
+                qmax = dq1;
+            }
+            else {
+                qmin = dq1;
+                qmax = 0.0;
+            }
+
+            // Limit the gradient
+            limit_gradient(dqv, qmin, qmax, beta_w);
+
+            //for (i=0;i<3;i++)
+            //ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0];
+            //ymom_vertex_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
+            //ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
+
+            for (i = 0; i < 3; i++) {
+                ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
+            }
+            //ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0];
+            //ymom_vertex_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
+            //ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
+        } // else [number_of_boundaries==2]
+    } // for k=0 to number_of_elements-1
+
+    if (extrapolate_velocity_second_order == 1) {
+        // Convert back from velocity to momentum
+        for (k = 0; k < number_of_elements; k++) {
+            k3 = 3 * k;
+            //dv0 = max(stage_vertex_values[k3]-elevation_vertex_values[k3],minimum_allowed_height);
+            //dv1 = max(stage_vertex_values[k3+1]-elevation_vertex_values[k3+1],minimum_allowed_height);
+            //dv2 = max(stage_vertex_values[k3+2]-elevation_vertex_values[k3+2],minimum_allowed_height);
+            dv0 = max(stage_vertex_values[k3] - elevation_vertex_values[k3], 0.);
+            dv1 = max(stage_vertex_values[k3 + 1] - elevation_vertex_values[k3 + 1], 0.);
+            dv2 = max(stage_vertex_values[k3 + 2] - elevation_vertex_values[k3 + 2], 0.);
+
+            //Correct centroid and vertex values
+            xmom_centroid_values[k] = xmom_centroid_store[k];
+            xmom_vertex_values[k3] = xmom_vertex_values[k3] * dv0;
+            xmom_vertex_values[k3 + 1] = xmom_vertex_values[k3 + 1] * dv1;
+            xmom_vertex_values[k3 + 2] = xmom_vertex_values[k3 + 2] * dv2;
+
+            ymom_centroid_values[k] = ymom_centroid_store[k];
+            ymom_vertex_values[k3] = ymom_vertex_values[k3] * dv0;
+            ymom_vertex_values[k3 + 1] = ymom_vertex_values[k3 + 1] * dv1;
+            ymom_vertex_values[k3 + 2] = ymom_vertex_values[k3 + 2] * dv2;
+
+        }
+    }
+
+    return 0;
+}
+
+
+
+// Computational routine
+int _extrapolate_second_order_sw(struct domain *D) {
+
+
+  // Domain Variables
+    int number_of_elements;
+    double epsilon;
+    double minimum_allowed_height;
+    double beta_w;
+    double beta_w_dry;
+    double beta_uh;
+    double beta_uh_dry;
+    double beta_vh;
+    double beta_vh_dry;
+    long* surrogate_neighbours;
+    long* number_of_boundaries;
+    double* centroid_coordinates;
+    double* stage_centroid_values;
+    double* xmom_centroid_values;
+    double* ymom_centroid_values;
+    double* bed_centroid_values;
+    double* edge_coordinates;
+    double* vertex_coordinates;
+    double* stage_edge_values;
+    double* xmom_edge_values;
+    double* ymom_edge_values;
+    double* bed_edge_values;
+    double* stage_vertex_values;
+    double* xmom_vertex_values;
+    double* ymom_vertex_values;
+    double* bed_vertex_values;
+    int optimise_dry_cells;
+    int extrapolate_velocity_second_order;
+
+    // Local variables
+    double a, b; // Gradient vector used to calculate edge values from centroids
+    int k, k0, k1, k2, k3, k6, coord_index, i;
+    double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
+    double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
+    double dqv[3], qmin, qmax, hmin, hmax;
+    double hc, h0, h1, h2, beta_tmp, hfactor;
+    //double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale;
+    double dk, dv0, dv1, dv2;
+
+    double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store;
+
+
+    // Associate memory location of Domain varialbe with local aliases
+    number_of_elements     = D->number_of_elements;
+    epsilon                = D->epsilon;
+    minimum_allowed_height = D->minimum_allowed_height;
+    beta_w                 = D->beta_w;
+    beta_w_dry             = D->beta_w_dry;
+    beta_uh                = D->beta_uh;
+    beta_uh_dry            = D->beta_uh_dry;
+    beta_vh                = D->beta_vh;
+    beta_vh_dry            = D->beta_vh_dry;
+    optimise_dry_cells     = D->optimise_dry_cells;
+    extrapolate_velocity_second_order = D->extrapolate_velocity_second_order;
+
+    surrogate_neighbours      = D->surrogate_neighbours;
+    number_of_boundaries      = D->number_of_boundaries;
+    centroid_coordinates      = D->centroid_coordinates;
+    stage_centroid_values     = D->stage_centroid_values;
+    xmom_centroid_values      = D->xmom_centroid_values;
+    ymom_centroid_values      = D->ymom_centroid_values;
+    bed_centroid_values       = D->bed_centroid_values;
+    edge_coordinates          = D->edge_coordinates;
+    vertex_coordinates        = D->vertex_coordinates;
+    stage_edge_values         = D->stage_edge_values;
+    xmom_edge_values          = D->xmom_edge_values;
+    ymom_edge_values          = D->ymom_edge_values;
+    bed_edge_values           = D->bed_edge_values;
+    stage_vertex_values       = D->stage_vertex_values;
+    xmom_vertex_values        = D->xmom_vertex_values;
+    ymom_vertex_values        = D->ymom_vertex_values;
+    bed_vertex_values         = D->bed_vertex_values;
+
+
+
+
+/*
 int _extrapolate_second_order_sw(int number_of_elements,
         double epsilon,
@@ -2094 +2688 @@
     double hc, h0, h1, h2, beta_tmp, hfactor;
     double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements], dk, dv0, dv1, dv2;
+*/
+
+   // Use malloc to avoid putting these variables on the stack, which can cause
+   // segfaults in large model runs
+    xmom_centroid_store = malloc(number_of_elements*sizeof(double));
+    ymom_centroid_store = malloc(number_of_elements*sizeof(double));
+    stage_centroid_store = malloc(number_of_elements*sizeof(double));
 
     if (extrapolate_velocity_second_order == 1) {
@@ -2100 +2701 @@
         for (k = 0; k < number_of_elements; k++) {
 
-            dk = max(stage_centroid_values[k] - elevation_centroid_values[k], minimum_allowed_height);
+            dk = max(stage_centroid_values[k] - bed_centroid_values[k], minimum_allowed_height);
             xmom_centroid_store[k] = xmom_centroid_values[k];
             xmom_centroid_values[k] = xmom_centroid_values[k] / dk;
@@ -2219 +2820 @@
 
             // Calculate heights of neighbouring cells
-            hc = stage_centroid_values[k] - elevation_centroid_values[k];
-            h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
-            h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
-            h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
+            hc = stage_centroid_values[k] - bed_centroid_values[k];
+            h0 = stage_centroid_values[k0] - bed_centroid_values[k0];
+            h1 = stage_centroid_values[k1] - bed_centroid_values[k1];
+            h2 = stage_centroid_values[k2] - bed_centroid_values[k2];
             hmin = min(min(h0, min(h1, h2)), hc);
             //hfactor = hc/(hc + 1.0);
@@ -2546 +3147 @@
         for (k = 0; k < number_of_elements; k++) {
             k3 = 3 * k;
-            //dv0 = max(stage_vertex_values[k3]-elevation_vertex_values[k3],minimum_allowed_height);
-            //dv1 = max(stage_vertex_values[k3+1]-elevation_vertex_values[k3+1],minimum_allowed_height);
-            //dv2 = max(stage_vertex_values[k3+2]-elevation_vertex_values[k3+2],minimum_allowed_height);
-            dv0 = max(stage_vertex_values[k3] - elevation_vertex_values[k3], 0.);
-            dv1 = max(stage_vertex_values[k3 + 1] - elevation_vertex_values[k3 + 1], 0.);
-            dv2 = max(stage_vertex_values[k3 + 2] - elevation_vertex_values[k3 + 2], 0.);
+            //dv0 = max(stage_vertex_values[k3]-bed_vertex_values[k3],minimum_allowed_height);
+            //dv1 = max(stage_vertex_values[k3+1]-bed_vertex_values[k3+1],minimum_allowed_height);
+            //dv2 = max(stage_vertex_values[k3+2]-bed_vertex_values[k3+2],minimum_allowed_height);
+            dv0 = max(stage_vertex_values[k3] - bed_vertex_values[k3], 0.);
+            dv1 = max(stage_vertex_values[k3 + 1] - bed_vertex_values[k3 + 1], 0.);
+            dv2 = max(stage_vertex_values[k3 + 2] - bed_vertex_values[k3 + 2], 0.);
 
             //Correct centroid and vertex values
@@ -2566 +3167 @@
         }
     }
+    
+    free(xmom_centroid_store);
+    free(ymom_centroid_store);
+    free(stage_centroid_store);
+
 
     return 0;
 }
 
-PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
+
+// Computational routine
+int _extrapolate_second_order_edge_sw(struct domain *D) {
+
+
+  // Domain Variables
+    int number_of_elements;
+    double epsilon;
+    double minimum_allowed_height;
+    double beta_w;
+    double beta_w_dry;
+    double beta_uh;
+    double beta_uh_dry;
+    double beta_vh;
+    double beta_vh_dry;
+    long* surrogate_neighbours;
+    long* number_of_boundaries;
+    double* centroid_coordinates;
+    double* stage_centroid_values;
+    double* xmom_centroid_values;
+    double* ymom_centroid_values;
+    double* bed_centroid_values;
+    double* edge_coordinates;
+    double* stage_edge_values;
+    double* xmom_edge_values;
+    double* ymom_edge_values;
+    double* bed_edge_values;
+    double* stage_vertex_values;
+    double* xmom_vertex_values;
+    double* ymom_vertex_values;
+    double* bed_vertex_values;
+    int optimise_dry_cells;
+    int extrapolate_velocity_second_order;
+
+    // Local variables
+    double a, b; // Gradient vector used to calculate edge values from centroids
+    int k, k0, k1, k2, k3, k6, coord_index, i;
+    double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
+    double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
+    double dqv[3], qmin, qmax, hmin, hmax;
+    double hc, h0, h1, h2, beta_tmp, hfactor;
+    //double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale;
+    double dk, de[3];
+
+    double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store;
+
+
+    // Associate memory location of Domain varialbe with local aliases
+    number_of_elements     = D->number_of_elements;
+    epsilon                = D->epsilon;
+    minimum_allowed_height = D->minimum_allowed_height;
+    beta_w                 = D->beta_w;
+    beta_w_dry             = D->beta_w_dry;
+    beta_uh                = D->beta_uh;
+    beta_uh_dry            = D->beta_uh_dry;
+    beta_vh                = D->beta_vh;
+    beta_vh_dry            = D->beta_vh_dry;
+    optimise_dry_cells     = D->optimise_dry_cells;
+    extrapolate_velocity_second_order = D->extrapolate_velocity_second_order;
+
+    surrogate_neighbours      = D->surrogate_neighbours;
+    number_of_boundaries      = D->number_of_boundaries;
+    centroid_coordinates      = D->centroid_coordinates;
+    stage_centroid_values     = D->stage_centroid_values;
+    xmom_centroid_values      = D->xmom_centroid_values;
+    ymom_centroid_values      = D->ymom_centroid_values;
+    bed_centroid_values       = D->bed_centroid_values;
+    edge_coordinates          = D->edge_coordinates;
+    stage_edge_values         = D->stage_edge_values;
+    xmom_edge_values          = D->xmom_edge_values;
+    ymom_edge_values          = D->ymom_edge_values;
+    bed_edge_values           = D->bed_edge_values;
+    stage_vertex_values       = D->stage_vertex_values;
+    xmom_vertex_values        = D->xmom_vertex_values;
+    ymom_vertex_values        = D->ymom_vertex_values;
+    bed_vertex_values         = D->bed_vertex_values;
+
+
+
+  // Use malloc to avoid putting these variables on the stack, which can cause
+  // segfaults in large model runs
+  xmom_centroid_store = malloc(number_of_elements*sizeof(double));
+  ymom_centroid_store = malloc(number_of_elements*sizeof(double));
+  stage_centroid_store = malloc(number_of_elements*sizeof(double));
+
+  if(extrapolate_velocity_second_order==1){
+      // Replace momentum centroid with velocity centroid to allow velocity
+      // extrapolation This will be changed back at the end of the routine
+      for (k=0; k<number_of_elements; k++){
+
+          dk = max(stage_centroid_values[k]-bed_centroid_values[k],minimum_allowed_height);
+          xmom_centroid_store[k] = xmom_centroid_values[k];
+          xmom_centroid_values[k] = xmom_centroid_values[k]/dk;
+
+          ymom_centroid_store[k] = ymom_centroid_values[k];
+          ymom_centroid_values[k] = ymom_centroid_values[k]/dk;
+
+      }
+  }
+
+  // Begin extrapolation routine
+  for (k = 0; k < number_of_elements; k++)
+  {
+    k3=k*3;
+    k6=k*6;
+
+    if (number_of_boundaries[k]==3)
+    //if (0==0)
+    {
+      // No neighbours, set gradient on the triangle to zero
+
+      stage_edge_values[k3]   = stage_centroid_values[k];
+      stage_edge_values[k3+1] = stage_centroid_values[k];
+      stage_edge_values[k3+2] = stage_centroid_values[k];
+      xmom_edge_values[k3]    = xmom_centroid_values[k];
+      xmom_edge_values[k3+1]  = xmom_centroid_values[k];
+      xmom_edge_values[k3+2]  = xmom_centroid_values[k];
+      ymom_edge_values[k3]    = ymom_centroid_values[k];
+      ymom_edge_values[k3+1]  = ymom_centroid_values[k];
+      ymom_edge_values[k3+2]  = ymom_centroid_values[k];
+
+      continue;
+    }
+    else
+    {
+      // Triangle k has one or more neighbours.
+      // Get centroid and edge coordinates of the triangle
+
+      // Get the edge coordinates
+      xv0 = edge_coordinates[k6];
+      yv0 = edge_coordinates[k6+1];
+      xv1 = edge_coordinates[k6+2];
+      yv1 = edge_coordinates[k6+3];
+      xv2 = edge_coordinates[k6+4];
+      yv2 = edge_coordinates[k6+5];
+
+      // Get the centroid coordinates
+      coord_index = 2*k;
+      x = centroid_coordinates[coord_index];
+      y = centroid_coordinates[coord_index+1];
+
+      // Store x- and y- differentials for the edges of
+      // triangle k relative to the centroid
+      dxv0 = xv0 - x;
+      dxv1 = xv1 - x;
+      dxv2 = xv2 - x;
+      dyv0 = yv0 - y;
+      dyv1 = yv1 - y;
+      dyv2 = yv2 - y;
+      // Compute the minimum distance from the centroid to an edge
+      //demin=min(dxv0*dxv0 +dyv0*dyv0, min(dxv1*dxv1+dyv1*dyv1, dxv2*dxv2+dyv2*dyv2));
+      //demin=sqrt(demin);
+    }
+
+
+
+    if (number_of_boundaries[k]<=1)
+    {
+      //==============================================
+      // Number of boundaries <= 1
+      //==============================================
+
+
+      // If no boundaries, auxiliary triangle is formed
+      // from the centroids of the three neighbours
+      // If one boundary, auxiliary triangle is formed
+      // from this centroid and its two neighbours
+
+      k0 = surrogate_neighbours[k3];
+      k1 = surrogate_neighbours[k3 + 1];
+      k2 = surrogate_neighbours[k3 + 2];
+
+      // Test to see whether we accept the surrogate neighbours
+      // Note that if ki is replaced with k in more than 1 neighbour, then the
+      // triangle area will be zero, and a first order extrapolation will be
+      // used
+      if(stage_centroid_values[k2]<=bed_centroid_values[k2]){
+          k2 = k ;
+      }
+      if(stage_centroid_values[k0]<=bed_centroid_values[k0]){
+          k0 = k ;
+      }
+      if(stage_centroid_values[k1]<=bed_centroid_values[k1]){
+          k1 = k ;
+      }
+      // Alternative approach (BRUTAL) -- if the max neighbour bed elevation is greater
+      // than the min neighbour stage, then we use first order extrapolation
+      //bedmax = max(bed_centroid_values[k],
+      //             max(bed_centroid_values[k0],
+      //                 max(bed_centroid_values[k1], bed_centroid_values[k2])));
+      //stagemin = min(stage_centroid_values[k],
+      //             min(stage_centroid_values[k0],
+      //                 min(stage_centroid_values[k1], stage_centroid_values[k2])));
+      //
+      //if(stagemin < bedmax){
+      //   // This will cause first order extrapolation
+      //   k2 = k;
+      //   k0 = k;
+      //   k1 = k;
+      //}
+
+      // Get the auxiliary triangle's vertex coordinates
+      // (really the centroids of neighbouring triangles)
+      coord_index = 2*k0;
+      x0 = centroid_coordinates[coord_index];
+      y0 = centroid_coordinates[coord_index+1];
+
+      coord_index = 2*k1;
+      x1 = centroid_coordinates[coord_index];
+      y1 = centroid_coordinates[coord_index+1];
+
+      coord_index = 2*k2;
+      x2 = centroid_coordinates[coord_index];
+      y2 = centroid_coordinates[coord_index+1];
+
+      // compute the maximum distance from the centroid to a neighbouring
+      // centroid
+      //dcmax=max( (x0-x)*(x0-x) + (y0-y)*(y0-y),
+      //           max((x1-x)*(x1-x) + (y1-y)*(y1-y),
+      //               (x2-x)*(x2-x) + (y2-y)*(y2-y)));
+      //dcmax=sqrt(dcmax);
+      //// Ratio of centroid to edge distance -- useful in attempting to adapt limiter
+      //if(dcmax>0.){
+      //    r0scale=demin/dcmax;
+      //    //printf("%f \n", r0scale);
+      //}else{
+      //    r0scale=0.5;
+      //}
+
+      // Store x- and y- differentials for the vertices
+      // of the auxiliary triangle
+      dx1 = x1 - x0;
+      dx2 = x2 - x0;
+      dy1 = y1 - y0;
+      dy2 = y2 - y0;
+
+      // Calculate 2*area of the auxiliary triangle
+      // The triangle is guaranteed to be counter-clockwise
+      area2 = dy2*dx1 - dy1*dx2;
+
+      // If the mesh is 'weird' near the boundary,
+      // the triangle might be flat or clockwise
+      // Default to zero gradient
+      if (area2 <= 0)
+      {
+          //printf("Error negative triangle area \n");
+          //report_python_error(AT, "Negative triangle area");
+          //return -1;
+
+          stage_edge_values[k3]   = stage_centroid_values[k];
+          stage_edge_values[k3+1] = stage_centroid_values[k];
+          stage_edge_values[k3+2] = stage_centroid_values[k];
+          xmom_edge_values[k3]    = xmom_centroid_values[k];
+          xmom_edge_values[k3+1]  = xmom_centroid_values[k];
+          xmom_edge_values[k3+2]  = xmom_centroid_values[k];
+          ymom_edge_values[k3]    = ymom_centroid_values[k];
+          ymom_edge_values[k3+1]  = ymom_centroid_values[k];
+          ymom_edge_values[k3+2]  = ymom_centroid_values[k];
+
+          continue;
+      }
+
+      // Calculate heights of neighbouring cells
+      hc = stage_centroid_values[k]  - bed_centroid_values[k];
+      h0 = stage_centroid_values[k0] - bed_centroid_values[k0];
+      h1 = stage_centroid_values[k1] - bed_centroid_values[k1];
+      h2 = stage_centroid_values[k2] - bed_centroid_values[k2];
+      hmin = min(min(h0, min(h1, h2)), hc);
+      //hmin = min(h0, min(h1, h2));
+      //hmin = max(hmin, 0.0);
+      //hfactor = hc/(hc + 1.0);
+
+      hfactor = 0.0;
+      //if (hmin > 0.001)
+      if (hmin > 0.)
+      //if (hc>0.0)
+      {
+        hfactor = 1.0 ;//hmin/(hmin + 0.004);
+        //hfactor=hmin/(hmin + 0.004);
+      }
+
+      if (optimise_dry_cells)
+      {
+        // Check if linear reconstruction is necessary for triangle k
+        // This check will exclude dry cells.
+
+        //hmax = max(h0, max(h1, max(h2, hc)));
+        hmax = max(h0, max(h1, h2));
+        if (hmax < epsilon)
+        {
+            continue;
+        }
+      }
+
+      //-----------------------------------
+      // stage
+      //-----------------------------------
+
+      // Calculate the difference between vertex 0 of the auxiliary
+      // triangle and the centroid of triangle k
+      dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
+
+      // Calculate differentials between the vertices
+      // of the auxiliary triangle (centroids of neighbouring triangles)
+      dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
+      dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
+
+      inv_area2 = 1.0/area2;
+      // Calculate the gradient of stage on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a *= inv_area2;
+      b = dx1*dq2 - dx2*dq1;
+      b *= inv_area2;
+
+      // Calculate provisional jumps in stage from the centroid
+      // of triangle k to its vertices, to be limited
+      dqv[0] = a*dxv0 + b*dyv0;
+      dqv[1] = a*dxv1 + b*dyv1;
+      dqv[2] = a*dxv2 + b*dyv2;
+
+      // Now we want to find min and max of the centroid and the
+      // vertices of the auxiliary triangle and compute jumps
+      // from the centroid to the min and max
+      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+
+      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
+
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_tmp);
+      //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
+
+      //for (i=0;i<3;i++)
+      stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
+      stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1];
+      stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
+
+      //-----------------------------------
+      // xmomentum
+      //-----------------------------------
+
+      // Calculate the difference between vertex 0 of the auxiliary
+      // triangle and the centroid of triangle k
+      dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
+
+      // Calculate differentials between the vertices
+      // of the auxiliary triangle
+      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
+      dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
+
+      // Calculate the gradient of xmom on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a *= inv_area2;
+      b = dx1*dq2 - dx2*dq1;
+      b *= inv_area2;
+
+      // Calculate provisional jumps in stage from the centroid
+      // of triangle k to its vertices, to be limited
+      dqv[0] = a*dxv0+b*dyv0;
+      dqv[1] = a*dxv1+b*dyv1;
+      dqv[2] = a*dxv2+b*dyv2;
+
+      // Now we want to find min and max of the centroid and the
+      // vertices of the auxiliary triangle and compute jumps
+      // from the centroid to the min and max
+      //
+      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+
+      beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
+
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_tmp);
+      //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
+
+
+      for (i=0; i < 3; i++)
+      {
+        xmom_edge_values[k3+i] = xmom_centroid_values[k] + dqv[i];
+      }
+
+      //-----------------------------------
+      // ymomentum
+      //-----------------------------------
+
+      // Calculate the difference between vertex 0 of the auxiliary
+      // triangle and the centroid of triangle k
+      dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
+
+      // Calculate differentials between the vertices
+      // of the auxiliary triangle
+      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
+      dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
+
+      // Calculate the gradient of xmom on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a *= inv_area2;
+      b = dx1*dq2 - dx2*dq1;
+      b *= inv_area2;
+
+      // Calculate provisional jumps in stage from the centroid
+      // of triangle k to its vertices, to be limited
+      dqv[0] = a*dxv0 + b*dyv0;
+      dqv[1] = a*dxv1 + b*dyv1;
+      dqv[2] = a*dxv2 + b*dyv2;
+
+      // Now we want to find min and max of the centroid and the
+      // vertices of the auxiliary triangle and compute jumps
+      // from the centroid to the min and max
+      //
+      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+
+      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;
+
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_tmp);
+      //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
+
+      for (i=0;i<3;i++)
+      {
+        ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
+      }
+
+    } // End number_of_boundaries <=1
+    else
+    {
+
+      //==============================================
+      // Number of boundaries == 2
+      //==============================================
+
+      // One internal neighbour and gradient is in direction of the neighbour's centroid
+
+      // Find the only internal neighbour (k1?)
+      for (k2 = k3; k2 < k3 + 3; k2++)
+      {
+      // Find internal neighbour of triangle k
+      // k2 indexes the edges of triangle k
+
+          if (surrogate_neighbours[k2] != k)
+          {
+             break;
+          }
+      }
+
+      if ((k2 == k3 + 3))
+      {
+        // If we didn't find an internal neighbour
+        report_python_error(AT, "Internal neighbour not found");
+        return -1;
+      }
+
+      k1 = surrogate_neighbours[k2];
+
+      // The coordinates of the triangle are already (x,y).
+      // Get centroid of the neighbour (x1,y1)
+      coord_index = 2*k1;
+      x1 = centroid_coordinates[coord_index];
+      y1 = centroid_coordinates[coord_index + 1];
+
+      // Compute x- and y- distances between the centroid of
+      // triangle k and that of its neighbour
+      dx1 = x1 - x;
+      dy1 = y1 - y;
+
+      // Set area2 as the square of the distance
+      area2 = dx1*dx1 + dy1*dy1;
+
+      // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2)
+      // and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
+      // respectively correspond to the x- and y- gradients
+      // of the conserved quantities
+      dx2 = 1.0/area2;
+      dy2 = dx2*dy1;
+      dx2 *= dx1;
+
+
+      //-----------------------------------
+      // stage
+      //-----------------------------------
+
+      // Compute differentials
+      dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
+
+      // Calculate the gradient between the centroid of triangle k
+      // and that of its neighbour
+      a = dq1*dx2;
+      b = dq1*dy2;
+
+      // Calculate provisional edge jumps, to be limited
+      dqv[0] = a*dxv0 + b*dyv0;
+      dqv[1] = a*dxv1 + b*dyv1;
+      dqv[2] = a*dxv2 + b*dyv2;
+
+      // Now limit the jumps
+      if (dq1>=0.0)
+      {
+        qmin=0.0;
+        qmax=dq1;
+      }
+      else
+      {
+        qmin = dq1;
+        qmax = 0.0;
+      }
+
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_w);
+
+      //for (i=0; i < 3; i++)
+      //{
+      stage_edge_values[k3] = stage_centroid_values[k] + dqv[0];
+      stage_edge_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
+      stage_edge_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
+      //}
+
+      //-----------------------------------
+      // xmomentum
+      //-----------------------------------
+
+      // Compute differentials
+      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
+
+      // Calculate the gradient between the centroid of triangle k
+      // and that of its neighbour
+      a = dq1*dx2;
+      b = dq1*dy2;
+
+      // Calculate provisional edge jumps, to be limited
+      dqv[0] = a*dxv0+b*dyv0;
+      dqv[1] = a*dxv1+b*dyv1;
+      dqv[2] = a*dxv2+b*dyv2;
+
+      // Now limit the jumps
+      if (dq1 >= 0.0)
+      {
+        qmin = 0.0;
+        qmax = dq1;
+      }
+      else
+      {
+        qmin = dq1;
+        qmax = 0.0;
+      }
+
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_w);
+
+      //for (i=0;i<3;i++)
+      //xmom_edge_values[k3] = xmom_centroid_values[k] + dqv[0];
+      //xmom_edge_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
+      //xmom_edge_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
+
+      for (i = 0; i < 3;i++)
+      {
+          xmom_edge_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
+      }
+
+      //-----------------------------------
+      // ymomentum
+      //-----------------------------------
+
+      // Compute differentials
+      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
+
+      // Calculate the gradient between the centroid of triangle k
+      // and that of its neighbour
+      a = dq1*dx2;
+      b = dq1*dy2;
+
+      // Calculate provisional edge jumps, to be limited
+      dqv[0] = a*dxv0 + b*dyv0;
+      dqv[1] = a*dxv1 + b*dyv1;
+      dqv[2] = a*dxv2 + b*dyv2;
+
+      // Now limit the jumps
+      if (dq1>=0.0)
+      {
+        qmin = 0.0;
+        qmax = dq1;
+      }
+      else
+      {
+        qmin = dq1;
+        qmax = 0.0;
+      }
+
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_w);
+
+      for (i=0;i<3;i++)
+              {
+              ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
+              }
+    } // else [number_of_boundaries==2]
+  } // for k=0 to number_of_elements-1
+
+  // Compute vertex values of quantities
+  for (k=0; k<number_of_elements; k++){
+      k3=3*k;
+
+      // Compute stage vertex values
+      stage_vertex_values[k3] = stage_edge_values[k3+1] + stage_edge_values[k3+2] -stage_edge_values[k3] ;
+      stage_vertex_values[k3+1] =  stage_edge_values[k3] + stage_edge_values[k3+2]-stage_edge_values[k3+1];
+      stage_vertex_values[k3+2] =  stage_edge_values[k3] + stage_edge_values[k3+1]-stage_edge_values[k3+2];
+
+     // Compute xmom vertex values
+      xmom_vertex_values[k3] = xmom_edge_values[k3+1] + xmom_edge_values[k3+2] -xmom_edge_values[k3] ;
+      xmom_vertex_values[k3+1] =  xmom_edge_values[k3] + xmom_edge_values[k3+2]-xmom_edge_values[k3+1];
+      xmom_vertex_values[k3+2] =  xmom_edge_values[k3] + xmom_edge_values[k3+1]-xmom_edge_values[k3+2];
+
+      // Compute ymom vertex values
+      ymom_vertex_values[k3] = ymom_edge_values[k3+1] + ymom_edge_values[k3+2] -ymom_edge_values[k3] ;
+      ymom_vertex_values[k3+1] =  ymom_edge_values[k3] + ymom_edge_values[k3+2]-ymom_edge_values[k3+1];
+      ymom_vertex_values[k3+2] =  ymom_edge_values[k3] + ymom_edge_values[k3+1]-ymom_edge_values[k3+2];
+
+      // If needed, convert from velocity to momenta
+      if(extrapolate_velocity_second_order==1){
+          //Convert velocity back to momenta at centroids
+          xmom_centroid_values[k] = xmom_centroid_store[k];
+          ymom_centroid_values[k] = ymom_centroid_store[k];
+
+          // Re-compute momenta at edges
+          for (i=0; i<3; i++){
+              de[i] = max(stage_edge_values[k3+i]-bed_edge_values[k3+i],0.0);
+              xmom_edge_values[k3+i]=xmom_edge_values[k3+i]*de[i];
+              ymom_edge_values[k3+i]=ymom_edge_values[k3+i]*de[i];
+          }
+
+          // Re-compute momenta at vertices
+          for (i=0; i<3; i++){
+              de[i] = max(stage_vertex_values[k3+i]-bed_vertex_values[k3+i],0.0);
+              xmom_vertex_values[k3+i]=xmom_vertex_values[k3+i]*de[i];
+              ymom_vertex_values[k3+i]=ymom_vertex_values[k3+i]*de[i];
+          }
+
+      }
+
+
+  }
+
+  free(xmom_centroid_store);
+  free(ymom_centroid_store);
+  free(stage_centroid_store);
+
+  return 0;
+}
+
+
+
+
+
+
+
+PyObject *extrapolate_second_order_sw_old(PyObject *self, PyObject *args) {
     /*Compute the vertex values based on a linear reconstruction
       on each triangle
@@ -2681 +3939 @@
 
     // Call underlying computational routine
-    e = _extrapolate_second_order_sw(number_of_elements,
+    e = _extrapolate_second_order_sw_old(number_of_elements,
             epsilon,
             minimum_allowed_height,
@@ -2713 +3971 @@
     return Py_BuildValue("");
 
+}// extrapolate_second-order_sw_old
+
+PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
+  /*Compute the edge values based on a linear reconstruction
+    on each triangle
+
+    Python call:
+    extrapolate_second_order_edge_sw(domain)
+
+    Post conditions:
+        The edges of each triangle have values from a
+        limited linear reconstruction
+        based on centroid values
+
+  */
+
+    struct domain D;
+    PyObject *domain;
+    int err;
+
+    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);
+
+    print_domain_struct(&D);
+
+    // Call underlying flux computation routine and update
+    // the explicit update arrays
+    err = _extrapolate_second_order_sw(&D);
+
+    if (err == -1) {
+        // Use error string set inside computational routine
+        return NULL;
+    }
+
+  return Py_BuildValue("");
+
 }// extrapolate_second-order_sw
+
+PyObject *extrapolate_second_order_edge_sw(PyObject *self, PyObject *args) {
+  /*Compute the edge values based on a linear reconstruction
+    on each triangle
+
+    Python call:
+    extrapolate_second_order_edge_sw(domain)
+
+    Post conditions:
+        The edges of each triangle have values from a
+        limited linear reconstruction
+        based on centroid values
+
+  */
+    
+    struct domain D;
+    PyObject *domain;
+    int err;
+
+    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
+    err = _extrapolate_second_order_edge_sw(&D);
+
+    if (err == -1) {
+        // Use error string set inside computational routine
+        return NULL;
+    }
+
+  return Py_BuildValue("");
+
+}// extrapolate_second-order_edge_sw
+
+
 
 PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
@@ -4416 +5758 @@
     {"rotate", (PyCFunction) rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
     {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
+    {"extrapolate_second_order_sw_old", extrapolate_second_order_sw_old, METH_VARARGS, "Print out"},
+    {"extrapolate_second_order_edge_sw", extrapolate_second_order_edge_sw, METH_VARARGS, "Print out"},
     {"compute_fluxes_ext_wb", compute_fluxes_ext_wb, METH_VARARGS, "Print out"},
     {"compute_fluxes_ext_wb_3", compute_fluxes_ext_wb_3, METH_VARARGS, "Print out"},

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

@@ -18 +18 @@
     long    optimise_dry_cells;
     double  evolve_max_timestep;
-
-    // The values in the python object will be changed
+    long    extrapolate_velocity_second_order;
+    double  minimum_allowed_height;
+
+    double beta_w;
+    double beta_w_dry;
+    double beta_uh;
+    double beta_uh_dry;
+    double beta_vh;
+    double beta_vh_dry;
+
+
+    // Changing values in these arrays will change the values in the python object
     long*   neighbours;
     long*   neighbour_edges;
+    long*   surrogate_neighbours;
     double* normals;
     double* edgelengths;
@@ -27 +38 @@
     double* areas;
 
+
+
     long*   tri_full_flag;
     long*   already_computed_flux;
@@ -32 +45 @@
 
     double* vertex_coordinates;
-
+    double* edge_coordinates;
+    double* centroid_coordinates;
+
+    long*   number_of_boundaries;
     double* stage_edge_values;
     double* xmom_edge_values;
@@ -139 +155 @@
             *already_computed_flux,
             *vertex_coordinates,
+            *edge_coordinates,
+            *centroid_coordinates,
+            *number_of_boundaries,
+            *surrogate_neighbours,
             *max_speed;
 
     PyObject *quantities;
 
-    D->number_of_elements = get_python_integer(domain, "number_of_elements");
-    D->epsilon = get_python_double(domain, "epsilon");
-    D->H0 = get_python_double(domain, "H0");
-    D->g = get_python_double(domain, "g");
-    D->optimise_dry_cells = get_python_integer(domain, "optimise_dry_cells");
-    D->evolve_max_timestep = get_python_double(domain, "evolve_max_timestep");
-
+    D->number_of_elements   = get_python_integer(domain, "number_of_elements");
+    D->epsilon              = get_python_double(domain, "epsilon");
+    D->H0                   = get_python_double(domain, "H0");
+    D->g                    = get_python_double(domain, "g");
+    D->optimise_dry_cells   = get_python_integer(domain, "optimise_dry_cells");
+    D->evolve_max_timestep  = get_python_double(domain, "evolve_max_timestep");
+    D->minimum_allowed_height = get_python_double(domain, "minimum_allowed_height");
+
+    D->extrapolate_velocity_second_order  = get_python_integer(domain, "extrapolate_velocity_second_order");
+
+    D->beta_w      = get_python_double(domain, "beta_w");;
+    D->beta_w_dry  = get_python_double(domain, "beta_w_dry");
+    D->beta_uh     = get_python_double(domain, "beta_uh");
+    D->beta_uh_dry = get_python_double(domain, "beta_uh_dry");
+    D->beta_vh     = get_python_double(domain, "beta_vh");
+    D->beta_vh_dry = get_python_double(domain, "beta_vh_dry");
+
+
+    
     neighbours = get_consecutive_array(domain, "neighbours");
     D->neighbours = (long *) neighbours->data;
 
+    surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
+    D->surrogate_neighbours = (long *) surrogate_neighbours->data;
+
     neighbour_edges = get_consecutive_array(domain, "neighbour_edges");
     D->neighbour_edges = (long *) neighbour_edges->data;
@@ -176 +211 @@
     vertex_coordinates = get_consecutive_array(domain, "vertex_coordinates");
     D->vertex_coordinates = (double *) vertex_coordinates->data;
+
+    edge_coordinates = get_consecutive_array(domain, "edge_coordinates");
+    D->edge_coordinates = (double *) edge_coordinates->data;
+
+
+    centroid_coordinates = get_consecutive_array(domain, "centroid_coordinates");
+    D->edge_coordinates = (double *) centroid_coordinates->data;
     
     max_speed = get_consecutive_array(domain, "max_speed");
     D->max_speed = (double *) max_speed->data;
+
+    number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
+    D->number_of_boundaries = (long *) number_of_boundaries->data;
+
 
 
@@ -210 +256 @@
     return D;
 }
+
+
+
+int print_domain_struct(struct domain *D) {
+
+
+    printf("D->number_of_elements %ld  \n", D->number_of_elements);
+    printf("D->epsilon %g \n", D->epsilon);
+    printf("D->H0 %g \n", D->H0);
+    printf("D->g %g \n", D->g);
+    printf("D->optimise_dry_cells %ld \n", D->optimise_dry_cells);
+    printf("D->evolve_max_timestep %g \n", D->evolve_max_timestep);
+    printf("D->minimum_allowed_height %g \n", D->minimum_allowed_height);
+    printf("D->extrapolate_velocity_second_order %ld \n", D->extrapolate_velocity_second_order);
+    printf("D->beta_w %g \n", D->beta_w);
+    printf("D->beta_w_dry %g \n", D->beta_w_dry);
+    printf("D->beta_uh %g \n", D->beta_uh);
+    printf("D->beta_uh_dry %g \n", D->beta_uh_dry);
+    printf("D->beta_vh %g \n", D->beta_vh);
+    printf("D->beta_vh_dry %g \n", D->beta_vh_dry);
+
+
+
+    printf("D->neighbours            %p \n", D->neighbours);
+    printf("D->surrogate_neighbours  %p \n", D->surrogate_neighbours);
+    printf("D->neighbour_edges       %p \n", D->neighbour_edges);
+    printf("D->normals               %p \n", D->normals);
+    printf("D->edgelengths           %p \n", D->edgelengths);
+    printf("D->radii                 %p \n", D->radii);
+    printf("D->areas                 %p \n", D->areas);
+    printf("D->tri_full_flag         %p \n", D->tri_full_flag);
+    printf("D->already_computed_flux %p \n", D->already_computed_flux);
+    printf("D->vertex_coordinates    %p \n", D->vertex_coordinates);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+    printf("D->neighbours %p \n", D->neighbours);
+
+
+
+//
+//    edge_coordinates = get_consecutive_array(domain, "edge_coordinates");
+//    D->edge_coordinates = (double *) edge_coordinates->data;
+//
+//
+//    centroid_coordinates = get_consecutive_array(domain, "centroid_coordinates");
+//    D->edge_coordinates = (double *) centroid_coordinates->data;
+//
+//    max_speed = get_consecutive_array(domain, "max_speed");
+//    D->max_speed = (double *) max_speed->data;
+//
+//    number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
+//    D->number_of_boundaries = (long *) number_of_boundaries->data;
+//
+//
+//
+//    quantities = get_python_object(domain, "quantities");
+//
+//    D->stage_edge_values     = get_python_array_data_from_dict(quantities, "stage",     "edge_values");
+//    D->xmom_edge_values      = get_python_array_data_from_dict(quantities, "xmomentum", "edge_values");
+//    D->ymom_edge_values      = get_python_array_data_from_dict(quantities, "ymomentum", "edge_values");
+//    D->bed_edge_values       = get_python_array_data_from_dict(quantities, "elevation", "edge_values");
+//
+//    D->stage_centroid_values     = get_python_array_data_from_dict(quantities, "stage",     "centroid_values");
+//    D->xmom_centroid_values      = get_python_array_data_from_dict(quantities, "xmomentum", "centroid_values");
+//    D->ymom_centroid_values      = get_python_array_data_from_dict(quantities, "ymomentum", "centroid_values");
+//    D->bed_centroid_values       = get_python_array_data_from_dict(quantities, "elevation", "centroid_values");
+//
+//    D->stage_vertex_values     = get_python_array_data_from_dict(quantities, "stage",     "vertex_values");
+//    D->xmom_vertex_values      = get_python_array_data_from_dict(quantities, "xmomentum", "vertex_values");
+//    D->ymom_vertex_values      = get_python_array_data_from_dict(quantities, "ymomentum", "vertex_values");
+//    D->bed_vertex_values       = get_python_array_data_from_dict(quantities, "elevation", "vertex_values");
+//
+//    D->stage_boundary_values = get_python_array_data_from_dict(quantities, "stage",     "boundary_values");
+//    D->xmom_boundary_values  = get_python_array_data_from_dict(quantities, "xmomentum", "boundary_values");
+//    D->ymom_boundary_values  = get_python_array_data_from_dict(quantities, "ymomentum", "boundary_values");
+//    D->bed_boundary_values   = get_python_array_data_from_dict(quantities, "elevation", "boundary_values");
+//
+//    D->stage_explicit_update = get_python_array_data_from_dict(quantities, "stage",     "explicit_update");
+//    D->xmom_explicit_update  = get_python_array_data_from_dict(quantities, "xmomentum", "explicit_update");
+//    D->ymom_explicit_update  = get_python_array_data_from_dict(quantities, "ymomentum", "explicit_update");
+//
+
+    return 0;
+}
+

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

@@ -3836 +3836 @@
         assert num.allclose(Y[1,:], Y_EX)
 
+
+
+
+    def Xtest_extrapolate_second_order_sw_real_data(self):
+        #Using test data generated by abstract_2d_finite_volumes-2
+        #Assuming no friction and flat bed (0.0)
+
+        a = [0.0, 0.0]
+        b = [0.0, 1.0/5]
+        c = [0.0, 2.0/5]
+        d = [1.0/5, 0.0]
+        e = [1.0/5, 1.0/5]
+        f = [1.0/5, 2.0/5]
+        g = [2.0/5, 2.0/5]
+
+        points = [a, b, c, d, e, f, g]
+        #             bae,     efb,     cbf,     feg
+        vertices = [[1,0,4], [4,5,1], [2,1,5], [5,4,6]]
+
+        domain = Domain(points, vertices)
+
+        def slope(x, y):
+            return -x/3
+
+        domain.set_quantity('elevation', slope)
+        domain.set_quantity('stage',
+                            [0.01298164, 0.00365611,
+                             0.01440365, -0.0381856437096],
+                            location='centroids')
+        domain.set_quantity('xmomentum',
+                            [0.00670439, 0.01263789,
+                             0.00647805, 0.0178180740668],
+                            location='centroids')
+        domain.set_quantity('ymomentum',
+                            [-7.23510980e-004, -6.30413883e-005,
+                             6.30413883e-005, 0.000200907255866],
+                            location='centroids')
+
+        E = domain.quantities['elevation'].vertex_values
+        L = domain.quantities['stage'].vertex_values
+        X = domain.quantities['xmomentum'].vertex_values
+        Y = domain.quantities['ymomentum'].vertex_values
+
+        domain._order_ = 2
+        domain.beta_w = 0.9
+        domain.beta_w_dry = 0.9
+        domain.beta_uh = 0.9
+        domain.beta_uh_dry = 0.9
+        domain.beta_vh = 0.9
+        domain.beta_vh_dry = 0.9
+
+        # FIXME (Ole): Need tests where this is commented out
+        domain.tight_slope_limiters = 0       # Backwards compatibility (14/4/7)
+        domain.use_centroid_velocities = 0    # Backwards compatibility (7/5/8)
+        domain.optimised_gradient_limiter = True
+
+        #domain.extrapolate_second_order_sw()
+
+
+        from  anuga.shallow_water.shallow_water_domain import extrapolate_second_order_sw
+
+        extrapolate_second_order_sw(domain)
+
+        #domain.distribute_to_vertices_and_edges()
+
+
+        print L[1,:]
+        print X[1,:]
+        print Y[1,:]
+
+        L_EX = [-0.00137913, -0.00098143,  0.0133289 ]
+        X_EX = [ 0.01998251,  0.01948814,  0.00247195]
+        Y_EX = [ -2.18223005e-04,   2.25636080e-04,  -5.36417393e-05]
+
+        
+        #L_EX = [-0.00825736, -0.00801881,  0.0272445 ]
+        #X_EX = [ 0.0170432 ,  0.01674667,  0.00654035]
+        #Y_EX = [ -1.56119371e-04,   1.10107449e-04,  -5.74034758e-05]
+
+        assert num.allclose(L[1,:], L_EX)
+        assert num.allclose(X[1,:], X_EX)
+        assert num.allclose(Y[1,:], Y_EX)
+
     def test_balance_deep_and_shallow(self):
         """Test that balanced limiters preserve conserved quantites.
@@ -8348 +8431 @@
 
 if __name__ == "__main__":
-    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_rainfall_forcing_with_evolve')
+    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_extrapolate_second_order_sw')
     suite = unittest.makeSuite(Test_Shallow_Water, 'test_')
     runner = unittest.TextTestRunner(verbosity=1)