Changeset 6902 for branches/numpy/anuga/shallow_water

Timestamp:

Apr 24, 2009, 5:22:14 PM (16 years ago)

Author:

rwilson

Message:

Back-merge from Numeric trunk to numpy branch.

Location:

branches/numpy/anuga/shallow_water

Files:

• data_manager.py (modified) (10 diffs)
• shallow_water_ext.c (modified) (86 diffs)
• test_data_manager.py (modified) (33 diffs)
• test_shallow_water_domain.py (modified) (11 diffs)
• urs_ext.c (modified) (6 diffs)

branches/numpy/anuga/shallow_water/data_manager.py

@@ -3575 +3575 @@
     other_quantities.remove('z')
     other_quantities.remove('volumes')
-    #try:
-    #    other_quantities.remove('stage_range')
-    #    other_quantities.remove('xmomentum_range')
-    #    other_quantities.remove('ymomentum_range')
-    #    other_quantities.remove('elevation_range')
-    #except:
-    #    pass
+    try:
+        other_quantities.remove('stage_range')
+        other_quantities.remove('xmomentum_range')
+        other_quantities.remove('ymomentum_range')
+        other_quantities.remove('elevation_range')
+    except:
+        pass
 
     conserved_quantities.remove('time')
@@ -5338 +5338 @@
         j += 1
 
-    #if verbose: sww.verbose_quantities(outfile)
+    if verbose: sww.verbose_quantities(outfile)
 
     outfile.close()
@@ -5855 +5855 @@
     # @param smoothing True if smoothing is to be used.
     # @param order 
-    # @param sww_precision Data type of the quantitiy written (netcdf constant)
+    # @param sww_precision Data type of the quantity written (netcdf constant)
     # @param verbose True if this function is to be verbose.
     # @note If 'times' is a list, the info will be made relative.
@@ -5936 +5936 @@
                                ('number_of_points',))
         q = 'elevation'
-        #outfile.createVariable(q + Write_sww.RANGE, sww_precision,
-        #                       ('numbers_in_range',))
-
-
-        ## Initialise ranges with small and large sentinels.
-        ## If this was in pure Python we could have used None sensibly
-        #outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
-        #outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
+        outfile.createVariable(q + Write_sww.RANGE, sww_precision,
+                               ('numbers_in_range',))
+
+        # Initialise ranges with small and large sentinels.
+        # If this was in pure Python we could have used None sensibly
+        outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
+        outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
 
         # FIXME: Backwards compatibility
@@ -5958 +5957 @@
             outfile.createVariable(q, sww_precision, ('number_of_timesteps',
                                                       'number_of_points'))
-            #outfile.createVariable(q + Write_sww.RANGE, sww_precision,
-            #                       ('numbers_in_range',))
+            outfile.createVariable(q + Write_sww.RANGE, sww_precision,
+                                   ('numbers_in_range',))
 
             # Initialise ranges with small and large sentinels.
             # If this was in pure Python we could have used None sensibly
-            #outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
-            #outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
+            outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
+            outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
 
         if isinstance(times, (list, num.ndarray)):
@@ -6073 +6072 @@
         q = 'elevation'
         # This updates the _range values
-        #outfile.variables[q + Write_sww.RANGE][0] = min(elevation)
-        #outfile.variables[q + Write_sww.RANGE][1] = max(elevation)
+        outfile.variables[q + Write_sww.RANGE][0] = num.min(elevation)
+        outfile.variables[q + Write_sww.RANGE][1] = num.max(elevation)
 
 
@@ -6126 +6125 @@
                                 q_values.astype(sww_precision)
         
-        #        # This updates the _range values
-        #        q_range = outfile.variables[q + Write_sww.RANGE][:]
-        #        q_values_min = min(q_values)
-        #        if q_values_min < q_range[0]:
-        #            outfile.variables[q + Write_sww.RANGE][0] = q_values_min
-        #        q_values_max = max(q_values)
-        #        if q_values_max > q_range[1]:
-        #            outfile.variables[q + Write_sww.RANGE][1] = q_values_max
+                # This updates the _range values
+                q_range = outfile.variables[q + Write_sww.RANGE][:]
+                q_values_min = num.min(q_values)
+                if q_values_min < q_range[0]:
+                    outfile.variables[q + Write_sww.RANGE][0] = q_values_min
+                q_values_max = num.max(q_values)
+                if q_values_max > q_range[1]:
+                    outfile.variables[q + Write_sww.RANGE][1] = q_values_max
 
     ##
     # @brief Print the quantities in the underlying file.
     # @param outfile UNUSED.
-    #def verbose_quantities(self, outfile):
-    #    print '------------------------------------------------'
-    #    print 'More Statistics:'
-    #    for q in Write_sww.sww_quantities:
-    #        print '  %s in [%f, %f]' % (q,
-    #                                    outfile.variables[q+Write_sww.RANGE][0],
-    #                                    outfile.variables[q+Write_sww.RANGE][1])
-    #    print '------------------------------------------------'
+    def verbose_quantities(self, outfile):
+        print '------------------------------------------------'
+        print 'More Statistics:'
+        for q in Write_sww.sww_quantities:
+            print '  %s in [%f, %f]' % (q,
+                                        outfile.variables[q+Write_sww.RANGE][0],
+                                        outfile.variables[q+Write_sww.RANGE][1])
+        print '------------------------------------------------'
 
 
@@ -6496 +6495 @@
 
                 # This updates the _range values
-                #q_range = outfile.variables[q + Write_sts.RANGE][:]
-                #q_values_min = min(q_values)
-                #if q_values_min < q_range[0]:
-                #    outfile.variables[q + Write_sts.RANGE][0] = q_values_min
-                #q_values_max = max(q_values)
-                #if q_values_max > q_range[1]:
-                #    outfile.variables[q + Write_sts.RANGE][1] = q_values_max
+                q_range = outfile.variables[q + Write_sts.RANGE][:]
+                q_values_min = num.min(q_values)
+                if q_values_min < q_range[0]:
+                    outfile.variables[q + Write_sts.RANGE][0] = q_values_min
+                q_values_max = num.max(q_values)
+                if q_values_max > q_range[1]:
+                    outfile.variables[q + Write_sts.RANGE][1] = q_values_max
 
 
@@ -6673 +6672 @@
 
 
-# FIXME (DSG): Add unit test, make general, not just 2 files,
-# but any number of files.
 ##
 # @brief Copy a file to a directory, and optionally append another file to it.
@@ -6680 +6677 @@
 # @param filename Path to file to copy to directory 'dir_name'.
 # @param filename2 Optional path to file to append to copied file.
