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

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r6653 r6654 781 781 782 782 783 def compute_ volumetric_balance(self):784 """Compute volumetric balanceat current timestep.783 def compute_boundary_flows(self): 784 """Compute boundary flows at current timestep. 785 785 786 786 Quantities computed are: 787 787 Total inflow across boundary 788 788 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 794 790 """ 795 791 … … 797 793 # the normal momentum ((uh, vh) dot normal) times segment length. 798 794 # 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 = stageelevation 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 799 892 # The go through explicit forcing update and record the rate of change for stage and 800 893 # record into forcing_inflow and forcing_outflow. Finally compute integral 801 894 # 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 823 903 824 904 #=============== End of class Shallow Water Domain ===============================
Note: See TracChangeset
for help on using the changeset viewer.