Changeset 215

Timestamp:

Aug 24, 2004, 6:11:26 PM (21 years ago)

Author:

ole

Message:

More testing of gravity, forcing terms, distributors

Location:

inundation/ga/storm_surge/pyvolution

Files:

• domain.py (modified) (1 diff)
• quantity.py (modified) (1 diff)
• shallow_water.py (modified) (3 diffs)
• test_domain.py (modified) (3 diffs)
• test_shallow_water.py (modified) (4 diffs)

inundation/ga/storm_surge/pyvolution/domain.py

@@ -438 +438 @@
         timestep = self.timestep
 
-        #Reset semi_implicit_updates
-        #(explicit updates are already set by compute_fluyzes)
-        for i, name in enumerate(self.conserved_quantities):
-            Q = self.quantities[name]
-            Q.semi_implicit_update[:] = 0.
-
         #Compute forcing terms    
         self.compute_forcing_terms()    
 
         #Update conserved_quantities from explicit updates
-        for i, name in enumerate(self.conserved_quantities):
+        for name in self.conserved_quantities:
             Q = self.quantities[name]
             Q.update(timestep)
+
+            #Clean up
+            Q.semi_implicit_update[:] = 0.0
+            Q.explicit_update[:] = 0.0 #Not necessary as fluxes will set it
             
 

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -328 +328 @@
             msg = 'Invalid location: %s' %location
             raise msg
+
+        if X is None:
+            msg = 'Given values are None'
+            raise msg            
         
         import types

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -586 +586 @@
     for k in range(domain.number_of_elements):    
         avg_h = sum( h[k,:] )/3
-
+        
         #Compute bed slope
         x0, y0, x1, y1, x2, y2 = V[k,:]    
         z0, z1, z2 = Elevation.vertex_values[k,:]
+
         zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
 
@@ -595 +596 @@
         Xmom.explicit_update[k] += -g*zx*avg_h
         Ymom.explicit_update[k] += -g*zy*avg_h        
-    
+
 
 def manning_friction(domain):
@@ -627 +628 @@
 
 ###########################
+###########################
+#Geometries
+
+
+#FIXME: Rethink this way of creating values.
+
+
+class Weir:
+    """Set a bathymetry for weir with a hole and a downstream gutter
+    x,y are assumed to be in the unit square
+    """
+    
+    def __init__(self, stage):
+        self.inflow_stage = stage
+
+    def __call__(self, x, y):    
+        from Numeric import zeros, Float
+        from math import sqrt
+    
+        N = len(x)
+        assert N == len(y)    
+
+        z = zeros(N, Float)
+        for i in range(N):
+            z[i] = -x[i]/2  #General slope
+
+            #Flattish bit to the left
+            if x[i] < 0.3:
+                z[i] = -x[i]/10 
+            
+            #Weir
+            if x[i] >= 0.3 and x[i] < 0.4:
+                z[i] = -x[i]+0.9
+
+            #Dip
+            x0 = 0.6
+            #depth = -1.3
+            depth = -1.0
+            #plateaux = -0.9
+            plateaux = -0.6            
+            if y[i] < 0.7:
+                if x[i] > x0 and x[i] < 0.9:
+                    z[i] = depth
+
+                #RHS plateaux
+                if x[i] >= 0.9:
+                    z[i] = plateaux
+                    
+                    
+            elif y[i] >= 0.7 and y[i] < 1.5:
+                #Restrict and deepen
+                if x[i] >= x0 and x[i] < 0.8:
+                    z[i] = depth-(y[i]/3-0.3)
+                    #z[i] = depth-y[i]/5
+                    #z[i] = depth
+                elif x[i] >= 0.8:    
+                    #RHS plateaux
+                    z[i] = plateaux
+                    
+            elif y[i] >= 1.5:
+                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
+                    #Widen up and stay at constant depth
+                    z[i] = depth-1.5/5
+                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:                        
+                    #RHS plateaux
+                    z[i] = plateaux                                    
+                    
+
+            #Hole in weir (slightly higher than inflow condition)
+            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
+                z[i] = -x[i]+self.inflow_stage + 0.02
+
+            #Channel behind weir
+            x0 = 0.5
+            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
+                z[i] = -x[i]+self.inflow_stage + 0.02
+                
+            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
+                #Flatten it out towards the end
+                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
+
+            #Hole to the east
+            x0 = 1.1; y0 = 0.35
+            #if x[i] < -0.2 and y < 0.5:
+            if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
+                z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
+
+            #Tiny channel draining hole
+            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
+                z[i] = -0.9 #North south
+                
+            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
+                z[i] = -1.0 + (x[i]-0.9)/3 #East west
+            
+            
+
+            #Stuff not in use
+            
+            #Upward slope at inlet to the north west
+            #if x[i] < 0.0: # and y[i] > 0.5:
+            #    #z[i] = -y[i]+0.5  #-x[i]/2
+            #    z[i] = x[i]/4 - y[i]**2 + 0.5
+
+            #Hole to the west
+            #x0 = -0.4; y0 = 0.35 # center
+            #if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
+            #    z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2
+
+
+
+
+
+        return z/2
+
+class Weir_simple:
+    """Set a bathymetry for weir with a hole and a downstream gutter
+    x,y are assumed to be in the unit square
+    """
+    
+    def __init__(self, stage):
+        self.inflow_stage = stage
+
+    def __call__(self, x, y):    
+        from Numeric import zeros, Float 
+    
+        N = len(x)
+        assert N == len(y)    
+
+        z = zeros(N, Float)
+        for i in range(N):
+            z[i] = -x[i]  #General slope
+
+            #Flat bit to the left
+            if x[i] < 0.3:
+                z[i] = -x[i]/10  #General slope
+            
+            #Weir
+            if x[i] > 0.3 and x[i] < 0.4:
+                z[i] = -x[i]+0.9
+
+            #Dip
+            if x[i] > 0.6 and x[i] < 0.9:
+                z[i] = -x[i]-0.5  #-y[i]/5
+
+            #Hole in weir (slightly higher than inflow condition)
+            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
+                z[i] = -x[i]+self.inflow_stage + 0.05       
+            
+
+        return z/2
+        
+
+            
 class Constant_height:
     """Set an initial condition with constant water height, e.g
     stage s = z+h
     """
