Changeset 1387 for inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

Timestamp:

May 13, 2005, 6:15:08 PM (20 years ago)

Author:

steve

Message:

Need to look at uint test for inscribed circle

File:

• inundation/ga/storm_surge/pyvolution/shallow_water_ext.c (modified) (40 diffs)

inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

@@ -2 +2 @@
 //
 // To compile (Python2.3):
-//  gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O 
-//  gcc -shared domain_ext.o  -o domain_ext.so  
+//  gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O
+//  gcc -shared domain_ext.o  -o domain_ext.so
 //
 // or use python compile.py
@@ -11 +11 @@
 //
 // Ole Nielsen, GA 2004
-    
-    
+
+
 #include "Python.h"
 #include "Numeric/arrayobject.h"
 #include "math.h"
+#include <stdio.h>
 
 //Shared code snippets
@@ -26 +27 @@
   /*Rotate the momentum component q (q[1], q[2])
     from x,y coordinates to coordinates based on normal vector (n1, n2).
-    
-    Result is returned in array 3x1 r     
+
+    Result is returned in array 3x1 r
     To rotate in opposite direction, call rotate with (q, n1, -n2)
-    
-    Contents of q are changed by this function */    
+
+    Contents of q are changed by this function */
 
 
   double q1, q2;
-  
+
   //Shorthands
   q1 = q[1];  //uh momentum
@@ -41 +42 @@
   //Rotate
   q[1] =  n1*q1 + n2*q2;
-  q[2] = -n2*q1 + n1*q2; 
+  q[2] = -n2*q1 + n1*q2;
 
   return 0;
-}  
-
-
-        
+}
+
+
+
 // Computational function for flux computation (using stage w=z+h)
-int flux_function(double *q_left, double *q_right, 
-                  double z_left, double z_right, 
-                  double n1, double n2, 
-                  double epsilon, double g, 
+int flux_function(double *q_left, double *q_right,
+                  double z_left, double z_right,
+                  double n1, double n2,
+                  double epsilon, 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 
+    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 
+
+    The implemented formula is given in equation (3.15) on page 714
   */
-       
+
   int i;
-       
+
   double w_left, h_left, uh_left, vh_left, u_left;
   double w_right, h_right, uh_right, vh_right, u_right;
   double s_min, s_max, soundspeed_left, soundspeed_right;
   double denom, z;
-  double q_left_copy[3], q_right_copy[3];    
+  double q_left_copy[3], q_right_copy[3];
   double flux_right[3], flux_left[3];
-  
+
   //Copy conserved quantities to protect from modification
   for (i=0; i<3; i++) {
     q_left_copy[i] = q_left[i];
     q_right_copy[i] = q_right[i];
-  } 
-    
+  }
+
   //Align x- and y-momentum with x-axis
   _rotate(q_left_copy, n1, n2);
-  _rotate(q_right_copy, n1, n2);    
+  _rotate(q_right_copy, n1, n2);
 
   z = (z_left+z_right)/2; //Take average of field values
@@ -95 +96 @@
     h_left = 0.0;  //Could have been negative
     u_left = 0.0;
-  } else { 
+  } else {
     u_left = uh_left/h_left;
   }
-  
+
   w_right = q_right_copy[0];
   h_right = w_right-z;
@@ -106 +107 @@
     h_right = 0.0; //Could have been negative
     u_right = 0.0;
-  } else { 
+  } else {
     u_right = uh_right/h_right;
   }
 
-  //Momentum in y-direction              
+  //Momentum in y-direction
   vh_left  = q_left_copy[2];
-  vh_right = q_right_copy[2];   
-       
+  vh_right = q_right_copy[2];
+
 
   //Maximal and minimal wave speeds
-  soundspeed_left  = sqrt(g*h_left); 
+  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; 
-  
+  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;   
-  
-  //Flux formulas  
+  if (s_min > 0.0) s_min = 0.0;
+
+  //Flux formulas
   flux_left[0] = u_left*h_left;
   flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
   flux_left[2] = u_left*vh_left;
