Changeset 7090
- Timestamp:
- May 26, 2009, 5:41:49 PM (16 years ago)
- Location:
- anuga_core/source/anuga/shallow_water
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r7043 r7090 880 880 """Create volumetric balance report suitable for printing or logging. 881 881 """ 882 883 boundary_flows, total_boundary_inflow, total_boundary_outflow = self.compute_boundary_flows() 882 883 X = self.compute_boundary_flows() 884 boundary_flows, total_boundary_inflow, total_boundary_outflow = X 884 885 885 886 s = '---------------------------\n' 886 887 s += 'Volumetric balance report:\n' 887 888 s += '--------------------------\n' 888 s += 'Total boundary inflow [m^3/s]: %.2 f\n' % total_boundary_inflow889 s += 'Total boundary outflow [m^3/s]: %.2 f\n' % total_boundary_outflow889 s += 'Total boundary inflow [m^3/s]: %.2e\n' % total_boundary_inflow 890 s += 'Total boundary outflow [m^3/s]: %.2e\n' % total_boundary_outflow 890 891 s += 'Net boundary flow by tags [m^3/s]\n' 891 892 for tag in boundary_flows: 892 s += ' %s [m^3/s]: %.2 f\n' % (tag, boundary_flows[tag])893 894 s += 'Total net boundary flow [m^3/s]: %.2 f\n' % (total_boundary_inflow + total_boundary_outflow)895 s += 'Total volume in domain [m^3]: %.2 f\n' % self.compute_total_volume()893 s += ' %s [m^3/s]: %.2e\n' % (tag, boundary_flows[tag]) 894 895 s += 'Total net boundary flow [m^3/s]: %.2e\n' % (total_boundary_inflow + total_boundary_outflow) 896 s += 'Total volume in domain [m^3]: %.2e\n' % self.compute_total_volume() 896 897 897 898 # 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 2279 2280 stage quantity. 2281 2282 Note also that any reference to the internal attribute 'rate' 2283 later will use the one converted to m/s. 2280 2284 2281 2285 polygon: Specifies a polygon to restrict the rainfall. -
anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
r7043 r7090 6745 6745 6746 6746 6747 #---------------------------------------------------------------------- --------6747 #---------------------------------------------------------------------- 6748 6748 # Import necessary modules 6749 #---------------------------------------------------------------------- --------6749 #---------------------------------------------------------------------- 6750 6750 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 6751 6751 from anuga.shallow_water import Domain 6752 6752 6753 6753 6754 #---------------------------------------------------------------------- --------6754 #---------------------------------------------------------------------- 6755 6755 # Setup computational domain 6756 #---------------------------------------------------------------------- --------6756 #---------------------------------------------------------------------- 6757 6757 6758 6758 length = 100. … … 6765 6765 6766 6766 6767 #---------------------------------------------------------------------- --------6767 #---------------------------------------------------------------------- 6768 6768 # Simple flat bottom bathtub 6769 #---------------------------------------------------------------------- --------6769 #---------------------------------------------------------------------- 6770 6770 6771 6771 d = 1.0 … … 6776 6776 6777 6777 6778 #---------------------------------------------------------------------- --------6778 #---------------------------------------------------------------------- 6779 6779 # Slope 6780 #---------------------------------------------------------------------- --------6780 #---------------------------------------------------------------------- 6781 6781 6782 6782 slope = 1.0/10 # RHS drops to -10m … … 6812 6812 6813 6813 6814 #--------------------------------------------------------------------- ---------6814 #--------------------------------------------------------------------- 6815 6815 # Import necessary modules 6816 #--------------------------------------------------------------------- ---------6816 #--------------------------------------------------------------------- 6817 6817 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 6818 6818 from anuga.shallow_water import Domain … … 6822 6822 from anuga.shallow_water.data_manager import get_flow_through_cross_section 6823 6823 6824 #---------------------------------------------------------------------- --------6824 #---------------------------------------------------------------------- 6825 6825 # Setup computational domain 6826 #---------------------------------------------------------------------- --------6826 #---------------------------------------------------------------------- 6827 6827 finaltime = 300.0 6828 6828 … … 6846 6846 6847 6847 6848 #---------------------------------------------------------------------- --------6848 #---------------------------------------------------------------------- 6849 6849 # Setup initial conditions 6850 #---------------------------------------------------------------------- --------6850 #---------------------------------------------------------------------- 6851 6851 slope = 0.0 6852 6852 def topography(x, y): … … 6860 6860 expression='elevation') 6861 6861 6862 #---------------------------------------------------------------------- --------6862 #---------------------------------------------------------------------- 6863 6863 # Setup boundary conditions 6864 #---------------------------------------------------------------------- --------6864 #---------------------------------------------------------------------- 6865 6865 6866 6866 Br = Reflective_boundary(domain) # Solid reflective wall … … 6875 6875 6876 6876 6877 #---------------------------------------------------------------------- --------6877 #---------------------------------------------------------------------- 6878 6878 # Evolve system through time 6879 #---------------------------------------------------------------------- --------6879 #---------------------------------------------------------------------- 6880 6880 for t in domain.evolve(yieldstep=100.0, finaltime=finaltime): 6881 6881 … … 6887 6887 6888 6888 6889 6889 6890 def test_volume_conservation_inflow(self): 6891 """test_volume_conservation 6892 6893 Test that total volume in domain is as expected, based on questions 6894 raised by Petar Milevski in May 2009. 6895 6896 This test adds inflow at a known rate and verifies that the total 6897 terminal volume is as expected. 6898 6899 """ 6900 6901 verbose = False 6902 6903 6904 #--------------------------------------------------------------------- 6905 # Import necessary modules 6906 #--------------------------------------------------------------------- 6907 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 6908 from anuga.shallow_water import Domain 6909 from anuga.shallow_water.shallow_water_domain import Reflective_boundary 6910 from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary 6911 from anuga.shallow_water.shallow_water_domain import Inflow 6912 from anuga.shallow_water.data_manager import get_flow_through_cross_section 6913 6914 #---------------------------------------------------------------------- 6915 # Setup computational domain 6916 #---------------------------------------------------------------------- 6917 finaltime = 200.0 6918 6919 length = 300. 6920 width = 20. 6921 dx = dy = 5 # Resolution: of grid on both axes 6922 6923 6924 points, vertices, boundary = rectangular_cross(int(length/dx), 6925 int(width/dy), 6926 len1=length, len2=width) 6927 6928 6929 domain = Domain(points, vertices, boundary) 6930 domain.set_name('Inflow_volume_test') # Output name 6931 6932 6933 #---------------------------------------------------------------------- 6934 # Setup initial conditions 6935 #---------------------------------------------------------------------- 6936 slope = 0.0 6937 def topography(x, y): 6938 z=-x * slope 6939 return z 6940 6941 domain.set_quantity('elevation', topography) # Use function for elevation 6942 domain.set_quantity('friction', 0.0) # Constant friction 6943 6944 domain.set_quantity('stage', 6945 expression='elevation') # Dry initially 6946 6947 6948 #-------------------------------------------------------------- 6949 # Setup Inflow 6950 #-------------------------------------------------------------- 6951 6952 # Fixed Flowrate onto Area 6953 fixed_inflow = Inflow(domain, 6954 center=(10.0, 10.0), 6955 radius=5.00, 6956 rate=10.00) 6957 6958 domain.forcing_terms.append(fixed_inflow) 6959 6960 #---------------------------------------------------------------------- 6961 # Setup boundary conditions 6962 #---------------------------------------------------------------------- 6963 6964 Br = Reflective_boundary(domain) # Solid reflective wall 6965 domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) 6966 6967 6968 #---------------------------------------------------------------------- 6969 # Evolve system through time 6970 #---------------------------------------------------------------------- 6971 ref_volume = 0.0 6972 ys = 10.0 # Yieldstep 6973 for t in domain.evolve(yieldstep=ys, finaltime=finaltime): 6974 6975 # Check volume 6976 assert num.allclose(domain.compute_total_volume(), ref_volume) 6977 6978 if verbose : 6979 print domain.timestepping_statistics() 6980 print domain.volumetric_balance_statistics() 6981 print 'reference volume', ref_volume 6982 6983 6984 # Update reference volume 6985 ref_volume += ys * fixed_inflow.rate 6986 6987 6988 6989 6990 6991 def test_volume_conservation_rain(self): 6992 """test_volume_conservation 6993 6994 Test that total volume in domain is as expected, based on questions 6995 raised by Petar Milevski in May 2009. 6996 6997 This test adds rain at a known rate and verifies that the total 6998 terminal volume is as expected. 6999 7000 """ 7001 7002 verbose = False 7003 7004 7005 #--------------------------------------------------------------------- 7006 # Import necessary modules 7007 #--------------------------------------------------------------------- 7008 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 7009 from anuga.shallow_water import Domain 7010 from anuga.shallow_water.shallow_water_domain import Reflective_boundary 7011 from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary 7012 from anuga.shallow_water.shallow_water_domain import Rainfall 7013 from anuga.shallow_water.data_manager import get_flow_through_cross_section 7014 7015 #---------------------------------------------------------------------- 7016 # Setup computational domain 7017 #---------------------------------------------------------------------- 7018 finaltime = 200.0 7019 7020 length = 300. 7021 width = 20. 7022 dx = dy = 5 # Resolution: of grid on both axes 7023 7024 7025 points, vertices, boundary = rectangular_cross(int(length/dx), 7026 int(width/dy), 7027 len1=length, len2=width) 7028 7029 7030 domain = Domain(points, vertices, boundary) 7031 domain.set_name('Rain_volume_test') # Output name 7032 7033 7034 #---------------------------------------------------------------------- 7035 # Setup initial conditions 7036 #---------------------------------------------------------------------- 7037 slope = 0.0 7038 def topography(x, y): 7039 z=-x * slope 7040 return z 7041 7042 domain.set_quantity('elevation', topography) # Use function for elevation 7043 domain.set_quantity('friction', 0.0) # Constant friction 7044 7045 domain.set_quantity('stage', 7046 expression='elevation') # Dry initially 7047 7048 7049 #-------------------------------------------------------------- 7050 # Setup rain 7051 #-------------------------------------------------------------- 7052 7053 # Fixed rain onto small circular area 7054 fixed_rain = Rainfall(domain, 7055 center=(10.0, 10.0), 7056 radius=5.00, 7057 rate=10.00) # 10 mm/s 7058 7059 domain.forcing_terms.append(fixed_rain) 7060 7061 #---------------------------------------------------------------------- 7062 # Setup boundary conditions 7063 #---------------------------------------------------------------------- 7064 7065 Br = Reflective_boundary(domain) # Solid reflective wall 7066 domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) 7067 7068 7069 #---------------------------------------------------------------------- 7070 # Evolve system through time 7071 #---------------------------------------------------------------------- 7072 ref_volume = 0.0 7073 ys = 10.0 # Yieldstep 7074 for t in domain.evolve(yieldstep=ys, finaltime=finaltime): 7075 7076 # Check volume 7077 V = domain.compute_total_volume() 7078 msg = 'V = %e, Ref = %e' % (V, ref_volume) 7079 assert num.allclose(V, ref_volume), msg 7080 7081 if verbose : 7082 print domain.timestepping_statistics() 7083 print domain.volumetric_balance_statistics() 7084 print 'reference volume', ref_volume 7085 print V 7086 7087 7088 # Update reference volume. 7089 # FIXME: Note that rate has now been redefined 7090 # as m/s internally. This is a little confusing 7091 # when it was specfied as mm/s. 7092 7093 delta_V = fixed_rain.rate*fixed_rain.exchange_area 7094 ref_volume += ys * delta_V 7095 7096 7097 7098 7099 def Xtest_rain_conservation_and_runoff(self): 7100 """test_rain_conservation_and_runoff 7101 7102 Test that total volume in domain is as expected, based on questions 7103 raised by Petar Milevski in May 2009. 7104 7105 This test adds rain at a known rate and verifies that the total 7106 volume and outflows are as expected. 7107 7108 """ 7109 7110 # FIXME (Ole): Does not work yet. Investigate boundary flows 7111 7112 verbose = True #False 7113 7114 7115 #--------------------------------------------------------------------- 7116 # Import necessary modules 7117 #--------------------------------------------------------------------- 7118 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 7119 from anuga.shallow_water import Domain 7120 from anuga.shallow_water.shallow_water_domain import Reflective_boundary 7121 from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary 7122 from anuga.shallow_water.shallow_water_domain import Rainfall 7123 from anuga.shallow_water.data_manager import get_flow_through_cross_section 7124 7125 #---------------------------------------------------------------------- 7126 # Setup computational domain 7127 #---------------------------------------------------------------------- 7128 finaltime = 500.0 7129 7130 length = 300. 7131 width = 20. 7132 dx = dy = 5 # Resolution: of grid on both axes 7133 7134 7135 points, vertices, boundary = rectangular_cross(int(length/dx), 7136 int(width/dy), 7137 len1=length, len2=width) 7138 7139 7140 domain = Domain(points, vertices, boundary) 7141 domain.set_name('Rain_volume_runoff_test') # Output name 7142 7143 7144 #---------------------------------------------------------------------- 7145 # Setup initial conditions 7146 #---------------------------------------------------------------------- 7147 slope = 0.0 7148 def topography(x, y): 7149 z=-x * slope 7150 return z 7151 7152 domain.set_quantity('elevation', topography) # Use function for elevation 7153 domain.set_quantity('friction', 0.0) # Constant friction 7154 7155 domain.set_quantity('stage', 7156 expression='elevation') # Dry initially 7157 7158 7159 #-------------------------------------------------------------- 7160 # Setup rain 7161 #-------------------------------------------------------------- 7162 7163 # Fixed rain onto small circular area 7164 fixed_rain = Rainfall(domain, 7165 center=(10.0, 10.0), 7166 radius=5.00, 7167 rate=10.00) # 10 mm/s 7168 7169 domain.forcing_terms.append(fixed_rain) 7170 7171 #---------------------------------------------------------------------- 7172 # Setup boundary conditions 7173 #---------------------------------------------------------------------- 7174 7175 Br = Reflective_boundary(domain) # Solid reflective wall 7176 Bt = Transmissive_stage_zero_momentum_boundary(domain) 7177 Bd = Dirichlet_boundary([-10, 0, 0]) 7178 domain.set_boundary({'left': Bt, 'right': Bd, 'top': Bt, 'bottom': Bt}) 7179 7180 7181 #---------------------------------------------------------------------- 7182 # Evolve system through time 7183 #---------------------------------------------------------------------- 7184 ref_volume = 0.0 7185 ys = 10.0 # Yieldstep 7186 for t in domain.evolve(yieldstep=ys, finaltime=finaltime): 7187 7188 # Check volume 7189 V = domain.compute_total_volume() 7190 msg = 'V = %e, Ref = %e' % (V, ref_volume) 7191 #assert num.allclose(V, ref_volume) or V < ref_volume, msg 7192 7193 if verbose: 7194 print domain.timestepping_statistics() 7195 print domain.volumetric_balance_statistics() 7196 print 'reference volume', ref_volume 7197 print V 7198 7199 7200 # Update reference volume. 7201 # FIXME: Note that rate has now been redefined 7202 # as m/s internally. This is a little confusing 7203 # when it was specfied as mm/s. 7204 7205 delta_V = fixed_rain.rate*fixed_rain.exchange_area 7206 ref_volume += ys * delta_V 7207 7208 # Compute outflow at right hand downstream boundary 7209 boundary_flows, inflow , outflow = domain.compute_boundary_flows() 7210 net_outflow = outflow - inflow 7211 7212 outflow = boundary_flows['right'] 7213 if verbose: 7214 print 'Outflow', outflow 7215 print 'Net outflow', net_outflow 7216 7217 # Update reference volume 7218 ref_volume += ys * outflow 7219 7220 7221 6890 7222 6891 7223 def test_inflow_using_flowline(self): … … 6929 7261 dx = dy = 5 # Resolution: of grid on both axes 6930 7262 6931 points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy), 7263 points, vertices, boundary = rectangular_cross(int(length/dx), 7264 int(width/dy), 6932 7265 len1=length, len2=width) 6933 7266 … … 7013 7346 #if verbose : 7014 7347 # print domain.timestepping_statistics() 7015 # print domain.volumetric_balance_statistics() 7016 7348 # print domain.volumetric_balance_statistics() 7017 7349 7018 7350 #-------------------------------------------------------------- … … 7021 7353 7022 7354 # Square on flowline at 200m 7023 q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,20.0]]) 7024 msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow) 7355 q=domain.get_flow_through_cross_section([[200.0,0.0], 7356 [200.0,20.0]]) 7357 msg = 'Predicted flow was %f, should have been %f' %\ 7358 (q, ref_flow) 7025 7359 if verbose: 7026 print '90 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow) 7360 print '90 degree flowline: ANUGA = %f, Ref = %f' %\ 7361 (q, ref_flow) 7027 7362 assert num.allclose(q, ref_flow, rtol=1.0e-2), msg 7028 7363 7029 7364 7030 7365 # 45 degree flowline at 200m 7031 q=domain.get_flow_through_cross_section([[200.0,0.0],[220.0,20.0]]) 7032 msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow) 7366 q=domain.get_flow_through_cross_section([[200.0,0.0], 7367 [220.0,20.0]]) 7368 msg = 'Predicted flow was %f, should have been %f' %\ 7369 (q, ref_flow) 7033 7370 if verbose: 7034 print '45 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow) 7371 print '45 degree flowline: ANUGA = %f, Ref = %f' %\ 7372 (q, ref_flow) 7035 7373 7036 7374 assert num.allclose(q, ref_flow, rtol=1.0e-2), msg … … 7155 7493 #if verbose : 7156 7494 # print domain.timestepping_statistics() 7157 # print domain.volumetric_balance_statistics() 7495 # print domain.volumetric_balance_statistics() 7496 7158 7497 7159 7498
Note: See TracChangeset
for help on using the changeset viewer.