Changeset 907

Timestamp:

Feb 17, 2005, 5:36:29 PM (20 years ago)

Author:

ole

Message:

Addressed problem with artificial momentum generated by discontinuous water depths in the presence of steep slopes.
Now very shallow water is limited with a separate h-limiter controlled by beta_h
(see config.py) for details.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• config.py (modified) (1 diff)
• domain.py (modified) (1 diff)
• netherlands.py (modified) (2 diffs)
• quantity.py (modified) (2 diffs)
• quantity_ext.c (modified) (4 diffs)
• quantity_ext.h (modified) (1 diff)
• shallow_water.py (modified) (13 diffs)
• shallow_water_ext.c (modified) (10 diffs)
• test_data_manager.py (modified) (1 diff)
• test_interpolate_sww.py (modified) (2 diffs)
• test_shallow_water.py (modified) (6 diffs)
• util_ext.h (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/config.py

@@ -45 +45 @@
 
 
+#Betas [0;1] control the allowed steepness of gradient for second order
+#extrapolations. Values of 1 allow the steepes gradients while
+#lower values are more conservative. Values of 0 correspond to
+#1'st order extrapolations.
+#
+# Large values of beta_h may cause simulations to require more timesteps
+# as surface will 'hug' closer to the bed.
+# Small values of beta_h will make code faster, but one may experience
+# artificial momenta caused by discontinuities in water depths in
+# the presence of steep slopes. One example of this would be
+# stationary water 'lapping' upwards to a higher point on the coast.
+# 
+#
+#
+#There are separate betas for the w-limiter and the h-limiter
+#
+#
+#
+#
+#Good values are:
+#beta_w = 0.9
+#beta_h = 0.2
 
-beta = 0.9  #]0;1] (e.g. 0.9). 1 allows the steepest gradients,
-            #lower values are the more conservative
-            #the limiter being a first order method.
+
+
+beta_w = 0.9
+beta_h = 0.2 
+
 
 pmesh_filename = '.\\pmesh'

inundation/ga/storm_surge/pyvolution/domain.py

@@ -52 +52 @@
 
         #Defaults
-        from config import max_smallsteps, beta, epsilon
-        self.beta = beta
+        from config import max_smallsteps, beta_w, beta_h, epsilon
+        self.beta_w = beta_w
+        self.beta_h = beta_h        
         self.epsilon = epsilon
+
+        #FIXME: Maybe have separate orders for h-limiter and w-limiter?
+        #Or maybe get rid of order altogether and use beta_w and beta_h
         self.default_order = 1
         self.order = self.default_order

inundation/ga/storm_surge/pyvolution/netherlands.py

@@ -99 +99 @@
 #N = 140
 
-N = 15
+#N = 15
 
 print 'Creating domain'
@@ -110 +110 @@
 domain.check_integrity()
 domain.default_order = 2
+#domain.beta_h=0
 
 #Output params

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -813 +813 @@
     N = quantity.domain.number_of_elements
     
-    beta = quantity.domain.beta
+    beta_w = quantity.domain.beta_w
         
     qc = quantity.centroid_values
@@ -852 +852 @@
             if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
             
-            phi = min( min(r*beta, 1), phi )    
+            phi = min( min(r*beta_w, 1), phi )    
 
         #Then update using phi limiter

inundation/ga/storm_surge/pyvolution/quantity_ext.c

@@ -15 +15 @@
 
 //Shared code snippets
-#include "util_ext.h"
-#include "quantity_ext.h"
+#include "util_ext.h" 
+//#include "quantity_ext.h" //Obsolete
 
 
@@ -450 +450 @@
     
   int k, i, n, N, k3;
-  double beta; //Safety factor
+  double beta_w; //Safety factor
   double *qmin, *qmax, qn;
   
@@ -464 +464 @@
   neighbours = get_consecutive_array(domain, "neighbours");
 
-  //Get safety factor beta
-  Tmp = PyObject_GetAttrString(domain, "beta"); 
+  //Get safety factor beta_w
+  Tmp = PyObject_GetAttrString(domain, "beta_w"); 
   if (!Tmp) 
     return NULL;
       
-  beta = PyFloat_AsDouble(Tmp);  
+  beta_w = PyFloat_AsDouble(Tmp);  
   
   Py_DECREF(Tmp);    
@@ -502 +502 @@
   
   // Call underlying routine
-  _limit(N, beta, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
+  _limit(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
          
   free(qmin);

inundation/ga/storm_surge/pyvolution/quantity_ext.h

@@ -1 @@
-
-void _limit(int N, double beta, double* qc, double* qv, 
-            double* qmin, double* qmax) { 
-
-  int k, i, k3;
-  double dq, dqa[3], phi, r;
-  
-  for (k=0; k<N; k++) {
-    k3 = k*3;
-    
-    //Find the gradient limiter (phi) across vertices  
-    phi = 1.0;
-    for (i=0; i<3; i++) {    
-      r = 1.0;
-      
-      dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
-      dqa[i] = dq;              //Save dq for use in the next loop
-      
-      if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
-      if (dq < 0.0) r = (qmin[k] - qc[k])/dq;      
-  
-  
-      phi = min( min(r*beta, 1.0), phi);    
-    }
-    
-    //Then update using phi limiter
-    for (i=0; i<3; i++) {    
-      qv[k3+i] = qc[k] + phi*dqa[i];
-    }
-  }
-}

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -521 +521 @@
             raise 'Unknown order'
 
-    #Take bed elevation into account when water heights are small    
+    #Take bed elevation into account when water heights are small
     balance_deep_and_shallow(domain)
 
@@ -594 +594 @@
     protect(domain.minimum_allowed_height, wc, zc, xmomc, ymomc)
     
-    
-def balance_deep_and_shallow_old(domain):
-    """Compute linear combination between stage as computed by
-    gradient-limiters and stage computed as constant height above bed.
-    The former takes precedence when heights are large compared to the
-    bed slope while the latter takes precedence when heights are
-    relatively small.  Anything in between is computed as a balanced
-    linear combination in order to avoid numerical disturbances which
-    would otherwise appear as a result of hard switching between
-    modes.
-    """
-    
-    #Shortcuts
-    wc = domain.quantities['stage'].centroid_values
-    zc = domain.quantities['elevation'].centroid_values
-    hc = wc - zc
-    
-    wv = domain.quantities['stage'].vertex_values
-    zv = domain.quantities['elevation'].vertex_values
-    hv = wv-zv
-
-
-    #Computed linear combination between constant stages and and
-    #stages parallel to the bed elevation.     
-    for k in range(domain.number_of_elements):
-        #Compute maximal variation in bed elevation
-        #  This quantitiy is
-        #    dz = max_i abs(z_i - z_c)
-        #  and it is independent of dimension
-        #  In the 1d case zc = (z0+z1)/2
-        #  In the 2d case zc = (z0+z1+z2)/3
-
-        dz = max(abs(zv[k,0]-zc[k]),
-                 abs(zv[k,1]-zc[k]),
-                 abs(zv[k,2]-zc[k]))
-
-        
-        hmin = min( hv[k,:] )
-
-        #Create alpha in [0,1], where alpha==0 means using shallow 
-        #first order scheme and alpha==1 means using the stage w as
-        #computed by the gradient limiter (1st or 2nd order)
-        #
-        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
-        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
-      
-        if dz > 0.0:
-            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
-        else:
-            #Flat bed 
-            alpha = 1.0
-
-        #Weighted balance between stage parallel to bed elevation 
-        #(wvi = zvi + hc) and stage as computed by 1st or 2nd 
-        #order gradient limiter
-        #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
-        #
-        #It follows that the updated wvi is
-        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = 
-        #  zvi + hc + alpha*(hvi - hc)
-        #
-        #Note that hvi = zc+hc-zvi in the first order case (constant).
-
-        if alpha < 1:         
-            for i in range(3):
-                wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
-            
-
-            #Momentums at centroids
-            xmomc = domain.quantities['xmomentum'].centroid_values
-            ymomc = domain.quantities['ymomentum'].centroid_values        
-
-            #Momentums at vertices
-            xmomv = domain.quantities['xmomentum'].vertex_values
-            ymomv = domain.quantities['ymomentum'].vertex_values         
-
-            # Update momentum as a linear combination of 
-            # xmomc and ymomc (shallow) and momentum
-            # from extrapolator xmomv and ymomv (deep).
-            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]            
-            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
-            
-
-def limit_on_h(domain):        
+
+
+def h_limiter(domain):        
     """Limit slopes for each volume to eliminate artificial variance
     introduced by e.g. second order extrapolator
     
     limit on h = w-z
+
+    This limiter depends on two quantities (w,z) so it resides within
+    this module rather than within quantity.py
     """
 
@@ -687 +609 @@
 
     N = domain.number_of_elements
-    beta = domain.beta    
+    beta_h = domain.beta_h
     
     #Shortcuts
@@ -734 +656 @@
             if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
             
-            phi = min( min(r*beta, 1), phi )    
+            phi = min( min(r*beta_h, 1), phi )    
 
         #Then update using phi limiter
@@ -741 +663 @@
 
     return hvbar
-    
+
+
+
+def h_limiter_c(domain):        
+    """Limit slopes for each volume to eliminate artificial variance
+    introduced by e.g. second order extrapolator
+    
+    limit on h = w-z
+
+    This limiter depends on two quantities (w,z) so it resides within
+    this module rather than within quantity.py
+
+    Wrapper for c-extension
+    """
+
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+    beta_h = domain.beta_h
+    
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+    
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv-zv
+
+    #Call C-extension
+    from shallow_water_ext import h_limiter
+    hvbar = h_limiter(domain, hc, hv)
+
+    return hvbar
+
                                             
 def balance_deep_and_shallow(domain):
@@ -776 +732 @@
     hv = wv-zv
 
+    #Limit h
+    hvbar = h_limiter(domain)
+
     for k in range(domain.number_of_elements):
         #Compute maximal variation in bed elevation
@@ -791 +750 @@
         hmin = min( hv[k,:] )
 
-        #Create alpha in [0,1], where alpha==0 means using shallow 
-        #first order scheme and alpha==1 means using the stage w as
-        #computed by the gradient limiter (1st or 2nd order)
-        #
+        #Create alpha in [0,1], where alpha==0 means using the h-limited 
+        #stage and alpha==1 means using the w-limited stage as
+        #computed by the gradient limiter (both 1st or 2nd order)
+    
         #If hmin > dz/2 then alpha = 1 and the bed will have no effect
         #If hmin < 0 then alpha = 0 reverting to constant height above bed.
-      
+
         if dz > 0.0:
             alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
@@ -805 +764 @@
 
         #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~)
         #  
         #    
@@ -814 +773 @@
         #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~
+        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
         #
-        # Momentum is balanced between constant and limited ???
-
-        #print 'Alpha = ', alpha
-        #FIXME: if alpha == 0 compute only h-limited slope, 
-        #if alpha == 1 compute only w-limited slope
-        #
+        # Momentum is balanced between constant and limited
+
         if alpha < 1:         
         
-            #Limit h
-            hvbar = limit_on_h(domain)
-    
             for i in range(3):
                 wv[k,i] = zv[k,i] + (1-alpha)*hvbar[k,i] + alpha*hv[k,i]
@@ -848 +800 @@
                                     
             
-
 def balance_deep_and_shallow_c(domain):
     """Wrapper for C implementation
@@ -870 +821 @@
     ymomv = domain.quantities['ymomentum'].vertex_values         
 
-    
+    #Limit h
+    hvbar = h_limiter(domain)    
 
     from shallow_water_ext import balance_deep_and_shallow
-    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
-                             xmomc, ymomc, xmomv, ymomv)
+    #balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
+    #                         xmomc, ymomc, xmomv, ymomv)
+
+    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv, hvbar,
+                             xmomc, ymomc, xmomv, ymomv)    
 
     
@@ -1272 +1227 @@
         ymom_update[k] += S*v        
             
+
+
+##############################
+#OBSOLETE STUFF
+    
+def balance_deep_and_shallow_old(domain):
+    """Compute linear combination between stage as computed by
+    gradient-limiters and stage computed as constant height above bed.
+    The former takes precedence when heights are large compared to the
+    bed slope while the latter takes precedence when heights are
+    relatively small.  Anything in between is computed as a balanced
+    linear combination in order to avoid numerical disturbances which
+    would otherwise appear as a result of hard switching between
+    modes.
+    """
+
+    #OBSOLETE
+    
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+    
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv-zv
+
+
+    #Computed linear combination between constant stages and and
+    #stages parallel to the bed elevation.     
+    for k in range(domain.number_of_elements):
+        #Compute maximal variation in bed elevation
+        #  This quantitiy is
+        #    dz = max_i abs(z_i - z_c)
+        #  and it is independent of dimension
+        #  In the 1d case zc = (z0+z1)/2
+        #  In the 2d case zc = (z0+z1+z2)/3
+
+        dz = max(abs(zv[k,0]-zc[k]),
+                 abs(zv[k,1]-zc[k]),
+                 abs(zv[k,2]-zc[k]))
+
+        
+        hmin = min( hv[k,:] )
+
+        #Create alpha in [0,1], where alpha==0 means using shallow 
+        #first order scheme and alpha==1 means using the stage w as
+        #computed by the gradient limiter (1st or 2nd order)
+        #
+        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
+        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
+      
+        if dz > 0.0:
+            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
+        else:
+            #Flat bed 
+            alpha = 1.0
+
+        #Weighted balance between stage parallel to bed elevation 
+        #(wvi = zvi + hc) and stage as computed by 1st or 2nd 
+        #order gradient limiter
+        #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
+        #
+        #It follows that the updated wvi is
+        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = 
+        #  zvi + hc + alpha*(hvi - hc)
+        #
+        #Note that hvi = zc+hc-zvi in the first order case (constant).
+
+        if alpha < 1:         
+            for i in range(3):
+                wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
+            
+
+            #Momentums at centroids
+            xmomc = domain.quantities['xmomentum'].centroid_values
+            ymomc = domain.quantities['ymomentum'].centroid_values        
+
+            #Momentums at vertices
+            xmomv = domain.quantities['xmomentum'].vertex_values
+            ymomv = domain.quantities['ymomentum'].vertex_values         
+
+            # Update momentum as a linear combination of 
+            # xmomc and ymomc (shallow) and momentum
+            # from extrapolator xmomv and ymomv (deep).
+            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]            
+            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
+            
 
 
@@ -1457 +1500 @@
     gravity = gravity_c
     manning_friction = manning_friction_c
+    h_limiter = h_limiter_c
     balance_deep_and_shallow = balance_deep_and_shallow_c
-    #balance_deep_and_shallow = balance_deep_and_shallow_old
     protect_against_infinitesimal_and_negative_heights = protect_against_infinitesimal_and_negative_heights_c
 

inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

@@ -190 +190 @@
                               double* zv, 
                               double* hv,
+                              double* hvbar,                          
                               double* xmomc, 
                               double* ymomc, 
@@ -198 +199 @@
   double dz, hmin, alpha;
   
-  //Compute linear combination between constant stages and and
-  //stages parallel to the bed elevation.     
+  //Compute linear combination between w-limited stages and
+  //h-limited stages close to the bed elevation.     
   
   for (k=0; k<N; k++) {
@@ -220 +221 @@
 
     
-    //Create alpha in [0,1], where alpha==0 means using shallow 
-    //first order scheme and alpha==1 means using the stage w as
-    //computed by the gradient limiter (1st or 2nd order)
+    //Create alpha in [0,1], where alpha==0 means using the h-limited 
+    //stage and alpha==1 means using the w-limited stage as
+    //computed by the gradient limiter (both 1st or 2nd order)
     //
     //If hmin > dz/2 then alpha = 1 and the bed will have no effect
@@ -232 +233 @@
       alpha = 1.0;  //Flat bed
 
-      
-    //Weighted balance between stage parallel to bed elevation 
-    //(wvi = zvi + hc) and stage as computed by 1st or 2nd 
-    //order gradient limiter
-    //(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
-    //
-    //It follows that the updated wvi is
-    //  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = 
-    //  zvi + hc + alpha*(hvi - hc)
-    //
-    //Note that hvi = zc+hc-zvi in the first order case (constant).
+
+    //  Let 
+    //  
+    //    wvi be the w-limited stage (wvi = zvi + hvi)        
+    //    wvi- be the h-limited state (wvi- = zvi + hvi-)
+    //      
+    //      
+    //  where i=0,1,2 denotes the vertex ids
+    //  
+    //  Weighted balance between w-limited and h-limited stage is  
+    //    
+    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)   
+    //  
+    //  It follows that the updated wvi is
+    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
+    //   
+    //   Momentum is balanced between constant and limited
 
     if (alpha < 1) {         
       for (i=0; i<3; i++) {
-        wv[k3+i] = zv[k3+i] + hc[k] + alpha*(hv[k3+i]-hc[k]);
+        wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i];
             
       
@@ -286 +293 @@
   return 0;
 }
+
+
+
+
+
 
 
@@ -683 +695 @@
     *wv,            //Stage at vertices
     *zv,            //Elevation at vertices
-    *hv,            //Heights at vertices    
+    *hv,            //Depths at vertices    
+    *hvbar,         //h-Limited depths at vertices        
     *xmomc,         //Momentums at centroids and vertices
     *ymomc, 
@@ -692 +705 @@
   
   // Convert Python arguments to C  
-  if (!PyArg_ParseTuple(args, "OOOOOOOOOO",  
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",  
                         &wc, &zc, &hc, 
-                        &wv, &zv, &hv,
+                        &wv, &zv, &hv, &hvbar,
                         &xmomc, &ymomc, &xmomv, &ymomv))
     return NULL;
@@ -707 +720 @@
                             (double*) zv -> data, 
                             (double*) hv -> data,
+                            (double*) hvbar -> data,                        
                             (double*) xmomc -> data, 
                             (double*) ymomc -> data, 
@@ -718 +732 @@
 
 
+PyObject *h_limiter(PyObject *self, PyObject *args) {
+  
+  PyObject *domain, *Tmp;
+  PyArrayObject 
+    *hv, *hc, //Depth at vertices and centroids        
+    *hvbar,   //Limited depth at vertices (return values)
+    *neighbours;
+    
+  int k, i, n, N, k3;
+  int dimensions[2];
+  double beta_h; //Safety factor (see config.py)
+  double *hmin, *hmax, hn;
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) 
+    return NULL;
+
+  neighbours = get_consecutive_array(domain, "neighbours");
+
+  //Get safety factor beta_h
+  Tmp = PyObject_GetAttrString(domain, "beta_h"); 
+  if (!Tmp) 
+    return NULL;
+      
+  beta_h = PyFloat_AsDouble(Tmp);  
+  
+  Py_DECREF(Tmp);    
+
+  N = hc -> dimensions[0];
+  
+  //Create hvbar
+  dimensions[0] = N;
+  dimensions[1] = 3;  
+  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
+
+      
+  //Find min and max of this and neighbour's centroid values
+  hmin = malloc(N * sizeof(double));
+  hmax = malloc(N * sizeof(double));  
+  for (k=0; k<N; k++) {
+    k3=k*3;
+    
+    hmin[k] = ((double*) hc -> data)[k];
+    hmax[k] = hmin[k];
+    
+    for (i=0; i<3; i++) {
+      n = ((long*) neighbours -> data)[k3+i];
+      if (n >= 0) {
+        hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
+          
+        hmin[k] = min(hmin[k], hn);
+        hmax[k] = max(hmax[k], hn);
+      }
+    }
+  }
+  
+  
+  // Call underlying routine
+  _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
+         
+  // // //Py_DECREF(domain); //FIXME: NEcessary?         
+  free(hmin);
+  free(hmax);  
+  
+  //return result using PyArray to avoid memory leak  
+  return PyArray_Return(hvbar);
+}
 
 
@@ -776 +857 @@
   {"balance_deep_and_shallow", balance_deep_and_shallow, 
    METH_VARARGS, "Print out"},          
+  {"h_limiter", h_limiter, 
+   METH_VARARGS, "Print out"},             
   {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},   
   {"assign_windfield_values", assign_windfield_values, 

inundation/ga/storm_surge/pyvolution/test_data_manager.py

@@ -414 +414 @@
         self.domain.smooth = False
         self.domain.store = True
+        self.domain.beta_h = 0
 
         #Evolution

inundation/ga/storm_surge/pyvolution/test_interpolate_sww.py

@@ -28 +28 @@
         domain = Domain(points, vertices, boundary)
         domain.default_order=2
+        domain.beta_h = 0
 
         
@@ -43 +44 @@
         ###################### 
         #Initial condition - with jumps
-
 
         bed = domain.quantities['elevation'].vertex_values

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -1434 +1434 @@
 
         #print E
-        domain.order = 2        
+        domain.order = 2
+        domain.beta_h = 0.0 #Use first order in h-limiter        
         domain.distribute_to_vertices_and_edges()
 
@@ -1837 +1838 @@
         domain.smooth = False
         domain.default_order=1
+        domain.beta_h = 0.0 #Use first order in h-limiter
 
         #Bed-slope and friction
@@ -1986 +1988 @@
         domain.smooth = False
         domain.default_order=2
+        domain.beta_h = 0.0 #Use first order in h-limiter        
 
         #Bed-slope and friction at vertices (and interpolated elsewhere)
@@ -2078 +2081 @@
         domain.smooth = False
         domain.default_order=2
+        domain.beta_h = 0.0 #Use first order in h-limiter        
 
         #Bed-slope and friction at vertices (and interpolated elsewhere)
@@ -2166 +2170 @@
         domain.smooth = False
         domain.default_order=2
+        domain.beta_h = 0.0 #Use first order in h-limiter        
 
         #Bed-slope and friction at vertices (and interpolated elsewhere)
@@ -2282 +2287 @@
         domain.smooth = False
         domain.default_order=2
+        domain.beta_h = 0.0 #Use first order in h-limiter        
 
         #Bed-slope and friction at vertices (and interpolated elsewhere)

inundation/ga/storm_surge/pyvolution/util_ext.h

@@ -51 +51 @@
 
 
+void _limit(int N, double beta, double* qc, double* qv, 
+            double* qmin, double* qmax) { 
+
+  //N are the number of elements
+  int k, i, k3;
+  double dq, dqa[3], phi, r;
+  
+  //printf("INSIDE\n");
+  for (k=0; k<N; k++) {
+    k3 = k*3;
+    
+    //Find the gradient limiter (phi) across vertices  
+    phi = 1.0;
+    for (i=0; i<3; i++) {    
+      r = 1.0;
+      
+      dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
+      dqa[i] = dq;              //Save dq for use in the next loop
+      
+      if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
+      if (dq < 0.0) r = (qmin[k] - qc[k])/dq;      
+  
+  
+      phi = min( min(r*beta, 1.0), phi);    
+    }
+    
+    //Then update using phi limiter
+    for (i=0; i<3; i++) {    
+      qv[k3+i] = qc[k] + phi*dqa[i];
+    }
+  }
+}
+
+
 void print_numeric_array(PyArrayObject *x) {  
   int i, j;