Changeset 7090

Timestamp:

May 26, 2009, 5:41:49 PM (14 years ago)

Author:

ole

Message:

Added tests for volumetric conservation in the presence of inflow or rainfall.

Location:

anuga_core/source/anuga/shallow_water

Files:

• shallow_water_domain.py (modified) (2 diffs)
• test_shallow_water_domain.py (modified) (13 diffs)

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -880 +880 @@
         """Create volumetric balance report suitable for printing or logging.
         """
-        
-        boundary_flows, total_boundary_inflow, total_boundary_outflow = self.compute_boundary_flows() 
+
+        X = self.compute_boundary_flows() 
+        boundary_flows, total_boundary_inflow, total_boundary_outflow = X
         
         s = '---------------------------\n'        
         s += 'Volumetric balance report:\n'
         s += '--------------------------\n'
-        s += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
-        s += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow        
+        s += 'Total boundary inflow [m^3/s]: %.2e\n' % total_boundary_inflow
+        s += 'Total boundary outflow [m^3/s]: %.2e\n' % total_boundary_outflow        
         s += 'Net boundary flow by tags [m^3/s]\n'
         for tag in boundary_flows:
-            s += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
-        
-        s += 'Total net boundary flow [m^3/s]: %.2f\n' % (total_boundary_inflow + total_boundary_outflow) 
-        s += 'Total volume in domain [m^3]: %.2f\n' % self.compute_total_volume()
+            s += '    %s [m^3/s]: %.2e\n' % (tag, boundary_flows[tag])
+        
+        s += 'Total net boundary flow [m^3/s]: %.2e\n' % (total_boundary_inflow + total_boundary_outflow) 
+        s += 'Total volume in domain [m^3]: %.2e\n' % self.compute_total_volume()
         
         # The go through explicit forcing update and record the rate of change for stage and 
@@ -2278 +2279 @@
                   the inflow region and then applied to update the
                   stage quantity.
+                  
+                  Note also that any reference to the internal attribute 'rate'
+                  later will use the one converted to m/s.
 
     polygon: Specifies a polygon to restrict the rainfall.

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -6745 +6745 @@
         
 
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
         # Import necessary modules
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
         from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
         from anuga.shallow_water import Domain
 
 
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
         # Setup computational domain
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
 
         length = 100.
@@ -6765 +6765 @@
                 
 
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
         # Simple flat bottom bathtub
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
 
         d = 1.0
@@ -6776 +6776 @@
 
 
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
         # Slope
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
                 
         slope = 1.0/10  # RHS drops to -10m
@@ -6812 +6812 @@
         
 
-        #------------------------------------------------------------------------------
+        #---------------------------------------------------------------------
         # Import necessary modules
-        #------------------------------------------------------------------------------
+        #---------------------------------------------------------------------
         from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
         from anuga.shallow_water import Domain
@@ -6822 +6822 @@
         from anuga.shallow_water.data_manager import get_flow_through_cross_section
 
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
         # Setup computational domain
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
         finaltime = 300.0
 
@@ -6846 +6846 @@
                 
 
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
         # Setup initial conditions
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
         slope = 0.0
         def topography(x, y):
@@ -6860 +6860 @@
                             expression='elevation')
 
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
         # Setup boundary conditions
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
 
         Br = Reflective_boundary(domain)      # Solid reflective wall
@@ -6875 +6875 @@
 
         
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
         # Evolve system through time
-        #------------------------------------------------------------------------------
+        #----------------------------------------------------------------------
         for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
         
