Changeset 1011

Timestamp:

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

Author:

ole

Message:

Added test for conservation of stage with various beds

Location:

inundation/ga/storm_surge/pyvolution

Files:

• general_mesh.py (modified) (1 diff)
• mesh.py (modified) (1 diff)
• quantity.py (modified) (1 diff)
• test_general_mesh.py (modified) (1 diff)
• test_quantity.py (modified) (3 diffs)
• test_shallow_water.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/general_mesh.py

@@ -267 +267 @@
 
 
-
-
+    def get_area(self):
+        """Return total area of mesh
+        """
+
+        return sum(self.areas)          
+        

inundation/ga/storm_surge/pyvolution/mesh.py

@@ -414 +414 @@
         return polygon
         
-            
-
-        
 
     def check_integrity(self):

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -568 +568 @@
             
 
+    def get_integral(self):
+        """Compute the integral of quantity across entire domain
+        """
+        integral = 0
+        for k in range(self.domain.number_of_elements):
+            area = self.domain.areas[k]
+            qc = self.centroid_values[k]     
+            integral += qc*area
+            
+        return integral    
+        
 
 class Conserved_quantity(Quantity):

inundation/ga/storm_surge/pyvolution/test_general_mesh.py

@@ -34 +34 @@
         assert domain.get_vertices([0,4]) == [domain.triangles[0],
                                               domain.triangles[4]]
+    def test_areas(self):
+        from mesh_factory import rectangular
+        from shallow_water import Domain
+        from Numeric import zeros, Float
+        
+        #Create basic mesh
+        points, vertices, boundary = rectangular(1, 3)
+        domain = Domain(points, vertices, boundary)
+
+        assert domain.get_area() == 1.0         
+
 
     def test_get_unique_vertex_values(self):

inundation/ga/storm_surge/pyvolution/test_quantity.py

@@ -16 +16 @@
         a = [0.0, 0.0]
         b = [0.0, 2.0]
-        c = [2.0,0.0]
+        c = [2.0, 0.0]
         d = [0.0, 4.0]
         e = [2.0, 2.0]
-        f = [4.0,0.0]
+        f = [4.0, 0.0]
 
         points = [a, b, c, d, e, f]
@@ -196 +196 @@
         assert allclose(quantity.centroid_values,
                         [4.0/3, 8.0/3, 10.0/3, 10.0/3])
+
+
+    def test_integral(self):
+        quantity = Quantity(self.mesh4)
+
+        #Try constants first    
+        const = 5
+        quantity.set_values(const, location = 'vertices')       
+        #print 'Q', quantity.get_integral()
+                
+        assert allclose(quantity.get_integral(), 
+                        self.mesh4.get_area() * const)
+        
+        #Try with a linear function
+        def f(x, y):
+            return x+y
+
+        quantity.set_values(f, location = 'vertices')
+        
+        
+        ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
+        
+        assert allclose (quantity.get_integral(), ref_integral) 
+
 
 
@@ -900 +924 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test_limiter2')
+    suite = unittest.makeSuite(TestCase,'test')
     #print "restricted test"
     #suite = unittest.makeSuite(TestCase,'test_set_vertex_values_subset')

inundation/ga/storm_surge/pyvolution/test_shallow_water.py

@@ -1522 +1522 @@
 
         
-        
+    def test_conservation_1(self):
+        """Test that stage is conserved globally
+        
+        This one uses a flat 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=1
+
+        #IC
+        def x_slope(x, y):
+            return x/3
+
+        domain.set_quantity('elevation', 0)
+        domain.set_quantity('friction', 0)      
+        domain.set_quantity('stage', x_slope)   
+        
+        # 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()      
+        
+        #print initial_xmom
+
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
+            volume =  domain.quantities['stage'].get_integral()         
+            assert allclose (volume, initial_volume) 
+         
+            #I don't believe that the total momentum should be the same
+            #It starts with zero and ends with zero though
+            #xmom = domain.quantities['xmomentum'].get_integral()   
+            #print xmom
+            #assert allclose (xmom, initial_xmom) 
+        
+
+        
+    def test_conservation_2(self):
+        """Test that stage is conserved globally
+        
+        This one uses a slopy 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=1
+
+        #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)
+        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()      
+        
+        #print initial_xmom
+
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
+            volume =  domain.quantities['stage'].get_integral()         
+            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_3(self):
+        """Test that stage is conserved globally
+        
+        This one uses a 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=1
+
+        #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.4
+                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 
+        #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?
+        domain.distribute_to_vertices_and_edges()
+        initial_volume = domain.quantities['stage'].get_integral()         
+                
+        #print initial_volume
+        
+        #print initial_xmom
+
+        #Evolution
+        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
+            volume =  domain.quantities['stage'].get_integral()         
+            #print t, volume, initial_volume
+            
+            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_second_order_flat_bed_onestep(self):