-def copy_code_files(dir_name, filename1, filename2=None):
+# @param verbose True if this function is to be verbose.
+# @note Allow filenames to be either a string or sequence of strings.
+def copy_code_files(dir_name, filename1, filename2=None, verbose=False):
     """Copies "filename1" and "filename2" to "dir_name".
-    Very useful for information management
-    filename1 and filename2 are both absolute pathnames
-    """
-
+
+    Each 'filename' may be a string or list of filename strings.
+
+    Filenames must be absolute pathnames
+    """
+
+    ##
+    # @brief copies a file or sequence to destination directory.
+    # @param dest The destination directory to copy to.
+    # @param file A filename string or sequence of filename strings.
+    def copy_file_or_sequence(dest, file):
+        if hasattr(file, '__iter__'):
+            for f in file:
+                shutil.copy(f, dir_name)
+                if verbose:
+                    print 'File %s copied' % (f)
+        else:
+            shutil.copy(file, dir_name)
+            if verbose:
+                print 'File %s copied' % (file)
+
+    # check we have a destination directory, create if necessary
     if access(dir_name, F_OK) == 0:
-        print 'Make directory %s' % dir_name
-        mkdir (dir_name,0777)
-
-    shutil.copy(filename1, dir_name + sep + basename(filename1))
-
-    if filename2 != None:
-        shutil.copy(filename2, dir_name + sep + basename(filename2))
-        print 'Files %s and %s copied' %(filename1, filename2)
-    else:
-        print 'File %s copied' %(filename1)
+        if verbose:
+            print 'Make directory %s' % dir_name
+        mkdir(dir_name, 0777)
+
+    copy_file_or_sequence(dir_name, filename1)
+
+    if not filename2 is None:
+        copy_file_or_sequence(dir_name, filename2)
 
 

branches/numpy/anuga/shallow_water/shallow_water_ext.c

@@ -28 +28 @@
 
 // Computational function for rotation
+// FIXME: Perhaps inline this and profile
 int _rotate(double *q, double n1, double n2) {
   /*Rotate the momentum component q (q[1], q[2])
@@ -52 +53 @@
 
 int find_qmin_and_qmax(double dq0, double dq1, double dq2, 
-                       double *qmin, double *qmax){
+               double *qmin, double *qmax){
   // Considering the centroid of an FV triangle and the vertices of its 
   // auxiliary triangle, find 
@@ -67 +68 @@
     if (dq1>=dq2){
       if (dq1>=0.0)
-        *qmax=dq0+dq1;
+    *qmax=dq0+dq1;
       else
-        *qmax=dq0;
-        
+    *qmax=dq0;
+    
       *qmin=dq0+dq2;
       if (*qmin>=0.0) *qmin = 0.0;
@@ -76 +77 @@
     else{// dq1<dq2
       if (dq2>0)
-        *qmax=dq0+dq2;
+    *qmax=dq0+dq2;
       else
-        *qmax=dq0;
-        
-      *qmin=dq0+dq1;    
+    *qmax=dq0;
+    
+      *qmin=dq0+dq1;    
       if (*qmin>=0.0) *qmin=0.0;
     }
@@ -87 +88 @@
     if (dq1<=dq2){
       if (dq1<0.0)
-        *qmin=dq0+dq1;
+    *qmin=dq0+dq1;
       else
-        *qmin=dq0;
-        
-      *qmax=dq0+dq2;    
+    *qmin=dq0;
+    
+      *qmax=dq0+dq2;    
       if (*qmax<=0.0) *qmax=0.0;
     }
     else{// dq1>dq2
       if (dq2<0.0)
-        *qmin=dq0+dq2;
+    *qmin=dq0+dq2;
       else
-        *qmin=dq0;
-        
+    *qmin=dq0;
+    
       *qmax=dq0+dq1;
       if (*qmax<=0.0) *qmax=0.0;
@@ -133 +134 @@
   
   phi=min(r*beta_w,1.0);
-  for (i=0;i<3;i++)
-    dqv[i]=dqv[i]*phi;
+  //for (i=0;i<3;i++)
+  dqv[0]=dqv[0]*phi;
+  dqv[1]=dqv[1]*phi;
+  dqv[2]=dqv[2]*phi;
     
   return 0;
@@ -141 +144 @@
 
 double compute_froude_number(double uh,
-                             double h, 
-                             double g,
-                             double epsilon) {
-                          
+                 double h, 
+                 double g,
+                 double epsilon) {
+              
   // Compute Froude number; v/sqrt(gh)
   // FIXME (Ole): Not currently in use
@@ -166 +169 @@
 // This is used by flux functions
 // Input parameters uh and h may be modified by this function.
+
+// FIXME: Perhaps inline this and profile
 double _compute_speed(double *uh, 
-                      double *h, 
-                      double epsilon, 
-                      double h0) {
+              double *h, 
+              double epsilon, 
+              double h0) {
   
   double u;
@@ -188 +193 @@
 }
 
+
+// Optimised squareroot computation
+//
+//See 
+//http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
+//and http://mail.opensolaris.org/pipermail/tools-gcc/2005-August/000047.html
+float fast_squareroot_approximation(float number) {
+  float x; 
+  const float f = 1.5;
+
+  // Allow i and y to occupy the same memory
+  union
+  {
+    long i;
+    float y;
+  } u;
+  
+  // Good initial guess  
+  x = number * 0.5;  
+  u.y  = number;
+  u.i  = 0x5f3759df - ( u.i >> 1 );
+  
+  // Take a few iterations
+  u.y  = u.y * ( f - ( x * u.y * u.y ) );
+  u.y  = u.y * ( f - ( x * u.y * u.y ) );
+    
+  return number * u.y;
+}
+
+
+
+// Optimised squareroot computation (double version, slower)
+double Xfast_squareroot_approximation(double number) {
+  double x;
+  const double f = 1.5;
+    
+  // Allow i and y to occupy the same memory    
+  union
+  {
+    long long i;
+    double y;
+  } u;
+
+  // Good initial guess
+  x = number * 0.5;
+  u.y  = number;
+  u.i  = 0x5fe6ec85e7de30daLL - ( u.i >> 1 );
+  
+  // Take a few iterations  
+  u.y  = u.y * ( f - ( x * u.y * u.y ) );
+  u.y  = u.y * ( f - ( x * u.y * u.y ) );
+
+  return number * u.y;
+}
+
+
 // Innermost flux function (using stage w=z+h)
 int _flux_function_central(double *q_left, double *q_right,
-                           double z_left, double z_right,
-                           double n1, double n2,
-                           double epsilon, double H0, double g,
-                           double *edgeflux, double *max_speed) {
+                           double z_left, double z_right,
+                           double n1, double n2,
+                           double epsilon, double H0, double g,
+                           double *edgeflux, double *max_speed) 
+{
 
   /*Compute fluxes between volumes for the shallow water wave equation
@@ -212 +274 @@
   double v_left, v_right;  
   double s_min, s_max, soundspeed_left, soundspeed_right;
-  double denom, z;
+  double denom, inverse_denominator, z;
 
   // Workspace (allocate once, use many)
@@ -219 +281 @@
   double h0 = H0*H0; // This ensures a good balance when h approaches H0.
                      // But evidence suggests that h0 can be as little as
-                     // epsilon!
+             // epsilon!
   
   // Copy conserved quantities to protect from modification
-  for (i=0; i<3; i++) {
-    q_left_rotated[i] = q_left[i];
-    q_right_rotated[i] = q_right[i];
-  }
+  q_left_rotated[0] = q_left[0];
+  q_right_rotated[0] = q_right[0];
+  q_left_rotated[1] = q_left[1];
+  q_right_rotated[1] = q_right[1];
+  q_left_rotated[2] = q_left[2];
+  q_right_rotated[2] = q_right[2];
 
   // Align x- and y-momentum with x-axis
@@ -231 +295 @@
   _rotate(q_right_rotated, n1, n2);
 
-  z = (z_left+z_right)/2; // Average elevation values. 
-                          // Even though this will nominally allow for discontinuities 
-                          // in the elevation data, there is currently no numerical 
-                          // support for this so results may be strange near jumps in the bed.
+  z = 0.5*(z_left + z_right); // Average elevation values. 
+                            // Even though this will nominally allow for discontinuities 
+                            // in the elevation data, there is currently no numerical 
+                            // support for this so results may be strange near jumps in the bed.
 
   // Compute speeds in x-direction
   w_left = q_left_rotated[0];          
-  h_left = w_left-z;
+  h_left = w_left - z;
   uh_left = q_left_rotated[1];
   u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
 
   w_right = q_right_rotated[0];
-  h_right = w_right-z;
+  h_right = w_right - z;
   uh_right = q_right_rotated[1];
   u_right = _compute_speed(&uh_right, &h_right, epsilon, h0);  
@@ -257 +321 @@
   // Maximal and minimal wave speeds
   soundspeed_left  = sqrt(g*h_left);
-  soundspeed_right = sqrt(g*h_right);
-
-  s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
-  if (s_max < 0.0) s_max = 0.0;
-
-  s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
-  if (s_min > 0.0) s_min = 0.0;
-
+  soundspeed_right = sqrt(g*h_right);  
+  
+  // Code to use fast square root optimisation if desired.
+  //soundspeed_left  = fast_squareroot_approximation(g*h_left);
+  //soundspeed_right = fast_squareroot_approximation(g*h_right);
+
+  s_max = max(u_left + soundspeed_left, u_right + soundspeed_right);
+  if (s_max < 0.0) 
+  {
+    s_max = 0.0;
+  }
+
+  s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
+  if (s_min > 0.0)
+  {
+    s_min = 0.0;
+  }
   
   // Flux formulas
@@ -275 +348 @@
   flux_right[2] = u_right*vh_right;
 
-  
   // Flux computation
-  denom = s_max-s_min;
-  if (denom < epsilon) { // FIXME (Ole): Try using H0 here
-    for (i=0; i<3; i++) edgeflux[i] = 0.0;
+  denom = s_max - s_min;
+  if (denom < epsilon) 
+  { // FIXME (Ole): Try using H0 here
+    memset(edgeflux, 0, 3*sizeof(double));
     *max_speed = 0.0;
-  } else {
-    for (i=0; i<3; i++) {
+  } 
+  else 
+  {
+    inverse_denominator = 1.0/denom;
+    for (i = 0; i < 3; i++) 
+    {
       edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
-      edgeflux[i] += s_max*s_min*(q_right_rotated[i]-q_left_rotated[i]);
-      edgeflux[i] /= denom;
+      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i]);
+      edgeflux[i] *= inverse_denominator;
     }
 
@@ -315 +392 @@
 // Computational function for flux computation (using stage w=z+h)
 int flux_function_kinetic(double *q_left, double *q_right,
-                  double z_left, double z_right,
-                  double n1, double n2,
-                  double epsilon, double H0, double g,
-                  double *edgeflux, double *max_speed) {
+          double z_left, double z_right,
+          double n1, double n2,
+          double epsilon, double H0, double g,
+          double *edgeflux, double *max_speed) {
 
   /*Compute fluxes between volumes for the shallow water wave equation
@@ -413 +490 @@
 
 void _manning_friction(double g, double eps, int N,
-                       double* w, double* z,
-                       double* uh, double* vh,
-                       double* eta, double* xmom, double* ymom) {
+               double* w, double* z,
+               double* uh, double* vh,
+               double* eta, double* xmom, double* ymom) {
 
   int k;
@@ -426 +503 @@
         S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
         S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
-        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
+        //S /= exp((7.0/3.0)*log(h));      //seems to save about 15% over manning_friction
         //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
 
@@ -441 +518 @@
 /*
 void _manning_friction_explicit(double g, double eps, int N,
-                       double* w, double* z,
-                       double* uh, double* vh,
-                       double* eta, double* xmom, double* ymom) {
+               double* w, double* z,
+               double* uh, double* vh,
+               double* eta, double* xmom, double* ymom) {
 
   int k;
@@ -452 +529 @@
       h = w[k]-z[k];
       if (h >= eps) {
-        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
-        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
-        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
-        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
-
-
-        //Update momentum
-        xmom[k] += S*uh[k];
-        ymom[k] += S*vh[k];
+    S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
+    S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
+    //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
+    //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
+
+
+    //Update momentum
+    xmom[k] += S*uh[k];
+    ymom[k] += S*vh[k];
       }
     }
@@ -471 +548 @@
 
 double velocity_balance(double uh_i, 
-                        double uh,
-                        double h_i, 
-                        double h, 
-                        double alpha,
-                        double epsilon) {
+            double uh,
+            double h_i, 
+            double h, 
+            double alpha,
+            double epsilon) {
   // Find alpha such that speed at the vertex is within one
-  // order of magnitude of the centroid speed   
+  // order of magnitude of the centroid speed   
 
   // FIXME(Ole): Work in progress
@@ -486 +563 @@
   
   printf("alpha = %f, uh_i=%f, uh=%f, h_i=%f, h=%f\n",
-         alpha, uh_i, uh, h_i, h);
+     alpha, uh_i, uh, h_i, h);
       
   
@@ -500 +577 @@
   if (h < epsilon) {
     b = 1.0e10; // Limit
-  } else {       
+  } else {   
     b = m*fabs(h_i - h)/h;
   }
@@ -507 +584 @@
 
   if (a > b) {
-    estimate = (m-1)/(a-b);             
+    estimate = (m-1)/(a-b);     
     
     printf("Alpha %f, estimate %f\n", 
-           alpha, estimate);    
-           
+       alpha, estimate);    
+       
     if (alpha < estimate) {
       printf("Adjusting alpha from %f to %f\n", 
-             alpha, estimate);
+         alpha, estimate);
       alpha = estimate;
     }
@@ -520 +597 @@
   
     if (h < h_i) {
-      estimate = (m-1)/(a-b);                 
+      estimate = (m-1)/(a-b);             
     
       printf("Alpha %f, estimate %f\n", 
-             alpha, estimate);    
-           
+         alpha, estimate);    
+       
       if (alpha < estimate) {
-        printf("Adjusting alpha from %f to %f\n", 
-               alpha, estimate);
-        alpha = estimate;
+    printf("Adjusting alpha from %f to %f\n", 
+           alpha, estimate);
+    alpha = estimate;
       }    
     }
@@ -540 +617 @@
 
 int _balance_deep_and_shallow(int N,
-                              double* wc,
-                              double* zc,
-                              double* wv,
-                              double* zv,
-                              double* hvbar, // Retire this
-                              double* xmomc,
-                              double* ymomc,
-                              double* xmomv,
-                              double* ymomv,
-                              double H0,
-                              int tight_slope_limiters,
-                              int use_centroid_velocities,
-                              double alpha_balance) {
+                  double* wc,
+                  double* zc,
+                  double* wv,
+                  double* zv,
+                  double* hvbar, // Retire this
+                  double* xmomc,
+                  double* ymomc,
+                  double* xmomv,
+                  double* ymomv,
+                  double H0,
+                  int tight_slope_limiters,
+                  int use_centroid_velocities,
+                  double alpha_balance) {
 
   int k, k3, i; 
@@ -579 +656 @@
       // FIXME: Try with this one precomputed
       for (i=0; i<3; i++) {
-        dz = max(dz, fabs(zv[k3+i]-zc[k]));
+    dz = max(dz, fabs(zv[k3+i]-zc[k]));
       }
     }
@@ -592 +669 @@
 
     //if (hmin < 0.0 ) {
-    //  printf("hmin = %f\n",hmin);
+    //  printf("hmin = %f\n",hmin);
     //}
     
@@ -607 +684 @@
       
       if (dz > 0.0) {
-        alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );      
+    alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );      
       } else {
-        alpha = 1.0;  // Flat bed
+    alpha = 1.0;  // Flat bed
       }
       //printf("Using old style limiter\n");
@@ -622 +699 @@
     
       if (hmin < H0) {
-        alpha = 1.0;
-        for (i=0; i<3; i++) {
-          
-          h_diff = hc_k - hv[i];          
-          if (h_diff <= 0) {
-            // Deep water triangle is further away from bed than 
-            // shallow water (hbar < h). Any alpha will do
-          
-          } else {  
-            // Denominator is positive which means that we need some of the 
-            // h-limited stage.
-            
-            alpha = min(alpha, (hc_k - H0)/h_diff);         
-          }
-        }
-
-        // Ensure alpha in [0,1]
-        if (alpha>1.0) alpha=1.0;
-        if (alpha<0.0) alpha=0.0;
-        
+    alpha = 1.0;
+    for (i=0; i<3; i++) {
+      
+      h_diff = hc_k - hv[i];      
+      if (h_diff <= 0) {
+        // Deep water triangle is further away from bed than 
+        // shallow water (hbar < h). Any alpha will do
+      
+      } else {  
+        // Denominator is positive which means that we need some of the 
+        // h-limited stage.
+        
+        alpha = min(alpha, (hc_k - H0)/h_diff);     
+      }
+    }
+
+    // Ensure alpha in [0,1]
+    if (alpha>1.0) alpha=1.0;
+    if (alpha<0.0) alpha=0.0;
+    
       } else {
-        // Use w-limited stage exclusively in deeper water.
-        alpha = 1.0;       
+    // Use w-limited stage exclusively in deeper water.
+    alpha = 1.0;       
       }
     }
 
 
-    //  Let
+    //  Let
     //
-    //    wvi be the w-limited stage (wvi = zvi + hvi)
-    //    wvi- be the h-limited state (wvi- = zvi + hvi-)
+    //    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
+    //  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)
+    //    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
+    //   Momentum is balanced between constant and limited
 
 
@@ -670 +747 @@
       for (i=0; i<3; i++) {
       
-        wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];     
-
-        // Update momentum at vertices
-        if (use_centroid_velocities == 1) {
-          // This is a simple, efficient and robust option
-          // It uses first order approximation of velocities, but retains
-          // the order used by stage.
-        
-          // Speeds at centroids
-          if (hc_k > epsilon) {
-            uc = xmomc[k]/hc_k;
-            vc = ymomc[k]/hc_k;
-          } else {
-            uc = 0.0;
-            vc = 0.0;
-          }
-          
-          // Vertex momenta guaranteed to be consistent with depth guaranteeing
-          // controlled speed
-          hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
-          xmomv[k3+i] = uc*hv[i];
-          ymomv[k3+i] = vc*hv[i];
-          
-        } else {
-          // Update momentum as a linear combination of
-          // xmomc and ymomc (shallow) and momentum
-          // from extrapolator xmomv and ymomv (deep).
-          // This assumes that values from xmomv and ymomv have
-          // been established e.g. by the gradient limiter.
-
-          // FIXME (Ole): I think this should be used with vertex momenta
-          // computed above using centroid_velocities instead of xmomc 
-          // and ymomc as they'll be more representative first order
-          // values.
-          
-          xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
-          ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
-        
-        }
+    wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i]; 
+
+    // Update momentum at vertices
+    if (use_centroid_velocities == 1) {
+      // This is a simple, efficient and robust option
+      // It uses first order approximation of velocities, but retains
+      // the order used by stage.
+    
+      // Speeds at centroids
+      if (hc_k > epsilon) {
+        uc = xmomc[k]/hc_k;
+        vc = ymomc[k]/hc_k;
+      } else {
+        uc = 0.0;
+        vc = 0.0;
+      }
+      
+      // Vertex momenta guaranteed to be consistent with depth guaranteeing
+      // controlled speed
+      hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
+      xmomv[k3+i] = uc*hv[i];
+      ymomv[k3+i] = vc*hv[i];
+      
+    } else {
+      // Update momentum as a linear combination of
+      // xmomc and ymomc (shallow) and momentum
+      // from extrapolator xmomv and ymomv (deep).
+      // This assumes that values from xmomv and ymomv have
+      // been established e.g. by the gradient limiter.
+
+      // FIXME (Ole): I think this should be used with vertex momenta
+      // computed above using centroid_velocities instead of xmomc 
+      // and ymomc as they'll be more representative first order
+      // values.
+      
+      xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
+      ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
+    
+    }
       }
     }
@@ -719 +796 @@
 
 int _protect(int N,
-             double minimum_allowed_height,
-             double maximum_allowed_speed,
-             double epsilon,
-             double* wc,
-             double* zc,
-             double* xmomc,
-             double* ymomc) {
+         double minimum_allowed_height,
+         double maximum_allowed_speed,
+         double epsilon,
+         double* wc,
+         double* zc,
+         double* xmomc,
+         double* ymomc) {
 
   int k;
@@ -738 +815 @@
 
       if (hc < minimum_allowed_height) {
-        
-        // Set momentum to zero and ensure h is non negative
-        xmomc[k] = 0.0;
-        ymomc[k] = 0.0;
-        if (hc <= 0.0) wc[k] = zc[k];
+        
+    // Set momentum to zero and ensure h is non negative
+    xmomc[k] = 0.0;
+    ymomc[k] = 0.0;
+    if (hc <= 0.0) wc[k] = zc[k];
       }
     }
@@ -757 +834 @@
               
         if (hc <= 0.0) {
-                wc[k] = zc[k];
-        xmomc[k] = 0.0;
-        ymomc[k] = 0.0;
+            wc[k] = zc[k];
+        xmomc[k] = 0.0;
+        ymomc[k] = 0.0;
         } else {
           //Reduce excessive speeds derived from division by small hc
-        //FIXME (Ole): This may be unnecessary with new slope limiters 
-        //in effect.
+        //FIXME (Ole): This may be unnecessary with new slope limiters 
+        //in effect.
           
           u = xmomc[k]/hc;
-          if (fabs(u) > maximum_allowed_speed) {
-            reduced_speed = maximum_allowed_speed * u/fabs(u);
-            //printf("Speed (u) has been reduced from %.3f to %.3f\n",
-            //   u, reduced_speed);
-            xmomc[k] = reduced_speed * hc;
-          }
+      if (fabs(u) > maximum_allowed_speed) {
+        reduced_speed = maximum_allowed_speed * u/fabs(u);
+        //printf("Speed (u) has been reduced from %.3f to %.3f\n",
+        //   u, reduced_speed);
+        xmomc[k] = reduced_speed * hc;
+      }
 
           v = ymomc[k]/hc;
-          if (fabs(v) > maximum_allowed_speed) {
-            reduced_speed = maximum_allowed_speed * v/fabs(v);
-            //printf("Speed (v) has been reduced from %.3f to %.3f\n",
-            //   v, reduced_speed);
-            ymomc[k] = reduced_speed * hc;
-          }
+      if (fabs(v) > maximum_allowed_speed) {
+        reduced_speed = maximum_allowed_speed * v/fabs(v);
+        //printf("Speed (v) has been reduced from %.3f to %.3f\n",
+        //   v, reduced_speed);
+        ymomc[k] = reduced_speed * hc;
+      }
         }
       }
@@ -790 +867 @@
 
 int _assign_wind_field_values(int N,
-                              double* xmom_update,
-                              double* ymom_update,
-                              double* s_vec,
-                              double* phi_vec,
-                              double cw) {
+                  double* xmom_update,
+                  double* ymom_update,
+                  double* s_vec,
+                  double* phi_vec,
+                  double cw) {
 
   // Assign windfield values to momentum updates
@@ -839 +916 @@
 
   if (!PyArg_ParseTuple(args, "OOOddOddd",
-                        &normal, &ql, &qr, &zl, &zr, &edgeflux,
-                        &epsilon, &g, &H0)) {
+            &normal, &ql, &qr, &zl, &zr, &edgeflux,
+            &epsilon, &g, &H0)) {
     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: flux_function_central could not parse input arguments");
     return NULL;
@@ -847 +924 @@
   
   _flux_function_central((double*) ql -> data, 
-                         (double*) qr -> data, 
-                         zl, 
-                         zr,                                             
-                         ((double*) normal -> data)[0],
-                         ((double*) normal -> data)[1],                  
-                         epsilon, H0, g,
-                         (double*) edgeflux -> data, 
-                         &max_speed);
+             (double*) qr -> data, 
+             zl, 
+             zr,                         
+             ((double*) normal -> data)[0],
+             ((double*) normal -> data)[1],          
+             epsilon, H0, g,
+             (double*) edgeflux -> data, 
+             &max_speed);
   
   return Py_BuildValue("d", max_speed);  
@@ -874 +951 @@
 
   if (!PyArg_ParseTuple(args, "dOOOOO",
-                        &g, &h, &v, &x,
-                        &xmom, &ymom)) {
+            &g, &h, &v, &x,
+            &xmom, &ymom)) {
     //&epsilon)) {
     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
@@ -940 +1017 @@
 
   if (!PyArg_ParseTuple(args, "ddOOOOOOO",
-                        &g, &eps, &w, &z, &uh, &vh, &eta,
-                        &xmom, &ymom)) {
+            &g, &eps, &w, &z, &uh, &vh, &eta,
+            &xmom, &ymom)) {
     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction could not parse input arguments");
     return NULL;
@@ -957 +1034 @@
   N = w -> dimensions[0];
   _manning_friction(g, eps, N,
-                    (double*) w -> data,
-                    (double*) z -> data,
-                    (double*) uh -> data,
-                    (double*) vh -> data,
-                    (double*) eta -> data,
-                    (double*) xmom -> data,
-                    (double*) ymom -> data);
+            (double*) w -> data,
+            (double*) z -> data,
+            (double*) uh -> data,
+            (double*) vh -> data,
+            (double*) eta -> data,
+            (double*) xmom -> data,
+            (double*) ymom -> data);
 
   return Py_BuildValue("");
@@ -981 +1058 @@
 
   if (!PyArg_ParseTuple(args, "ddOOOOOOO",
-                        &g, &eps, &w, &z, &uh, &vh, &eta,
-                        &xmom, &ymom))
+            &g, &eps, &w, &z, &uh, &vh, &eta,
+            &xmom, &ymom))
     return NULL;
 
@@ -996 +1073 @@
   N = w -> dimensions[0];
   _manning_friction_explicit(g, eps, N,
-                    (double*) w -> data,
-                    (double*) z -> data,
-                    (double*) uh -> data,
-                    (double*) vh -> data,
-                    (double*) eta -> data,
-                    (double*) xmom -> data,
-                    (double*) ymom -> data);
+            (double*) w -> data,
+            (double*) z -> data,
+            (double*) uh -> data,
+            (double*) vh -> data,
+            (double*) eta -> data,
+            (double*) xmom -> data,
+            (double*) ymom -> data);
 
   return Py_BuildValue("");
@@ -1012 +1089 @@
 // Computational routine
 int _extrapolate_second_order_sw(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) {
-                                  
-                                  
+                                 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) {
+                  
+                  
 
   // 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;
+  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;
 
-        
-  for (k=0; k<number_of_elements; k++) {
+    
+  for (k = 0; k < number_of_elements; k++) 
+  {
     k3=k*3;
     k6=k*6;
 
-    
-    if (number_of_boundaries[k]==3){
+    if (number_of_boundaries[k]==3)
+    {
       // No neighbours, set gradient on the triangle to zero
       
@@ -1065 +1143 @@
       continue;
     }
-    else {
+    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];
+      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];
+      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;
+      dxv0 = xv0 - x; 
+      dxv1 = xv1 - x; 
+      dxv2 = xv2 - x;
+      dyv0 = yv0 - y; 
+      dyv1 = yv1 - y; 
+      dyv2 = yv2 - y;
     }
 
 
             
-    if (number_of_boundaries[k]<=1){
-    
+    if (number_of_boundaries[k]<=1)
+    {
       //==============================================
       // Number of boundaries <= 1
@@ -1099 +1185 @@
       // from this centroid and its two neighbours
       
-      k0=surrogate_neighbours[k3];
-      k1=surrogate_neighbours[k3+1];
-      k2=surrogate_neighbours[k3+2];
+      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];
+      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;
+      dx1 = x1 - x0; 
+      dx2 = x2 - x0;
+      dy1 = y1 - y0; 
+      dy2 = y2 - y0;
       
       // Calculate 2*area of the auxiliary triangle
@@ -1128 +1216 @@
       // If the mesh is 'weird' near the boundary, 
       // the triangle might be flat or clockwise:
-      if (area2<=0) {
-        PyErr_SetString(PyExc_RuntimeError, 
-                        "shallow_water_ext.c: negative triangle area encountered");
-        return -1;
+      if (area2 <= 0) 
+      {
+        PyErr_SetString(PyExc_RuntimeError, 
+                        "shallow_water_ext.c: negative triangle area encountered");
+    
+        return -1;
       }  
       
@@ -1139 +1229 @@
       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);
+      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;
-        }
-      }
-
-            
+      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
@@ -1164 +1256 @@
       // 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];
+      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];
-      
+      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 /= area2;
+      a *= inv_area2;
       b = dx1*dq2 - dx2*dq1;
-      b /= area2;
+      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;
+      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);
+      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
       
       // Playing with dry wet interface
@@ -1193 +1286 @@
       //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);
@@ -1199 +1292 @@
       //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+i]=stage_centroid_values[k]+dqv[i];
+      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];
       
       
@@ -1211 +1306 @@
       // 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];
+      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];
+      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 /= area2;
+      a *= inv_area2;
       b = dx1*dq2 - dx2*dq1;
-      b /= area2;
+      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;
+      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);
+      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
       //beta_tmp = beta_uh;
       //if (hmin<minimum_allowed_height)
@@ -1240 +1335 @@
 
       // 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];
-      
+      limit_gradient(dqv, qmin, qmax, beta_tmp);
+
+      for (i=0; i < 3; i++)
+      {
+        xmom_vertex_values[k3+i] = xmom_centroid_values[k] + dqv[i];
+      }
       
       //-----------------------------------
@@ -1252 +1348 @@
       // 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];
+      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];
+      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 /= area2;
+      a *= inv_area2;
       b = dx1*dq2 - dx2*dq1;
-      b /= area2;
+      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;
+      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);
+      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
       
       //beta_tmp = beta_vh;
@@ -1280 +1376 @@
       //if (hmin<minimum_allowed_height)
       //beta_tmp = beta_vh_dry;
-      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;       
+      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
 
       // Limit the gradient
-      limit_gradient(dqv,qmin,qmax,beta_tmp);
+      limit_gradient(dqv, qmin, qmax, beta_tmp);
       
       for (i=0;i<3;i++)
-        ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
-        
+      {
+        ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
+      }
     } // End number_of_boundaries <=1 
-    else{
+    else
+    {
 
       //==============================================
@@ -1298 +1396 @@
       
       // 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
-        PyErr_SetString(PyExc_RuntimeError, 
-                        "shallow_water_ext.c: Internal neighbour not found");      
-        return -1;
-      }
-      
-      k1=surrogate_neighbours[k2];
+      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
+        PyErr_SetString(PyExc_RuntimeError, 
+                        "shallow_water_ext.c: 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];
+      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;
+      dx1 = x1 - x; 
+      dy1 = y1 - y;
       
       // Set area2 as the square of the distance
-      area2=dx1*dx1+dy1*dy1;
+      area2 = dx1*dx1 + dy1*dy1;
       
       // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) 
@@ -1332 +1435 @@
       // respectively correspond to the x- and y- gradients 
       // of the conserved quantities
-      dx2=1.0/area2;
-      dy2=dx2*dy1;
-      dx2*=dx1;
+      dx2 = 1.0/area2;
+      dy2 = dx2*dy1;
+      dx2 *= dx1;
       
       
@@ -1342 +1445 @@
 
       // Compute differentials
-      dq1=stage_centroid_values[k1]-stage_centroid_values[k];
+      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;
+      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;
+      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;
+      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+i]=stage_centroid_values[k]+dqv[i];
-      
+      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];
+      //}
       
       //-----------------------------------
@@ -1376 +1484 @@
       
       // Compute differentials
-      dq1=xmom_centroid_values[k1]-xmom_centroid_values[k];
+      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;
+      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;
+      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;
+      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+i]=xmom_centroid_values[k]+dqv[i];
-      
+      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];
+      }
       
       //-----------------------------------
@@ -1410 +1526 @@
 
       // Compute differentials
-      dq1=ymom_centroid_values[k1]-ymom_centroid_values[k];
+      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;
+      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;
+      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;
+      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+i]=ymom_centroid_values[k]+dqv[i];
-        
+      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
   
   return 0;
-}                       
+}           
 
 
@@ -1477 +1604 @@
     Post conditions:
             The vertices of each triangle have values from a 
-            limited linear reconstruction
-            based on centroid values
+        limited linear reconstruction
+        based on centroid values
 
   */
@@ -1505 +1632 @@
   // Convert Python arguments to C
   if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
-                        &domain,
-                        &surrogate_neighbours,
-                        &number_of_boundaries,
-                        &centroid_coordinates,
-                        &stage_centroid_values,
-                        &xmom_centroid_values,
-                        &ymom_centroid_values,
-                        &elevation_centroid_values,
-                        &vertex_coordinates,
-                        &stage_vertex_values,
-                        &xmom_vertex_values,
-                        &ymom_vertex_values,
-                        &elevation_vertex_values,
-                        &optimise_dry_cells)) {                 
-                        
+            &domain,
+            &surrogate_neighbours,
+            &number_of_boundaries,
+            &centroid_coordinates,
+            &stage_centroid_values,
+            &xmom_centroid_values,
+            &ymom_centroid_values,
+            &elevation_centroid_values,
+            &vertex_coordinates,
+            &stage_vertex_values,
+            &xmom_vertex_values,
+            &ymom_vertex_values,
+            &elevation_vertex_values,
+            &optimise_dry_cells)) {         
+            
     PyErr_SetString(PyExc_RuntimeError, 
-                    "Input arguments to extrapolate_second_order_sw failed");
+            "Input arguments to extrapolate_second_order_sw failed");
     return NULL;
   }
@@ -1557 +1684 @@
   // Call underlying computational routine
   e = _extrapolate_second_order_sw(number_of_elements,
-                                   epsilon,
-                                   minimum_allowed_height,
-                                   beta_w,
-                                   beta_w_dry,
-                                   beta_uh,
-                                   beta_uh_dry,
-                                   beta_vh,
-                                   beta_vh_dry,
-                                   (long*) surrogate_neighbours -> data,
-                                   (long*) number_of_boundaries -> data,
-                                   (double*) centroid_coordinates -> data,
-                                   (double*) stage_centroid_values -> data,
-                                   (double*) xmom_centroid_values -> data,
-                                   (double*) ymom_centroid_values -> data,
-                                   (double*) elevation_centroid_values -> data,
-                                   (double*) vertex_coordinates -> data,
-                                   (double*) stage_vertex_values -> data,
-                                   (double*) xmom_vertex_values -> data,
-                                   (double*) ymom_vertex_values -> data,
-                                   (double*) elevation_vertex_values -> data,
-                                   optimise_dry_cells);
+                   epsilon,
+                   minimum_allowed_height,
+                   beta_w,
+                   beta_w_dry,
+                   beta_uh,
+                   beta_uh_dry,
+                   beta_vh,
+                   beta_vh_dry,
+                   (long*) surrogate_neighbours -> data,
+                   (long*) number_of_boundaries -> data,
+                   (double*) centroid_coordinates -> data,
+                   (double*) stage_centroid_values -> data,
+                   (double*) xmom_centroid_values -> data,
+                   (double*) ymom_centroid_values -> data,
+                   (double*) elevation_centroid_values -> data,
+                   (double*) vertex_coordinates -> data,
+                   (double*) stage_vertex_values -> data,
+                   (double*) xmom_vertex_values -> data,
+                   (double*) ymom_vertex_values -> data,
+                   (double*) elevation_vertex_values -> data,
+                   optimise_dry_cells);
   if (e == -1) {
     // Use error string set inside computational routine
     return NULL;
-  }                               
+  }                   
   
   
@@ -1609 +1736 @@
   // Convert Python arguments to C
   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
-                                   &Q, &Normal, &direction)) {
+                   &Q, &Normal, &direction)) {
     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
     return NULL;
@@ -1654 +1781 @@
 // Computational function for flux computation
 double _compute_fluxes_central(int number_of_elements,
-                               double timestep,
-                               double epsilon,
-                               double H0,
-                               double g,
-                               long* neighbours,
-                               long* neighbour_edges,
-                               double* normals,
-                               double* edgelengths, 
-                               double* radii, 
-                               double* areas,
-                               long* tri_full_flag,
-                               double* stage_edge_values,
-                               double* xmom_edge_values,
-                               double* ymom_edge_values,
-                               double* bed_edge_values,
-                               double* stage_boundary_values,
-                               double* xmom_boundary_values,
-                               double* ymom_boundary_values,
-                               double* stage_explicit_update,
-                               double* xmom_explicit_update,
-                               double* ymom_explicit_update,
-                               long* already_computed_flux,
-                               double* max_speed_array,
-                               int optimise_dry_cells) {
-                               
+                               double timestep,
+                               double epsilon,
+                               double H0,
+                               double g,
+                               long* neighbours,
+                               long* neighbour_edges,
+                               double* normals,
+                               double* edgelengths, 
+                               double* radii, 
+                               double* areas,
+                               long* tri_full_flag,
+                               double* stage_edge_values,
+                               double* xmom_edge_values,
+                               double* ymom_edge_values,
+                               double* bed_edge_values,
+                               double* stage_boundary_values,
+                               double* xmom_boundary_values,
+                               double* ymom_boundary_values,
+                               double* stage_explicit_update,
+                               double* xmom_explicit_update,
+                               double* ymom_explicit_update,
+                               long* already_computed_flux,
+                               double* max_speed_array,
+                               int optimise_dry_cells) 
+{                   
   // Local variables
-  double max_speed, length, area, zl, zr;
+  double max_speed, length, inv_area, zl, zr;
   int k, i, m, n;
   int ki, nm=0, ki2; // Index shorthands
@@ -1687 +1814 @@
   double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
 
-  static long call=1; // Static local variable flagging already computed flux
-                               
-        
+  static long call = 1; // Static local variable flagging already computed flux
+                   
   // Start computation
   call++; // Flag 'id' of flux calculation for this timestep 
-  
+
   // Set explicit_update to zero for all conserved_quantities.
   // This assumes compute_fluxes called before forcing terms
-  for (k=0; k<number_of_elements; k++) {
-    stage_explicit_update[k]=0.0;
-    xmom_explicit_update[k]=0.0;
-    ymom_explicit_update[k]=0.0;  
-  }
-
+  memset((char*) stage_explicit_update, 0, number_of_elements*sizeof(double));
+  memset((char*) xmom_explicit_update, 0, number_of_elements*sizeof(double));
+  memset((char*) ymom_explicit_update, 0, number_of_elements*sizeof(double));    
+  
   // For all triangles
-  for (k=0; k<number_of_elements; k++) {
-    
+  for (k = 0; k < number_of_elements; k++) 
+  {  
     // Loop through neighbours and compute edge flux for each  
-    for (i=0; i<3; i++) {
-      ki = k*3+i; // Linear index (triangle k, edge i)
+    for (i = 0; i < 3; i++) 
+    {
+      ki = k*3 + i; // Linear index to edge i of triangle k
       
       if (already_computed_flux[ki] == call)
+      {
         // We've already computed the flux across this edge
-        continue;
-        
-        
+        continue;
+      }
+    
+      // Get left hand side values from triangle k, edge i 
       ql[0] = stage_edge_values[ki];
       ql[1] = xmom_edge_values[ki];
@@ -1718 +1845 @@
       zl = bed_edge_values[ki];
 
-      // Quantities at neighbour on nearest face
+      // Get right hand side values either from neighbouring triangle
+      // or from boundary array (Quantities at neighbour on nearest face).
       n = neighbours[ki];
-      if (n < 0) {
+      if (n < 0) 
+      {
         // Neighbour is a boundary condition
-        m = -n-1; // Convert negative flag to boundary index
-        
-        qr[0] = stage_boundary_values[m];
-        qr[1] = xmom_boundary_values[m];
-        qr[2] = ymom_boundary_values[m];
-        zr = zl; // Extend bed elevation to boundary
-      } else {
-        // Neighbour is a real element
-        m = neighbour_edges[ki];
-        nm = n*3+m; // Linear index (triangle n, edge m)
-        
-        qr[0] = stage_edge_values[nm];
-        qr[1] = xmom_edge_values[nm];
-        qr[2] = ymom_edge_values[nm];
-        zr = bed_edge_values[nm];
-      }
-      
-      
-      if (optimise_dry_cells) {      
-        // Check if flux calculation is necessary across this edge
-        // This check will exclude dry cells.
-        // This will also optimise cases where zl != zr as 
-        // long as both are dry
-
-        if ( fabs(ql[0] - zl) < epsilon && 
-             fabs(qr[0] - zr) < epsilon ) {
-          // Cell boundary is dry
-          
-          already_computed_flux[ki] = call; // #k Done  
-          if (n>=0)
-            already_computed_flux[nm] = call; // #n Done
-        
-          max_speed = 0.0;
-          continue;
-        }
-      }
-    
+        m = -n - 1; // Convert negative flag to boundary index
+    
+        qr[0] = stage_boundary_values[m];
+        qr[1] = xmom_boundary_values[m];
+        qr[2] = ymom_boundary_values[m];
+        zr = zl; // Extend bed elevation to boundary
+      } 
+      else 
+      {
+        // Neighbour is a real triangle
+        m = neighbour_edges[ki];
+        nm = n*3 + m; // Linear index (triangle n, edge m)
+        
+        qr[0] = stage_edge_values[nm];
+        qr[1] = xmom_edge_values[nm];
+        qr[2] = ymom_edge_values[nm];
+        zr = bed_edge_values[nm];
+      }
+      
+      // Now we have values for this edge - both from left and right side.
+      
+      if (optimise_dry_cells) 
+      {     
+        // Check if flux calculation is necessary across this edge
+        // This check will exclude dry cells.
+        // This will also optimise cases where zl != zr as 
+        // long as both are dry
+
+        if (fabs(ql[0] - zl) < epsilon && 
+            fabs(qr[0] - zr) < epsilon) 
+        {
+          // Cell boundary is dry
+          
+          already_computed_flux[ki] = call; // #k Done  
+          if (n >= 0)
+          {
+            already_computed_flux[nm] = call; // #n Done
+          }
+        
+          max_speed = 0.0;
+          continue;
+        }
+      }
       
       // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
@@ -1765 +1900 @@
       // Edge flux computation (triangle k, edge i)
       _flux_function_central(ql, qr, zl, zr,
-                             normals[ki2], normals[ki2+1],
-                             epsilon, H0, g,
-                             edgeflux, &max_speed);
+                             normals[ki2], normals[ki2+1],
+                             epsilon, H0, g,
+                             edgeflux, &max_speed);
       
       
@@ -1786 +1921 @@
       
       // Update neighbour n with same flux but reversed sign
-      if (n>=0) {
-        stage_explicit_update[n] += edgeflux[0];
-        xmom_explicit_update[n] += edgeflux[1];
-        ymom_explicit_update[n] += edgeflux[2];
-        
-        already_computed_flux[nm] = call; // #n Done
-      }
-
+      if (n >= 0) 
+      {
+        stage_explicit_update[n] += edgeflux[0];
+        xmom_explicit_update[n] += edgeflux[1];
+        ymom_explicit_update[n] += edgeflux[2];
+        
+        already_computed_flux[nm] = call; // #n Done
+      }
 
       // Update timestep based on edge i and possibly neighbour n
-      if (tri_full_flag[k] == 1) {
-        if (max_speed > epsilon) {
-        
-          // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
-          
-          // CFL for triangle k
-          timestep = min(timestep, radii[k]/max_speed);
-          
-          if (n>=0) 
-            // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
-            timestep = min(timestep, radii[n]/max_speed);
-          
-          // Ted Rigby's suggested less conservative version
-          //if (n>=0) {             
-          //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
-          //} else {
-          //  timestep = min(timestep, radii[k]/max_speed);
-          // }
-        }
+      if (tri_full_flag[k] == 1) 
+      {
+        if (max_speed > epsilon) 
+        {
+          // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
+          
+          // CFL for triangle k
+          timestep = min(timestep, radii[k]/max_speed);
+          
+          if (n >= 0)
+          {
+            // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
+            timestep = min(timestep, radii[n]/max_speed);
+          }
+
+          // Ted Rigby's suggested less conservative version
+          //if (n>=0) {         
+          //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
+          //} else {
+          //  timestep = min(timestep, radii[k]/max_speed);
+          // }
+        }
       }
       
@@ -1822 +1960 @@
     // Normalise triangle k by area and store for when all conserved
     // quantities get updated
-    area = areas[k];
-    stage_explicit_update[k] /= area;
-    xmom_explicit_update[k] /= area;
-    ymom_explicit_update[k] /= area;
+    inv_area = 1.0/areas[k];
+    stage_explicit_update[k] *= inv_area;
+    xmom_explicit_update[k] *= inv_area;
+    ymom_explicit_update[k] *= inv_area;
     
    
@@ -1832 +1970 @@
     
   } // End triangle k
-        
-                               
-                               
+             
   return timestep;
 }
@@ -1869 +2005 @@
 
     PyObject 
-        *domain,
-        *stage, 
-        *xmom, 
-        *ymom, 
-        *bed;
+    *domain,
+    *stage, 
+    *xmom, 
+    *ymom, 
+    *bed;
 
     PyArrayObject 
-        *neighbours, 
-        *neighbour_edges,
-        *normals, 
-        *edgelengths, 
-        *radii, 
-        *areas,
-        *tri_full_flag,
-        *stage_edge_values,
-        *xmom_edge_values,
-        *ymom_edge_values,
-        *bed_edge_values,
-        *stage_boundary_values,
-        *xmom_boundary_values,
-        *ymom_boundary_values,
-        *stage_explicit_update,
-        *xmom_explicit_update,
-        *ymom_explicit_update,
-        *already_computed_flux, //Tracks whether the flux across an edge has already been computed
-        *max_speed_array; //Keeps track of max speeds for each triangle
+    *neighbours, 
+    *neighbour_edges,
+    *normals, 
+    *edgelengths, 
+    *radii, 
+    *areas,
+    *tri_full_flag,
+    *stage_edge_values,
+    *xmom_edge_values,
+    *ymom_edge_values,
+    *bed_edge_values,
+    *stage_boundary_values,
+    *xmom_boundary_values,
+    *ymom_boundary_values,
+    *stage_explicit_update,
+    *xmom_explicit_update,
+    *ymom_explicit_update,
+    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
+    *max_speed_array; //Keeps track of max speeds for each triangle
 
     
@@ -1902 +2038 @@
     // Convert Python arguments to C
     if (!PyArg_ParseTuple(args, "dOOOO", &timestep, &domain, &stage, &xmom, &ymom, &bed )) {
-        PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
-        return NULL;
+    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
+    return NULL;
     }
 
@@ -1940 +2076 @@
   // the explicit update arrays 
   timestep = _compute_fluxes_central(number_of_elements,
-                                     timestep,
-                                     epsilon,
-                                     H0,
-                                     g,
-                                     (long*) neighbours -> data,
-                                     (long*) neighbour_edges -> data,
-                                     (double*) normals -> data,
-                                     (double*) edgelengths -> data, 
-                                     (double*) radii -> data, 
-                                     (double*) areas -> data,
-                                     (long*) tri_full_flag -> data,
-                                     (double*) stage_edge_values -> data,
-                                     (double*) xmom_edge_values -> data,
-                                     (double*) ymom_edge_values -> data,
-                                     (double*) bed_edge_values -> data,
-                                     (double*) stage_boundary_values -> data,
-                                     (double*) xmom_boundary_values -> data,
-                                     (double*) ymom_boundary_values -> data,
-                                     (double*) stage_explicit_update -> data,
-                                     (double*) xmom_explicit_update -> data,
-                                     (double*) ymom_explicit_update -> data,
-                                     (long*) already_computed_flux -> data,
-                                     (double*) max_speed_array -> data,
-                                     optimise_dry_cells);
+                     timestep,
+                     epsilon,
+                     H0,
+                     g,
+                     (long*) neighbours -> data,
+                     (long*) neighbour_edges -> data,
+                     (double*) normals -> data,
+                     (double*) edgelengths -> data, 
+                     (double*) radii -> data, 
+                     (double*) areas -> data,
+                     (long*) tri_full_flag -> data,
+                     (double*) stage_edge_values -> data,
+                     (double*) xmom_edge_values -> data,
+                     (double*) ymom_edge_values -> data,
+                     (double*) bed_edge_values -> data,
+                     (double*) stage_boundary_values -> data,
+                     (double*) xmom_boundary_values -> data,
+                     (double*) ymom_boundary_values -> data,
+                     (double*) stage_explicit_update -> data,
+                     (double*) xmom_explicit_update -> data,
+                     (double*) ymom_explicit_update -> data,
+                     (long*) already_computed_flux -> data,
+                     (double*) max_speed_array -> data,
+                     optimise_dry_cells);
 
   Py_DECREF(neighbours);
@@ -2013 +2149 @@
     domain.timestep = compute_fluxes(timestep,
                                      domain.epsilon,
-                                     domain.H0,
+                     domain.H0,
                                      domain.g,
                                      domain.neighbours,
@@ -2033 +2169 @@
                                      Ymom.explicit_update,
                                      already_computed_flux,
-                                     optimise_dry_cells)                                     
+                     optimise_dry_cells)                     
 
 
@@ -2066 +2202 @@
   // Convert Python arguments to C
   if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
-                        &timestep,
-                        &epsilon,
-                        &H0,
-                        &g,
-                        &neighbours,
-                        &neighbour_edges,
-                        &normals,
-                        &edgelengths, &radii, &areas,
-                        &tri_full_flag,
-                        &stage_edge_values,
-                        &xmom_edge_values,
-                        &ymom_edge_values,
-                        &bed_edge_values,
-                        &stage_boundary_values,
-                        &xmom_boundary_values,
-                        &ymom_boundary_values,
-                        &stage_explicit_update,
-                        &xmom_explicit_update,
-                        &ymom_explicit_update,
-                        &already_computed_flux,
-                        &max_speed_array,
-                        &optimise_dry_cells)) {
+            &timestep,
+            &epsilon,
+            &H0,
+            &g,
+            &neighbours,
+            &neighbour_edges,
+            &normals,
+            &edgelengths, &radii, &areas,
+            &tri_full_flag,
+            &stage_edge_values,
+            &xmom_edge_values,
+            &ymom_edge_values,
+            &bed_edge_values,
+            &stage_boundary_values,
+            &xmom_boundary_values,
+            &ymom_boundary_values,
+            &stage_explicit_update,
+            &xmom_explicit_update,
+            &ymom_explicit_update,
+            &already_computed_flux,
+            &max_speed_array,
+            &optimise_dry_cells)) {
     PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     return NULL;
@@ -2118 +2254 @@
   // the explicit update arrays 
   timestep = _compute_fluxes_central(number_of_elements,
-                                     timestep,
-                                     epsilon,
-                                     H0,
-                                     g,
-                                     (long*) neighbours -> data,
-                                     (long*) neighbour_edges -> data,
-                                     (double*) normals -> data,
-                                     (double*) edgelengths -> data, 
-                                     (double*) radii -> data, 
-                                     (double*) areas -> data,
-                                     (long*) tri_full_flag -> data,
-                                     (double*) stage_edge_values -> data,
-                                     (double*) xmom_edge_values -> data,
-                                     (double*) ymom_edge_values -> data,
-                                     (double*) bed_edge_values -> data,
-                                     (double*) stage_boundary_values -> data,
-                                     (double*) xmom_boundary_values -> data,
-                                     (double*) ymom_boundary_values -> data,
-                                     (double*) stage_explicit_update -> data,
-                                     (double*) xmom_explicit_update -> data,
-                                     (double*) ymom_explicit_update -> data,
-                                     (long*) already_computed_flux -> data,
-                                     (double*) max_speed_array -> data,
-                                     optimise_dry_cells);
+                     timestep,
+                     epsilon,
+                     H0,
+                     g,
+                     (long*) neighbours -> data,
+                     (long*) neighbour_edges -> data,
+                     (double*) normals -> data,
+                     (double*) edgelengths -> data, 
+                     (double*) radii -> data, 
+                     (double*) areas -> data,
+                     (long*) tri_full_flag -> data,
+                     (double*) stage_edge_values -> data,
+                     (double*) xmom_edge_values -> data,
+                     (double*) ymom_edge_values -> data,
+                     (double*) bed_edge_values -> data,
+                     (double*) stage_boundary_values -> data,
+                     (double*) xmom_boundary_values -> data,
+                     (double*) ymom_boundary_values -> data,
+                     (double*) stage_explicit_update -> data,
+                     (double*) xmom_explicit_update -> data,
+                     (double*) ymom_explicit_update -> data,
+                     (long*) already_computed_flux -> data,
+                     (double*) max_speed_array -> data,
+                     optimise_dry_cells);
   
   // Return updated flux timestep
@@ -2228 +2364 @@
   // Convert Python arguments to C
   if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
-                        &timestep,
-                        &epsilon,
-                        &H0,
-                        &g,
-                        &neighbours,
-                        &neighbour_edges,
-                        &normals,
-                        &edgelengths, &radii, &areas,
-                        &tri_full_flag,
-                        &stage_edge_values,
-                        &xmom_edge_values,
-                        &ymom_edge_values,
-                        &bed_edge_values,
-                        &stage_boundary_values,
-                        &xmom_boundary_values,
-                        &ymom_boundary_values,
-                        &stage_explicit_update,
-                        &xmom_explicit_update,
-                        &ymom_explicit_update,
-                        &already_computed_flux)) {
+            &timestep,
+            &epsilon,
+            &H0,
+            &g,
+            &neighbours,
+            &neighbour_edges,
+            &normals,
+            &edgelengths, &radii, &areas,
+            &tri_full_flag,
+            &stage_edge_values,
+            &xmom_edge_values,
+            &ymom_edge_values,
+            &bed_edge_values,
+            &stage_boundary_values,
+            &xmom_boundary_values,
+            &ymom_boundary_values,
+            &stage_explicit_update,
+            &xmom_explicit_update,
+            &ymom_explicit_update,
+            &already_computed_flux)) {
     PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     return NULL;
@@ -2265 +2401 @@
       ki = k*3+i;
       if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
-        continue;
+    continue;
       ql[0] = ((double *) stage_edge_values -> data)[ki];
       ql[1] = ((double *) xmom_edge_values -> data)[ki];
@@ -2274 +2410 @@
       n = ((long *) neighbours -> data)[ki];
       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];
-        zr = zl; //Extend bed elevation to boundary
+    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];
+    zr = zl; //Extend bed elevation to boundary
       } else {
-        m = ((long *) neighbour_edges -> data)[ki];
-        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];
+    m = ((long *) neighbour_edges -> data)[ki];
+    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];
       }
       // Outward pointing normal vector
@@ -2294 +2430 @@
       //Edge flux computation
       flux_function_kinetic(ql, qr, zl, zr,
-                    normal[0], normal[1],
-                    epsilon, H0, g,
-                    edgeflux, &max_speed);
+            normal[0], normal[1],
+            epsilon, H0, g,
+            edgeflux, &max_speed);
       //update triangle k
       ((long *) already_computed_flux->data)[ki]=call;
@@ -2304 +2440 @@
       //update the neighbour n
       if (n>=0){
-        ((long *) already_computed_flux->data)[nm]=call;
-        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
-        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
-        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
+    ((long *) already_computed_flux->data)[nm]=call;
+    ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
+    ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
+    ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
       }
       ///for (j=0; j<3; j++) {
-        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
-        ///}
-        //Update timestep
-        //timestep = min(timestep, domain.radii[k]/max_speed)
-        //FIXME: SR Add parameter for CFL condition
+    ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
+    ///}
+    //Update timestep
+    //timestep = min(timestep, domain.radii[k]/max_speed)
+    //FIXME: SR Add parameter for CFL condition
     if ( ((long *) tri_full_flag -> data)[k] == 1) {
-            if (max_speed > epsilon) {
-                timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
-                //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
-                if (n>=0)
-                    timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
-            }
+        if (max_speed > epsilon) {
+            timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
+            //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
+            if (n>=0)
+                timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
+        }
     }
     } // end for i
@@ -2350 +2486 @@
   // Convert Python arguments to C
   if (!PyArg_ParseTuple(args, "dddOOOO",
-                        &minimum_allowed_height,
-                        &maximum_allowed_speed,
-                        &epsilon,
-                        &wc, &zc, &xmomc, &ymomc)) {
+            &minimum_allowed_height,
+            &maximum_allowed_speed,
+            &epsilon,
+            &wc, &zc, &xmomc, &ymomc)) {
     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
     return NULL;
@@ -2361 +2497 @@
 
   _protect(N,
-           minimum_allowed_height,
-           maximum_allowed_speed,
-           epsilon,
-           (double*) wc -> data,
-           (double*) zc -> data,
-           (double*) xmomc -> data,
-           (double*) ymomc -> data);
+       minimum_allowed_height,
+       maximum_allowed_speed,
+       epsilon,
+       (double*) wc -> data,
+       (double*) zc -> data,
+       (double*) xmomc -> data,
+       (double*) ymomc -> data);
 
   return Py_BuildValue("");
@@ -2409 +2545 @@
   // Convert Python arguments to C
   if (!PyArg_ParseTuple(args, "OOOOOOOOOO",
-                        &domain,
-                        &wc, &zc, 
-                        &wv, &zv, &hvbar,
-                        &xmomc, &ymomc, &xmomv, &ymomv)) {
+            &domain,
+            &wc, &zc, 
+            &wv, &zv, &hvbar,
+            &xmomc, &ymomc, &xmomv, &ymomv)) {
     PyErr_SetString(PyExc_RuntimeError, 
-                    "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
+            "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
     return NULL;
   }  
-          
-          
+      
+      
   // FIXME (Ole): I tested this without GetAttrString and got time down
   // marginally from 4.0s to 3.8s. Consider passing everything in 
@@ -2463 +2599 @@
   N = wc -> dimensions[0];
   _balance_deep_and_shallow(N,
-                            (double*) wc -> data,
-                            (double*) zc -> data,
-                            (double*) wv -> data,
-                            (double*) zv -> data,
-                            (double*) hvbar -> data,
-                            (double*) xmomc -> data,
-                            (double*) ymomc -> data,
-                            (double*) xmomv -> data,
-                            (double*) ymomv -> data,
-                            H0,
-                            (int) tight_slope_limiters,
-                            (int) use_centroid_velocities,                          
-                            alpha_balance);
+                (double*) wc -> data,
+                (double*) zc -> data,
+                (double*) wv -> data,
+                (double*) zv -> data,
+                (double*) hvbar -> data,
+                (double*) xmomc -> data,
+                (double*) ymomc -> data,
+                (double*) xmomv -> data,
+                (double*) ymomv -> data,
+                H0,
+                (int) tight_slope_limiters,
+                (int) use_centroid_velocities,              
+                alpha_balance);
 
 
@@ -2503 +2639 @@
   // Convert Python arguments to C
   if (!PyArg_ParseTuple(args, "OOOOd",
-                        &xmom_update,
-                        &ymom_update,
-                        &s_vec, &phi_vec,
-                        &cw)) {
+            &xmom_update,
+            &ymom_update,
+            &s_vec, &phi_vec,
+            &cw)) {
     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
     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*) phi_vec -> data,
-           cw);
+       (double*) xmom_update -> data,
+       (double*) ymom_update -> data,
+       (double*) s_vec -> data,
+       (double*) phi_vec -> data,
+       cw);
 
   return Py_BuildValue("");

branches/numpy/anuga/shallow_water/test_data_manager.py

@@ -15 +15 @@
 from struct import pack, unpack
 from sets import ImmutableSet
+import shutil
 
 from anuga.shallow_water import *
@@ -47 +48 @@
         
         self.verbose = Test_Data_Manager.verbose
-        #Create basic mesh
+        # Create basic mesh
         points, vertices, boundary = rectangular(2, 2)
 
-        #Create shallow water domain
+        # Create shallow water domain
         domain = Domain(points, vertices, boundary)
         domain.default_order = 2
 
-        #Set some field values
+        # Set some field values
         domain.set_quantity('elevation', lambda x,y: -x)
         domain.set_quantity('friction', 0.03)
@@ -252 +253 @@
 
         # Get the variables
-        #range = fid.variables['stage_range'][:]
+        range = fid.variables['stage_range'][:]
         #print range
-        #assert num.allclose(range,[-0.93519, 0.15]) or\
-        #       num.allclose(range,[-0.9352743, 0.15]) or\
-        #       num.allclose(range,[-0.93522203, 0.15000001]) # Old slope limiters
-        
-        #range = fid.variables['xmomentum_range'][:]
+        assert num.allclose(range,[-0.93519, 0.15]) or\
+               num.allclose(range,[-0.9352743, 0.15]) or\
+               num.allclose(range,[-0.93522203, 0.15000001]) # Old slope limiters
+        
+        range = fid.variables['xmomentum_range'][:]
         ##print range
-        #assert num.allclose(range,[0,0.4695096]) or \
-        #       num.allclose(range,[0,0.47790655]) or\
-        #       num.allclose(range,[0,0.46941957]) or\
-        #       num.allclose(range,[0,0.47769409])
-
-        
-        #range = fid.variables['ymomentum_range'][:]
+        assert num.allclose(range,[0,0.4695096]) or \
+               num.allclose(range,[0,0.47790655]) or\
+               num.allclose(range,[0,0.46941957]) or\
+               num.allclose(range,[0,0.47769409])
+
+        
+        range = fid.variables['ymomentum_range'][:]
         ##print range
-        #assert num.allclose(range,[0,0.02174380]) or\
-        #       num.allclose(range,[0,0.02174439]) or\
-        #       num.allclose(range,[0,0.02283983]) or\
-        #       num.allclose(range,[0,0.0217342]) or\
-        #       num.allclose(range,[0,0.0227564]) # Old slope limiters 
+        assert num.allclose(range,[0,0.02174380]) or\
+               num.allclose(range,[0,0.02174439]) or\
+               num.allclose(range,[0,0.02283983]) or\
+               num.allclose(range,[0,0.0217342]) or\
+               num.allclose(range,[0,0.0227564]) # Old slope limiters 
         
         fid.close()
@@ -1232 +1233 @@
         from Scientific.IO.NetCDF import NetCDFFile
 
-        #Setup
+        # Setup
         self.domain.set_name('datatest')
 
@@ -1245 +1246 @@
         self.domain.set_quantity('stage', 1.0)
 
-        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+        self.domain.geo_reference = Geo_reference(56, 308500, 6189000)
 
         sww = get_dataobject(self.domain)
@@ -5563 +5564 @@
             quantities_init[0].append(this_ha) # HA
             quantities_init[1].append(this_ua) # UA
-            quantities_init[2].append(this_va) # VA
+            quantities_init[2].append(this_va) # VA 
                 
         file_handle, base_name = tempfile.mkstemp("")
@@ -6104 +6105 @@
                         #print 'writing', time, point_i, q_time[time, point_i]
                         f.write(pack('f', q_time[time, point_i]))
-
             f.close()
 
@@ -6215 +6215 @@
 
         msg='Incorrect gauge depths returned'
-        assert num.allclose(elevation,-depth),msg
+        assert num.allclose(elevation, -depth),msg
         msg='incorrect gauge height time series returned'
-        assert num.allclose(stage,ha)
+        assert num.allclose(stage, ha)
         msg='incorrect gauge ua time series returned'
-        assert num.allclose(xvelocity,ua)
+        assert num.allclose(xvelocity, ua)
         msg='incorrect gauge va time series returned'
         assert num.allclose(yvelocity, -va) # South is positive in MUX
@@ -6397 +6397 @@
                 if j == 2: assert num.allclose(data[i][:parameters_index], -va0[permutation[i], :])
         
+        self.delete_mux(filesI)
 
 
@@ -6407 +6408 @@
         wrong values Win32
 
-        This test does not pass on Windows but test_read_mux_platform_problem1 does
+        This test does not pass on Windows but test_read_mux_platform_problem1
+        does
         """
         