-
-    #FIXME: Rethink this way of creating values.
-
     def __init__(self, W, h):
         self.W = W

inundation/ga/storm_surge/pyvolution/test_domain.py

@@ -165 +165 @@
 
 
-    def test_distribute(self):
+    def test_distribute_first_order(self):
         """Domain implements a default first order gradient limiter
         """
@@ -206 +206 @@
 
 
-        #print domain.quantities['level'].centroid_values
         domain.distribute_to_vertices_and_edges()
 
@@ -244 +243 @@
 
 
-        domain.set_quantity('level', [[1,2,3], [5,5,5],
-                                      [0,0,9], [-6, 3, 3]])
-
-        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
-                                          [3,3,3], [4, 4, 4]])
-
-        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
-                                          [30,30,30], [40, 40, 40]])        
-
+        domain.set_quantity('level', [1,2,3,4], 'centroids')
+        domain.set_quantity('xmomentum', [1,2,3,4], 'centroids')
+        domain.set_quantity('ymomentum', [1,2,3,4], 'centroids')
+                                          
 
         #Assign some values to update vectors
-        #domain.explicit_updates[]
+        #Set explicit_update
+
+        for name in domain.conserved_quantities:
+            domain.quantities[name].explicit_update = array([4.,3.,2.,1.])
+            domain.quantities[name].semi_implicit_update = array([1.,1.,1.,1.])
+            
+
+        #Update with given timestep (assuming no other forcing terms)
+        domain.timestep = 0.1
+        domain.update_conserved_quantities()
+
+        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
+        x /= array( [.9,.9,.9,.9] )
+
+        for name in domain.conserved_quantities:
+            assert allclose(domain.quantities[name].centroid_values, x)        
+
         
         

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -534 +534 @@
             pass
 