-  
+
   flux_right[0] = u_right*h_right;
   flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
   flux_right[2] = u_right*vh_right;
-    
-
-  //Flux computation    
+
+
+  //Flux computation
   denom = s_max-s_min;
   if (denom == 0.0) {
     for (i=0; i<3; i++) edgeflux[i] = 0.0;
     *max_speed = 0.0;
-  } else {    
+  } else {
     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_copy[i]-q_left_copy[i]);
       edgeflux[i] /= denom;
-    }  
-        
+    }
+
     //Maximal wavespeed
     *max_speed = max(fabs(s_max), fabs(s_min));
-    
-    //Rotate back       
+
+    //Rotate back
     _rotate(edgeflux, n1, -n2);
   }
   return 0;
 }
-        
-void _manning_friction(double g, double eps, int N, 
-                       double* w, double* z, 
-                       double* uh, double* vh, 
-                       double* eta, double* xmom, double* ymom) {     
+
+void _manning_friction(double g, double eps, int N,
+                       double* w, double* z,
+                       double* uh, double* vh,
+                       double* eta, double* xmom, double* ymom) {
 
   int k;
   double S, h;
-  
+
   for (k=0; k<N; k++) {
     if (eta[k] > eps) {
@@ -180 +181 @@
   }
 }
-            
+
 
 
 int _balance_deep_and_shallow(int N,
                               double* wc,
-                              double* zc, 
-                              double* hc,                             
-                              double* wv, 
-                              double* zv, 
+                              double* zc,
+                              double* hc,
+                              double* wv,
+                              double* zv,
                               double* hv,
-                              double* hvbar,                          
-                              double* xmomc, 
-                              double* ymomc, 
-                              double* xmomv, 
-                              double* ymomv) { 
-  
+                              double* hvbar,
+                              double* xmomc,
+                              double* ymomc,
+                              double* xmomv,
+                              double* ymomv) {
+
   int k, k3, i;
   double dz, hmin, alpha;
-  
+
   //Compute linear combination between w-limited stages and
-  //h-limited stages close to the bed elevation.     
-  
+  //h-limited stages close to the bed elevation.
+
   for (k=0; k<N; k++) {
     // Compute maximal variation in bed elevation
@@ -211 +212 @@
 
     k3 = 3*k;
-    
+
     //FIXME: Try with this one precomputed
     dz = 0.0;
@@ -220 +221 @@
     }
 
-    
-    //Create alpha in [0,1], where alpha==0 means using the h-limited 
+
+    //Create alpha in [0,1], where alpha==0 means using the h-limited
     //stage and alpha==1 means using the w-limited stage as
     //computed by the gradient limiter (both 1st or 2nd order)
@@ -227 +228 @@
     //If hmin > dz/2 then alpha = 1 and the bed will have no effect
     //If hmin < 0 then alpha = 0 reverting to constant height above bed.
-   
+
 
     if (dz > 0.0)
-      //if (hmin<0.0) 
+      //if (hmin<0.0)
       //        alpha = 0.0;
       //else
       //  alpha = max( min( hc[k]/dz, 1.0), 0.0 );
       alpha = max( min( 2.0*hmin/dz, 1.0), 0.0 );
-   
     else
       alpha = 1.0;  //Flat bed
 
+    //alpha = 1.0;
+
     //printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
 
-    //  Let 
-    //  
-    //    wvi be the w-limited stage (wvi = zvi + hvi)        
+    //  Let
+    //
+    //    wvi be the w-limited stage (wvi = zvi + hvi)
     //    wvi- be the h-limited state (wvi- = zvi + hvi-)
-    //      
-    //      
+    //
+    //
     //  where i=0,1,2 denotes the vertex ids
-    //  
-    //  Weighted balance between w-limited and h-limited stage is  
-    //    
-    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)   
-    //  
+    //
+    //  Weighted balance between w-limited and h-limited stage is
+    //
+    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
+    //
     //  It follows that the updated wvi is
     //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
-    //   
+    //
     //   Momentum is balanced between constant and limited
 
-    if (alpha < 1) {         
+    if (alpha < 1) {
       for (i=0; i<3; i++) {
         wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i];
-            
-        //Update momentum as a linear combination of 
+
+        //Update momentum as a linear combination of
         //xmomc and ymomc (shallow) and momentum
         //from extrapolator xmomv and ymomv (deep).
-        xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];            
+        xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
         ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
-      } 
+      }
     }