@@ -6422 +6424 @@
         times_ref = num.arange(0, time_step_count*time_step, time_step)
         
-        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115), (-21.,115.), (-22., 117.)]
+        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115),
+                           (-21.,115.), (-22., 117.)]
         n = len(lat_long_points)
         
@@ -6466 +6469 @@
              for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
                  # For timesteps before and after recording range
-                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0                                  
+                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0 
 
 
@@ -6474 +6477 @@
         #print 'va', va1[0,:]
         
-        # Write second mux file to be combined by urs2sts                                             
+        # Write second mux file to be combined by urs2sts 
         base_nameII, filesII = self.write_mux2(lat_long_points,
                                                time_step_count, time_step,
@@ -6605 +6608 @@
 
             f.close()
-
-
-
-
-
-
-        
                                                
         # Create ordering file
         permutation = ensure_numeric([4,0,2])
 
-        _, ordering_filename = tempfile.mkstemp('')
-        order_fid = open(ordering_filename, 'w')  
-        order_fid.write('index, longitude, latitude\n')
-        for index in permutation:
-            order_fid.write('%d, %f, %f\n' %(index, 
-                                             lat_long_points[index][1], 
-                                             lat_long_points[index][0]))
-        order_fid.close()
-        
-        
-
+       #  _, ordering_filename = tempfile.mkstemp('')
+#         order_fid = open(ordering_filename, 'w')  
+#         order_fid.write('index, longitude, latitude\n')
+#         for index in permutation:
+#             order_fid.write('%d, %f, %f\n' %(index, 
+#                                              lat_long_points[index][1], 
+#                                              lat_long_points[index][0]))
+#         order_fid.close()
+        
         # -------------------------------------
         # Now read files back and check values
@@ -6638 +6632 @@
         for j, file in enumerate(filesII):
             # Read stage, u, v enumerated as j
-
-            
             #print 'Reading', j, file
             data = read_mux2(1, [file], weights, file_params, permutation, verbose)
@@ -6645 +6637 @@
             #print 'Data received by Python'
             #print data[1][8]
-
-            
             number_of_selected_stations = data.shape[0]
 
@@ -6659 +6649 @@
                 #print i, parameters_index
                 #print quantity[i][:]
-
-                
                 if j == 0: assert num.allclose(data[i][:parameters_index], ha1[permutation[i], :])
                 if j == 1: assert num.allclose(data[i][:parameters_index], ua1[permutation[i], :])
@@ -6669 +6657 @@
                     #print j, i
                     #print 'Input'
-                    #print 'u', ua1[permutation[i], 8]                                        
+                    #print 'u', ua1[permutation[i], 8]       
                     #print 'v', va1[permutation[i], 8]
                 
                     #print 'Output'
-                    #print 'v ', data[i][:parameters_index][8]                    
+                    #print 'v ', data[i][:parameters_index][8]  
+
+                    # South is positive in MUX
+                    #print "data[i][:parameters_index]", data[i][:parameters_index]
+                    #print "-va1[permutation[i], :]", -va1[permutation[i], :]
+                    assert num.allclose(data[i][:parameters_index], -va1[permutation[i], :])
+        
+        self.delete_mux(filesII)
+           
+    def test_read_mux_platform_problem3(self):
+        
+        # This is to test a situation where read_mux returned 
+        # wrong values Win32
+
+        
+        from urs_ext import read_mux2 
+        
+        from anuga.config import single_precision as epsilon        
+        
+        verbose = False
                 
-                    # South is positive in MUX
+        tide = 1.5
+        time_step_count = 10
+        time_step = 0.02
+
+        '''
+        Win results
+        time_step = 0.2000001
+        This is OK        
+        '''
+        
+        '''
+        Win results
+        time_step = 0.20000001
+
+        ======================================================================
+ERROR: test_read_mux_platform_problem3 (__main__.Test_Data_Manager)
+----------------------------------------------------------------------
+Traceback (most recent call last):
+  File "test_data_manager.py", line 6718, in test_read_mux_platform_problem3
+    ha1[0]=num.sin(times_ref)
+ValueError: matrices are not aligned for copy
+
+        '''
+
+        '''
+        Win results
+        time_step = 0.200000001
+        FAIL
+         assert num.allclose(data[i][:parameters_index],
+         -va1[permutation[i], :])
+        '''
+        times_ref = num.arange(0, time_step_count*time_step, time_step)
+        #print "times_ref", times_ref
+        
+        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115),
+                           (-21.,115.), (-22., 117.)]
+        stations = len(lat_long_points)
+        
+        # Create different timeseries starting and ending at different times 
+        first_tstep=num.ones(stations, num.int)
+        first_tstep[0]+=2   # Point 0 starts at 2
+        first_tstep[1]+=4   # Point 1 starts at 4        
+        first_tstep[2]+=3   # Point 2 starts at 3
+        
+        last_tstep=(time_step_count)*num.ones(stations, num.int)
+        last_tstep[0]-=1    # Point 0 ends 1 step early
+        last_tstep[1]-=2    # Point 1 ends 2 steps early                
+        last_tstep[4]-=3    # Point 4 ends 3 steps early        
+        
+        # Create varying elevation data (positive values for seafloor)
+        gauge_depth=20*num.ones(stations, num.float)
+        for i in range(stations):
+            gauge_depth[i] += i**2
+            
+        # Create data to be written to second mux file        
+        ha1=num.ones((stations,time_step_count), num.float)
+        ha1[0]=num.sin(times_ref)
+        ha1[1]=2*num.sin(times_ref - 3)
+        ha1[2]=5*num.sin(4*times_ref)
+        ha1[3]=num.sin(times_ref)
+        ha1[4]=num.sin(2*times_ref-0.7)
+                
+        ua1=num.zeros((stations,time_step_count),num.float)
+        ua1[0]=3*num.cos(times_ref)        
+        ua1[1]=2*num.sin(times_ref-0.7)   
+        ua1[2]=num.arange(3*time_step_count,4*time_step_count)
+        ua1[4]=2*num.ones(time_step_count)
+        
+        va1=num.zeros((stations,time_step_count),num.float)
+        va1[0]=2*num.cos(times_ref-0.87)        
+        va1[1]=3*num.ones(time_step_count)
+        va1[3]=2*num.sin(times_ref-0.71)        
+        #print "va1[0]", va1[0]  # The 8th element is what will go bad.
+        # Ensure data used to write mux file to be zero when gauges are
+        # not recording
+        for i in range(stations):
+             # For each point
+             for j in range(0, first_tstep[i]-1) + range(last_tstep[i],
+                                                         time_step_count):
+                 # For timesteps before and after recording range
+                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0 
+
+
+        #print 'Second station to be written to MUX'
+        #print 'ha', ha1[0,:]
+        #print 'ua', ua1[0,:]
+        #print 'va', va1[0,:]
+        
+        # Write second mux file to be combined by urs2sts 
+        base_nameII, filesII = self.write_mux2(lat_long_points,
+                                               time_step_count, time_step,
+                                               first_tstep, last_tstep,
+                                               depth=gauge_depth,
+                                               ha=ha1,
+                                               ua=ua1,
+                                               va=va1)
+        #print "filesII", filesII
+
+
+
+
+        # Read mux file back and verify it's correcness
+
+        ####################################################
+        # FIXME (Ole): This is where the test should
+        # verify that the MUX files are correct.
+
+        #JJ: It appears as though
+        #that certain quantities are not being stored with enough precision
+        #inn muxfile or more likely that they are being cast into a
+        #lower precision when read in using read_mux2 Time step and q_time
+        # are equal but only to approx 1e-7
+        ####################################################
+
+        #define information as it should be stored in mus2 files
+        points_num=len(lat_long_points)
+        depth=gauge_depth
+        ha=ha1
+        ua=ua1
+        va=va1
+        
+        quantities = ['HA','UA','VA']
+        mux_names = [WAVEHEIGHT_MUX2_LABEL,
+                     EAST_VELOCITY_MUX2_LABEL,
+                     NORTH_VELOCITY_MUX2_LABEL]
+        quantities_init = [[],[],[]]
+        latlondeps = []
+        #irrelevant header information
+        ig=ilon=ilat=0
+        mcolat=mcolon=centerlat=centerlon=offset=az=baz=id=0.0
+        # urs binary is latitude fastest
+        for i,point in enumerate(lat_long_points):
+            lat = point[0]
+            lon = point[1]
+            _ , e, n = redfearn(lat, lon)
+            if depth is None:
+                this_depth = n
+            else:
+                this_depth = depth[i]
+            latlondeps.append([lat, lon, this_depth])
+
+            if ha is None:
+                this_ha = e
+                quantities_init[0].append(num.ones(time_step_count,
+                                                   num.float)*this_ha) # HA
+            else:
+                quantities_init[0].append(ha[i])
+            if ua is None:
+                this_ua = n
+                quantities_init[1].append(num.ones(time_step_count,
+                                                   num.float)*this_ua) # UA
+            else:
+                quantities_init[1].append(ua[i])
+            if va is None:
+                this_va = e
+                quantities_init[2].append(num.ones(time_step_count,
+                                                   num.float)*this_va) #
+            else:
+                quantities_init[2].append(va[i])
+
+        for i, q in enumerate(quantities):
+            #print
+            #print i, q
+            
+            q_time = num.zeros((time_step_count, points_num), num.float)
+            quantities_init[i] = ensure_numeric(quantities_init[i])
+            for time in range(time_step_count):
+                #print i, q, time, quantities_init[i][:,time]
+                q_time[time,:] = quantities_init[i][:,time]
+                #print i, q, time, q_time[time, :]
+
+            
+            filename = base_nameII + mux_names[i]
+            f = open(filename, 'rb')
+            if self.verbose: print 'Reading' + filename
+            assert abs(points_num-unpack('i',f.read(4))[0])<epsilon
+            #write mux 2 header
+            for latlondep in latlondeps:
+                assert abs(latlondep[0]-unpack('f',f.read(4))[0])<epsilon
+                assert abs(latlondep[1]-unpack('f',f.read(4))[0])<epsilon
+                assert abs(mcolat-unpack('f',f.read(4))[0])<epsilon
+                assert abs(mcolon-unpack('f',f.read(4))[0])<epsilon
+                assert abs(ig-unpack('i',f.read(4))[0])<epsilon
+                assert abs(ilon-unpack('i',f.read(4))[0])<epsilon
+                assert abs(ilat-unpack('i',f.read(4))[0])<epsilon
+                assert abs(latlondep[2]-unpack('f',f.read(4))[0])<epsilon
+                assert abs(centerlat-unpack('f',f.read(4))[0])<epsilon
+                assert abs(centerlon-unpack('f',f.read(4))[0])<epsilon
+                assert abs(offset-unpack('f',f.read(4))[0])<epsilon
+                assert abs(az-unpack('f',f.read(4))[0])<epsilon
+                assert abs(baz-unpack('f',f.read(4))[0])<epsilon
+                
+                x = unpack('f', f.read(4))[0]
+                #print time_step
+                #print x
+                assert abs(time_step-x)<epsilon
+                assert abs(time_step_count-unpack('i',f.read(4))[0])<epsilon
+                for j in range(4): # identifier
+                    assert abs(id-unpack('i',f.read(4))[0])<epsilon 
+
+            #first_tstep=1
+            #last_tstep=time_step_count
+            for i,latlondep in enumerate(latlondeps):
+                assert abs(first_tstep[i]-unpack('i',f.read(4))[0])<epsilon
+            for i,latlondep in enumerate(latlondeps):
+                assert abs(last_tstep[i]-unpack('i',f.read(4))[0])<epsilon
+
+            # Find when first station starts recording
+            min_tstep = min(first_tstep)
+            # Find when all stations have stopped recording
+            max_tstep = max(last_tstep)
+
+            #for time in  range(time_step_count):
+            for time in range(min_tstep-1,max_tstep):
+                assert abs(time*time_step-unpack('f',f.read(4))[0])<epsilon
+                for point_i in range(points_num):
+                    if time+1>=first_tstep[point_i] and time+1<=last_tstep[point_i]:
+                        x = unpack('f',f.read(4))[0]
+                        #print time, x, q_time[time, point_i]
+                        if q == 'VA': x = -x # South is positive in MUX
+                        #print q+" q_time[%d, %d] = %f" %(time, point_i, 
+                                                      #q_time[time, point_i])
+                        assert abs(q_time[time, point_i]-x)<epsilon
+
+            f.close()
+                            
+        permutation = ensure_numeric([4,0,2])
+                   
+        # Create ordering file
+#         _, ordering_filename = tempfile.mkstemp('')
+#         order_fid = open(ordering_filename, 'w')  
+#         order_fid.write('index, longitude, latitude\n')
+#         for index in permutation:
+#             order_fid.write('%d, %f, %f\n' %(index, 
+#                                              lat_long_points[index][1], 
+#                                              lat_long_points[index][0]))
+#         order_fid.close()
+        
+        # -------------------------------------
+        # Now read files back and check values
+        weights = ensure_numeric([1.0])
+
+        # For each quantity read the associated list of source mux2 file with 
+        # extention associated with that quantity
+        file_params=-1*num.ones(3,num.float) # [nsta,dt,nt]
+        OFFSET = 5
+
+        for j, file in enumerate(filesII):
+            # Read stage, u, v enumerated as j
+            #print 'Reading', j, file
+            #print "file", file
+            data = read_mux2(1, [file], weights, file_params,
+                             permutation, verbose)
+            #print str(j) + "data", data
+
+            #print 'Data received by Python'
+            #print data[1][8]
+            number_of_selected_stations = data.shape[0]
+            #print "number_of_selected_stations", number_of_selected_stations
+            #print "stations", stations
+
+            # Index where data ends and parameters begin
+            parameters_index = data.shape[1]-OFFSET          
+                 
+            for i in range(number_of_selected_stations):
+        
+                #print i, parameters_index
+                if j == 0:
+                    assert num.allclose(data[i][:parameters_index],
+                                        ha1[permutation[i], :])
+                    
+                if j == 1: assert num.allclose(data[i][:parameters_index], ua1[permutation[i], :])
+                if j == 2:
                     assert num.allclose(data[i][:parameters_index], -va1[permutation[i], :])