-
-
+    #####################################################    
+    def test_gravity(self):
+        #Assuming no friction
+
+        from config import g
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+        
+        domain = Domain(points, vertices)
+
+        #Set up for a gradient of (3,0) at mid triangle         
+        def slope(x, y):
+            return 3*x
+
+        h = 0.1
+        def level(x,y):
+            return slope(x,y)+h
+
+        domain.set_quantity('elevation', slope)
+        domain.set_quantity('level', level)
+
+        for name in domain.conserved_quantities:
+            assert allclose(domain.quantities[name].explicit_update, 0)
+            assert allclose(domain.quantities[name].semi_implicit_update, 0)
+            
+        domain.compute_forcing_terms()                                    
+
+        assert allclose(domain.quantities['level'].explicit_update, 0)
+        assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
+        assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
+        
+        
+
+    #####################################################
     def test_first_order_extrapolator_const_z(self):
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+        
+        domain = Domain(points, vertices)
+        val0 = 2.+2.0/3
+        val1 = 4.+4.0/3
+        val2 = 8.+2.0/3
+        val3 = 2.+8.0/3
+
+        zl=zr=-3.75 #Assume constant bed (must be less than level)        
+        domain.set_quantity('elevation', zl*ones( (4,3) ))
+        domain.set_quantity('level', [[val0, val0-1, val0-2],
+                                      [val1, val1+1, val1],
+                                      [val2, val2-2, val2],
+                                      [val3-0.5, val3, val3]])
+
+
+
+        domain.order = 1
+        domain.distribute_to_vertices_and_edges()
+
+        #Check that centroid values were distributed to vertices
+        C = domain.quantities['level'].centroid_values
+        for i in range(3): 
+            assert allclose( domain.quantities['level'].vertex_values[:,i], C)
+        
+
+    def test_first_order_limiter_variable_z(self):
+        #Check that first order limiter follows bed_slope
+        from Numeric import alltrue, greater_equal
+        from config import epsilon
         
         a = [0.0, 0.0]
@@ -555 +636 @@
         val3 = 2.+8.0/3
 
-        zl=zr=-3.75 #Assume constant bed (must be less than level)        
-        domain.set_quantity('elevation', zl*ones( (4,3) ))
-        domain.set_quantity('level', [[val0, val0-1, val0-2],
-                                      [val1, val1+1, val1],
-                                      [val2, val2-2, val2],
-                                      [val3-0.5, val3, val3]])
-
-
-
-        domain.order = 1
-        domain.distribute_to_vertices_and_edges()
-
-        #Check that centroid values were distributed to vertices
-        C = domain.quantities['level'].centroid_values
-        for i in range(3): 
-            assert allclose( domain.quantities['level'].vertex_values[:,i], C)
-        
-
-    def test_first_order_limiter_variable_z(self):
-        #Check that first order limiter follows bed_slope
-        from Numeric import alltrue, greater_equal
-        from config import epsilon
-        
-        a = [0.0, 0.0]
-        b = [0.0, 2.0]
-        c = [2.0,0.0]
-        d = [0.0, 4.0]
-        e = [2.0, 2.0]
-        f = [4.0,0.0]
-
-        points = [a, b, c, d, e, f]
-        #bac, bce, ecf, dbe
-        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
-        
-        domain = Domain(points, vertices)
-        val0 = 2.+2.0/3
-        val1 = 4.+4.0/3
-        val2 = 8.+2.0/3
-        val3 = 2.+8.0/3
-
         domain.set_quantity('elevation', [[0,0,0], [6,0,0],
                                           [6,6,6], [6,6,6]])
@@ -616 +657 @@
         assert alltrue(alltrue(greater_equal(L,E-epsilon)))
 
+
+    #####################################################    
+    def test_distribute_near_bed(self):
+        #Assuming no friction
+
+        from config import g
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+        
+        domain = Domain(points, vertices)
+        domain.order = 2
+
+        #Set up for a gradient of (10,0) at mid triangle         
+        def slope(x, y):
+            return 10*x
+
+        h = 0.1
+        def level(x,y):
+            return slope(x,y)+h
+
+        domain.set_quantity('elevation', slope)
+        domain.set_quantity('level', level, 'centroids')
+
+        E = domain.quantities['elevation'].vertex_values
+        L = domain.quantities['level'].vertex_values
+
+        #print E
+        domain.distribute_to_vertices_and_edges()
+        #print L
+
+        #FIXME: HERTIL 24/8/4. Write similar test using pyvolution2 and
+        #get test data from that
 
         
@@ -850 +932 @@
         domain.default_order=2
 
-        #Bed-slope and friction
+        #Bed-slope and friction at vertices (and interpolated elsewhere)
         def x_slope(x, y):
             return -x/3