-  }         
+  }
   return 0;
 }
@@ -275 +277 @@
 
 
-int _protect(int N, 
-             double minimum_allowed_height,        
+int _protect(int N,
+             double minimum_allowed_height,
              double* wc,
-             double* zc, 
-             double* xmomc, 
+             double* zc,
+             double* xmomc,
              double* ymomc) {
-  
-  int k; 
+
+  int k;
   double hc;
-  
+
   //Protect against initesimal and negative heights
-  
+
   for (k=0; k<N; k++) {
     hc = wc[k] - zc[k];
-    
-    if (hc < minimum_allowed_height) {    
-      wc[k] = zc[k];       
+
+    if (hc < minimum_allowed_height) {
+      wc[k] = zc[k];
       xmomc[k] = 0.0;
-      ymomc[k] = 0.0;      
+      ymomc[k] = 0.0;
     }
-    
+
   }
   return 0;
@@ -309 +311 @@
 int _assign_wind_field_values(int N,
                               double* xmom_update,
-                              double* ymom_update, 
+                              double* ymom_update,
                               double* s_vec,
                               double* phi_vec,
@@ -318 +320 @@
   int k;
   double S, s, phi, u, v;
-  
+
   for (k=0; k<N; k++) {
-    
+
     s = s_vec[k];
     phi = phi_vec[k];
@@ -330 +332 @@
     u = s*cos(phi);
     v = s*sin(phi);
-        
+
     //Compute wind stress
     S = cw * sqrt(u*u + v*v);
     xmom_update[k] += S*u;
-    ymom_update[k] += S*v;        
-  }
-  return 0;  
-}            
+    ymom_update[k] += S*v;
+  }
+  return 0;
+}
 
 
@@ -348 +350 @@
   //  gravity(g, h, v, x, xmom, ymom)
   //
-  
-  
+
+
   PyArrayObject *h, *v, *x, *xmom, *ymom;
   int k, i, N, k3, k6;
   double g, avg_h, zx, zy;
   double x0, y0, x1, y1, x2, y2, z0, z1, z2;
-    
+
   if (!PyArg_ParseTuple(args, "dOOOOO",
-                        &g, &h, &v, &x, 
-                        &xmom, &ymom))  
-    return NULL;  
+                        &g, &h, &v, &x,
+                        &xmom, &ymom))
+    return NULL;
 
   N = h -> dimensions[0];
   for (k=0; k<N; k++) {
     k3 = 3*k;  // base index
-    k6 = 6*k;  // base index    
-    
+    k6 = 6*k;  // base index
+
     avg_h = 0.0;
     for (i=0; i<3; i++) {
       avg_h += ((double *) h -> data)[k3+i];
-    }   
+    }
     avg_h /= 3;
-        
-    
+
+
     //Compute bed slope
     x0 = ((double*) x -> data)[k6 + 0];
-    y0 = ((double*) x -> data)[k6 + 1];    
+    y0 = ((double*) x -> data)[k6 + 1];
     x1 = ((double*) x -> data)[k6 + 2];
-    y1 = ((double*) x -> data)[k6 + 3];        
+    y1 = ((double*) x -> data)[k6 + 3];
     x2 = ((double*) x -> data)[k6 + 4];
-    y2 = ((double*) x -> data)[k6 + 5];           
+    y2 = ((double*) x -> data)[k6 + 5];
 
 
     z0 = ((double*) v -> data)[k3 + 0];
     z1 = ((double*) v -> data)[k3 + 1];
