Changeset 6654

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

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

 r6191 """ 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

 r6191 #Check neighbour structure # Check neighbour structure for i in range(N): # For each triangle 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

 r6541 # @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 # 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

 r5897 int _interpolate_from_vertices_to_edges(int N, double* vertex_values, double* edge_values) { double* vertex_values, double* edge_values) { int k, k3; 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

 r6653 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 """ # 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

 r6653 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 #------------------------------------------------------------------------------ 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 def test_friction_dependent_flow_using_flowline(self): def Xtest_friction_dependent_flow_using_flowline(self): """test_friction_dependent_flow_using_flowline if verbose : print domain.timestepping_statistics() print domain.volumetric_balance_statistics() # 90 degree flowline at 200m 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)
Note: See TracChangeset for help on using the changeset viewer.