Changeset 2632

Timestamp:

Mar 29, 2006, 4:33:30 PM (17 years ago)

Author:

ole

Message:

Work on creep demo

Location:

inundation

Files:

• examples/island.py (modified) (6 diffs)
• examples/island_timeseries.py (modified) (6 diffs)
• pyvolution/shallow_water_ext.c (modified) (4 diffs)

inundation/examples/island.py

@@ -7 +7 @@
 
 
-#-------------------------------------------------------------------------------
+#------------------------------------------------------------------------------
 # Import necessary modules
-#-------------------------------------------------------------------------------
+#------------------------------------------------------------------------------
 
 # Standard modules
@@ -21 +21 @@
 
 
-#-------------------------------------------------------------------------------   
+#------------------------------------------------------------------------------
 # Setup computational domain
-#-------------------------------------------------------------------------------
+#------------------------------------------------------------------------------
 
 #Create basic mesh
@@ -37 +37 @@
                           maximum_triangle_area = 25,
                           filename = 'island.msh',
-                          interior_regions=[ ([[55,25], [75,25],
-                                               [75,75], [55,75]], 3)]
+                          interior_regions=[ ([[50,25], [70,25],
+                                               [70,75], [50,75]], 3)]
                           )
                           
@@ -63 +63 @@
 
 
-#-------------------------------------------------------------------------------
+#------------------------------------------------------------------------------
 # Setup initial conditions
-#-------------------------------------------------------------------------------
+#------------------------------------------------------------------------------
 
 def island(x, y):
@@ -73 +73 @@
     return z
 
-domain.set_quantity('friction', 0.0)
+domain.set_quantity('friction', 2.0)
 domain.set_quantity('elevation', island)
 domain.set_quantity('stage', 1)
 
 
-#-------------------------------------------------------------------------------
+#------------------------------------------------------------------------------
 # Setup boundary conditions (all reflective)
-#-------------------------------------------------------------------------------
+#------------------------------------------------------------------------------
 
 Br = Reflective_boundary(domain)
@@ -88 +88 @@
 
 
-#-------------------------------------------------------------------------------
+#------------------------------------------------------------------------------
 # Evolve system through time
-#-------------------------------------------------------------------------------
+#------------------------------------------------------------------------------
 
 import time
-for t in domain.evolve(yieldstep = 0.5, finaltime = 100):
+for t in domain.evolve(yieldstep = 1, finaltime = 100):
     domain.write_time()
+#for t in domain.evolve(yieldstep = 0.05, finaltime = 5):
+#    domain.write_time()
+#    #print '    Volume: ', domain.get_quantity('stage').get_integral()
 
 print 'Done'

inundation/examples/island_timeseries.py

@@ -10 +10 @@
 
 swwfile = 'island.sww' 
-#gauges = [[0,0.5], [0.2,0.5], [0.4,0.5], [0.6,0.5], [0.8,0.5], [1,0.5]]
-#gauges = [[20, 50], [30, 50], [40, 50], [50, 50]]
-gauges = [[42, 48], [58, 48]]     
+#gauges = [[56, 48], [58, 48], [60, 48], [62, 48], [64, 48]]
+gauges = [[60, 48], [62, 48], [64, 48]]    
 
+
+store = False
 
 #Read model output
@@ -103 +104 @@
     hold(False)
 
-    if elevations[0] < -10:
+    if elevations[0] < 0:
         plot(model_time, stages, '-b')
     else:    
         plot(model_time, stages, '-b',
              model_time, elevations, '-k')
-    axis([0, 100, z-0.01, z+0.005])
+    #axis([0, 100, z-0.01, z+0.005])
         
     #name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1])
@@ -120 +121 @@
            shadow=True,
            loc='upper right')
-    #savefig('Gauge_%d_stage' %k)
-    savefig('Gauge_%s_stage' %myloc)
+
+    if store is True: savefig('Gauge_%s_stage' %myloc)
 
     raw_input('Next')
@@ -136 +137 @@
     ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')    
     #savefig('Gauge_%d_momentum' %k)
-    savefig('Gauge_%s_momentum' %myloc)
+    if store is True: savefig('Gauge_%s_momentum' %myloc)
     
     raw_input('Next')
@@ -160 +161 @@
     ylabel(' atan(vh/uh) [degrees from North]')    
     #savefig('Gauge_%d_bearing' %k)
-    savefig('Gauge_%s_bearing' %myloc)
+    if store is True: savefig('Gauge_%s_bearing' %myloc)
     
     raw_input('Next')
@@ -176 +177 @@
     ylabel('sqrt( uh^2 + vh^2 ) / depth [m/s]')    
     #savefig('Gauge_%d_speed' %k)
-    savefig('Gauge_%s_speed' %myloc)
+    if store is True: savefig('Gauge_%s_speed' %myloc)
     
     raw_input('Next')

inundation/pyvolution/shallow_water_ext.c

@@ -246 +246 @@
       if (h >= eps) {
         S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
-        //S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
-        S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
+        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
+        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
         //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
 
@@ -372 +372 @@
 
     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;      
+      //wc[k] = zc[k]; 
+      if (hc <= 0.0) wc[k] = zc[k]; 
+      */
+
+      //New code: Adjust momentum to guarantee speeds are physical
+      //          ensure h is non negative      
+   
       if (hc <= 0.0) {
-        wc[k] = zc[k]; //Contain 'lost mass' error
+        wc[k] = zc[k]; 
         xmomc[k] = 0.0;
         ymomc[k] = 0.0;
@@ -383 +395 @@
           reduced_speed = maximum_allowed_speed * u/fabs(u);
           //printf("Speed (u) has been reduced from %.3f to %.3f\n", 
-          //        u, reduced_speed);  
+          //     u, reduced_speed);     
           xmomc[k] = reduced_speed * hc;
         }
@@ -391 +403 @@
           reduced_speed = maximum_allowed_speed * v/fabs(v);    
           //printf("Speed (v) has been reduced from %.3f to %.3f\n",   
-          //        v, reduced_speed);  
+          //     v, reduced_speed);     
           ymomc[k] = reduced_speed * hc;          
         }       
       }
-      
-      //Old code
-      //wc[k] = zc[k]; //Contain 'lost mass' error
-      //printf("Speed has been reduced to %.3f\n", 0.0);              
-      //xmomc[k] = 0.0;
-      //ymomc[k] = 0.0;      
     }
   }