-    z2 = ((double*) v -> data)[k3 + 2];        
+    z2 = ((double*) v -> data)[k3 + 2];
 
     _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
@@ -389 +391 @@
     //Update momentum
     ((double*) xmom -> data)[k] += -g*zx*avg_h;
-    ((double*) ymom -> data)[k] += -g*zy*avg_h;        
-  }
-    
+    ((double*) ymom -> data)[k] += -g*zy*avg_h;
+  }
+
   return Py_BuildValue("");
 }
@@ -400 +402 @@
   // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
   //
-  
-  
+
+
   PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
   int N;
   double g, eps;
-    
+
   if (!PyArg_ParseTuple(args, "ddOOOOOOO",
-                        &g, &eps, &w, &z, &uh, &vh, &eta, 
-                        &xmom, &ymom))  
-    return NULL;  
-
-  N = w -> dimensions[0];    
+                        &g, &eps, &w, &z, &uh, &vh, &eta,
+                        &xmom, &ymom))
+    return NULL;
+
+  N = w -> dimensions[0];
   _manning_friction(g, eps, N,
                     (double*) w -> data,
-                    (double*) z -> data,                    
-                    (double*) uh -> data, 
-                    (double*) vh -> data, 
+                    (double*) z -> data,
+                    (double*) uh -> data,
+                    (double*) vh -> data,
                     (double*) eta -> data,
-                    (double*) xmom -> data, 
+                    (double*) xmom -> data,
                     (double*) ymom -> data);
 
   return Py_BuildValue("");
-}                   
+}
 
 PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
@@ -431 +433 @@
   // normal a Float numeric array of length 2.
 
-  
+
   PyObject *Q, *Normal;
   PyArrayObject *q, *r, *normal;
-  
+
   static char *argnames[] = {"q", "normal", "direction", NULL};
   int dimensions[1], i, direction=1;
   double n1, n2;
 
-  // Convert Python arguments to C  
-  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames, 
-                                   &Q, &Normal, &direction))  
-    return NULL;  
+  // Convert Python arguments to C
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
+                                   &Q, &Normal, &direction))
+    return NULL;
 
   //Input checks (convert sequences into numeric arrays)
-  q = (PyArrayObject *) 
+  q = (PyArrayObject *)
     PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
-  normal = (PyArrayObject *) 
+  normal = (PyArrayObject *)
     PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
 
-  
+
   if (normal -> dimensions[0] != 2) {
     PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
@@ -456 +458 @@
   }
 
-  //Allocate space for return vector r (don't DECREF) 
+  //Allocate space for return vector r (don't DECREF)
   dimensions[0] = 3;
   r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
@@ -462 +464 @@
   //Copy
   for (i=0; i<3; i++) {
-    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];  
-  }
-  
+    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
+  }
+
   //Get normal and direction
-  n1 = ((double *) normal -> data)[0]; 
-  n2 = ((double *) normal -> data)[1];   
+  n1 = ((double *) normal -> data)[0];
+  n2 = ((double *) normal -> data)[1];
   if (direction == -1) n2 = -n2;
 
@@ -474 +476 @@
 
   //Release numeric arrays
-  Py_DECREF(q);    
+  Py_DECREF(q);
   Py_DECREF(normal);
-        
+
   //return result using PyArray to avoid memory leak
   return PyArray_Return(r);
-}    
+}
 
 
@@ -490 +492 @@
     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 
+    explicit_update for each of the three conserved quantities
     stage, xmomentum and ymomentum
 
@@ -496 +498 @@
     is converted to a timestep that must not be exceeded. The minimum of
     those is computed as the next overall timestep.
-    
+
     Python call:
-    domain.timestep = compute_fluxes(timestep, 
-                                     domain.epsilon, 
+    domain.timestep = compute_fluxes(timestep,
+                                     domain.epsilon,
                                      domain.g,
                                      domain.neighbours,
                                      domain.neighbour_edges,
                                      domain.normals,
-                                     domain.edgelengths,                       
+                                     domain.edgelengths,
                                      domain.radii,
                                      domain.areas,
-                                     Stage.edge_values, 
-                                     Xmom.edge_values, 
-                                     Ymom.edge_values,  
-                                     Bed.edge_values,   
+                                     Stage.edge_values,
+                                     Xmom.edge_values,
+                                     Ymom.edge_values,
+                                     Bed.edge_values,
                                      Stage.boundary_values,
                                      Xmom.boundary_values,
@@ -517 +519 @@
                                      Xmom.explicit_update,
                                      Ymom.explicit_update)
-        
+
 
     Post conditions:
       domain.explicit_update is reset to computed flux values
-      domain.timestep is set to the largest step satisfying all volumes. 
-
-            
+      domain.timestep is set to the largest step satisfying all volumes.
+
+
   */
 
-  
+
   PyArrayObject *neighbours, *neighbour_edges,
     *normals, *edgelengths, *radii, *areas,
-    *stage_edge_values, 
-    *xmom_edge_values, 
-    *ymom_edge_values,  
-    *bed_edge_values,   
+    *stage_edge_values,
+    *xmom_edge_values,
+    *ymom_edge_values,
+    *bed_edge_values,
     *stage_boundary_values,
     *xmom_boundary_values,
@@ -540 +542 @@
     *ymom_explicit_update;
 
-    
-  //Local variables  
+
+  //Local variables
   double timestep, max_speed, epsilon, g;
   double normal[2], ql[3], qr[3], zl, zr;
@@ -548 +550 @@
   int number_of_elements, k, i, j, m, n;
   int ki, nm, ki2; //Index shorthands
-  
- 
-  // Convert Python arguments to C  
+
+
+  // Convert Python arguments to C
   if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOO",
                         &timestep,
                         &epsilon,
                         &g,
-                        &neighbours, 
+                        &neighbours,
                         &neighbour_edges,
-                        &normals, 
+                        &normals,
                         &edgelengths, &radii, &areas,
-                        &stage_edge_values, 
-                        &xmom_edge_values, 
-                        &ymom_edge_values,  
-                        &bed_edge_values,   
+                        &stage_edge_values,
+                        &xmom_edge_values,
+                        &ymom_edge_values,
+                        &bed_edge_values,
                         &stage_boundary_values,
                         &xmom_boundary_values,
@@ -574 +576 @@
 
   number_of_elements = stage_edge_values -> dimensions[0];
-  
-    
+
+
   for (k=0; k<number_of_elements; k++) {
 
     //Reset work array
     for (j=0; j<3; j++) flux[j] = 0.0;
-                          
+
     //Loop through neighbours and compute edge flux for each
     for (i=0; i<3; i++) {
@@ -586 +588 @@
       ql[0] = ((double *) stage_edge_values -> data)[ki];
       ql[1] = ((double *) xmom_edge_values -> data)[ki];
-      ql[2] = ((double *) ymom_edge_values -> data)[ki];            
-      zl =    ((double *) bed_edge_values -> data)[ki];                  
+      ql[2] = ((double *) ymom_edge_values -> data)[ki];
+      zl =    ((double *) bed_edge_values -> data)[ki];
 
       //Quantities at neighbour on nearest face
@@ -593 +595 @@
       if (n < 0) {
         m = -n-1; //Convert negative flag to index
-        qr[0] = ((double *) stage_boundary_values -> data)[m];  
-        qr[1] = ((double *) xmom_boundary_values -> data)[m];   
-        qr[2] = ((double *) ymom_boundary_values -> data)[m];   
+        qr[0] = ((double *) stage_boundary_values -> data)[m];
+        qr[1] = ((double *) xmom_boundary_values -> data)[m];
+        qr[2] = ((double *) ymom_boundary_values -> data)[m];
         zr = zl; //Extend bed elevation to boundary
-      } else {    
+      } else {
         m = ((long *) neighbour_edges -> data)[ki];
-        
-        nm = n*3+m;     
+
+        nm = n*3+m;
         qr[0] = ((double *) stage_edge_values -> data)[nm];
         qr[1] = ((double *) xmom_edge_values -> data)[nm];
-        qr[2] = ((double *) ymom_edge_values -> data)[nm];            
-        zr =    ((double *) bed_edge_values -> data)[nm];                  
+        qr[2] = ((double *) ymom_edge_values -> data)[nm];
+        zr =    ((double *) bed_edge_values -> data)[nm];
       }
 
       //printf("%d %d [%d] (%d, %d): %.2f %.2f %.2f\n", k, i, ki, n, m,
-      //             ql[0], ql[1], ql[2]);       
-
-      
-      // Outward pointing normal vector   
+      //             ql[0], ql[1], ql[2]);
+
+
+      // Outward pointing normal vector
       // normal = domain.normals[k, 2*i:2*i+2]
       ki2 = 2*ki; //k*6 + i*2
       normal[0] = ((double *) normals -> data)[ki2];
-      normal[1] = ((double *) normals -> data)[ki2+1];      
+      normal[1] = ((double *) normals -> data)[ki2+1];
 
       //Edge flux computation
-      flux_function(ql, qr, zl, zr, 
+      flux_function(ql, qr, zl, zr,
                     normal[0], normal[1],
-                    epsilon, g, 
+                    epsilon, g,
                     edgeflux, &max_speed);
 
-                    
+
       //flux -= edgeflux * edgelengths[k,i]
-      for (j=0; j<3; j++) { 
+      for (j=0; j<3; j++) {
         flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
       }
-      
+
       //Update timestep
       //timestep = min(timestep, domain.radii[k]/max_speed)
@@ -634 +636 @@
       if (max_speed > epsilon) {
         timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
-      }   
+      }
     } // end for i
-    
+
     //Normalise by area and store for when all conserved
     //quantities get updated
     // flux /= areas[k]
-    for (j=0; j<3; j++) { 
+    for (j=0; j<3; j++) {
       flux[j] /= ((double *) areas -> data)[k];
     }
@@ -646 +648 @@
     ((double *) stage_explicit_update -> data)[k] = flux[0];
     ((double *) xmom_explicit_update -> data)[k] = flux[1];
-    ((double *) ymom_explicit_update -> data)[k] = flux[2];        
+    ((double *) ymom_explicit_update -> data)[k] = flux[2];
 
   } //end for k
 
   return Py_BuildValue("d", timestep);
-}    
+}
 
 
@@ -658 +660 @@
   //
   //    protect(minimum_allowed_height, wc, zc, xmomc, ymomc)
-  
-
-  PyArrayObject 
+
+
+  PyArrayObject
   *wc,            //Stage at centroids
-  *zc,            //Elevation at centroids    
+  *zc,            //Elevation at centroids
   *xmomc,         //Momentums at centroids
-  *ymomc; 
-
-    
-  int N; 
+  *ymomc;
+
+
+  int N;
   double minimum_allowed_height;
-  
-  // Convert Python arguments to C  
-  if (!PyArg_ParseTuple(args, "dOOOO",  
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dOOOO",
                         &minimum_allowed_height,
                         &wc, &zc, &xmomc, &ymomc))
@@ -677 +679 @@
 
   N = wc -> dimensions[0];
-    
+
   _protect(N,
            minimum_allowed_height,
            (double*) wc -> data,
-           (double*) zc -> data, 
-           (double*) xmomc -> data, 
+           (double*) zc -> data,
+           (double*) xmomc -> data,
            (double*) ymomc -> data);
- 
-  return Py_BuildValue("");  
+
+  return Py_BuildValue("");
 }
 
