Changeset 5544
 Timestamp:
 Jul 22, 2008, 4:44:03 AM (15 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/shallow_water/test_data_manager.py
r5541 r5544 5977 5977 os.remove(sww_file) 5978 5978 5979 def write_mux2(self, lat_long_points, time_step_count, time_step,5979 def write_mux2(self, lat_long_points, time_step_count, time_step, 5980 5980 first_tstep, last_tstep, 5981 5981 depth=None, ha=None, ua=None, va=None): … … 6547 6547 """ 6548 6548 6549 tide =06550 time_step_count = 36549 tide = 0.35 6550 time_step_count = 6 # I made this a little longer (Ole) 6551 6551 time_step = 2 6552 6552 lat_long_points =[(21.5,114.5),(21,114.5),(21.5,115), (21.,115.)] … … 6682 6682 6683 6683 6684 assert allclose(2.0*transpose(ha_permutation) , stage) # Meters6684 assert allclose(2.0*transpose(ha_permutation)+tide, stage) # Meters 6685 6685 6686 6686 #Check the momentums  ua … … 6794 6794 6795 6795 6796 6797 6798 def test_urs2sts_file_boundary_stsI(self): 6796 6797 6798 def test_urs2sts_ordering_different_sources(self): 6799 """Test multiple sources with ordering file, different source files and weights. 6800 This test also has more variable values than the previous ones 6801 """ 6802 6803 from Numeric import sin, cos 6804 6805 tide = 1.5 6806 time_step_count = 10 6807 time_step = 0.2 6808 6809 times_ref = arange(0, time_step_count*time_step, time_step) 6810 #print 'time vector', times_ref 6811 6812 lat_long_points = [(21.5,114.5), (21,114.5), (21.5,115), (21.,115.), (22., 117.)] 6813 n = len(lat_long_points) 6814 6815 # Create nontrivial weights 6816 weights = [0.8, 1.5] 6817 6818 # Create different timeseries starting and ending at different times 6819 first_tstep=ones(n,Int) 6820 first_tstep[0]+=2 # Point 0 starts at 2 6821 first_tstep[1]+=4 # Point 1 starts at 4 6822 first_tstep[2]+=3 # Point 2 starts at 3 6823 6824 last_tstep=(time_step_count)*ones(n,Int) 6825 last_tstep[0]=1 # Point 0 ends 1 step early 6826 last_tstep[1]=2 # Point 1 ends 2 steps early 6827 last_tstep[4]=3 # Point 4 ends 3 steps early 6828 6829 #print 6830 #print 'time_step_count', time_step_count 6831 #print 'time_step', time_step 6832 #print 'first_tstep', first_tstep 6833 #print 'last_tstep', last_tstep 6834 6835 6836 # Create varying elevation data (positive values for seafloor) 6837 gauge_depth=20*ones(n,Float) 6838 for i in range(n): 6839 gauge_depth[i] += i**2 6840 6841 #print 'gauge_depth', gauge_depth 6842 6843 # Create data to be written to first mux file 6844 ha0=2*ones((n,time_step_count),Float) 6845 ha0[0]=arange(0,time_step_count) 6846 ha0[1]=arange(time_step_count,2*time_step_count) 6847 ha0[2]=arange(2*time_step_count,3*time_step_count) 6848 ha0[3]=arange(3*time_step_count,4*time_step_count) 6849 ua0=5*ones((n,time_step_count),Float) 6850 va0=10*ones((n,time_step_count),Float) 6851 6852 # Ensure data used to write mux file to be zero when gauges are 6853 # not recording 6854 for i in range(n): 6855 # For each point 6856 6857 for j in range(0, first_tstep[i]1) + range(last_tstep[i], time_step_count): 6858 # For timesteps before and after recording range 6859 ha0[i][j] = ua0[i][j] = va0[i][j] = 0.0 6860 6861 6862 #print 6863 #print 'using varying start and end time' 6864 #print 'ha0', ha0 6865 #print 'ua0', ua0 6866 #print 'va0', va0 6867 6868 # Write first mux file to be combined by urs2sts 6869 base_nameI, filesI = self.write_mux2(lat_long_points, 6870 time_step_count, time_step, 6871 first_tstep, last_tstep, 6872 depth=gauge_depth, 6873 ha=ha0, 6874 ua=ua0, 6875 va=va0) 6876 6877 6878 6879 # Create data to be written to second mux file 6880 ha1=zeros((n,time_step_count),Float) 6881 ha1[0]=sin(times_ref) 6882 ha1[1]=2*sin(times_ref  3) 6883 ha1[2]=5*sin(4*times_ref) 6884 ha1[3]=sin(times_ref) 6885 ha1[4]=sin(2*times_ref0.7) 6886 6887 ua1=zeros((n,time_step_count),Float) 6888 ua1[0]=3*cos(times_ref) 6889 ua1[1]=2*sin(times_ref0.7) 6890 ua1[2]=arange(3*time_step_count,4*time_step_count) 6891 ua1[4]=2*ones(time_step_count) 6892 6893 va1=zeros((n,time_step_count),Float) 6894 va1[0]=2*cos(times_ref0.87) 6895 va1[1]=3*ones(time_step_count) 6896 va1[3]=2*sin(times_ref0.71) 6897 6898 6899 # Ensure data used to write mux file to be zero when gauges are 6900 # not recording 6901 for i in range(n): 6902 # For each point 6903 6904 for j in range(0, first_tstep[i]1) + range(last_tstep[i], time_step_count): 6905 # For timesteps before and after recording range 6906 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0 6907 6908 6909 #print 6910 #print 'using varying start and end time' 6911 #print 'ha1', ha1 6912 #print 'ua1', ua1 6913 #print 'va1', va1 6914 6915 6916 # Write second mux file to be combined by urs2sts 6917 base_nameII, filesII = self.write_mux2(lat_long_points, 6918 time_step_count, time_step, 6919 first_tstep, last_tstep, 6920 depth=gauge_depth, 6921 ha=ha1, 6922 ua=ua1, 6923 va=va1) 6924 6925 6926 # Create ordering file 6927 permutation = [4,0,2] 6928 6929 _, ordering_filename = tempfile.mkstemp('') 6930 order_fid = open(ordering_filename, 'w') 6931 order_fid.write('index, longitude, latitude\n') 6932 for index in permutation: 6933 order_fid.write('%d, %f, %f\n' %(index, 6934 lat_long_points[index][1], 6935 lat_long_points[index][0])) 6936 order_fid.close() 6937 6938 6939 6940 6941 # Call urs2sts with multiple mux files 6942 urs2sts([base_nameI, base_nameII], 6943 basename_out=base_nameI, 6944 ordering_filename=ordering_filename, 6945 weights=weights, 6946 mean_stage=tide, 6947 verbose=False) 6948 6949 # Now read the sts file and check that values have been stored correctly. 6950 sts_file = base_nameI + '.sts' 6951 fid = NetCDFFile(sts_file) 6952 6953 # Make x and y absolute 6954 x = fid.variables['x'][:] 6955 y = fid.variables['y'][:] 6956 6957 geo_reference = Geo_reference(NetCDFObject=fid) 6958 points = geo_reference.get_absolute(map(None, x, y)) 6959 points = ensure_numeric(points) 6960 6961 x = points[:,0] 6962 y = points[:,1] 6963 6964 for i, index in enumerate(permutation): 6965 # Check that STS points are stored in the correct order 6966 6967 # Work out the UTM coordinates sts point i 6968 zone, e, n = redfearn(lat_long_points[index][0], 6969 lat_long_points[index][1]) 6970 6971 #print i, [x[i],y[i]], [e,n] 6972 assert allclose([x[i],y[i]], [e,n]) 6973 6974 6975 # Check the time vector 6976 times = fid.variables['time'][:] 6977 assert allclose(ensure_numeric(times), 6978 ensure_numeric(times_ref)) 6979 6980 6981 # Check sts values 6982 stage = fid.variables['stage'][:] 6983 xmomentum = fid.variables['xmomentum'][:] 6984 ymomentum = fid.variables['ymomentum'][:] 6985 elevation = fid.variables['elevation'][:] 6986 6987 6988 # The quantities stored in the .sts file should be the weighted sum of the 6989 # quantities written to the mux2 files subject to the permutation vector. 6990 6991 ha_ref = weights[0]*take(ha0, permutation) + weights[1]*take(ha1, permutation) 6992 ua_ref = weights[0]*take(ua0, permutation) + weights[1]*take(ua1, permutation) 6993 va_ref = weights[0]*take(va0, permutation) + weights[1]*take(va1, permutation) 6994 6995 gauge_depth_ref = take(gauge_depth, permutation) 6996 6997 6998 assert allclose(transpose(ha_ref)+tide, stage) # Meters 6999 7000 #Check the momentums  ua 7001 #momentum = velocity*(stageelevation) 7002 # elevation =  depth 7003 #momentum = velocity_ua *(stage+depth) 7004 7005 depth_ref = zeros((len(permutation), time_step_count), Float) 7006 for i in range(len(permutation)): 7007 depth_ref[i]=gauge_depth_ref[i]+tide+ha_ref[i] 7008 7009 7010 7011 7012 # The xmomentum stored in the .sts file should be the sum of the ua 7013 # in the two mux2 files multiplied by the depth. 7014 assert allclose(transpose(ua_ref*depth_ref), xmomentum) 7015 7016 #Check the momentums  va 7017 #momentum = velocity*(stageelevation) 7018 # elevation =  depth 7019 #momentum = velocity_va *(stage+depth) 7020 7021 # The ymomentum stored in the .sts file should be the sum of the va 7022 # in the two mux2 files multiplied by the depth. 7023 assert allclose(transpose(va_ref*depth_ref), ymomentum) 7024 7025 # check the elevation values. 7026 # ve since urs measures depth, sww meshers height, 7027 assert allclose(gauge_depth_ref, elevation) #Meters 7028 7029 fid.close() 7030 self.delete_mux(filesI) 7031 self.delete_mux(filesII) 7032 os.remove(sts_file) 7033 7034 7035 7036 7037 7038 7039 7040 7041 def test_file_boundary_stsI(self): 7042 """test_file_boundary_stsI(self): 7043 """ 7044 6799 7045 from anuga.shallow_water import Domain 6800 7046 from anuga.shallow_water import Reflective_boundary … … 6802 7048 from anuga.shallow_water import File_boundary 6803 7049 from anuga.pmesh.mesh_interface import create_mesh_from_regions 6804 from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance 7050 6805 7051 bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]] 6806 tide =0.7052 tide = 0.37 6807 7053 time_step_count = 5 6808 7054 time_step = 2 … … 6811 7057 first_tstep=ones(n,Int) 6812 7058 last_tstep=(time_step_count)*ones(n,Int) 6813 gauge_depth=20*ones(n,Float) 6814 ha=2*ones((n,time_step_count),Float) 6815 ua=10*ones((n,time_step_count),Float) 6816 va=10*ones((n,time_step_count),Float) 7059 7060 h = 20 7061 w = 2 7062 u = 10 7063 v = 10 7064 gauge_depth=h*ones(n,Float) 7065 ha=w*ones((n,time_step_count),Float) 7066 ua=u*ones((n,time_step_count),Float) 7067 va=v*ones((n,time_step_count),Float) 6817 7068 base_name, files = self.write_mux2(lat_long_points, 6818 time_step_count, time_step,6819 first_tstep, last_tstep,6820 depth=gauge_depth,6821 ha=ha,6822 ua=ua,6823 va=va)7069 time_step_count, time_step, 7070 first_tstep, last_tstep, 7071 depth=gauge_depth, 7072 ha=ha, 7073 ua=ua, 7074 va=va) 6824 7075 6825 7076 sts_file=base_name 6826 7077 urs2sts(base_name, 6827 7078 sts_file, 6828 mean_stage=tide,verbose=False) 7079 mean_stage=tide, 7080 verbose=False) 6829 7081 self.delete_mux(files) 6830 7082 … … 6840 7092 interior_regions=interior_regions,verbose=False) 6841 7093 6842 domain_fbound = pmesh_to_domain_instance(meshname, Domain)7094 domain_fbound = Domain(meshname) 6843 7095 domain_fbound.set_quantity('stage', tide) 6844 7096 Bf = File_boundary(sts_file+'.sts', domain_fbound, boundary_polygon=bounding_polygon) 6845 7097 Br = Reflective_boundary(domain_fbound) 6846 Bd=Dirichlet_boundary([2.0,220,220]) 7098 6847 7099 domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br}) 6848 7100 finaltime=time_step*(time_step_count1) 6849 7101 yieldstep=time_step 6850 7102 temp_fbound=zeros(int(finaltime/yieldstep)+1,Float) 6851 i=0 6852 for t indomain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime,6853 skip_initial_step = False):7103 7104 for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime, 7105 skip_initial_step = False)): 6854 7106 temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2] 6855 i+=1 6856 6857 domain_drchlt = pmesh_to_domain_instance(meshname, Domain)7107 7108 7109 domain_drchlt = Domain(meshname) 6858 7110 domain_drchlt.set_quantity('stage', tide) 6859 7111 Br = Reflective_boundary(domain_drchlt) 6860 Bd =Dirichlet_boundary([2.0,220,220])7112 Bd = Dirichlet_boundary([w+tide, u*(w+h+tide), v*(w+h+tide)]) 6861 7113 domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br}) 6862 7114 temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float) 6863 i=0 6864 for t indomain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime,6865 skip_initial_step = False):7115 7116 for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime, 7117 skip_initial_step = False)): 6866 7118 temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2] 6867 i+=1 6868 7119 7120 #print domain_fbound.quantities['stage'].vertex_values 7121 #print domain_drchlt.quantities['stage'].vertex_values 7122 6869 7123 assert allclose(temp_fbound,temp_drchlt) 7124 7125 assert allclose(domain_fbound.quantities['stage'].vertex_values, 7126 domain_drchlt.quantities['stage'].vertex_values) 7127 7128 assert allclose(domain_fbound.quantities['xmomentum'].vertex_values, 7129 domain_drchlt.quantities['xmomentum'].vertex_values) 7130 7131 assert allclose(domain_fbound.quantities['ymomentum'].vertex_values, 7132 domain_drchlt.quantities['ymomentum'].vertex_values) 7133 7134 6870 7135 os.remove(sts_file+'.sts') 6871 7136 os.remove(meshname) 6872 7137 6873 def test_urs2sts_file_boundary_stsII(self): 6874 """ mux2 file has points not included in boundary 7138 def test_file_boundary_stsII(self): 7139 """test_file_boundary_stsII(self): 7140 7141 mux2 file has points not included in boundary 6875 7142 mux2 gauges are not stored with the same order as they are 6876 7143 found in bounding_polygon. This does not matter as long as bounding … … 6878 7145 the correct order). 6879 7146 """ 7147 6880 7148 from anuga.shallow_water import Domain 6881 7149 from anuga.shallow_water import Reflective_boundary … … 6883 7151 from anuga.shallow_water import File_boundary 6884 7152 from anuga.pmesh.mesh_interface import create_mesh_from_regions 6885 from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance 7153 6886 7154 bounding_polygon=[[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02],[6.0,97.0]] 6887 tide =0.6888 time_step_count = 57155 tide = 0.0 # FIXME (Ole): For some reason, this one has to be zero 7156 time_step_count = 20 6889 7157 time_step = 2 6890 7158 lat_long_points=bounding_polygon[0:2] … … 6917 7185 interior_regions=None 6918 7186 boundary_tags={'ocean': [0,4], 'otherocean': [1,2,3]} 6919 # have to change boundary tags from last example because now bounding6920 # polygon starts in different place.7187 # have to change boundary tags from last example because now bounding 7188 # polygon starts in different place. 6921 7189 create_mesh_from_regions(bounding_polygon,boundary_tags=boundary_tags, 6922 7190 maximum_triangle_area=extent_res,filename=meshname, 6923 7191 interior_regions=interior_regions,verbose=False) 6924 7192 6925 domain_fbound = pmesh_to_domain_instance(meshname, Domain)7193 domain_fbound = Domain(meshname) 6926 7194 domain_fbound.set_quantity('stage', tide) 6927 7195 Bf = File_boundary(sts_file+'.sts', domain_fbound, boundary_polygon=bounding_polygon) 6928 7196 Br = Reflective_boundary(domain_fbound) 6929 Bd=Dirichlet_boundary([2.0,220,220]) 7197 6930 7198 domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br}) 6931 7199 finaltime=time_step*(time_step_count1) 6932 7200 yieldstep=time_step 6933 7201 temp_fbound=zeros(int(finaltime/yieldstep)+1,Float) 6934 i=06935 for t indomain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime,6936 skip_initial_step = False):7202 7203 for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime, 7204 skip_initial_step = False)): 6937 7205 temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2] 6938 i+=1 6939 6940 domain_drchlt = pmesh_to_domain_instance(meshname, Domain) 7206 7207 domain_drchlt = Domain(meshname) 6941 7208 domain_drchlt.set_quantity('stage', tide) 6942 7209 Br = Reflective_boundary(domain_drchlt) 6943 Bd =Dirichlet_boundary([2.0,220,220])7210 Bd = Dirichlet_boundary([2.0+tide,220+10*tide,22010*tide]) 6944 7211 domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br}) 6945 7212 temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float) 6946 i=0 6947 for t indomain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime,6948 skip_initial_step = False):7213 7214 for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime, 7215 skip_initial_step = False)): 6949 7216 temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2] 6950 i+=1 6951 6952 assert allclose(temp_fbound,temp_drchlt) 7217 7218 7219 assert allclose(temp_fbound,temp_drchlt) 7220 7221 #print domain_fbound.quantities['stage'].vertex_values 7222 #print domain_drchlt.quantities['stage'].vertex_values 7223 7224 7225 assert allclose(domain_fbound.quantities['stage'].vertex_values, 7226 domain_drchlt.quantities['stage'].vertex_values) 7227 7228 assert allclose(domain_fbound.quantities['xmomentum'].vertex_values, 7229 domain_drchlt.quantities['xmomentum'].vertex_values) 7230 7231 assert allclose(domain_fbound.quantities['ymomentum'].vertex_values, 7232 domain_drchlt.quantities['ymomentum'].vertex_values) 7233 7234 7235 7236 7237 6953 7238 os.remove(sts_file+'.sts') 6954 7239 os.remove(meshname) 6955 7240 6956 def test_urs2sts_file_boundary_stsIII_ordering(self): 6957 """Read correct points from ordering file 7241 7242 7243 def test_file_boundary_stsIII_ordering(self): 7244 """test_file_boundary_stsIII_ordering(self): 7245 Read correct points from ordering file and apply sts to boundary 6958 7246 """ 6959 7247 from anuga.shallow_water import Domain … … 6962 7250 from anuga.shallow_water import File_boundary 6963 7251 from anuga.pmesh.mesh_interface import create_mesh_from_regions 6964 from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance6965 7252 6966 7253 lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]] 6967 7254 bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]] 6968 tide =0.6969 time_step_count = 5 7255 tide = 0.0 # FIXME (Ole): For some reason, this one has to be zero 7256 time_step_count = 50 6970 7257 time_step = 2 6971 7258 n=len(lat_long_points) … … 6984 7271 va=va) 6985 7272 6986 # Write order file7273 # Write order file 6987 7274 file_handle, order_base_name = tempfile.mkstemp("") 6988 7275 os.close(file_handle) … … 6991 7278 order_file=order_base_name+'order.txt' 6992 7279 fid=open(order_file,'w') 6993 #Write Header 7280 7281 # Write Header 6994 7282 header='index, longitude, latitude\n' 6995 7283 fid.write(header) … … 7012 7300 os.remove(order_file) 7013 7301 7014 # Append the remaining part of the boundary polygon to be defined by7015 # the user7302 # Append the remaining part of the boundary polygon to be defined by 7303 # the user 7016 7304 bounding_polygon_utm=[] 7017 7305 for point in bounding_polygon: … … 7039 7327 interior_regions=None 7040 7328 boundary_tags={'ocean': [0,4], 'otherocean': [1,2,3]} 7041 #have to change boundary tags from last example because now bounding 7042 #polygon starts in different place. 7329 7330 # have to change boundary tags from last example because now bounding 7331 # polygon starts in different place. 7043 7332 create_mesh_from_regions(boundary_polygon,boundary_tags=boundary_tags, 7044 7333 maximum_triangle_area=extent_res,filename=meshname, 7045 7334 interior_regions=interior_regions,verbose=False) 7046 7335 7047 domain_fbound = pmesh_to_domain_instance(meshname, Domain)7336 domain_fbound = Domain(meshname) 7048 7337 domain_fbound.set_quantity('stage', tide) 7049 7338 Bf = File_boundary(sts_file+'.sts', domain_fbound, boundary_polygon=boundary_polygon) 7050 7339 Br = Reflective_boundary(domain_fbound) 7051 Bd=Dirichlet_boundary([2.0,220,220]) 7340 7052 7341 domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br}) 7053 7342 finaltime=time_step*(time_step_count1) 7054 7343 yieldstep=time_step 7055 7344 temp_fbound=zeros(int(finaltime/yieldstep)+1,Float) 7056 i=07057 for t indomain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime,7058 skip_initial_step = False):7345 7346 for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime, 7347 skip_initial_step = False)): 7059 7348 temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2] 7060 i+=17061 7062 domain_drchlt = pmesh_to_domain_instance(meshname, Domain)7349 7350 7351 domain_drchlt = Domain(meshname) 7063 7352 domain_drchlt.set_quantity('stage', tide) 7064 7353 Br = Reflective_boundary(domain_drchlt) 7065 Bd =Dirichlet_boundary([2.0,220,220])7354 Bd = Dirichlet_boundary([2.0+tide,220+10*tide,22010*tide]) 7066 7355 domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br}) 7067 7356 temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float) 7068 i=07069 for t indomain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime,7070 skip_initial_step = False):7357 7358 for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime, 7359 skip_initial_step = False)): 7071 7360 temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2] 7072 i+=1 7073 7361 7362 7363 #print domain_fbound.quantities['stage'].vertex_values 7364 #print domain_drchlt.quantities['stage'].vertex_values 7365 7074 7366 assert allclose(temp_fbound,temp_drchlt) 7367 7368 7369 assert allclose(domain_fbound.quantities['stage'].vertex_values, 7370 domain_drchlt.quantities['stage'].vertex_values) 7371 7372 assert allclose(domain_fbound.quantities['xmomentum'].vertex_values, 7373 domain_drchlt.quantities['xmomentum'].vertex_values) 7374 7375 assert allclose(domain_fbound.quantities['ymomentum'].vertex_values, 7376 domain_drchlt.quantities['ymomentum'].vertex_values) 7377 7378 # Use known Dirichlet condition (if sufficient timesteps have been taken) 7379 assert allclose(domain_drchlt.quantities['stage'].vertex_values[6], 2) 7380 assert allclose(domain_fbound.quantities['stage'].vertex_values[6], 2) 7381 7382 7075 7383 7076 7384 try: … … 7081 7389 7082 7390 os.remove(meshname) 7391 7392 7393 7394 def Xtest_file_boundary_stsIV_sinewave_ordering(self): 7395 """test_file_boundary_stsIV_sinewave_ordering(self): 7396 Read correct points from ordering file and apply sts to boundary 7397 This one uses a sine wave and compares to time boundary 7398 """ 7399 #FIXME (Ole): Under construction 7400 7401 from anuga.shallow_water import Domain 7402 from anuga.shallow_water import Reflective_boundary 7403 from anuga.shallow_water import Dirichlet_boundary 7404 from anuga.shallow_water import File_boundary 7405 from anuga.pmesh.mesh_interface import create_mesh_from_regions 7406 7407 lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]] 7408 bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]] 7409 tide = 0.35 7410 time_step_count = 50 7411 time_step = 2 7412 7413 n=len(lat_long_points) 7414 first_tstep=ones(n,Int) 7415 last_tstep=(time_step_count)*ones(n,Int) 7416 7417 gauge_depth=20*ones(n,Float) 7418 ha=2*ones((n,time_step_count),Float) 7419 ua=10*ones((n,time_step_count),Float) 7420 va=10*ones((n,time_step_count),Float) 7421 7422 7423 base_name, files = self.write_mux2(lat_long_points, 7424 time_step_count, time_step, 7425 first_tstep, last_tstep, 7426 depth=gauge_depth, 7427 ha=ha, 7428 ua=ua, 7429 va=va) 7430 7431 # Write order file 7432 file_handle, order_base_name = tempfile.mkstemp("") 7433 os.close(file_handle) 7434 os.remove(order_base_name) 7435 d="," 7436 order_file=order_base_name+'order.txt' 7437 fid=open(order_file,'w') 7438 7439 # Write Header 7440 header='index, longitude, latitude\n' 7441 fid.write(header) 7442 indices=[3,0,1] 7443 for i in indices: 7444 line=str(i)+d+str(lat_long_points[i][1])+d+\ 7445 str(lat_long_points[i][0])+"\n" 7446 fid.write(line) 7447 fid.close() 7448 7449 sts_file=base_name 7450 urs2sts(base_name, basename_out=sts_file, 7451 ordering_filename=order_file, 7452 mean_stage=tide, 7453 verbose=False) 7454 self.delete_mux(files) 7455 7456 7457 7458 # Now read the sts file and check that values have been stored correctly. 7459 fid = NetCDFFile(sts_file + '.sts') 7460 7461 # Check the time vector 7462 times = fid.variables['time'][:] 7463 7464 #print times 7465 7466 # Check sts quantities 7467 stage = fid.variables['stage'][:] 7468 xmomentum = fid.variables['xmomentum'][:] 7469 ymomentum = fid.variables['ymomentum'][:] 7470 elevation = fid.variables['elevation'][:] 7471 7472 #print stage 7473 #print xmomentum 7474 #print ymomentum 7475 #print elevation 7476 7477 7478 7479 # Create beginnings of boundary polygon based on sts_boundary 7480 boundary_polygon = create_sts_boundary(base_name) 7481 7482 os.remove(order_file) 7483 7484 # Append the remaining part of the boundary polygon to be defined by 7485 # the user 7486 bounding_polygon_utm=[] 7487 for point in bounding_polygon: 7488 zone,easting,northing=redfearn(point[0],point[1]) 7489 bounding_polygon_utm.append([easting,northing]) 7490 7491 boundary_polygon.append(bounding_polygon_utm[3]) 7492 boundary_polygon.append(bounding_polygon_utm[4]) 7493 7494 #print 'boundary_polygon', boundary_polygon 7495 7496 plot=False 7497 if plot: 7498 from pylab import plot,show,axis 7499 boundary_polygon=ensure_numeric(boundary_polygon) 7500 bounding_polygon_utm=ensure_numeric(bounding_polygon_utm) 7501 #plot(lat_long_points[:,0],lat_long_points[:,1],'o') 7502 plot(boundary_polygon[:,0], boundary_polygon[:,1]) 7503 plot(bounding_polygon_utm[:,0],bounding_polygon_utm[:,1]) 7504 show() 7505 7506 assert allclose(bounding_polygon_utm,boundary_polygon) 7507 7508 7509 extent_res=1000000 7510 meshname = 'urs_test_mesh' + '.tsh' 7511 interior_regions=None 7512 boundary_tags={'ocean': [0,4], 'otherocean': [1,2,3]} 7513 7514 # have to change boundary tags from last example because now bounding 7515 # polygon starts in different place. 7516 create_mesh_from_regions(boundary_polygon, 7517 boundary_tags=boundary_tags, 7518 maximum_triangle_area=extent_res, 7519 filename=meshname, 7520 interior_regions=interior_regions, 7521 verbose=False) 7522 7523 domain_fbound = Domain(meshname) 7524 domain_fbound.set_quantity('stage', tide) 7525 Bf = File_boundary(sts_file+'.sts', 7526 domain_fbound, 7527 boundary_polygon=boundary_polygon) 7528 Br = Reflective_boundary(domain_fbound) 7529 7530 domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br}) 7531 finaltime=time_step*(time_step_count1) 7532 yieldstep=time_step 7533 temp_fbound=zeros(int(finaltime/yieldstep)+1,Float) 7534 7535 for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep, 7536 finaltime=finaltime, 7537 skip_initial_step=False)): 7538 temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2] 7539 7540 7541 domain_drchlt = Domain(meshname) 7542 domain_drchlt.set_quantity('stage', tide) 7543 Br = Reflective_boundary(domain_drchlt) 7544 w = 2.0+tide 7545 h = 20+w 7546 Bd = Dirichlet_boundary([w, 10*h,10*h]) 7547 domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br}) 7548 7549 temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float) 7550 7551 for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime, 7552 skip_initial_step=False)): 7553 temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2] 7554 7555 7556 7557 print temp_fbound 7558 print temp_drchlt 7559 7560 print domain_fbound.quantities['stage'].vertex_values 7561 print domain_drchlt.quantities['stage'].vertex_values 7562 7563 assert allclose(temp_fbound, temp_drchlt) 7564 assert allclose(domain_fbound.quantities['stage'].vertex_values, 7565 domain_drchlt.quantities['stage'].vertex_values) 7566 7567 assert allclose(domain_fbound.quantities['xmomentum'].vertex_values, 7568 domain_drchlt.quantities['xmomentum'].vertex_values) 7569 7570 assert allclose(domain_fbound.quantities['ymomentum'].vertex_values, 7571 domain_drchlt.quantities['ymomentum'].vertex_values) 7572 7573 7574 try: 7575 os.remove(sts_file+'.sts') 7576 except: 7577 # Windoze can't remove this file for some reason 7578 pass 7579 7580 os.remove(meshname) 7581 7582 7583 7083 7584 7084 7585 def test_lon_lat2grid(self): … … 8944 9445 if __name__ == "__main__": 8945 9446 8946 #suite = unittest.makeSuite(Test_Data_Manager,'test')8947 #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_read_mux2_pyI')8948 9447 suite = unittest.makeSuite(Test_Data_Manager,'test') 9448 #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_sts') 9449 #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_ordering_different_sources') 8949 9450 #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo') 8950 9451 #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
Note: See TracChangeset
for help on using the changeset viewer.