Changeset 212

Timestamp:

Aug 24, 2004, 4:11:25 PM (21 years ago)

Author:

ole

Message:

Testing

Location:

inundation/ga/storm_surge/pyvolution

Files:

• domain.py (modified) (5 diffs)
• shallow_water.py (modified) (9 diffs)
• test_quantity.py (modified) (1 diff)
• test_shallow_water.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/domain.py

@@ -38 +38 @@
 
 
+        #FIXME: Move these explanations elsewhere    
+            
         #Create an empty list for explicit forcing terms
         # 
@@ -51 +53 @@
         # FIXME: How to call and how function should look
 
-        self.explicit_forcing_terms = []
+        self.forcing_terms = []
 
 
@@ -64 +66 @@
         #   q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
         
-        self.semi_implicit_forcing_terms = []
+        ###self.semi_implicit_forcing_terms = []
 
 
@@ -216 +218 @@
                   %(self.time, self.min_timestep, self.number_of_steps,
                     self.number_of_first_order_steps)
-        else:    
+        elif self.min_timestep > self.max_timestep:
+            print 'Time = %.4f, steps=%d (%d)'\
+                  %(self.time, self.number_of_steps,
+                    self.number_of_first_order_steps)            
+        else:
             print 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
                   %(self.time, self.min_timestep,
@@ -416 +422 @@
         """
 
-        for f in self.explicit_forcing_terms:
-            
-            pass
-
+        for f in self.forcing_terms:
+            f(self)
 
 

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -33 +33 @@
         self.minimum_allowed_height = minimum_allowed_height
 
-        
+        self.forcing_terms.append(gravity)
+        self.forcing_terms.append(manning_friction)
 
 
@@ -280 +281 @@
         raise 'Unknown order'
 
+    #Compute edge values
+    for name in domain.conserved_quantities:
+        Q = domain.quantities[name]
+        Q.interpolate_from_vertices_to_edges()            
+
+
 def protect_against_infinitesimal_heights_centroid(domain):
     """Adjust height and momentum at centroid if height is less than
@@ -333 +340 @@
         Q = domain.quantities[name]
         Q.extrapolate_first_order()
-        Q.interpolate_from_vertices_to_edges()        
+
      
   
@@ -419 +426 @@
     """
     
-    #FIXME: first and second order migh merge
-
-    from Numeric import minimum, maximum
-    
+    #FIXME: first and second order might merge
+
     #Update conserved quantities using straight second order
     for name in domain.conserved_quantities:
         Q = domain.quantities[name]
-
         Q.extrapolate_second_order()
-        Q.limiter()
-        Q.interpolate_from_vertices_to_edges()                
+        Q.limit()
 
   
@@ -469 +472 @@
 
 
-        hmin = minimum( hv[k, :] )
+        hmin = min( hv[k, :] )
 
         #Create alpha in [0,1], where alpha==0 means using shallow 
@@ -500 +503 @@
         # from extrapolator xmomv and ymomv (deep).
         
-
-        ##f##or i in range(3):
-        #####    xmomv[k,i] = (1-alpha)*xmomc[k] + alpha*xmomv[k,i];
-            
         xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:];            
             
@@ -512 +511 @@
     hc = wc-zc    
     for k in range(domain.number_of_elements):
-        hmax = maximum(hv[k,:])
-        
-        if hmax < minimum_allowed_height:        
+        hmax = max(hv[k,:])
+        
+        if hmax < domain.minimum_allowed_height:        
             #Reset negative heights to bed elevation        
             if hc[k] < 0.0:
                 wc[k] = zc[k]
-                ###hc[k] = 0.0
             for i in range(3):    
                 if hv[k,i] < 0.0:
                     wv[k,i] = zv[k,i]
-                    ##hv0 = 0.0;}         
 
     
@@ -569 +566 @@
 #########################
 #Standard forcing terms:
-
+#
 def gravity(domain):
     """Implement forcing function for bed slope working with
@@ -629 +626 @@
             
 
+###########################
+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
+        self.h = h
+
+    def __call__(self, x, y):
+        if self.W is None:
+            from Numeric import ones, Float
+            return self.h*ones(len(x), Float)
+        else:
+            return self.W(x,y) + self.h
+
+
+
 ##############################################
 #Initialise module

inundation/ga/storm_surge/pyvolution/test_quantity.py

@@ -217 +217 @@
 
 
-#     def test_second_order_extrapolation2(self):
-
-#         initialise_consecutive_datastructure(points=6+4, elements=4)                
-        
-#         a = Point (0.0, 0.0)
-#         b = Point (0.0, 2.0)
-#         c = Point (2.0, 0.0)
-#         d = Point (0.0, 4.0)
-#         e = Point (2.0, 2.0)
-#         f = Point (4.0, 0.0)
-
-#         #Set up for a gradient of (3,1), f(x) = 3x+y
-#         v1 = Volume(b,a,c,array([2.0+2.0/3,0,0]))        
-#         v2 = Volume(b,c,e,array([4.0+4.0/3,0,0]))
-#         v3 = Volume(e,c,f,array([8.0+2.0/3,0,0]))
-#         v4 = Volume(d,b,e,array([2.0+8.0/3,0,0]))
-
-#         #Setup neighbour structure
-#         domain = Domain([v1,v2,v3,v4])
-#         domain.precompute()
-#       domain.check_integrity()
-        
-#         #Lets's check first order first, hey
-#       domain.order = 1
-#       domain.limiter = dummy_limiter
-#         distribute_to_vertices_and_edges(domain)              
-        
-#         assert allclose(v2.conserved_quantities_vertex0,
-#                         v2.conserved_quantities_centroid)
-#         assert allclose(v2.conserved_quantities_vertex1,
-#                         v2.conserved_quantities_centroid)
-#         assert allclose(v2.conserved_quantities_vertex2,
-#                         v2.conserved_quantities_centroid)                
-
-
-#         #Flux across right edge of volume 1
-#         #Outward pointing normal vector
-#         from shallow_water import flux_using_stage as flux_function
-#         normals = Volume.normals
-        
-#         normal = Vector.coordinates[normals[1][2]]
-#         ql = Volume.conserved_quantities_face2[1]
-#         qr = Volume.conserved_quantities_face1[0]
-#         fl = array([0.,0.])
-#         fr = array([0.,0.])
-#         flux0, max_speed = flux_function(normal, ql, qr, fl, fr)
-
-#         #print flux0, max_speed        
-
-#         #print
-#         #print v1.conserved_quantities_face0,\
-#         #      v2.conserved_quantities_face0,\
-#         #      v3.conserved_quantities_face0,\
-#         #      v4.conserved_quantities_face0              
-#         #print
-#         #edgelengths = Volume.geometric[:,2:]
-#         #print
-#         #print
-        
-#         from python_versions import compute_flux
-#         compute_flux(domain, 100)
-#         F1 = Volume.explicit_update
-
-#         from domain import compute_flux
-#         compute_flux(domain, 100)
-#         F2 = Volume.explicit_update        
-
-#         assert allclose(F1, F2)
-        
-#         #print F1
-#         #print F2
-        
-
-
-#         #Gradient of fitted pwl surface    
-#       a, b = compute_gradient(v2.id)  
-
-#         assert abs(a[0] - 3.0) < epsilon
-#         assert abs(b[0] - 1.0) < epsilon
-#         #assert qmin[0] == 2.0 + 2.0/3
-#         #assert qmax[0] == 8.0 + 2.0/3                
-
-#         #And now for the second order stuff        
-#         # - the full second order extrapolation
-#       domain.order = 2
-#       domain.limiter = dummy_limiter
-#         distribute_to_vertices_and_edges(domain)              
-
-#         assert allclose(v2.conserved_quantities_vertex0[0], 2.0)
-#         assert allclose(v2.conserved_quantities_vertex1[0], 6.0)        
-#         assert allclose(v2.conserved_quantities_vertex2[0], 8.0)
-
+    def test_second_order_extrapolation2(self):
+        quantity = Quantity(self.mesh4)
+
+        #Set up for a gradient of (3,1), f(x) = 3x+y
+        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
+                            location = 'centroids')
+        
+        #Gradients
+        a, b = quantity.compute_gradients()
+
+        #gradient bewteen t0 and t1 is undefined as det==0
+        assert a[0] == 0.0
+        assert b[0] == 0.0
+        #The others are OK
+        for i in range(1,4):
+            assert allclose(a[i], 3.0)
+            assert allclose(b[i], 1.0)
+
+
+        quantity.extrapolate_second_order()
+
+        assert allclose(quantity.vertex_values[1,0], 2.0)
+        assert allclose(quantity.vertex_values[1,1], 6.0)        
+        assert allclose(quantity.vertex_values[1,2], 8.0)
 
 

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -614 +614 @@
 
         #Check that all levels are above elevation (within eps)