@@ -694 +696 @@
   //    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
   //                             xmomc, ymomc, xmomv, ymomv)
-  
-
-  PyArrayObject 
+
+
+  PyArrayObject
     *wc,            //Stage at centroids
-    *zc,            //Elevation at centroids    
-    *hc,            //Height at centroids        
+    *zc,            //Elevation at centroids
+    *hc,            //Height at centroids
     *wv,            //Stage at vertices
     *zv,            //Elevation at vertices
-    *hv,            //Depths at vertices    
-    *hvbar,         //h-Limited depths at vertices        
+    *hv,            //Depths at vertices
+    *hvbar,         //h-Limited depths at vertices
     *xmomc,         //Momentums at centroids and vertices
-    *ymomc, 
-    *xmomv, 
-    *ymomv;    
-    
+    *ymomc,
+    *xmomv,
+    *ymomv;
+
   int N; //, err;
-  
-  // Convert Python arguments to C  
-  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",  
-                        &wc, &zc, &hc, 
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",
+                        &wc, &zc, &hc,
                         &wv, &zv, &hv, &hvbar,
                         &xmomc, &ymomc, &xmomv, &ymomv))
@@ -719 +721 @@
 
   N = wc -> dimensions[0];
-    
+
   _balance_deep_and_shallow(N,
                             (double*) wc -> data,
-                            (double*) zc -> data, 
-                            (double*) hc -> data,                           
-                            (double*) wv -> data, 
-                            (double*) zv -> data, 
+                            (double*) zc -> data,
+                            (double*) hc -> data,
+                            (double*) wv -> data,
+                            (double*) zv -> data,
                             (double*) hv -> data,
-                            (double*) hvbar -> data,                        
-                            (double*) xmomc -> data, 
-                            (double*) ymomc -> data, 
-                            (double*) xmomv -> data, 
-                            (double*) ymomv -> data);  
-  
- 
-  return Py_BuildValue("");  
+                            (double*) hvbar -> data,
+                            (double*) xmomc -> data,
+                            (double*) ymomc -> data,
+                            (double*) xmomv -> data,
+                            (double*) ymomv -> data);
+
+
+  return Py_BuildValue("");
 }
 
