Changeset 244

Timestamp:

Aug 30, 2004, 5:38:16 PM (20 years ago)

Author:

ole

Message:

 

Location:

inundation/ga/storm_surge/pyvolution

Files:

• run_profile.py (modified) (1 diff)
• shallow_water.py (modified) (4 diffs)

inundation/ga/storm_surge/pyvolution/run_profile.py

@@ -31 +31 @@
 #domain.visualise = True
 domain.newstyle = False  #For comparison
+
 
 print 'Field values'

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -410 +410 @@
     if not domain.newstyle:
         #Protect against infinitesimal heights and high speeds at_centroid
-        #FIXME: Try without when using the second balance function
+        #FIXME: To be phased out
         
         for k in range(domain.number_of_elements):
@@ -468 +468 @@
         Q.limit()        
 
-
     balance_deep_and_shallow(domain)
 
 
 def balance_deep_and_shallow(domain):
-    #FIXME: Grap math comments from old_first_order_balance....(below)
+    """Compute linear combination between stage as computed by gradient-limiters and
+    stage computed as constant height above bed.
+    The former takes precedence when heights are large compared to the bed slope while the latter
+    takes precedence when heights are relatively small.
+    Anything in between is computed as a balanced linear combination in order to avoid numerical
+    disturbances which would otherwise appear as a result of hard switching between modes.
+    """
     
     #Shortcuts
@@ -484 +489 @@
     hv = wv-zv
 
+
     #Computed linear combination between constant levels and and
     #levels parallel to the bed elevation.     
-    for k in range(domain.number_of_elements):                
-                
+    for k in range(domain.number_of_elements):
         #Compute maximal variation in bed elevation
         #  This quantitiy is
@@ -494 +499 @@
         #  In the 1d case zc = (z0+z1)/2
         #  In the 2d case zc = (z0+z1+z2)/3
-        #  
-        
-        z_range = max(abs(zv[k,0]-zc[k]),
-                      abs(zv[k,1]-zc[k]),
-                      abs(zv[k,2]-zc[k]))
-
-
-        hmin = min( hv[k, :] )
+
+        dz = max(abs(zv[k,0]-zc[k]),
+                 abs(zv[k,1]-zc[k]),
+                 abs(zv[k,2]-zc[k]))
+
+        
+        hmin = min( hv[k,:] )
 
         #Create alpha in [0,1], where alpha==0 means using shallow 
-        #first order scheme and alpha==1 means using the limited
-        #second order scheme for the stage w
-        #If hmin < 0 then alpha=0 reverting to first order.
+        #first order scheme and alpha==1 means using the stage w as
+        #computed by the gradient limiter (1st or 2nd order)
+        #
+        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
+        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
       
-        if z_range > 0.0:
-            alpha = max( min( 2*hmin/z_range, 1.0), 0.0)
-        else: 
+        if dz > 0.0:
+            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
+        else:
+            #Flat bed 
             alpha = 1.0
-   
-        #Update water levels
-        #  (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = 
-        #   zvi + hc + alpha*(hvi - hc)
+
+
+        #Weighted balance between stage parallel to bed elevation 
+        #(wvi = zvi + hc) and stage as computed by 1st or 2nd order gradient limiter
+        #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
+        #
+        #It follows that the updated wvi is
+        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = 
+        #  zvi + hc + alpha*(hvi - hc)
+        #
+        #Note that hvi = zc+hc-zvi in the first order case (constant).
+
         if alpha < 1:         
             for i in range(3):