Changeset 4239

Timestamp:

Feb 7, 2007, 6:18:28 PM (18 years ago)

Author:

ole

Message:

Implemented new balanced slope limiter.
To enable, set

domain.limit2007 = 1

in the ANUGA scripts.
Later, we may make it default depending on how it goes.

Location:

anuga_core/source/anuga

Files:

• config.py (modified) (1 diff)
• shallow_water/shallow_water_domain.py (modified) (3 diffs)
• shallow_water/shallow_water_ext.c (modified) (8 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (2 diffs)

anuga_core/source/anuga/config.py

@@ -88 +88 @@
 alpha_balance = 2.0 
 
+# Flag use of new limiters.
+# limit2007 = 0 means use old limiters (e.g. for some tests)
+# limit2007 = 1 means use new limiters that hug the bathymetry closer
+limit2007 = 0
+
 
 CFL = 1.0  #FIXME (ole): Is this in use yet??

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -100 +100 @@
 from anuga.config import minimum_allowed_height, maximum_allowed_speed
 from anuga.config import g, beta_h, beta_w, beta_w_dry,\
-     beta_uh, beta_uh_dry, beta_vh, beta_vh_dry
+     beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, limit2007
 from anuga.config import alpha_balance
 
@@ -158 +158 @@
         self.beta_h      = beta_h
         self.alpha_balance = alpha_balance
+
+        self.limit2007 = limit2007
 
         self.flux_function = flux_function_central
@@ -1184 +1186 @@
 
     #Limit h
-    hvbar = h_limiter(domain)
-
-    #This is how one would make a first order h_limited value
-    #as in the old balancer (pre 17 Feb 2005):
-    #from Numeric import zeros, Float
-    #hvbar = zeros( (len(hc), 3), Float)
-    #for i in range(3):
-    #    hvbar[:,i] = hc[:]
+    if domain.beta_h > 0:
+        hvbar = h_limiter(domain)
+    else:
+        # print 'Using first order h-limiter'
+      
+        
+        #This is how one would make a first order h_limited value
+        #as in the old balancer (pre 17 Feb 2005):
+        # If we wish to hard wire this, one should modify the C-code
+        from Numeric import zeros, Float
+        hvbar = zeros( (len(hc), 3), Float)
+        for i in range(3):
+            hvbar[:,i] = hc[:]
 
     from shallow_water_ext import balance_deep_and_shallow

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -428 +428 @@
                               double* xmomv,
                               double* ymomv,
+                              double H0,
+                              int limit2007,
                               double alpha_balance) {
 
   int k, k3, i;
-  double dz, hmin, alpha;
+  double dz, hmin, alpha, h_diff;
 
   //Compute linear combination between w-limited stages and
@@ -450 +452 @@
     hmin = hv[k3];
     for (i=0; i<3; i++) {
-      dz = max(dz, fabs(zv[k3+i]-zc[k]));
+      if (limit2007 == 0) { 
+        dz = max(dz, fabs(zv[k3+i]-zc[k]));
+      }
+      
       hmin = min(hmin, hv[k3+i]);
     }
@@ -458 +463 @@
     //stage and alpha==1 means using the w-limited stage as
     //computed by the gradient limiter (both 1st or 2nd order)
-    //
-    //If hmin > dz/alpha_balance then alpha = 1 and the bed will have no effect
-    //If hmin < 0 then alpha = 0 reverting to constant height above bed.
-    //The parameter alpha_balance==2 by default 
-
-
-    if (dz > 0.0) {
-      alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );      
+    
+    
+    if (limit2007 == 0) {     
+      //If hmin > dz/alpha_balance then alpha = 1 and the bed will have no 
+      //effect
+      //If hmin < 0 then alpha = 0 reverting to constant height above bed.
+      //The parameter alpha_balance==2 by default 
+
+      
+      if (dz > 0.0) {
+        alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );      
+      } else {
+        alpha = 1.0;  //Flat bed
+      }
+      //printf("Using old style limiter\n");
+      
     } else {
-      alpha = 1.0;  //Flat bed
-    }
-
+
+      // 2007 Balanced Limiter
     
+      // Make alpha as large as possible but still ensure that 
+      // final depth is positive
+    
+      if (hmin < H0) {
+        alpha = 1.0;
+        for (i=0; i<3; i++) {
+
+          h_diff = hvbar[k3+i] - hv[k3+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, (hvbar[k3+i] - 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
+        alpha = 1.0;       
+      }
+    }
+        
+    
+        
     //printf("k=%d, hmin=%.2f, dz=%.2f, alpha=%.2f, alpha_balance=%.2f\n", 
     //     k, hmin, dz, alpha, alpha_balance);
-      
-          
-    //alpha = 1.0;  //Always deep FIXME: This actually looks good now
-    //alpha = 0.2;
 
     //printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
@@ -499 +540 @@
     if (alpha < 1) {
       for (i=0; i<3; i++) {
-         wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i];
+        wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i];
 
         //Update momentum as a linear combination of
         //xmomc and ymomc (shallow) and momentum
         //from extrapolator xmomv and ymomv (deep).
+        //FIXME (Ole): Is this really needed?
         xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
         ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
@@ -550 +592 @@
       } else {
         //Reduce excessive speeds derived from division by small hc
+        //FIXME (Ole): This may be unnecessary with new slope limiters 
+        //in effect.
         
         u = xmomc[k]/hc;
@@ -1594 +1638 @@
     
   double alpha_balance = 2.0;
-
-  int N; //, err;
+  double H0;
+
+  int N, limit2007; //, err;
 
   // Convert Python arguments to C
@@ -1616 +1661 @@
   Py_DECREF(Tmp);
 
-
+  
+  Tmp = PyObject_GetAttrString(domain, "H0");
+  if (!Tmp) {
+    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object H0 from domain");
+    return NULL;
+  }  
+  H0 = PyFloat_AsDouble(Tmp);
+  Py_DECREF(Tmp);
+
+  
+  Tmp = PyObject_GetAttrString(domain, "limit2007");
+  if (!Tmp) {
+    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object limit2007 from domain");
+    return NULL;
+  }  
+  limit2007 = PyInt_AsLong(Tmp);
+  Py_DECREF(Tmp);
+    
+
+  
+  
   N = wc -> dimensions[0];
 
@@ -1627 +1692 @@
                             (double*) hv -> data,
                             (double*) hvbar -> data,
-                                (double*) xmomc -> data,
+                            (double*) xmomc -> data,
                             (double*) ymomc -> data,
                             (double*) xmomv -> data,
                             (double*) ymomv -> data,
+                            H0,
+                            (int) limit2007,
                             alpha_balance);
 

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -1052 +1052 @@
 
         
-
-
         
         assert allclose(gauge_values[0], G0)
@@ -3418 +3416 @@
         domain.beta_vh_dry = 0.9
         domain.beta_h = 0.0 #Use first order in h-limiter
-        domain.H0 = 0 # Backwards compatibility (6/2/7)        
+        domain.H0 = 0 # Backwards compatibility (6/2/7)
+        
 
         #Bed-slope and friction at vertices (and interpolated elsewhere)