@@ -740 +742 @@
 
 PyObject *h_limiter(PyObject *self, PyObject *args) {
-  
+
   PyObject *domain, *Tmp;
-  PyArrayObject 
-    *hv, *hc, //Depth at vertices and centroids        
+  PyArrayObject
+    *hv, *hc, //Depth at vertices and centroids
     *hvbar,   //Limited depth at vertices (return values)
     *neighbours;
-    
+
   int k, i, n, N, k3;
   int dimensions[2];
   double beta_h; //Safety factor (see config.py)
   double *hmin, *hmax, hn;
-  
-  // Convert Python arguments to C  
-  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) 
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
     return NULL;
 
@@ -759 +761 @@
 
   //Get safety factor beta_h
-  Tmp = PyObject_GetAttrString(domain, "beta_h"); 
-  if (!Tmp) 
-    return NULL;
-      
-  beta_h = PyFloat_AsDouble(Tmp);  
-  
-  Py_DECREF(Tmp);    
+  Tmp = PyObject_GetAttrString(domain, "beta_h");
+  if (!Tmp)
+    return NULL;
+
+  beta_h = PyFloat_AsDouble(Tmp);
+
+  Py_DECREF(Tmp);
 
   N = hc -> dimensions[0];
-  
+
   //Create hvbar
   dimensions[0] = N;
-  dimensions[1] = 3;  
+  dimensions[1] = 3;
   hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
- 
+
 
   //Find min and max of this and neighbour's centroid values
   hmin = malloc(N * sizeof(double));
-  hmax = malloc(N * sizeof(double));  
+  hmax = malloc(N * sizeof(double));
   for (k=0; k<N; k++) {
     k3=k*3;
-    
+
     hmin[k] = ((double*) hc -> data)[k];
     hmax[k] = hmin[k];
-    
+
     for (i=0; i<3; i++) {
       n = ((long*) neighbours -> data)[k3+i];
-      
+
       //Initialise hvbar with values from hv
-      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];  
-      
+      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
+
       if (n >= 0) {
         hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
-          
+
         hmin[k] = min(hmin[k], hn);
         hmax[k] = max(hmax[k], hn);
@@ -798 +800 @@
     }
   }
