Changeset 1015

Timestamp:

Mar 5, 2005, 6:55:36 PM (20 years ago)

Author:

ole

Message:

Further testing of limiters and conservation

Location:

inundation/ga/storm_surge/pyvolution

Files:

• shallow_water.py (modified) (4 diffs)
• test_shallow_water.py (modified) (6 diffs)

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -526 +526 @@
         else:
             raise 'Unknown order'
-
+        
     #Take bed elevation into account when water heights are small
     balance_deep_and_shallow(domain)
@@ -567 +567 @@
     """
     
+    #FIXME: Look here for error
+    
     #Shortcuts
     wc = domain.quantities['stage'].centroid_values
@@ -574 +576 @@
     hc = wc - zc  #Water depths at centroids    
 
+    #print zc
+    #print '1', wc      
     #Update
     for k in range(domain.number_of_elements):
@@ -584 +588 @@
             xmomc[k] = ymomc[k] = 0.0
         
+    #print '2', wc              
 
 

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -1624 +1624 @@
             #print xmom
             #assert allclose (xmom, initial_xmom) 
-        
-        
+
+            
     def test_conservation_3(self):
         """Test that stage is conserved globally
         
-        This one uses a convoluted bed, reflective bdries and a suitable 
+        This one uses a larger grid, convoluted bed, reflective bdries and a suitable 
         initial condition
         """
@@ -1638 +1638 @@
 
         #Create basic mesh
-        points, vertices, boundary = rectangular(6, 6)
+        points, vertices, boundary = rectangular(2, 1)
 
         #Create shallow water domain
         domain = Domain(points, vertices, boundary)
         domain.smooth = False
-        domain.default_order=1  #FIXME: 2 will fail
+        domain.default_order=2
         domain.beta_h = 0
 
@@ -1655 +1655 @@
                     z[i] = -0.5
                 if 0.5 <= x[i] < 0.7:  
-                    z[i] = 0.4       
-                    #z[i] = 0.7                     
+                    z[i] = 0.39                     
                 if 0.7 <= x[i]:  
                     z[i] = x[i]/3                   
@@ -1678 +1677 @@
         initial_xmom = domain.quantities['xmomentum'].get_integral()      
         
+        import copy
+        ref_centroid_values =\
+                 copy.copy(domain.quantities['stage'].centroid_values)
+                                
+        #print 'ORG', domain.quantities['stage'].centroid_values                                
+        domain.distribute_to_vertices_and_edges()
+        
+        
+        #print domain.quantities['stage'].centroid_values
+        assert allclose(domain.quantities['stage'].centroid_values,
+                        ref_centroid_values)            
+        
+        
+        #Check that initial limiter doesn't violate cons quan
+        assert allclose (domain.quantities['stage'].get_integral(),
+                         initial_volume) 
+        
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 5):
+            volume =  domain.quantities['stage'].get_integral()         
+            #print t, volume, initial_volume
+            assert allclose (volume, initial_volume)        
+                    
+        
+    def test_conservation_4(self):
+        """Test that stage is conserved globally
+        
+        This one uses a larger grid, convoluted bed, reflective bdries and a suitable 
+        initial condition
+        """
+        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
+        domain.beta_h = 0.0
+
+        #IC
+        def x_slope(x, y):
+            z = 0*x
+            for i in range(len(x)):
+                if x[i] < 0.3:  
+                    z[i] = x[i]/3
+                if 0.3 <= x[i] < 0.5:  
+                    z[i] = -0.5
+                if 0.5 <= x[i] < 0.7:  
+                    z[i] = 0.34 #OK
+                    #z[i] = 0.35        #Fails      
+                if 0.7 <= x[i]:  
+                    z[i] = x[i]/3                   
+            return z                                
+                                    
+                    
+
+        domain.set_quantity('elevation', x_slope)
+        domain.set_quantity('friction', 0)      
+        domain.set_quantity('stage', 0.4) #Steady       
+        
+        # Boundary conditions (reflective everywhere)
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+        domain.check_integrity()
+
+        #domain.visualise = True #If you want to take a sticky beak  
+        
+        initial_volume = domain.quantities['stage'].get_integral()         
+        initial_xmom = domain.quantities['xmomentum'].get_integral()      
+        
         
         #Let surface settle according to limiter 
@@ -1684 +1759 @@
         #needs to respect the law of the limiters?
         domain.distribute_to_vertices_and_edges()
+        
+        
+        #FIXME: At this point study the centroid of individual triangles.
+        
+        #Check that initial limiter doesn't violate cons quan
+        assert allclose (domain.quantities['stage'].get_integral(),
+                         initial_volume) 
+        
+                         
+                                
         initial_volume = domain.quantities['stage'].get_integral()         
                 
@@ -2950 +3035 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(TestCase,'test')
-    #suite = unittest.makeSuite(TestCase,'test_balance_deep_and_shallow')
     runner = unittest.TextTestRunner()
     runner.run(suite)