Changeset 2561 for inundation/pyvolution

Timestamp:

Mar 21, 2006, 8:33:34 AM (19 years ago)

Author:

ole

Message:

Allowed particle speeds to be calculated even if depth i sless than minimum_allowed_height. A new protection, called maximum_allowed_speed, which is much less restrictive, has been suggested.

Location:

inundation/pyvolution

Files:

• island.py (modified) (4 diffs)
• shallow_water_ext.c (modified) (2 diffs)

inundation/pyvolution/island.py

@@ -13 +13 @@
      Constant_height
 from Numeric import array
-from util import Polygon_function, read_polygon
+from utilities.polygon import Polygon_function, read_polygon
 from math import exp
 
@@ -23 +23 @@
 #Create shallow water domain
 domain = Domain(points, vertices, boundary)
-domain.store = False
+domain.store = True
 domain.smooth = False
-domain.visualise = True
+domain.visualise = False
 domain.set_name('island')
 print 'Output being written to ' + domain.get_datadir() + sep + \
@@ -52 +52 @@
     return z
 
-domain.set_quantity('friction', 0.0)
+domain.set_quantity('friction', 0.1)
 domain.set_quantity('elevation', island)
 domain.set_quantity('stage', 1)
@@ -69 +69 @@
 #Evolution
 import time
-for t in domain.evolve(yieldstep = 1, finaltime = 40):
+for t in domain.evolve(yieldstep = 1, finaltime = 100):
     domain.write_time()
 

inundation/pyvolution/shallow_water_ext.c

@@ -364 +364 @@
   int k;
   double hc;
+  double vx, vy;  
+  double maximum_allowed_speed = 100.0;  //FIXME (Ole): Pass in
 
   //Protect against initesimal and negative heights
@@ -371 +373 @@
 
     if (hc < minimum_allowed_height) {
-      if (hc < epsilon) wc[k] = zc[k]; //Contain 'lost mass' error
-      xmomc[k] = 0.0;
-      ymomc[k] = 0.0;
+      if (hc < epsilon) {
+        wc[k] = zc[k]; //Contain 'lost mass' error
+        xmomc[k] = 0.0;
+        ymomc[k] = 0.0;
+      } else {
+        //Try to see if calculated speeds would blow up 
+        //FIXME (Ole): This has not been written in Python
+        vx = xmomc[k]/hc;
+        vy = ymomc[k]/hc;       
+        if (fabs(vx) > maximum_allowed_speed) {
+          printf("WARNING: Speed has %.3f has been reduced to %.3f\n", 
+                 vx, maximum_allowed_speed);    
+          xmomc[k] = maximum_allowed_speed * hc * vx/fabs(vx);
+        }
+        
+        if (fabs(vy) > maximum_allowed_speed) {
+          printf("WARNING: Speed has %.3f has been reduced to %.3f\n", 
+                 vy, maximum_allowed_speed);    
+          ymomc[k] = maximum_allowed_speed * hc * vy/fabs(vy);
+        }       
+        
+
+      }
     }