- Timestamp:
- May 26, 2009, 5:41:49 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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.
