Changeset 5315

Timestamp:

May 13, 2008, 6:08:40 PM (16 years ago)

Author:

ole

Message:

Rolled back changes from changeset:5297 regarding ticket:262.
All tests and validations passes.

Location:

anuga_core

Files:

• documentation/user_manual/demos/channel3.py (modified) (1 diff)
• source/anuga/config.py (modified) (1 diff)
• source/anuga/shallow_water/shallow_water_domain.py (modified) (1 diff)
• source/anuga/shallow_water/shallow_water_ext.c (modified) (11 diffs)

anuga_core/documentation/user_manual/demos/channel3.py

@@ -74 +74 @@
 # Evolve system through time
 #------------------------------------------------------------------------------
-for t in domain.evolve(yieldstep = 0.2, finaltime = 16.0):
+for t in domain.evolve(yieldstep = 0.1, finaltime = 16.0):
     print domain.timestepping_statistics()
 

anuga_core/source/anuga/config.py

@@ -88 +88 @@
 tight_slope_limiters = True
 
-# Use centroid velocities to reconstruct momentum at vertices in very shallow water
+# Use centroid velocities to reconstruct momentum at vertices in
+# very shallow water
 # This option has a first order flavour to it, but we still have second order
-# reconstruction of stage and this option only applies in balance_deep_and_shallow when
+# reconstruction of stage and this option only applies in
+# balance_deep_and_shallow when
 # alpha < 1 so in deeper water the full second order scheme is used.
 #

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -732 +732 @@
               Ymom.vertex_values,
               Elevation.vertex_values,
-              int(domain.optimise_dry_cells),
-              int(domain.use_centroid_velocities))
+              int(domain.optimise_dry_cells))
 
 

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -550 +550 @@
                               double H0,
                               int tight_slope_limiters,
-                              int use_centroid_velocities,                    
+                              int use_centroid_velocities,
                               double alpha_balance) {
 
-  int k, k3, i; //, excessive_froude_number=0;
+  int k, k3, i; 
 
   double dz, hmin, alpha, h_diff, hc_k;
   double epsilon = 1.0e-6; // FIXME: Temporary measure
   double hv[3]; // Depths at vertices
-  //double Fx, Fy; // Froude numbers
   double uc, vc; // Centroid speeds
 
@@ -651 +650 @@
         alpha = 1.0;       
       }
-
     }
-            
+
 
     //  Let
@@ -684 +682 @@
 
         // Update momentum at vertices
-        if (use_centroid_velocities == 0) {
+        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
@@ -691 +709 @@
           // 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];
@@ -697 +719 @@
         }
       }
-    } // If alpha == 1 use quantities as calculated by the gradient-limiter
-    
-    if (use_centroid_velocities == 1) {    
-      // This is a simple, efficient and robust option in shallow water
-      // It uses first order approximation of velocities, but retains
-      // the order used by stage.
-      // In this case the xmomv and ymomv calculations should be switched off
-      // in the gradient limiter.
-
-      for (i=0; i<3; i++) {    
-        
-        // 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];
-      } 
     }
   }
   return 0;
 }
+  
 
 
@@ -1020 +1017 @@
                                   double* ymom_vertex_values,
                                   double* elevation_vertex_values,
-                                  int optimise_dry_cells,
-                                  int use_centroid_velocities) {
+                                  int optimise_dry_cells) {
                                   
                                   
@@ -1027 +1023 @@
   // Local variables
   double a, b; // Gradient vector used to calculate vertex values from centroids
-  int k,k0,k1,k2,k3,k6,coord_index, i;
+  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;
@@ -1192 +1188 @@
       for (i=0;i<3;i++)
         stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
-      
-      
-      if (use_centroid_velocities == 1) {
-        // Use first order reconstruction using speeds only.
-        
-        // This happens in balance_deep_and_shallow so there
-        // is no need to do more here
-        
-        // FIXME (Ole): Optionally we could put the computation here but then
-        // it'd go into all other variants of the gradient limiter  
-        continue;
-      }
       
       
@@ -1500 +1484 @@
   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, use_centroid_velocities;
+  int optimise_dry_cells, number_of_elements, e;
   
   // 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, "OOOOOOOOOOOOOii",
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
                         &domain,
                         &surrogate_neighbours,
@@ -1519 +1503 @@
                         &ymom_vertex_values,
                         &elevation_vertex_values,
-                        &optimise_dry_cells,
-                        &use_centroid_velocities)) {                    
+                        &optimise_dry_cells)) {                 
                         
     PyErr_SetString(PyExc_RuntimeError, 
@@ -1565 +1548 @@
                                    (double*) ymom_vertex_values -> data,
                                    (double*) elevation_vertex_values -> data,
-                                   optimise_dry_cells,
-                                   use_centroid_velocities);
+                                   optimise_dry_cells);
   if (e == -1) {
     // Use error string set inside computational routine