@@ -6887 +6887 @@
             
         
-        
+            
+    def test_volume_conservation_inflow(self):
+        """test_volume_conservation
+        
+        Test that total volume in domain is as expected, based on questions
+        raised by Petar Milevski in May 2009.
+        
+        This test adds inflow at a known rate and verifies that the total 
+        terminal volume is as expected.
+        
+        """
+
+        verbose = False
+        
+
+        #---------------------------------------------------------------------
+        # Import necessary modules
+        #---------------------------------------------------------------------
+        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+        from anuga.shallow_water import Domain
+        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
+        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
+        from anuga.shallow_water.shallow_water_domain import Inflow
+        from anuga.shallow_water.data_manager import get_flow_through_cross_section
+
+        #----------------------------------------------------------------------
+        # Setup computational domain
+        #----------------------------------------------------------------------
+        finaltime = 200.0
+
+        length = 300.
+        width  = 20.
+        dx = dy = 5       # Resolution: of grid on both axes
+        
+
+        points, vertices, boundary = rectangular_cross(int(length/dx), 
+                                                       int(width/dy),
+                                                       len1=length, len2=width)
+
+
+        domain = Domain(points, vertices, boundary)   
+        domain.set_name('Inflow_volume_test')              # Output name
+                
+
+        #----------------------------------------------------------------------
+        # Setup initial conditions
+        #----------------------------------------------------------------------
+        slope = 0.0
+        def topography(x, y):
+            z=-x * slope
+            return z
+
+        domain.set_quantity('elevation', topography) # Use function for elevation
+        domain.set_quantity('friction', 0.0)         # Constant friction
+                
+        domain.set_quantity('stage',
+                            expression='elevation')  # Dry initially
+                            
+
+        #--------------------------------------------------------------
+        # Setup Inflow
+        #--------------------------------------------------------------
+
+        # Fixed Flowrate onto Area 
+        fixed_inflow = Inflow(domain,
+                              center=(10.0, 10.0),
+                              radius=5.00,
+                              rate=10.00)                               
+                            
+        domain.forcing_terms.append(fixed_inflow)                            
+        
+        #----------------------------------------------------------------------
+        # Setup boundary conditions
+        #----------------------------------------------------------------------
+
+        Br = Reflective_boundary(domain) # Solid reflective wall
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+        
+        #----------------------------------------------------------------------
+        # Evolve system through time
+        #----------------------------------------------------------------------
+        ref_volume = 0.0
+        ys = 10.0  # Yieldstep
+        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
+        
+            # Check volume
+            assert num.allclose(domain.compute_total_volume(), ref_volume)
+        
+            if verbose :
+                print domain.timestepping_statistics()
+                print domain.volumetric_balance_statistics()
+                print 'reference volume', ref_volume
+            
+            
+            # Update reference volume
+            ref_volume += ys * fixed_inflow.rate
+        
+            
+            
+
+            
+    def test_volume_conservation_rain(self):
+        """test_volume_conservation
+        
+        Test that total volume in domain is as expected, based on questions
+        raised by Petar Milevski in May 2009.
+        
+        This test adds rain at a known rate and verifies that the total 
+        terminal volume is as expected.
+        
+        """
+
+        verbose = False
+        
+
+        #---------------------------------------------------------------------
+        # Import necessary modules
+        #---------------------------------------------------------------------
+        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+        from anuga.shallow_water import Domain
+        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
+        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
+        from anuga.shallow_water.shallow_water_domain import Rainfall
+        from anuga.shallow_water.data_manager import get_flow_through_cross_section
+
+        #----------------------------------------------------------------------
+        # Setup computational domain
+        #----------------------------------------------------------------------
+        finaltime = 200.0
+
+        length = 300.
+        width  = 20.
+        dx = dy = 5       # Resolution: of grid on both axes
+        
+
+        points, vertices, boundary = rectangular_cross(int(length/dx), 
+                                                       int(width/dy),
+                                                       len1=length, len2=width)
+
+
+        domain = Domain(points, vertices, boundary)   
+        domain.set_name('Rain_volume_test')              # Output name
+                
+
+        #----------------------------------------------------------------------
+        # Setup initial conditions
+        #----------------------------------------------------------------------
+        slope = 0.0
+        def topography(x, y):
+            z=-x * slope
+            return z
+
+        domain.set_quantity('elevation', topography) # Use function for elevation
+        domain.set_quantity('friction', 0.0)         # Constant friction
+                
+        domain.set_quantity('stage',
+                            expression='elevation')  # Dry initially
+                            
+
+        #--------------------------------------------------------------
+        # Setup rain
+        #--------------------------------------------------------------
+
+        # Fixed rain onto small circular area 
+        fixed_rain = Rainfall(domain,
+                              center=(10.0, 10.0),
+                              radius=5.00,
+                              rate=10.00)   # 10 mm/s                            
+                            
+        domain.forcing_terms.append(fixed_rain)                            
+        
+        #----------------------------------------------------------------------
+        # Setup boundary conditions
+        #----------------------------------------------------------------------
+
+        Br = Reflective_boundary(domain) # Solid reflective wall
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+        
+        #----------------------------------------------------------------------
+        # Evolve system through time
+        #----------------------------------------------------------------------
+        ref_volume = 0.0
+        ys = 10.0  # Yieldstep
+        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
+        
+            # Check volume
+            V = domain.compute_total_volume()
+            msg = 'V = %e, Ref = %e' % (V, ref_volume)
+            assert num.allclose(V, ref_volume), msg
+        
+            if verbose :
+                print domain.timestepping_statistics()
+                print domain.volumetric_balance_statistics()
+                print 'reference volume', ref_volume
+                print V
+            
+            
+            # Update reference volume.
+            # FIXME: Note that rate has now been redefined
+            # as m/s internally. This is a little confusing
+            # when it was specfied as mm/s.
+            
+            delta_V = fixed_rain.rate*fixed_rain.exchange_area
+            ref_volume += ys * delta_V
+        
+            
+            
+            
+    def Xtest_rain_conservation_and_runoff(self):
+        """test_rain_conservation_and_runoff
+        
+        Test that total volume in domain is as expected, based on questions
+        raised by Petar Milevski in May 2009.
+        
+        This test adds rain at a known rate and verifies that the total 
+        volume and outflows are as expected.
+        
+        """
+
+        # FIXME (Ole): Does not work yet. Investigate boundary flows
+        
+        verbose = True #False
+        
+
+        #---------------------------------------------------------------------
+        # Import necessary modules
+        #---------------------------------------------------------------------
+        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+        from anuga.shallow_water import Domain
+        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
+        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
+        from anuga.shallow_water.shallow_water_domain import Rainfall
+        from anuga.shallow_water.data_manager import get_flow_through_cross_section
+
+        #----------------------------------------------------------------------
+        # Setup computational domain
+        #----------------------------------------------------------------------
+        finaltime = 500.0
+
+        length = 300.
+        width  = 20.
+        dx = dy = 5       # Resolution: of grid on both axes
+        
+
+        points, vertices, boundary = rectangular_cross(int(length/dx), 
+                                                       int(width/dy),
+                                                       len1=length, len2=width)
+
+
+        domain = Domain(points, vertices, boundary)   
+        domain.set_name('Rain_volume_runoff_test')         # Output name
+                
+
+        #----------------------------------------------------------------------
+        # Setup initial conditions
+        #----------------------------------------------------------------------
+        slope = 0.0
+        def topography(x, y):
+            z=-x * slope
+            return z
+
+        domain.set_quantity('elevation', topography) # Use function for elevation
+        domain.set_quantity('friction', 0.0)         # Constant friction
+                
+        domain.set_quantity('stage',
+                            expression='elevation')  # Dry initially
+                            
+
+        #--------------------------------------------------------------
+        # Setup rain
+        #--------------------------------------------------------------
+
+        # Fixed rain onto small circular area 
+        fixed_rain = Rainfall(domain,
+                              center=(10.0, 10.0),
+                              radius=5.00,
+                              rate=10.00)   # 10 mm/s                            
+                            
+        domain.forcing_terms.append(fixed_rain)                            
+        
+        #----------------------------------------------------------------------
+        # Setup boundary conditions
+        #----------------------------------------------------------------------
+
+        Br = Reflective_boundary(domain) # Solid reflective wall
+        Bt = Transmissive_stage_zero_momentum_boundary(domain)
+        Bd = Dirichlet_boundary([-10, 0, 0])
+        domain.set_boundary({'left': Bt, 'right': Bd, 'top': Bt, 'bottom': Bt})
+
+        
+        #----------------------------------------------------------------------
+        # Evolve system through time
+        #----------------------------------------------------------------------
+        ref_volume = 0.0
+        ys = 10.0  # Yieldstep
+        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
+        
+            # Check volume
+            V = domain.compute_total_volume()
+            msg = 'V = %e, Ref = %e' % (V, ref_volume)
+            #assert num.allclose(V, ref_volume) or V < ref_volume, msg
+        
+            if verbose:
+                print domain.timestepping_statistics()
+                print domain.volumetric_balance_statistics()
+                print 'reference volume', ref_volume
+                print V
+            
+            
+            # Update reference volume.
+            # FIXME: Note that rate has now been redefined
+            # as m/s internally. This is a little confusing
+            # when it was specfied as mm/s.
+            
+            delta_V = fixed_rain.rate*fixed_rain.exchange_area
+            ref_volume += ys * delta_V
+        
+            # Compute outflow at right hand downstream boundary
+            boundary_flows, inflow , outflow = domain.compute_boundary_flows()
+            net_outflow = outflow - inflow
+        
+            outflow = boundary_flows['right']
+            if verbose:
+                print 'Outflow', outflow
+                print 'Net outflow', net_outflow
+        
+            # Update reference volume
+            ref_volume += ys * outflow            
+
+        
+                
         
     def test_inflow_using_flowline(self):