-
-        
-        assert alltrue(alltrue(greater_equal(L,E-epsilon)))                
-        
-
+        assert alltrue(alltrue(greater_equal(L,E-epsilon)))
+
+
+        
+    def test_second_order_flat_bed_onestep(self):
+
+        from mesh_factory import rectangular
+        from shallow_water import Domain, Reflective_boundary,\
+             Dirichlet_boundary, Constant_height
+        from Numeric import array
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(6, 6)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.smooth = False
+        domain.default_order=2
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain)
+        Bd = Dirichlet_boundary([0.1, 0., 0.])
+        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
+
+        domain.check_integrity()
+
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
+            pass# domain.write_time()
+
+        #Data from earlier version of pyvolution
+        assert allclose(domain.min_timestep, 0.0396825396825)
+        assert allclose(domain.max_timestep, 0.0396825396825) 
+
+        assert allclose(domain.quantities['level'].centroid_values[:12],
+                        [0.00171396, 0.02656103, 0.00241523, 0.02656103,
+                        0.00241523, 0.02656103, 0.00241523, 0.02656103,
+                        0.00241523, 0.02656103, 0.00241523, 0.0272623])
+
+        domain.distribute_to_vertices_and_edges()
+        assert allclose(domain.quantities['level'].vertex_values[:12,0],
+                        [0.0001714, 0.02656103, 0.00024152,
+                        0.02656103, 0.00024152, 0.02656103,
+                        0.00024152, 0.02656103, 0.00024152,
+                        0.02656103, 0.00024152, 0.0272623])
+
+        assert allclose(domain.quantities['level'].vertex_values[:12,1],
+                        [0.00315012, 0.02656103, 0.00024152, 0.02656103,
+                         0.00024152, 0.02656103, 0.00024152, 0.02656103,
+                         0.00024152, 0.02656103, 0.00040506, 0.0272623])
+        
+        assert allclose(domain.quantities['level'].vertex_values[:12,2],
+                        [0.00182037, 0.02656103, 0.00676264,
+                         0.02656103, 0.00676264, 0.02656103,
+                         0.00676264, 0.02656103, 0.00676264,
+                         0.02656103, 0.0065991, 0.0272623])
+
+        assert allclose(domain.quantities['xmomentum'].centroid_values[:12],
+                        [0.00113961, 0.01302432, 0.00148672,
+                         0.01302432, 0.00148672, 0.01302432,
+                         0.00148672, 0.01302432, 0.00148672 ,
+                         0.01302432, 0.00148672, 0.01337143])
+
+        assert allclose(domain.quantities['ymomentum'].centroid_values[:12],
+                        [-2.91240050e-004 , 1.22721531e-004,
+                         -1.22721531e-004, 1.22721531e-004 ,
+                         -1.22721531e-004, 1.22721531e-004,
+                         -1.22721531e-004 , 1.22721531e-004,
+                         -1.22721531e-004, 1.22721531e-004,
+                         -1.22721531e-004, -4.57969873e-005])
+
+
+        
+    def test_second_order_flat_bed_moresteps(self):
+
+        from mesh_factory import rectangular
+        from shallow_water import Domain, Reflective_boundary,\
+             Dirichlet_boundary, Constant_height
+        from Numeric import array
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(6, 6)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.smooth = False
+        domain.default_order=2
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain)
+        Bd = Dirichlet_boundary([0.1, 0., 0.])
+        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
+
+        domain.check_integrity()
+
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
+            pass
+
+        #Data from earlier version of pyvolution
+        #assert allclose(domain.min_timestep, 0.0396825396825)
+        #assert allclose(domain.max_timestep, 0.0396825396825) 
+        #print domain.quantities['level'].centroid_values
+
+        
+        
+    def test_bedslope_problem_first_order(self):
+
+        from mesh_factory import rectangular
+        from shallow_water import Domain, Reflective_boundary, Constant_height
+        from Numeric import array
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(6, 6)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.smooth = False
+        domain.default_order=1
+
+        #Bed-slope and friction
+        def x_slope(x, y):
+            return -x/3
+
+        domain.set_quantity('elevation', x_slope)
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+        #Initial condition
+        domain.set_quantity('level', Constant_height(x_slope, 0.05))
+        domain.check_integrity()
+
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
+            pass# domain.write_time()
+
+        #Data from earlier version of pyvolution
+        assert allclose(domain.min_timestep, 0.0423354136479)
+        assert allclose(domain.max_timestep, 0.0479397182296)
+        assert allclose(domain.quantities['level'].centroid_values,
+                        [0.00974056, 0.02131643, 0.01373269, 0.02139684,
+                         0.01373269, 0.02139684, 0.01373269, 0.02139684,
+                         0.01373269, 0.02139684, 0.01386079, 0.0254346,
+                         -0.04562579, -0.025253, -0.04154658, -0.02512009,
+                         -0.04154658, -0.02512009, -0.04154658, -0.02512009,
+                         -0.04154658, -0.02512009, -0.04140825, -0.0210463,
+                         -0.10118135, -0.08080855, -0.09710213, -0.08067564,
+                         -0.09710213, -0.08067564, -0.09710213, -0.08067564,
+                         -0.09710213, -0.08067564, -0.0969638 , -0.07660185,
+                         -0.1567369, -0.13636411, -0.15265769, -0.1362312,
+                         -0.15265769, -0.1362312, -0.15265769, -0.1362312,
+                         -0.15265769, -0.1362312, -0.15251936, -0.13215741,
+                         -0.21229246, -0.19191967, -0.20821325, -0.19178675,
+                         -0.20821325, -0.19178675, -0.20821325, -0.19178675,
+                         -0.20821325, -0.19178675, -0.20807491, -0.18771296,
+                         -0.25875433, -0.24722502, -0.25470346, -0.24709274,
+                         -0.25470346, -0.24709274, -0.25470346, -0.24709274,
+                         -0.25470346, -0.24709274, -0.25460675, -0.24309962])
+
+
+        
+    def test_bedslope_problem_first_order_moresteps(self):
+
+        from mesh_factory import rectangular
+        from shallow_water import Domain, Reflective_boundary, Constant_height
+        from Numeric import array
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(6, 6)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.smooth = False
+        domain.default_order=1
+
+        #Bed-slope and friction
+        def x_slope(x, y):
+            return -x/3
+
+        domain.set_quantity('elevation', x_slope)
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+        #Initial condition
+        domain.set_quantity('level', Constant_height(x_slope, 0.05))
+        domain.check_integrity()
+
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
+            pass# domain.write_time()
+
+        #Data from earlier version of pyvolution
+        assert allclose(domain.quantities['level'].centroid_values,
+                        [-0.03558044, -0.01758186, -0.03457868,
+                         -0.01757463, -0.03400681, -0.01765817,
+                         -0.0334924, -0.01763917, -0.03304351,
+                         -0.01745604, -0.03259134, -0.01690249,
+                         -0.08343238, -0.06844054, -0.08114626,
+                         -0.06776272, -0.07985917, -0.0670156,
+                         -0.07857347, -0.06636824, -0.07751765,
+                         -0.06571005, -0.07694892, -0.06485356,
+                         -0.12816762, -0.11399292, -0.12388354,
+                         -0.11298543, -0.12146159, -0.11160487,
+                         -0.11890573, -0.11016723, -0.11720749,
+                         -0.10922083, -0.11748089, -0.10796284,
+                         -0.17324244, -0.15860786, -0.16784961,
+                         -0.15785693, -0.16501143, -0.1558837,
+                         -0.1618709,  -0.15384493, -0.15979744,
+                         -0.15251944, -0.15976661, -0.15094149,
+                         -0.17055387, -0.1934305, -0.17869999,
+                         -0.19402574, -0.17570318, -0.19030006,
+                         -0.17114482, -0.18720482, -0.16970896,
+                         -0.18836511, -0.17655148, -0.19499397,
+                         -0.14139494, -0.14492939, -0.14365011,
+                         -0.14969524, -0.14343543, -0.14754936,
+                         -0.13944213, -0.14337135, -0.13651514,
+                         -0.14183568, -0.13623283, -0.14529803])
+        
+    def test_bedslope_problem_second_order(self):
+
+        from mesh_factory import rectangular
+        from shallow_water import Domain, Reflective_boundary, Constant_height
+        from Numeric import array
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(6, 6)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.smooth = False
+        domain.default_order=2
+
+        #Bed-slope and friction
+        def x_slope(x, y):
+            return -x/3
+
+        domain.set_quantity('elevation', x_slope)
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+        #Initial condition
+        domain.set_quantity('level', Constant_height(x_slope, 0.05))
+        domain.check_integrity()
+
+        assert allclose(domain.quantities['level'].centroid_values,
+                        [0.01296296, 0.03148148, 0.01296296,
+                        0.03148148, 0.01296296, 0.03148148,
+                        0.01296296, 0.03148148, 0.01296296,
+                        0.03148148, 0.01296296, 0.03148148,
+                        -0.04259259, -0.02407407, -0.04259259,
+                        -0.02407407, -0.04259259, -0.02407407,
+                        -0.04259259, -0.02407407, -0.04259259,
+                        -0.02407407, -0.04259259, -0.02407407,
+                        -0.09814815, -0.07962963, -0.09814815,
+                        -0.07962963, -0.09814815, -0.07962963,
+                        -0.09814815, -0.07962963, -0.09814815,
+                        -0.07962963, -0.09814815, -0.07962963,
+                        -0.1537037 , -0.13518519, -0.1537037,
+                        -0.13518519, -0.1537037, -0.13518519,
+                        -0.1537037 , -0.13518519, -0.1537037,
+                        -0.13518519, -0.1537037, -0.13518519,
+                        -0.20925926, -0.19074074, -0.20925926,
+                        -0.19074074, -0.20925926, -0.19074074,
+                        -0.20925926, -0.19074074, -0.20925926,
+                        -0.19074074, -0.20925926, -0.19074074,
+                        -0.26481481, -0.2462963, -0.26481481,
+                        -0.2462963, -0.26481481, -0.2462963,
+                        -0.26481481, -0.2462963, -0.26481481,
+                        -0.2462963, -0.26481481, -0.2462963])
+
+
+        #print domain.quantities['level'].extrapolate_second_order()
+        #domain.distribute_to_vertices_and_edges()
+        #print domain.quantities['level'].vertex_values[:,0]        
+        
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
+            #domain.write_time()            
+            pass
+
+
+        #print domain.quantities['level'].centroid_values
+        
+        #Data from earlier version of pyvolution
+        assert allclose(domain.min_timestep, 0.0376895634803) 
+        assert allclose(domain.max_timestep, 0.0415635655309)
+
+
+        #FIXME: Fails
+        #assert allclose(domain.quantities['level'].centroid_values,
+        #                [0.00855788, 0.01575204, 0.00994606, 0.01717072,
+        #                 0.01005985, 0.01716362, 0.01005985, 0.01716299,
+        #                 0.01007098, 0.01736248, 0.01216452, 0.02026776,
+        #                 -0.04462374, -0.02479045, -0.04199789, -0.0229465,
+        #                 -0.04184033, -0.02295693, -0.04184013, -0.02295675,
+        #                 -0.04184486, -0.0228168, -0.04028876, -0.02036486,
+        #                 -0.10029444, -0.08170809, -0.09772846, -0.08021704,
+        #                 -0.09760006, -0.08022143, -0.09759984, -0.08022124,
+        #                 -0.09760261, -0.08008893, -0.09603914, -0.07758209,
+        #                 -0.15584152, -0.13723138, -0.15327266, -0.13572906,
+        #                 -0.15314427, -0.13573349, -0.15314405, -0.13573331,
+        #                 -0.15314679, -0.13560104, -0.15158523, -0.13310701,
+        #                 -0.21208605, -0.19283913, -0.20955631, -0.19134189,
+        #                 -0.20942821, -0.19134598, -0.20942799, -0.1913458,
+        #                 -0.20943005, -0.19120952, -0.20781177, -0.18869401,
+        #                 -0.25384082, -0.2463294, -0.25047649, -0.24464654,
+        #                 -0.25031159, -0.24464253, -0.25031112, -0.24464253,
+        #                 -0.25031463, -0.24454764, -0.24885323, -0.24286438])
+        
     
 #-------------------------------------------------------------