Changeset 4731

Timestamp:

Sep 12, 2007, 1:05:19 PM (18 years ago)

Author:

ole

Message:

More small optimisations

Location:

anuga_core/source/anuga/shallow_water

Files:

• shallow_water_ext.c (modified) (4 diffs)
• test_shallow_water_domain.py (modified) (4 diffs)

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -280 +280 @@
   // Flux computation
   denom = s_max-s_min;
-  if (denom == 0.0) {  // FIXME (Ole): Try using epsilon here
+  if (denom < epsilon) { // FIXME (Ole): Try using H0 here
     for (i=0; i<3; i++) edgeflux[i] = 0.0;
     *max_speed = 0.0;
@@ -623 +623 @@
   double u, v, reduced_speed;
 
-  // Protect against initesimal and negative heights
-  for (k=0; k<N; k++) {
-    hc = wc[k] - zc[k];
-
-    if (hc < minimum_allowed_height) {
+
+  // Protect against initesimal and negative heights  
+  if (maximum_allowed_speed < epsilon) {
+    for (k=0; k<N; k++) {
+      hc = wc[k] - zc[k];
+
+      if (hc < minimum_allowed_height) {
         
-      //Old code: Set momentum to zero and ensure h is non negative
-      //xmomc[k] = 0.0;
-      //ymomc[k] = 0.0;
-      //if (hc <= 0.0) wc[k] = zc[k];
-
-
-      //New code: Adjust momentum to guarantee speeds are physical
-      //          ensure h is non negative
-      //FIXME (Ole): This is only implemented in this C extension and
-      //             has no Python equivalent
-            
-      if (hc <= 0.0) {
-        wc[k] = zc[k];
+        // Set momentum to zero and ensure h is non negative
         xmomc[k] = 0.0;
         ymomc[k] = 0.0;
-      } 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;
-        if (fabs(u) > maximum_allowed_speed) {
-          reduced_speed = maximum_allowed_speed * u/fabs(u);
-          //printf("Speed (u) has been reduced from %.3f to %.3f\n",
-          //     u, reduced_speed);
-          xmomc[k] = reduced_speed * hc;
-        }
-
-        v = ymomc[k]/hc;
-        if (fabs(v) > maximum_allowed_speed) {
-          reduced_speed = maximum_allowed_speed * v/fabs(v);
-          //printf("Speed (v) has been reduced from %.3f to %.3f\n",
-          //     v, reduced_speed);
-          ymomc[k] = reduced_speed * hc;
-        }
+        if (hc <= 0.0) wc[k] = zc[k];
+      }
+    }
+  } else {
+   
+    // Protect against initesimal and negative heights
+    for (k=0; k<N; k++) {
+      hc = wc[k] - zc[k];
+
+      if (hc < minimum_allowed_height) {
+
+        //New code: Adjust momentum to guarantee speeds are physical
+        //          ensure h is non negative
+        //FIXME (Ole): This is only implemented in this C extension and
+        //             has no Python equivalent
+              
+        if (hc <= 0.0) {
+                wc[k] = zc[k];
+        xmomc[k] = 0.0;
+        ymomc[k] = 0.0;
+        } 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;
+          if (fabs(u) > maximum_allowed_speed) {
+            reduced_speed = maximum_allowed_speed * u/fabs(u);
+            //printf("Speed (u) has been reduced from %.3f to %.3f\n",
+            //   u, reduced_speed);
+            xmomc[k] = reduced_speed * hc;
+          }
+
+          v = ymomc[k]/hc;
+          if (fabs(v) > maximum_allowed_speed) {
+            reduced_speed = maximum_allowed_speed * v/fabs(v);
+            //printf("Speed (v) has been reduced from %.3f to %.3f\n",
+            //   v, reduced_speed);
+            ymomc[k] = reduced_speed * hc;
+          }
+        }
       }
     }
@@ -669 +680 @@
   return 0;
 }
-
 
 
@@ -680 +690 @@
                               double cw) {
 
-  //Assign windfield values to momentum updates
+  // Assign windfield values to momentum updates
 
   int k;

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -873 +873 @@
         N = 5
         points, vertices, boundary = rectangular_cross(N, N) 
-        domain = Domain(points, vertices, boundary)            
+        domain = Domain(points, vertices, boundary)
+        domain.set_maximum_allowed_speed(1.0)        
 
         #--------------------------------------------------------------
@@ -1020 +1021 @@
         domain = Domain(points, vertices, boundary)
         domain.set_name('runup_test')
+        domain.set_maximum_allowed_speed(1.0)
         #domain.tight_slope_limiters = 1 #FIXME: This works better with old limiters
 
@@ -2989 +2991 @@
         domain.beta_vh_dry = 0.9        
         #domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
-        domain.H0 = 0 # Backwards compatibility (6/2/7)        
+        domain.H0 = 0 # Backwards compatibility (6/2/7)
+        domain.set_maximum_allowed_speed(1.0)        
 
         # Boundary conditions
@@ -3104 +3107 @@
         domain.beta_vh     = 0.9
         domain.beta_vh_dry = 0.9
-        domain.H0 = 0 # Backwards compatibility (6/2/7)        
+        domain.H0 = 0 # Backwards compatibility (6/2/7)
+        domain.set_maximum_allowed_speed(1.0)        
 
         # Boundary conditions