@@ -6929 +7261 @@
         dx = dy = 5          # Resolution: of grid on both axes
 
-        points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
+        points, vertices, boundary = rectangular_cross(int(length/dx), 
+                                                       int(width/dy),
                                                        len1=length, len2=width)
 
@@ -7013 +7346 @@
                     #if verbose :
                     #    print domain.timestepping_statistics()
-                    #    print domain.volumetric_balance_statistics()                                    
-
+                    #    print domain.volumetric_balance_statistics()
 
                 #--------------------------------------------------------------
@@ -7021 +7353 @@
                     
                 # Square on flowline at 200m
-                q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,20.0]])
-                msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow)
+                q=domain.get_flow_through_cross_section([[200.0,0.0],
+                                                         [200.0,20.0]])
+                msg = 'Predicted flow was %f, should have been %f' %\
+                    (q, ref_flow)
                 if verbose:
-                    print '90 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow)
+                    print '90 degree flowline: ANUGA = %f, Ref = %f' %\
+                        (q, ref_flow)
                 assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
 
                            
                 # 45 degree flowline at 200m
-                q=domain.get_flow_through_cross_section([[200.0,0.0],[220.0,20.0]])
-                msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow)
+                q=domain.get_flow_through_cross_section([[200.0,0.0],
+                                                         [220.0,20.0]])
+                msg = 'Predicted flow was %f, should have been %f' %\
+                    (q, ref_flow)
                 if verbose:
-                    print '45 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow)
+                    print '45 degree flowline: ANUGA = %f, Ref = %f' %\
+                        (q, ref_flow)
                     
                 assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
@@ -7155 +7493 @@
                     #if verbose :
                     #    print domain.timestepping_statistics()
-                    #    print domain.volumetric_balance_statistics()                                    
+                    #    print domain.volumetric_balance_statistics()
+