-  
+
   // Call underlying standard routine
   _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
-         
-  // // //Py_DECREF(domain); //FIXME: NEcessary?         
+
+  // // //Py_DECREF(domain); //FIXME: NEcessary?
   free(hmin);
-  free(hmax);  
-  
-  //return result using PyArray to avoid memory leak  
+  free(hmax);
+
+  //return result using PyArray to avoid memory leak
   return PyArray_Return(hvbar);
-  //return Py_BuildValue("");    
+  //return Py_BuildValue("");
 }
 
@@ -819 +821 @@
   //                              s_vec, phi_vec, self.const)
 
-  
+
 
   PyArrayObject   //(one element per triangle)
-  *s_vec,         //Speeds 
-  *phi_vec,       //Bearings 
+  *s_vec,         //Speeds
+  *phi_vec,       //Bearings
   *xmom_update,   //Momentum updates
-  *ymom_update; 
-
-    
-  int N; 
+  *ymom_update;
+
+
+  int N;
   double cw;
-  
-  // Convert Python arguments to C  
-  if (!PyArg_ParseTuple(args, "OOOOd",  
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOOOd",
                         &xmom_update,
-                        &ymom_update,                   
-                        &s_vec, &phi_vec, 
+                        &ymom_update,
+                        &s_vec, &phi_vec,
                         &cw))
     return NULL;
 
   N = xmom_update -> dimensions[0];
-    
+
   _assign_wind_field_values(N,
            (double*) xmom_update -> data,
-           (double*) ymom_update -> data, 
-           (double*) s_vec -> data, 
+           (double*) ymom_update -> data,
+           (double*) s_vec -> data,
            (double*) phi_vec -> data,
            cw);
- 
-  return Py_BuildValue("");  
-}
-
-
-
-
-//////////////////////////////////////////      
+
+  return Py_BuildValue("");
+}
+
+
+
+
+//////////////////////////////////////////
 // Method table for python module
 static struct PyMethodDef MethodTable[] = {
@@ -861 +863 @@
    * three.
    */
-  
+
   {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
-  {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},    
-  {"gravity", gravity, METH_VARARGS, "Print out"},      
-  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},        
-  {"balance_deep_and_shallow", balance_deep_and_shallow, 
-   METH_VARARGS, "Print out"},          
-  {"h_limiter", h_limiter, 
-   METH_VARARGS, "Print out"},             
-  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},   
-  {"assign_windfield_values", assign_windfield_values, 
-   METH_VARARGS | METH_KEYWORDS, "Print out"},     
-  //{"distribute_to_vertices_and_edges", 
-  // distribute_to_vertices_and_edges, METH_VARARGS},    
-  //{"update_conserved_quantities", 
-  // update_conserved_quantities, METH_VARARGS},       
-  //{"set_initialcondition", 
-  // set_initialcondition, METH_VARARGS},    
+  {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},
+  {"gravity", gravity, METH_VARARGS, "Print out"},
+  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
+  {"balance_deep_and_shallow", balance_deep_and_shallow,
+   METH_VARARGS, "Print out"},
+  {"h_limiter", h_limiter,
+   METH_VARARGS, "Print out"},
+  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
+  {"assign_windfield_values", assign_windfield_values,
+   METH_VARARGS | METH_KEYWORDS, "Print out"},
+  //{"distribute_to_vertices_and_edges",
+  // distribute_to_vertices_and_edges, METH_VARARGS},
+  //{"update_conserved_quantities",
+  // update_conserved_quantities, METH_VARARGS},
+  //{"set_initialcondition",
+  // set_initialcondition, METH_VARARGS},
   {NULL, NULL}
 };
-        
-// Module initialisation   
+
+// Module initialisation
 void initshallow_water_ext(void){
   Py_InitModule("shallow_water_ext", MethodTable);
-  
-  import_array();     //Necessary for handling of NumPY structures  
-}
+
+  import_array();     //Necessary for handling of NumPY structures
+}