-                    
+        
+        self.delete_mux(filesII)            
         
     def test_urs2sts0(self):
@@ -6819 +7099 @@
 
         base_name, files = self.write_mux2(lat_long_points,
-                                      time_step_count, time_step,
-                                      first_tstep, last_tstep,
-                                      depth=gauge_depth,
-                                      ha=ha,
-                                      ua=ua,
-                                      va=va)
+                                           time_step_count, time_step,
+                                           first_tstep, last_tstep,
+                                           depth=gauge_depth,
+                                           ha=ha,
+                                           ua=ua,
+                                           va=va)
 
         urs2sts(base_name,
@@ -6850 +7130 @@
         y = points[:,1]
 
-        #Check that all coordinate are correctly represented       
-        #Using the non standard projection (50) 
+        # Check that all coordinate are correctly represented       
+        # Using the non standard projection (50) 
         for i in range(4):
             zone, e, n = redfearn(lat_long_points[i][0],
@@ -6858 +7138 @@
             assert num.allclose([x[i],y[i]], [e,n])
             assert zone==-1
+        
+        self.delete_mux(files)
 
             
@@ -6915 +7197 @@
         y = points[:,1]
 
-        #Check that all coordinate are correctly represented       
-        #Using the non standard projection (50) 
+        # Check that all coordinate are correctly represented       
+        # Using the non standard projection (50) 
         for i in range(4):
             zone, e, n = redfearn(lat_long_points[i][0], lat_long_points[i][1],
@@ -6922 +7204 @@
             assert num.allclose([x[i],y[i]], [e,n])
             assert zone==geo_reference.zone
+        
+        self.delete_mux(files)
 
             
@@ -7652 +7936 @@
             raise Exception, msg
 
+        
+        self.delete_mux(filesI)
+        self.delete_mux(filesII)
 
         
@@ -8034 +8321 @@
         os.remove(sts_file)
         
-        
-        
         #----------------------
         # Then read the mux files together and test
@@ -8148 +8433 @@
         os.remove(sts_file)
         
-
-        
         #---------------
         # "Manually" add the timeseries up with weights and test
@@ -8157 +8440 @@
         stage_man = weights[0]*(stage0-tide) + weights[1]*(stage1-tide) + tide
         assert num.allclose(stage_man, stage)
-        
-        
-        
-        
-        
                 
         
@@ -8288 +8566 @@
                         
         assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
-                            domain_drchlt.quantities['xmomentum'].vertex_values)                        
+                            domain_drchlt.quantities['xmomentum'].vertex_values) 
                         
         assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
-                            domain_drchlt.quantities['ymomentum'].vertex_values)                                                
+                            domain_drchlt.quantities['ymomentum'].vertex_values)
         
         
         os.remove(sts_file+'.sts')
         os.remove(meshname)
