Changeset 1016

Timestamp:

Mar 6, 2005, 4:18:34 PM (20 years ago)

Author:

ole

Message:

More testing of conservation.
Checked that momentum is conserved in a steady state scenario

Location:

inundation/ga/storm_surge/pyvolution

Files:

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

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -830 +830 @@
     wv = domain.quantities['stage'].vertex_values
     zv = domain.quantities['elevation'].vertex_values
-    hv = wv-zv
+    hv = wv - zv
 
     #Momentums at centroids
@@ -842 +842 @@
     #Limit h
     hvbar = h_limiter(domain)    
+    
+    #This is how one would make a first order h_limited value
+    #as in the old balancer (pre 17 Feb 2005):
+    #from Numeric import zeros, Float
+    #hvbar = zeros( (len(hc), 3), Float)
+    #for i in range(3):
+    #    hvbar[:,i] = hc[:]
 
     from shallow_water_ext import balance_deep_and_shallow
-    #balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
-    #                         xmomc, ymomc, xmomv, ymomv)
-
     balance_deep_and_shallow(wc, zc, hc, wv, zv, hv, hvbar,
                              xmomc, ymomc, xmomv, ymomv)    

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -1644 +1644 @@
         domain.smooth = False
         domain.default_order=2
-        domain.beta_h = 0
+        domain.beta_h = 0.2
 
         #IC
@@ -1695 +1695 @@
         
         #Evolution
-        for t in domain.evolve(yieldstep = 0.05, finaltime = 5):
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 10):
             volume =  domain.quantities['stage'].get_integral()         
             #print t, volume, initial_volume
@@ -1729 +1729 @@
                 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.5 <= x[i] < 0.7:
+                    #z[i] = 0.3 #OK with beta == 0.2              
+                    z[i] = 0.34 #OK with beta == 0.0
+                    #z[i] = 0.35#Fails after 80 timesteps with an error 
+                                #of the order 1.0e-5    
                 if 0.7 <= x[i]:  
                     z[i] = x[i]/3                   
@@ -1753 +1755 @@
         initial_xmom = domain.quantities['xmomentum'].get_integral()      
         
-        
-        #Let surface settle according to limiter 
-        #FIXME: What exactly is the significance of this?
-        #Is it a problem that an arbitrary initial condition 
-        #needs to respect the law of the limiters?
+        import copy
+        ref_centroid_values =\
+                 copy.copy(domain.quantities['stage'].centroid_values)
+                                
+        #Test limiter by itself
         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) 
-        
+        #NOTE: This would fail if any initial stage was less than the 
+        #corresponding bed elevation - but that is reasonable.
                          
-                                
-        initial_volume = domain.quantities['stage'].get_integral()         
                 
-        #print initial_volume
-        
-        #print initial_xmom
-
         #Evolution
         for t in domain.evolve(yieldstep = 0.05, finaltime = 10.0):
             volume =  domain.quantities['stage'].get_integral()         
+
             #print t, volume, initial_volume
+    
             
+            #if not allclose (volume, initial_volume): 
+            #    print 't==4.05'
+            #    for k in range(domain.number_of_elements):
+            #       pass
+            #       print domain.quantities['stage'].centroid_values[k] -\
+            #             ref_centroid_values[k] 
+            
             assert allclose (volume, initial_volume) 
-         
-            #FIXME: What would we expect from momentum
-            #xmom = domain.quantities['xmomentum'].get_integral()   
-            #print xmom
-            #assert allclose (xmom, initial_xmom) 
-        
-        
+        
+
+        
+        
+        
+    def test_conservation_5(self):
+        """Test that momentum is conserved globally in 
+        steady state scenario
+        
+        This one uses a slopy bed, dirichlet and reflective bdries
+        """
+        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
+
+        #IC
+        def x_slope(x, y):
+            return x/3
+
+        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)        
+        Bleft = Dirichlet_boundary([0.5,0,0])
+        Bright = Dirichlet_boundary([0.1,0,0])  
+        domain.set_boundary({'left': Bleft, 'right': Bright, 
+                             '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()      
+        
+
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 15.0):
+            stage =  domain.quantities['stage'].get_integral()         
+            xmom = domain.quantities['xmomentum'].get_integral()            
+            ymom = domain.quantities['ymomentum'].get_integral()                    
+            
+            if allclose(t, 6):  #Steady state reached
+                steady_xmom = domain.quantities['xmomentum'].get_integral()
+                steady_ymom = domain.quantities['ymomentum'].get_integral()
+                steady_stage = domain.quantities['stage'].get_integral()                
+                
+            if t > 6:
+                #print '%.2f %14.8f %14.8f' %(t, ymom, steady_ymom)
+                assert allclose(xmom, steady_xmom)              
+                assert allclose(ymom, steady_ymom)
+                assert allclose(stage, steady_stage)            
+                                                
+            
+                
                 
     def test_second_order_flat_bed_onestep(self):