Changeset 5290 for anuga_core/source/anuga/shallow_water/shallow_water_ext.c

Timestamp:

May 7, 2008, 5:40:11 PM (17 years ago)

Author:

ole

Message:

Refactored options, so that the use of centroid_velocities is a separate
option to tight_slope_limiters.

File:

• anuga_core/source/anuga/shallow_water/shallow_water_ext.c (modified) (16 diffs)

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -144 +144 @@
                           
   // Compute Froude number; v/sqrt(gh)
-  
+  // FIXME (Ole): Not currently in use
+    
   double froude_number;
   
@@ -549 +550 @@
                               double H0,
                               int tight_slope_limiters,
+                              int use_centroid_velocities,                            
                               double alpha_balance) {
 
-  int k, k3, i, excessive_froude_number=0;
+  int k, k3, i; //, excessive_froude_number=0;
 
   double dz, hmin, alpha, h_diff, hc_k;
   double epsilon = 1.0e-6; // FIXME: Temporary measure
-  double g = 9.81; // FIXME: Temporary measure
-  double hv[3], h; // Depths at vertices
-  double Fx, Fy; // Froude numbers
+  double hv[3]; // Depths at vertices
+  //double Fx, Fy; // Froude numbers
   double uc, vc; // Centroid speeds
 
@@ -701 +702 @@
     if (alpha < 1) {      
       for (i=0; i<3; i++) {
+      
         // FIXME (Ole): Simplify when (if) hvbar gets retired       
         if (beta_h > epsilon) {   
@@ -709 +711 @@
 
         // Update momentum at vertices
-
-
-        // Update momentum at vertices
-        // Update momentum as a linear combination of
-        // xmomc and ymomc (shallow) and momentum
-        // from extrapolator xmomv and ymomv (deep).
-
-
-        //xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
-        //ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
-
-        if (tight_slope_limiters == 1) {
-          // FIXME(Ole): Here's what I think (as of 17 Nov 2007)
-          // we need to do. Simple and efficient:
+        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
@@ -738 +730 @@
           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).
-          // FIXME (Ole): Is this really needed? Could we use the above
-          // instead?
+          // This assumes that values from xmomv and ymomv have
+          // been established e.g. by the gradient limiter.
+
           
           xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
@@ -749 +743 @@
         
         }
-      }
-    }
-        
-
-        
-    if (0) {  // FIXME(Ole): Disabled while testing balancing of velocities above
-      //if (tight_slope_limiters == 1) {                
-    
-      // Ensure that the Froude number is kept realistic at vertices
-      // FIXME (Ole): I think it could be used to adjust alpha down
-      // whenever Fr is too large. Possible make sure it doesn't deviate 
-      // too much from the value at the centroid. I like this idea!
-      
-      // FIXME (Ole): currently only used with tights_SL
-      
-      // FIXME (Ole): may not be necessary now
-
-      excessive_froude_number=0;    
-      for (i=0; i<3; i++) {    
-      
-        // Recalculate depth at vertex i
-        h = wv[k3+i] - zv[k3+i];
-        
-        Fx = compute_froude_number(xmomv[k3+i], h, g, epsilon);
-        Fy = compute_froude_number(ymomv[k3+i], h, g, epsilon);
-        
-        if ( (fabs(Fx) > 100.0) || (fabs(Fy) > 100.0)) {
-          // FIXME: use max_froude - or base it on centroid value of F
-          // printf("Excessive Froude number detected: %f or %f\n", Fx, Fy);
-          excessive_froude_number=1;          
-        }
-      }
-     
-      if (excessive_froude_number) { 
-    
-        // printf("Adjusting momentum to first order.\n"); 
-        // Go back to first order (alpha = 0) for this triangle
-        for (i=0; i<3; i++) {          
-          xmomv[k3+i] = xmomc[k];
-          ymomv[k3+i] = ymomc[k];
-          
-          if (beta_h > epsilon) {         
-            wv[k3+i] = zv[k3+i] + hvbar[k3+i];
-          } else {
-            wv[k3+i] = zv[k3+i] + hc_k; 
-          }
-        }       
       }
     }
@@ -1094 +1041 @@
                                   double* ymom_vertex_values,
                                   double* elevation_vertex_values,
-                                  int optimise_dry_cells) {
+                                  int optimise_dry_cells,
+                                  int use_centroid_velocities) {
                                   
                                   
@@ -1359 +1307 @@
       // One internal neighbour and gradient is in direction of the neighbour's centroid
       
-      // Find the only internal neighbour
+      // Find the only internal neighbour (k1?)
       for (k2=k3;k2<k3+3;k2++){
         // Find internal neighbour of triangle k      
@@ -1367 +1315 @@
           break;
       }
+      
       if ((k2==k3+3)) {
         // If we didn't find an internal neighbour
@@ -1396 +1345 @@
       dy2=dx2*dy1;
       dx2*=dx1;
-      
       
       
@@ -1561 +1509 @@
   double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry;    
   double minimum_allowed_height, epsilon;
-  int optimise_dry_cells, number_of_elements, e;
+  int optimise_dry_cells, number_of_elements, e, use_centroid_velocities;
   
   // Provisional jumps from centroids to v'tices and safety factor re limiting
   // by which these jumps are limited
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOii",
                         &domain,
                         &surrogate_neighbours,
@@ -1580 +1528 @@
                         &ymom_vertex_values,
                         &elevation_vertex_values,
-                        &optimise_dry_cells)) {                 
+                        &optimise_dry_cells,
+                        &use_centroid_velocities)) {                    
                         
     PyErr_SetString(PyExc_RuntimeError, 
@@ -1625 +1574 @@
                                    (double*) ymom_vertex_values -> data,
                                    (double*) elevation_vertex_values -> data,
-                                   optimise_dry_cells);
+                                   optimise_dry_cells,
+                                   use_centroid_velocities);
   if (e == -1) {
     // Use error string set inside computational routine
@@ -2271 +2221 @@
   double H0, beta_h;
 
-  int N, tight_slope_limiters; //, err;
+  int N, tight_slope_limiters, use_centroid_velocities; //, err;
 
   // Convert Python arguments to C
@@ -2318 +2268 @@
   
   
+  Tmp = PyObject_GetAttrString(domain, "use_centroid_velocities");
+  if (!Tmp) {
+    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object use_centroid_velocities from domain");
+    return NULL;
+  }  
+  use_centroid_velocities = PyInt_AsLong(Tmp);
+  Py_DECREF(Tmp);
+  
+  
+  
   N = wc -> dimensions[0];
   _balance_deep_and_shallow(N,
@@ -2332 +2292 @@
                             H0,
                             (int) tight_slope_limiters,
+                            (int) use_centroid_velocities,                          
                             alpha_balance);