-        
-        
-        
                 
         
@@ -8419 +8694 @@
 
 
-
-        
-        
-        
-            
-
-        
         domain_drchlt = Domain(meshname)
         domain_drchlt.set_quantity('stage', tide)
@@ -8447 +8715 @@
                         
         assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
-                            domain_drchlt.quantities['xmomentum'].vertex_values)                        
+                            domain_drchlt.quantities['xmomentum'].vertex_values) 
                         
         assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
-                            domain_drchlt.quantities['ymomentum'].vertex_values)                                                
-        
+                            domain_drchlt.quantities['ymomentum'].vertex_values) 
         
         os.remove(sts_file+'.sts')
         os.remove(meshname)
-
-
-        
                 
         
@@ -10934 +11198 @@
         domain.set_quantity('xmomentum', uh)
         domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
+
         for t in domain.evolve(yieldstep=1, finaltime = t_end):
             pass
@@ -11357 +11622 @@
 
             
+    def test_copy_code_files(self):
+        '''test that the copy_code_files() function is sane.'''
+
+        def create_file(f):
+            fd = open(f, 'w')
+            fd.write('%s\n' % f)
+            fd.close()
+
+        # create working directories and test files
+        work_dir = tempfile.mkdtemp()
+        dst_dir = tempfile.mkdtemp(dir=work_dir)
+        src_dir = tempfile.mkdtemp(dir=work_dir)
+
+        f1 = 'file1'        
+        filename1 = os.path.join(src_dir, f1)
+        create_file(filename1)
+        f2 = 'file2'        
+        filename2 = os.path.join(src_dir, f2)
+        create_file(filename2)
+        f3 = 'file3'        
+        filename3 = os.path.join(src_dir, f3)
+        create_file(filename3)
+        f4 = 'file4'        
+        filename4 = os.path.join(src_dir, f4)
+        create_file(filename4)
+        f5 = 'file5'        
+        filename5 = os.path.join(src_dir, f5)
+        create_file(filename5)
+
+        # exercise the copy function
+        copy_code_files(dst_dir, filename1)
+        copy_code_files(dst_dir, filename1, filename2)
+        copy_code_files(dst_dir, (filename4, filename5, filename3))
+
+        # test that files were actually copied
+        self.failUnless(access(os.path.join(dst_dir, f1), F_OK))
+        self.failUnless(access(os.path.join(dst_dir, f2), F_OK))
+        self.failUnless(access(os.path.join(dst_dir, f3), F_OK))
+        self.failUnless(access(os.path.join(dst_dir, f4), F_OK))
+        self.failUnless(access(os.path.join(dst_dir, f5), F_OK))
+
+        # clean up
+        shutil.rmtree(work_dir)
             
 #-------------------------------------------------------------
