Changeset 6654

Timestamp:

Mar 28, 2009, 6:54:15 PM (16 years ago)

Author:

ole

Message:

Added volumetric balance report for boundary flows and total volume.
Flows due to forcing terms have not yet been added.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/general_mesh.py (modified) (1 diff)
• abstract_2d_finite_volumes/neighbour_mesh.py (modified) (2 diffs)
• abstract_2d_finite_volumes/quantity.py (modified) (2 diffs)
• abstract_2d_finite_volumes/quantity_ext.c (modified) (2 diffs)
• shallow_water/shallow_water_domain.py (modified) (2 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (5 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -235 +235 @@
         """
         return self.normals[i, 2*j:2*j+2]
+        
+    def get_edgelength(self, i, j):
+        """Return length of j'th edge of the i'th triangle.
+        Return value is the numeric array slice [x, y]
+        """
+        return self.edgelengths[i, j]
+                
 
     def get_number_of_nodes(self):

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -687 +687 @@
 
 
-        #Check neighbour structure
+        # Check neighbour structure
         for i in range(N):
             # For each triangle
@@ -869 +869 @@
 
     def get_triangle_containing_point(self, point):
-        """Return triangle id for triangle containing specifiend point (x,y)
+        """Return triangle id for triangle containing specified point (x,y)
 
         If point isn't within mesh, raise exception

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -1111 +1111 @@
     # @param use_cache 
     # @param verbose True if this method is to be verbose.
-    def get_values(self, interpolation_points=None,
-                         location='vertices',
-                         indices=None,
-                         use_cache=False,
-                         verbose=False):
+    def get_values(self, 
+                   interpolation_points=None,
+                   location='vertices',
+                   indices=None,
+                   use_cache=False,
+                   verbose=False):
         """Get values for quantity
 
@@ -1168 +1169 @@
         # Edges have already been deprecated in set_values, see changeset:5521,
         # but *might* be useful in get_values. Any thoughts anyone?
+        # YES (Ole): Edge values are necessary for volumetric balance check
 
         if location not in ['vertices', 'centroids',

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

@@ -607 +607 @@
 
 int _interpolate_from_vertices_to_edges(int N,
-                 double* vertex_values,
-                 double* edge_values) {
+                                        double* vertex_values,
+                                        double* edge_values) {
 
         int k, k3;
@@ -630 +630 @@
 
 int _interpolate_from_edges_to_vertices(int N,
-                 double* vertex_values,
-                 double* edge_values) {
+                                        double* vertex_values,
+                                        double* edge_values) {
 
         int k, k3;

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -781 +781 @@
         
 
-    def compute_volumetric_balance(self):
-        """Compute volumetric balance at current timestep.
+    def compute_boundary_flows(self):
+        """Compute boundary flows at current timestep.
                         
         Quantities computed are:
            Total inflow across boundary
            Total outflow across boundary
-           Total inflow through forcing terms
-           Total outflow through forcing terms
-           Current total volume in domain
-        
-                        
+           Flow across each tagged boundary segment
         """
                 
@@ -797 +793 @@
         # the normal momentum ((uh, vh) dot normal) times segment length.
         # Based on sign accumulate this into boundary_inflow and boundary_outflow.
+                        
+        # Compute flows along boundary
+        
+        uh = self.get_quantity('xmomentum').get_values(location='edges')
+        vh = self.get_quantity('ymomentum').get_values(location='edges')        
+        
+        # Loop through edges that lie on the boundary and calculate 
+        # flows
+        boundary_flows = {}
+        total_boundary_inflow = 0.0
+        total_boundary_outflow = 0.0
+        for vol_id, edge_id in self.boundary:
+            # Compute normal flow across edge. Since normal vector points
+            # away from triangle, a positive sign means that water 
+            # flows *out* from this triangle.
+             
+            momentum = [uh[vol_id, edge_id], vh[vol_id, edge_id]]
+            normal = self.mesh.get_normal(vol_id, edge_id)
+            length = self.mesh.get_edgelength(vol_id, edge_id)            
+            normal_flow = num.dot(momentum, normal)*length
+            
+            # Reverse sign so that + is taken to mean inflow
+            # and - means outflow. This is more intuitive.
+            edge_flow = -normal_flow
+            
+            # Tally up inflows and outflows separately
+            if edge_flow > 0:
+                # Flow is inflow      
+                total_boundary_inflow += edge_flow                                  
+            else:
+                # Flow is outflow
+                total_boundary_outflow += edge_flow    
+
+            # Tally up flows by boundary tag
+            tag = self.boundary[(vol_id, edge_id)]
+            
+            if tag not in boundary_flows:
+                boundary_flows[tag] = 0.0
+            boundary_flows[tag] += edge_flow
+            
+                
+        return boundary_flows, total_boundary_inflow, total_boundary_outflow
+        
+
+    def compute_forcing_flows(self):
+        """Compute flows in and out of domain due to forcing terms.
+                        
+        Quantities computed are:
+                
+        
+           Total inflow through forcing terms
+           Total outflow through forcing terms
+           Current total volume in domain        
+
+        """
+
+        #FIXME(Ole): We need to separate what part of explicit_update was 
+        # due to the normal flux calculations and what is due to forcing terms.
+        
+        pass
+                        
+        
+    def compute_total_volume(self):
+        """Compute total volume (m^3) of water in entire domain
+        """
+        
+        area = self.mesh.get_areas()
+        volume = 0.0
+        
+        stage = self.get_quantity('stage').get_values(location='centroids')
+        elevation = self.get_quantity('elevation').get_values(location='centroids')        
+        depth = stage-elevation
+        
+        #print 'z', elevation
+        #print 'w', stage
+        #print 'h', depth
+        return num.sum(depth*area)
+        
+        
+    def volumetric_balance_statistics(self):                
+        """Create volumetric balance report suitable for printing or logging.
+        """
+        
+        boundary_flows, total_boundary_inflow, total_boundary_outflow = self.compute_boundary_flows() 
+        
+        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 += '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()
+        
         # The go through explicit forcing update and record the rate of change for stage and 
         # record into forcing_inflow and forcing_outflow. Finally compute integral 
         # of depth to obtain total volume of domain.
-                        
-                        
-        # Compute flows along boundary
-        
-        uh = self.get_quantity('xmomentum').get_values()
-        vh = self.get_quantity('ymomentum').get_values()        
-        
-        # Loop through edges that lie on the boundary and calculate 
-        # flows
-        inflow = 0.0
-        outflow = 0.0
-        for vol_id, edge_id in self.boundary:
-            print vol_id, edge_id, self.boundary[(vol_id, edge_id)]
-
-            # Pick edge and compute normal flow
-            print uh[vol_id, :]
-            print vh[vol_id, :]            
-         
-        
-        
-                
+        
+        # FIXME(Ole): This part is not yet done.                
+        
+        return s        
+           
+           
+            
+            
                 
 #=============== End of class Shallow Water Domain ===============================

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -6736 +6736 @@
 
             
+    def test_total_volume(self):        
+        """test_total_volume
+        
+        Test that total volume can be computed correctly
+        """            
+
+        
+
+        #------------------------------------------------------------------------------
+        # 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.
+        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)   
+                
+
+        #------------------------------------------------------------------------------
+        # Simple flat bottom bathtub
+        #------------------------------------------------------------------------------
+
+        d = 1.0
+        domain.set_quantity('elevation', 0.0)
+        domain.set_quantity('stage', d)
+        
+        assert num.allclose(domain.compute_total_volume(), length*width*d)
+
+
+        #------------------------------------------------------------------------------
+        # Slope
+        #------------------------------------------------------------------------------
+                
+        slope = 1.0/10  # RHS drops to -10m
+        def topography(x, y):
+            return -x * slope
+
+        domain.set_quantity('elevation', topography)
+        domain.set_quantity('stage', 0.0) # Domain full
+        
+        V = domain.compute_total_volume()
+        assert num.allclose(V, length*width*10/2)
+
+        
+        domain.set_quantity('stage', -5.0) # Domain 'half' full
+        
+        # IMPORTANT: Adjust stage to match elevation
+        domain.distribute_to_vertices_and_edges()
+        
+        V = domain.compute_total_volume()
+        assert num.allclose(V, width*(length/2)*5.0/2)
+        
+                
+        
             
-            
-    def Xtest_volumetric_balance_computation(self):
+    def test_volumetric_balance_computation(self):
         """test_volumetric_balance_computation
         
         Test that total in and out flows are computed correctly
+        FIXME(Ole): This test is more about looking at the printed report
         """
 
-        verbose = True
+        verbose = False
         
 
@@ -6814 +6878 @@
         #------------------------------------------------------------------------------
         for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
+        
+            S = domain.volumetric_balance_statistics()
             if verbose :
                 print domain.timestepping_statistics()
-                                    
-            print domain.compute_volumetric_balance()
+                print S
             
             
@@ -6962 +7027 @@
         
         
-    def test_friction_dependent_flow_using_flowline(self):
+    def Xtest_friction_dependent_flow_using_flowline(self):
         """test_friction_dependent_flow_using_flowline
         
@@ -7063 +7128 @@
                     if verbose :
                         print domain.timestepping_statistics()
+                        print domain.volumetric_balance_statistics()                        
                                     
                 # 90 degree flowline at 200m                                    
@@ -7102 +7168 @@
                                 
 if __name__ == "__main__":
-    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_friction_dependent_flow_using_flowline')
+    suite = unittest.makeSuite(Test_Shallow_Water, 'test_friction_dependent_flow_using_flowline')
     #suite = unittest.makeSuite(Test_Shallow_Water, 'test_volumetric_balance_computation') 
-    suite = unittest.makeSuite(Test_Shallow_Water, 'test')        
+    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_total_volume')            
+    #suite = unittest.makeSuite(Test_Shallow_Water, 'test')        
 
     runner = unittest.TextTestRunner(verbosity=1)