Changeset 6654


Ignore:
Timestamp:
Mar 28, 2009, 6:54:15 PM (15 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:
6 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

    r6191 r6654  
    235235        """
    236236        return self.normals[i, 2*j:2*j+2]
     237       
     238    def get_edgelength(self, i, j):
     239        """Return length of j'th edge of the i'th triangle.
     240        Return value is the numeric array slice [x, y]
     241        """
     242        return self.edgelengths[i, j]
     243               
    237244
    238245    def get_number_of_nodes(self):
  • anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

    r6191 r6654  
    687687
    688688
    689         #Check neighbour structure
     689        # Check neighbour structure
    690690        for i in range(N):
    691691            # For each triangle
     
    869869
    870870    def get_triangle_containing_point(self, point):
    871         """Return triangle id for triangle containing specifiend point (x,y)
     871        """Return triangle id for triangle containing specified point (x,y)
    872872
    873873        If point isn't within mesh, raise exception
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r6541 r6654  
    11111111    # @param use_cache
    11121112    # @param verbose True if this method is to be verbose.
    1113     def get_values(self, interpolation_points=None,
    1114                          location='vertices',
    1115                          indices=None,
    1116                          use_cache=False,
    1117                          verbose=False):
     1113    def get_values(self,
     1114                   interpolation_points=None,
     1115                   location='vertices',
     1116                   indices=None,
     1117                   use_cache=False,
     1118                   verbose=False):
    11181119        """Get values for quantity
    11191120
     
    11681169        # Edges have already been deprecated in set_values, see changeset:5521,
    11691170        # but *might* be useful in get_values. Any thoughts anyone?
     1171        # YES (Ole): Edge values are necessary for volumetric balance check
    11701172
    11711173        if location not in ['vertices', 'centroids',
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r5897 r6654  
    607607
    608608int _interpolate_from_vertices_to_edges(int N,
    609                  double* vertex_values,
    610                  double* edge_values) {
     609                                        double* vertex_values,
     610                                        double* edge_values) {
    611611
    612612        int k, k3;
     
    630630
    631631int _interpolate_from_edges_to_vertices(int N,
    632                  double* vertex_values,
    633                  double* edge_values) {
     632                                        double* vertex_values,
     633                                        double* edge_values) {
    634634
    635635        int k, k3;
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r6653 r6654  
    781781       
    782782
    783     def compute_volumetric_balance(self):
    784         """Compute volumetric balance at current timestep.
     783    def compute_boundary_flows(self):
     784        """Compute boundary flows at current timestep.
    785785                       
    786786        Quantities computed are:
    787787           Total inflow across boundary
    788788           Total outflow across boundary
    789            Total inflow through forcing terms
    790            Total outflow through forcing terms
    791            Current total volume in domain
    792        
    793                        
     789           Flow across each tagged boundary segment
    794790        """
    795791               
     
    797793        # the normal momentum ((uh, vh) dot normal) times segment length.
    798794        # Based on sign accumulate this into boundary_inflow and boundary_outflow.
     795                       
     796        # Compute flows along boundary
     797       
     798        uh = self.get_quantity('xmomentum').get_values(location='edges')
     799        vh = self.get_quantity('ymomentum').get_values(location='edges')       
     800       
     801        # Loop through edges that lie on the boundary and calculate
     802        # flows
     803        boundary_flows = {}
     804        total_boundary_inflow = 0.0
     805        total_boundary_outflow = 0.0
     806        for vol_id, edge_id in self.boundary:
     807            # Compute normal flow across edge. Since normal vector points
     808            # away from triangle, a positive sign means that water
     809            # flows *out* from this triangle.
     810             
     811            momentum = [uh[vol_id, edge_id], vh[vol_id, edge_id]]
     812            normal = self.mesh.get_normal(vol_id, edge_id)
     813            length = self.mesh.get_edgelength(vol_id, edge_id)           
     814            normal_flow = num.dot(momentum, normal)*length
     815           
     816            # Reverse sign so that + is taken to mean inflow
     817            # and - means outflow. This is more intuitive.
     818            edge_flow = -normal_flow
     819           
     820            # Tally up inflows and outflows separately
     821            if edge_flow > 0:
     822                # Flow is inflow     
     823                total_boundary_inflow += edge_flow                                 
     824            else:
     825                # Flow is outflow
     826                total_boundary_outflow += edge_flow   
     827
     828            # Tally up flows by boundary tag
     829            tag = self.boundary[(vol_id, edge_id)]
     830           
     831            if tag not in boundary_flows:
     832                boundary_flows[tag] = 0.0
     833            boundary_flows[tag] += edge_flow
     834           
     835               
     836        return boundary_flows, total_boundary_inflow, total_boundary_outflow
     837       
     838
     839    def compute_forcing_flows(self):
     840        """Compute flows in and out of domain due to forcing terms.
     841                       
     842        Quantities computed are:
     843               
     844       
     845           Total inflow through forcing terms
     846           Total outflow through forcing terms
     847           Current total volume in domain       
     848
     849        """
     850
     851        #FIXME(Ole): We need to separate what part of explicit_update was
     852        # due to the normal flux calculations and what is due to forcing terms.
     853       
     854        pass
     855                       
     856       
     857    def compute_total_volume(self):
     858        """Compute total volume (m^3) of water in entire domain
     859        """
     860       
     861        area = self.mesh.get_areas()
     862        volume = 0.0
     863       
     864        stage = self.get_quantity('stage').get_values(location='centroids')
     865        elevation = self.get_quantity('elevation').get_values(location='centroids')       
     866        depth = stage-elevation
     867       
     868        #print 'z', elevation
     869        #print 'w', stage
     870        #print 'h', depth
     871        return num.sum(depth*area)
     872       
     873       
     874    def volumetric_balance_statistics(self):               
     875        """Create volumetric balance report suitable for printing or logging.
     876        """
     877       
     878        boundary_flows, total_boundary_inflow, total_boundary_outflow = self.compute_boundary_flows()
     879       
     880        s = '---------------------------\n'       
     881        s += 'Volumetric balance report:\n'
     882        s += '--------------------------\n'
     883        s += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
     884        s += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow       
     885        s += 'Net boundary flow by tags [m^3/s]\n'
     886        for tag in boundary_flows:
     887            s += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
     888       
     889        s += 'Total net boundary flow [m^3/s]: %.2f\n' % (total_boundary_inflow + total_boundary_outflow)
     890        s += 'Total volume in domain [m^3]: %.2f\n' % self.compute_total_volume()
     891       
    799892        # The go through explicit forcing update and record the rate of change for stage and
    800893        # record into forcing_inflow and forcing_outflow. Finally compute integral
    801894        # of depth to obtain total volume of domain.
    802                        
    803                        
    804         # Compute flows along boundary
    805        
    806         uh = self.get_quantity('xmomentum').get_values()
    807         vh = self.get_quantity('ymomentum').get_values()       
    808        
    809         # Loop through edges that lie on the boundary and calculate
    810         # flows
    811         inflow = 0.0
    812         outflow = 0.0
    813         for vol_id, edge_id in self.boundary:
    814             print vol_id, edge_id, self.boundary[(vol_id, edge_id)]
    815 
    816             # Pick edge and compute normal flow
    817             print uh[vol_id, :]
    818             print vh[vol_id, :]           
    819          
    820        
    821        
    822                
     895       
     896        # FIXME(Ole): This part is not yet done.               
     897       
     898        return s       
     899           
     900           
     901           
     902           
    823903               
    824904#=============== End of class Shallow Water Domain ===============================
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r6653 r6654  
    67366736
    67376737           
     6738    def test_total_volume(self):       
     6739        """test_total_volume
     6740       
     6741        Test that total volume can be computed correctly
     6742        """           
     6743
     6744       
     6745
     6746        #------------------------------------------------------------------------------
     6747        # Import necessary modules
     6748        #------------------------------------------------------------------------------
     6749        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     6750        from anuga.shallow_water import Domain
     6751
     6752
     6753        #------------------------------------------------------------------------------
     6754        # Setup computational domain
     6755        #------------------------------------------------------------------------------
     6756
     6757        length = 100.
     6758        width  = 20.
     6759        dx = dy = 5       # Resolution: of grid on both axes
     6760       
     6761        points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
     6762                                                       len1=length, len2=width)
     6763        domain = Domain(points, vertices, boundary)   
     6764               
     6765
     6766        #------------------------------------------------------------------------------
     6767        # Simple flat bottom bathtub
     6768        #------------------------------------------------------------------------------
     6769
     6770        d = 1.0
     6771        domain.set_quantity('elevation', 0.0)
     6772        domain.set_quantity('stage', d)
     6773       
     6774        assert num.allclose(domain.compute_total_volume(), length*width*d)
     6775
     6776
     6777        #------------------------------------------------------------------------------
     6778        # Slope
     6779        #------------------------------------------------------------------------------
     6780               
     6781        slope = 1.0/10  # RHS drops to -10m
     6782        def topography(x, y):
     6783            return -x * slope
     6784
     6785        domain.set_quantity('elevation', topography)
     6786        domain.set_quantity('stage', 0.0) # Domain full
     6787       
     6788        V = domain.compute_total_volume()
     6789        assert num.allclose(V, length*width*10/2)
     6790
     6791       
     6792        domain.set_quantity('stage', -5.0) # Domain 'half' full
     6793       
     6794        # IMPORTANT: Adjust stage to match elevation
     6795        domain.distribute_to_vertices_and_edges()
     6796       
     6797        V = domain.compute_total_volume()
     6798        assert num.allclose(V, width*(length/2)*5.0/2)
     6799       
     6800               
     6801       
    67386802           
    6739            
    6740     def Xtest_volumetric_balance_computation(self):
     6803    def test_volumetric_balance_computation(self):
    67416804        """test_volumetric_balance_computation
    67426805       
    67436806        Test that total in and out flows are computed correctly
     6807        FIXME(Ole): This test is more about looking at the printed report
    67446808        """
    67456809
    6746         verbose = True
     6810        verbose = False
    67476811       
    67486812
     
    68146878        #------------------------------------------------------------------------------
    68156879        for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
     6880       
     6881            S = domain.volumetric_balance_statistics()
    68166882            if verbose :
    68176883                print domain.timestepping_statistics()
    6818                                    
    6819             print domain.compute_volumetric_balance()
     6884                print S
    68206885           
    68216886           
     
    69627027       
    69637028       
    6964     def test_friction_dependent_flow_using_flowline(self):
     7029    def Xtest_friction_dependent_flow_using_flowline(self):
    69657030        """test_friction_dependent_flow_using_flowline
    69667031       
     
    70637128                    if verbose :
    70647129                        print domain.timestepping_statistics()
     7130                        print domain.volumetric_balance_statistics()                       
    70657131                                   
    70667132                # 90 degree flowline at 200m                                   
     
    71027168                               
    71037169if __name__ == "__main__":
    7104     #suite = unittest.makeSuite(Test_Shallow_Water, 'test_friction_dependent_flow_using_flowline')
     7170    suite = unittest.makeSuite(Test_Shallow_Water, 'test_friction_dependent_flow_using_flowline')
    71057171    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_volumetric_balance_computation')
    7106     suite = unittest.makeSuite(Test_Shallow_Water, 'test')       
     7172    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_total_volume')           
     7173    #suite = unittest.makeSuite(Test_Shallow_Water, 'test')       
    71077174
    71087175    runner = unittest.TextTestRunner(verbosity=1)   
Note: See TracChangeset for help on using the changeset viewer.