@@ -11362 +11670 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_Data_Manager,'test')
-
-    # FIXME (Ole): This is the test that fails under Windows
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_read_mux_platform_problem2')
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_stsIV')
     
     if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':

branches/numpy/anuga/shallow_water/test_shallow_water_domain.py

@@ -27 +27 @@
 def linear_function(point):
     point = num.array(point)
-    return point[:,0] + point[:,1]
-
+    return point[:,0]+point[:,1]
 
 class Weir:
     """Set a bathymetry for weir with a hole and a downstream gutter
-
     x,y are assumed to be in the unit square
     """
@@ -45 +43 @@
         z = num.zeros(N, num.float)
         for i in range(N):
-            z[i] = -x[i] / 2    # General slope
-
-            # Flattish bit to the left
+            z[i] = -x[i]/2  #General slope
+
+            #Flattish bit to the left
             if x[i] < 0.3:
                 z[i] = -x[i]/10
 
-            # Weir
+            #Weir
             if x[i] >= 0.3 and x[i] < 0.4:
                 z[i] = -x[i]+0.9
 
-            # Dip
+            #Dip
             x0 = 0.6
             depth = -1.0
@@ -62 +60 @@
                 if x[i] > x0 and x[i] < 0.9:
                     z[i] = depth
-                # RHS plateaux
+                #RHS plateaux
                 if x[i] >= 0.9:
                     z[i] = plateaux
             elif y[i] >= 0.7 and y[i] < 1.5:
-                # Restrict and deepen
+                #Restrict and deepen
                 if x[i] >= x0 and x[i] < 0.8:
                     z[i] = depth - (y[i]/3 - 0.3)
                 elif x[i] >= 0.8:
-                    # RHS plateaux
+                    #RHS plateaux
                     z[i] = plateaux
             elif y[i] >= 1.5:
                 if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
-                    # Widen up and stay at constant depth
-                    z[i] = depth - 1.5/5
-                elif x[i] >= 0.8 + (y[i] - 1.5)/1.2:
-                    # RHS plateaux
+                    #Widen up and stay at constant depth
+                    z[i] = depth-1.5/5
+                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
+                    #RHS plateaux
                     z[i] = plateaux
 
-            # Hole in weir (slightly higher than inflow condition)
+            #Hole in weir (slightly higher than inflow condition)
             if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
-                z[i] = -x[i] + self.inflow_stage + 0.02
-
-            # Channel behind weir
+                z[i] = -x[i]+self.inflow_stage + 0.02
+
+            #Channel behind weir
             x0 = 0.5
             if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
-                z[i] = -x[i] + self.inflow_stage + 0.02
+                z[i] = -x[i]+self.inflow_stage + 0.02
 
             if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
-                # Flatten it out towards the end
-                z[i] = -x0 + self.inflow_stage + 0.02 + (x0 - x[i])/5
+                #Flatten it out towards the end
+                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
 
             # Hole to the east
