Changeset 2566

Timestamp:

Mar 21, 2006, 2:25:40 PM (19 years ago)

Author:

ole

Message:

Created maximum_allowed_speed in config and passed it through.
Also ditched epsilon and used h<=0 instead.

Location:

inundation

Files:

• pyvolution/config.py (modified) (2 diffs)
• pyvolution/shallow_water.py (modified) (3 diffs)
• pyvolution/shallow_water_ext.c (modified) (7 diffs)
• pyvolution/test_shallow_water.py (modified) (1 diff)
• test_all.py (modified) (1 diff)

inundation/pyvolution/config.py

@@ -33 +33 @@
 
 
-v_max = 100 #For use in domain_ext.c
+#v_max = 100 #For use in domain_ext.c
 sound_speed = 500
 
@@ -106 +106 @@
 #Specific to shallow water W.E.
 minimum_allowed_height = 1.0e-3 #Water depth below which it is considered to be 0
+maximum_allowed_speed = 100.0 #Maximal particle speed of water

inundation/pyvolution/shallow_water.py

@@ -79 +79 @@
                                 tagged_elements, geo_reference, use_inscribed_circle)
 
-        from config import minimum_allowed_height, g
+        from config import minimum_allowed_height, maximum_allowed_speed, g
         self.minimum_allowed_height = minimum_allowed_height
+        self.maximum_allowed_speed = maximum_allowed_speed      
         self.g = g
 
@@ -667 +668 @@
     hc = wc - zc  #Water depths at centroids
 
-    #Update
+    #Update 
+    #FIXME: Modify accroditg to c-version - or discard altogether.
     for k in range(domain.number_of_elements):
 
@@ -691 +693 @@
     from shallow_water_ext import protect
 
-    protect(domain.minimum_allowed_height, domain.epsilon,
-            wc, zc, xmomc, ymomc)
+    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed, 
+            domain.epsilon, wc, zc, xmomc, ymomc)
 
 

inundation/pyvolution/shallow_water_ext.c

@@ -356 +356 @@
 int _protect(int N,
              double minimum_allowed_height,
+             double maximum_allowed_speed,           
              double epsilon,
              double* wc,
@@ -364 +365 @@
   int k;
   double hc;
-  double vx, vy;  
-  double maximum_allowed_speed = 100.0;  //FIXME (Ole): Pass in
+  double u, v, tmp_speed;  
 
   //Protect against initesimal and negative heights
@@ -373 +373 @@
 
     if (hc < minimum_allowed_height) {
-      if (hc < epsilon) {
+    //if (1) {
+      //if (hc < epsilon) {
+      if (hc <= 0.0) {
         wc[k] = zc[k]; //Contain 'lost mass' error
         xmomc[k] = 0.0;
@@ -380 +382 @@
         //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);
+        u = xmomc[k]/hc;
+        if (fabs(u) > maximum_allowed_speed) {
+          tmp_speed = maximum_allowed_speed * u/fabs(u);
+          //printf("WARNING: Speed (u) has been reduced from %.3f to %.3f\n", 
+          //     u, tmp_speed); 
+          xmomc[k] = tmp_speed * hc;
         }
         
-        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);
+        v = ymomc[k]/hc;        
+        if (fabs(v) > maximum_allowed_speed) {
+          tmp_speed = maximum_allowed_speed * v/fabs(v);        
+          //printf("WARNING: Speed (v) has been reduced from %.3f to %.3f\n",   
+          //     v, maximum_allowed_speed);     
+          ymomc[k] = tmp_speed * hc;      
         }       
         
@@ -1043 +1047 @@
 PyObject *protect(PyObject *self, PyObject *args) {
   //
-  //    protect(minimum_allowed_height, wc, zc, xmomc, ymomc)
+  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
 
 
@@ -1054 +1058 @@
 
   int N;
-  double minimum_allowed_height, epsilon;
+  double minimum_allowed_height, maximum_allowed_speed, epsilon;
 
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "ddOOOO",
+  if (!PyArg_ParseTuple(args, "dddOOOO",
                         &minimum_allowed_height,
+                        &maximum_allowed_speed,                 
                         &epsilon,
                         &wc, &zc, &xmomc, &ymomc))
@@ -1067 +1072 @@
   _protect(N,
            minimum_allowed_height,
+           maximum_allowed_speed,
            epsilon,
            (double*) wc -> data,

inundation/pyvolution/test_shallow_water.py

@@ -2174 +2174 @@
         os.remove(domain.filename + '.sww')
 
+        
+    def test_flatbed_second_order_vmax_0(self):
+        from mesh_factory import rectangular
+        from shallow_water import Domain,\
+             Reflective_boundary, Dirichlet_boundary
+
+        from Numeric import array
+
+        #Create basic mesh
+        N = 8
+        points, vertices, boundary = rectangular(N, N)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.smooth = False
+        domain.visualise = False
+        domain.default_order=2
+        domain.maximum_allowed_speed = 0.0 #Makes it like the 'oldstyle'
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain)
+        Bd = Dirichlet_boundary([0.2,0.,0.])
+
+        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
+        domain.check_integrity()
+
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
+            pass
+
+
+        assert allclose(domain.min_timestep, 0.0210448446782)
+        assert allclose(domain.max_timestep, 0.0210448446782)
+
+        #FIXME: These numbers were from version before 21/3/6 - 
+        #could be recreated by setting maximum_allowed_speed to 0 maybe 
+        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
+                        [ 0.00064835, 0.03685719, 0.00085073, 0.03687313]) 
+        
+
+        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
+                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
+
+
+        os.remove(domain.filename + '.sww')
+
+        
 
     def test_flatbed_second_order_distribute(self):

inundation/test_all.py

@@ -23 +23 @@
 exclude_dirs = ['pypar_dist', #Special requirements
                 'props', 'wcprops', 'prop-base', 'text-base', '.svn', #Svn
-                'tmp'] 
-                #'pmesh']     #Name conflict on my home machine (Ole)       
+                'tmp', 
+                'pmesh']     #Name conflict on my home machine (Ole)