Changeset 901 for inundation/ga/storm_surge/pyvolution

Timestamp:

Feb 16, 2005, 10:18:23 PM (20 years ago)

Author:

ole

Message:

First experiment with new limiter (near bed)

Location:

inundation/ga/storm_surge/pyvolution

Files:

• flatbed.py (modified) (2 diffs)
• netherlands.py (modified) (1 diff)
• shallow_water.py (modified) (5 diffs)
• test_shallow_water.py (modified) (1 diff)
• twolevels.py (added)
• util.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/flatbed.py

@@ -14 +14 @@
 from Numeric import array
 from util import Polygon_function, read_polygon
-from config import data_dir
+
 
 #Create basic mesh
@@ -24 +24 @@
 domain.store = True
 domain.set_name('polygons')
-print "Output being written to " + data_dir + sep + \
+print "Output being written to " + domain.get_datadir() + sep + \
               domain.filename + "_size%d." %len(domain) + domain.format
         

inundation/ga/storm_surge/pyvolution/netherlands.py

@@ -97 +97 @@
 N = 60
 
-N = 140
+#N = 140
+
+N = 15
 
 print 'Creating domain'

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -595 +595 @@
     
     
-def balance_deep_and_shallow(domain):
+def balance_deep_and_shallow_old(domain):
     """Compute linear combination between stage as computed by
     gradient-limiters and stage computed as constant height above bed.
@@ -633 +633 @@
         hmin = min( hv[k,:] )
 
-        
         #Create alpha in [0,1], where alpha==0 means using shallow 
         #first order scheme and alpha==1 means using the stage w as
@@ -647 +646 @@
             alpha = 1.0
 
-
         #Weighted balance between stage parallel to bed elevation 
         #(wvi = zvi + hc) and stage as computed by 1st or 2nd 
@@ -664 +662 @@
             
 
+            #Momentums at centroids
+            xmomc = domain.quantities['xmomentum'].centroid_values
+            ymomc = domain.quantities['ymomentum'].centroid_values        
+
+            #Momentums at vertices
+            xmomv = domain.quantities['xmomentum'].vertex_values
+            ymomv = domain.quantities['ymomentum'].vertex_values         
+
+            # Update momentum as a linear combination of 
+            # xmomc and ymomc (shallow) and momentum
+            # from extrapolator xmomv and ymomv (deep).
+            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]            
+            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
+            
+
+def limit_on_h(domain):        
+    """Limit slopes for each volume to eliminate artificial variance
+    introduced by e.g. second order extrapolator
+    
+    limit on h = w-z
+    """
+
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+    beta = domain.beta    
+    
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+    
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv-zv
+
+    hvbar = zeros(hv.shape, Float) #h-limited values 
+    
+    #Find min and max of this and neighbour's centroid values
+    hmax = zeros(hc.shape, Float)
+    hmin = zeros(hc.shape, Float)        
+    
+    for k in range(N):
+        hmax[k] = hmin[k] = hc[k]
+        for i in range(3):
+            n = domain.neighbours[k,i]
+            if n >= 0:
+                hn = hc[n] #Neighbour's centroid value
+                
+                hmin[k] = min(hmin[k], hn)
+                hmax[k] = max(hmax[k], hn)
+      
+    
+    #Diffences between centroids and maxima/minima
+    dhmax = hmax - hc 
+    dhmin = hmin - hc
+        
+    #Deltas between vertex and centroid values
+    dh = zeros(hv.shape, Float)
+    for i in range(3):
+        dh[:,i] = hv[:,i] - hc
+
+    #Phi limiter    
+    for k in range(N):
+        
+        #Find the gradient limiter (phi) across vertices
+        phi = 1.0
+        for i in range(3):
+            r = 1.0
+            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
+            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
+            
+            phi = min( min(r*beta, 1), phi )    
+
+        #Then update using phi limiter
+        for i in range(3):            
+            hvbar[k,i] = hc[k] + phi*dh[k,i]
+
+    return hvbar
+    
+                                            
+def balance_deep_and_shallow(domain):
+    """Compute linear combination between stage as computed by
+    gradient-limiters limiting using w, and stage computed as 
+    constant height above bed and limited using h.
+    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.
+    """
+
+    #New idea.
+    #
+    # In the presence and near of bedslope it is necessary to 
+    # limit slopes based on differences in h rather than w 
+    # (which is what is needed away from the bed).
+    #
+    # So whether extrapolation was first order or second order,
+    # it will need to be balanced with a h-limited gradient.
+    #
+    # For this we will use the quantity alpha as before 
+    #
+        
+    #Shortcuts
+    wc = domain.quantities['stage'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+    
+    wv = domain.quantities['stage'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv-zv
+
+    for k in range(domain.number_of_elements):
+        #Compute maximal variation in bed elevation
+        #  This quantitiy is
+        #    dz = max_i abs(z_i - z_c)
+        #  and it is independent of dimension
+        #  In the 1d case zc = (z0+z1)/2
+        #  In the 2d case zc = (z0+z1+z2)/3
+
+        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 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 dz > 0.0:
+            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
+        else:
+            #Flat bed 
+            alpha = 1.0
+
+        #Let 
+        # 
+        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
+        #  wvi~ be the w-limited stage (wvi~ = zvi + hvi~)
+        #  
+        #    
+        #where i=0,1,2 denotes the vertex ids
+        #
+        #Weighted balance between w-limited and h-limited stage is  
+        #
+        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi~)   
+        #
+        #It follows that the updated wvi is
+        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi~
+        #
+        # Momentum is balanced between constant and limited ???
+
+        #print 'Alpha = ', alpha
+        #FIXME: if alpha == 0 compute only h-limited slope, 
+        #if alpha == 1 compute only w-limited slope
+        #
+        if alpha < 1:         
+        
+            #Limit h
+            hvbar = limit_on_h(domain)
+    
+            for i in range(3):
+                wv[k,i] = zv[k,i] + (1-alpha)*hvbar[k,i] + alpha*hv[k,i]
+            
             #Momentums at centroids
             xmomc = domain.quantities['xmomentum'].centroid_values
@@ -1289 +1458 @@
     manning_friction = manning_friction_c
     balance_deep_and_shallow = balance_deep_and_shallow_c
+    #balance_deep_and_shallow = balance_deep_and_shallow_old
     protect_against_infinitesimal_and_negative_heights = protect_against_infinitesimal_and_negative_heights_c
 

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -125 +125 @@
         assert allclose(flux, [0., 0.5*g*h**2, 0.])
         assert max_speed == sqrt(g*h)
+        
+    #def test_flux_slope(self):
+    #    #FIXME: TODO
+    #    w = 2.0
+    #   
+    #    normal = array([1.,0])        
+    #    ql = array([w, 0, 0])
+    #    qr = array([w, 0, 0])
+    #    zl = zr = 0.
+    #    h = w - (zl+zr)/2
+    #    
+    #    flux, max_speed = flux_function(normal, ql, qr, zl, zr)
+    #
+    #    assert allclose(flux, [0., 0.5*g*h**2, 0.])
+    #    assert max_speed == sqrt(g*h)
         
 

inundation/ga/storm_surge/pyvolution/util.py

@@ -263 +263 @@
             for i, t in enumerate(self.T):
                 #Interpolate quantities at this timestep
-                print ' time step %d of %d' %(i, len(self.T))
+                #print ' time step %d of %d' %(i, len(self.T))
                 for name in self.quantities:
                     self.values[name][i, :] =\