@@ -97 +95 @@
             y0 = 0.35
             if num.sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
-                z[i] = num.sqrt((x[i]-x0)**2 + (y[i]-y0)**2) - 1.0
-
-            # Tiny channel draining hole
+                z[i] = num.sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
+
+            #Tiny channel draining hole
             if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
-                z[i] = -0.9    # North south
+                z[i] = -0.9 #North south
 
             if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
-                z[i] = -1.0 + (x[i]-0.9)/3    # East west
+                z[i] = -1.0 + (x[i]-0.9)/3 #East west
 
             # Stuff not in use
@@ -136 +134 @@
         z = num.zeros(N, num.float)
         for i in range(N):
-            z[i] = -x[i]    # General slope
-
-            # Flat bit to the left
+            z[i] = -x[i]  #General slope
+
+            #Flat bit to the left
             if x[i] < 0.3:
-                z[i] = -x[i]/10    # General slope
-
-            # Weir
+                z[i] = -x[i]/10  #General slope
+
+            #Weir
             if x[i] > 0.3 and x[i] < 0.4:
-                z[i] = -x[i] + 0.9
-
-            # Dip
+                z[i] = -x[i]+0.9
+
+            #Dip
             if x[i] > 0.6 and x[i] < 0.9:
-                z[i] = -x[i] - 0.5    # -y[i]/5
-
-            # Hole in weir (slightly higher than inflow condition)
+                z[i] = -x[i]-0.5  #-y[i]/5
+
+            #Hole in weir (slightly higher than inflow condition)
             if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
-                z[i] = -x[i] + self.inflow_stage + 0.05
+                z[i] = -x[i]+self.inflow_stage + 0.05
+
 
         return z/2
@@ -170 +169 @@
 
     N = len(x)
-    s = 0 * x    # New array
+    s = 0*x  #New array
 
     for k in range(N):
@@ -263 +262 @@
         H0 = 0.0
 
-        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
-                                  epsilon, g, H0)
-
-        assert num.allclose(edgeflux, [0, 0, 0])
+        
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
+
+        assert num.allclose(edgeflux, [0,0,0])
         assert max_speed == 0.
 
@@ -272 +271 @@
         w = 2.0
 
-        normal = num.array([1., 0])
+        normal = num.array([1.,0])
         ql = num.array([w, 0, 0])
         qr = num.array([w, 0, 0])
@@ -280 +279 @@
         H0 = 0.0
 
-        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
-                                  epsilon, g, H0)
-
+        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)        
         assert num.allclose(edgeflux, [0., 0.5*g*h**2, 0.])
         assert max_speed == num.sqrt(g*h)
@@ -290 +287 @@
     #    w = 2.0
     #
-    #    normal = array([1., 0])
+    #    normal = array([1.,0])
     #    ql = array([w, 0, 0])
     #    qr = array([w, 0, 0])
@@ -1517 +1514 @@
         uh = u*h
 
-        Br = Reflective_boundary(domain)       # Side walls
-        Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
+        Br = Reflective_boundary(domain)     # Side walls
+        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet: 
+
 
         # Initial conditions

branches/numpy/anuga/shallow_water/urs_ext.c

@@ -104 +104 @@
                 memcpy(data + it, muxData + offset, sizeof(float));
                 
-                //printf("%d: muxdata=%f\n", it, muxData[offset]);                      
-                //printf("data[%d]=%f, offset=%d\n", it, data[it], offset);         
+                //printf("%d: muxdata=%f\n", it, muxData[offset]); 
+                //printf("data[%d]=%f, offset=%d\n", it, data[it], offset); 
                 offset++;
             }
@@ -221 +221 @@
 
         // Open the mux file
-        if((fp = fopen(muxFileName, "r")) == NULL)
+        if((fp = fopen(muxFileName, "rb")) == NULL)
         {
             char *err_msg = strerror(errno);
@@ -237 +237 @@
             }
         
-            fros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int)); 
+            fros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int));
             lros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int));
       
@@ -447 +447 @@
     temp_sts_data = (float*) calloc(len_sts_data, sizeof(float));
 
-    muxData = (float*) malloc(numDataMax*sizeof(float));
+    muxData = (float*) calloc(numDataMax, sizeof(float));
     
     // Loop over all sources
@@ -460 +460 @@
         // Read in data block from mux2 file
         muxFileName = muxFileNameArray[isrc];
-        if((fp = fopen(muxFileName, "r")) == NULL)
+        if((fp = fopen(muxFileName, "rb")) == NULL)
         {
             fprintf(stderr, "cannot open file %s\n", muxFileName);
@@ -470 +470 @@
         }
 
-        offset = sizeof(int) + total_number_of_stations*(sizeof(struct tgsrwg) + 2*sizeof(int));
-        fseek(fp, offset, 0);
+        offset = (long int)sizeof(int) + total_number_of_stations*(sizeof(struct tgsrwg) + 2*sizeof(int));
+        //printf("\n offset %i ", (long int)offset);
+        fseek(fp, offset, 0);
 
         numData = getNumData(fros_per_source, 
                              lros_per_source, 
                              total_number_of_stations);
-                             
-                             
-                                     
-        elements_read = fread(muxData, ((int) numData)*sizeof(float), 1, fp); 
+        // Note numData is larger than what it has to be.                    
+        //elements_read = fread(muxData, ((int) numData)*sizeof(float), 1, fp); 
+        elements_read = fread(muxData, (size_t) sizeof(float), (size_t) numData, fp); 
+        //printf("\n elements_read  %d, ", (int)elements_read);
+        //printf("\n ferror(fp)  %d, ", (int)ferror(fp));
         if ((int) elements_read == 0 && ferror(fp)) {
           fprintf(stderr, "Error reading mux data\n");
           return NULL;
-        }
-                
-        fclose(fp);
+        }       
         
-
-        // FIXME (Ole): This is where Nariman and Ole traced the platform dependent 
-        // difference on 11 November 2008. We don't think the problem lies in the 
-        // C code. Maybe it is a problem with the MUX files written by the unit test
-        // that fails on Windows but works OK on Linux. JJ's test on 17th November shows
-        // that as far as Python is concerned this file should be OK on both platforms.
-        
-        
-        
-        // These printouts are enough to show the problem when compared 
-        // on the two platforms 
-        //printf("\nRead %d elements, ", (int) numData);
-        //printf("muxdata[%d]=%f\n", 39, muxData[39]);          
-        
-        
-        /*
-        In Windows we get
-        
-        Read 85 elements, muxdata[39]=0.999574
-        Read 85 elements, muxdata[39]=-0.087599
-        Read 85 elements, muxdata[39]=-0.087599
-        
-        
-        In Linux we get (the correct)
-        
-        Read 85 elements, muxdata[39]=0.999574
-        Read 85 elements, muxdata[39]=-0.087599
-        Read 85 elements, muxdata[39]=1.490349
-        */
-        
-        
-        
-        
+        fclose(fp);  
         
         // loop over stations present in the permutation array