Changeset 8308

Timestamp:

Jan 18, 2012, 5:06:52 PM (12 years ago)

Author:

davies

Message:

Added extrapolate_velocity_second_order option in config.py
and getting various tests working in gareth/tests

Location:

trunk

Files:

• anuga_core/source/anuga/config.py (modified) (1 diff)
• anuga_core/source/anuga/shallow_water/shallow_water_domain.py (modified) (3 diffs)
• anuga_core/source/anuga/shallow_water/shallow_water_ext.c (modified) (7 diffs)
• anuga_work/development/gareth/tests (added)
• anuga_work/development/gareth/tests/channel_floodplain (added)
• anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py (added)
• anuga_work/development/gareth/tests/channel_floodplain/plotme.py (added)
• anuga_work/development/gareth/tests/parabolic (added)
• anuga_work/development/gareth/tests/parabolic/parabolaplot.py (added)
• anuga_work/development/gareth/tests/parabolic/parabolic.py (added)
• anuga_work/development/gareth/tests/runup (added)
• anuga_work/development/gareth/tests/runup/runup.py (added)
• anuga_work/development/gareth/tests/runup/runuplot.py (added)
• anuga_work/development/gareth/tests/runup_sinusoid (added)
• anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py (added)
• anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoidplot.py (added)
• anuga_work/development/gareth/tests/util.py (added)

trunk/anuga_core/source/anuga/config.py

@@ -95 +95 @@
 # FIXME (Ole) Maybe get rid of order altogether and use beta_w
 default_order = 2
+
+# Option to use velocity extrapolation instead of momentum extrapolation in the
+# routine domain.extrapolate_second_order_sw
+extrapolate_velocity_second_order=False
 
 ################################################################################

trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -210 +210 @@
         from anuga.config import g, beta_w, beta_w_dry, \
              beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
+        from anuga.config import extrapolate_velocity_second_order
         from anuga.config import alpha_balance
         from anuga.config import optimise_dry_cells
@@ -215 +216 @@
         from anuga.config import use_edge_limiter
         from anuga.config import use_centroid_velocities
+        
 
         self.set_minimum_allowed_height(minimum_allowed_height)
         self.maximum_allowed_speed = maximum_allowed_speed
+
+        self.extrapolate_velocity_second_order=extrapolate_velocity_second_order
 
         self.g = g
@@ -1136 +1140 @@
               Ymom.vertex_values,
               Elevation.vertex_values,
-              int(domain.optimise_dry_cells))
+              int(domain.optimise_dry_cells),
+              int(domain.extrapolate_velocity_second_order))
 
 def distribute_using_vertex_limiter(domain):

trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -1332 +1332 @@
                                  double* ymom_vertex_values,
                                  double* elevation_vertex_values,
-                                 int optimise_dry_cells) {
+                                 int optimise_dry_cells, 
+                                 int extrapolate_velocity_second_order) {
                   
                   
@@ -1343 +1344 @@
   double dqv[3], qmin, qmax, hmin, hmax;
   double hc, h0, h1, h2, beta_tmp, hfactor;
-
+  double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements], dk, dv0, dv1, dv2;
     
+  if(extrapolate_velocity_second_order==1){
+      // Replace momentum centroid with velocity centroid to allow velocity
+      // extrapolation This will be changed back at the end of the routine
+      for (k=0; k<number_of_elements; k++){
+
+          dk = max(stage_centroid_values[k]-elevation_centroid_values[k],minimum_allowed_height);
+          xmom_centroid_store[k] = xmom_centroid_values[k];
+          xmom_centroid_values[k] = xmom_centroid_values[k]/dk;
+
+          ymom_centroid_store[k] = ymom_centroid_values[k];
+          ymom_centroid_values[k] = ymom_centroid_values[k]/dk;
+      }
+  }
+
+  // Begin extrapolation routine
   for (k = 0; k < number_of_elements; k++) 
   {
@@ -1794 +1810 @@
 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];
         }
 //ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0];
@@ -1802 +1818 @@
   } // for k=0 to number_of_elements-1
   
+  if(extrapolate_velocity_second_order==1){
+      // Convert back from velocity to momentum
+      for (k=0; k<number_of_elements; k++){
+          k3=3*k;
+          //dv0 = max(stage_vertex_values[k3]-elevation_vertex_values[k3],minimum_allowed_height);
+          //dv1 = max(stage_vertex_values[k3+1]-elevation_vertex_values[k3+1],minimum_allowed_height);
+          //dv2 = max(stage_vertex_values[k3+2]-elevation_vertex_values[k3+2],minimum_allowed_height);
+          dv0 = max(stage_vertex_values[k3]-elevation_vertex_values[k3],0.);
+          dv1 = max(stage_vertex_values[k3+1]-elevation_vertex_values[k3+1],0.);
+          dv2 = max(stage_vertex_values[k3+2]-elevation_vertex_values[k3+2],0.);
+
+          //Correct centroid and vertex values
+          xmom_centroid_values[k] = xmom_centroid_store[k];
+          xmom_vertex_values[k3]=xmom_vertex_values[k3]*dv0;
+          xmom_vertex_values[k3+1]=xmom_vertex_values[k3+1]*dv1;
+          xmom_vertex_values[k3+2]=xmom_vertex_values[k3+2]*dv2;
+          
+          ymom_centroid_values[k] = ymom_centroid_store[k];
+          ymom_vertex_values[k3]=ymom_vertex_values[k3]*dv0;
+          ymom_vertex_values[k3+1]=ymom_vertex_values[k3+1]*dv1;
+          ymom_vertex_values[k3+2]=ymom_vertex_values[k3+2]*dv2;
+          
+      } 
+  }
+
   return 0;
 }           
@@ -1860 +1901 @@
   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, extrapolate_velocity_second_order;
   
   // 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,
@@ -1879 +1920 @@
             &ymom_vertex_values,
             &elevation_vertex_values,
-            &optimise_dry_cells)) {         
+            &optimise_dry_cells,
+            &extrapolate_velocity_second_order)) {         
 
     report_python_error(AT, "could not parse input arguments");
@@ -1937 +1979 @@
                    (double*) ymom_vertex_values -> data,
                    (double*) elevation_vertex_values -> data,
-                   optimise_dry_cells);
+                   optimise_dry_cells, 
+                   extrapolate_velocity_second_order);
+
   if (e == -1) {
     // Use error string set inside computational routine