Changeset 3940
- Timestamp:
- Nov 8, 2006, 4:04:55 PM (17 years ago)
- Location:
- anuga_work/production/dampier_2006
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/production/dampier_2006/build_dampier.py
r3905 r3940 32 32 from anuga.pmesh.mesh_interface import create_mesh_from_regions 33 33 from anuga.geospatial_data.geospatial_data import * 34 from anuga.abstract_2d_finite_volumes.util import Screen_Catcher34 from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files 35 35 36 36 # Application specific imports … … 42 42 #------------------------------------------------------------------------------ 43 43 44 if access(project.output_build_time_dir,F_OK) == 0: 45 mkdir (project.output_build_time_dir) 46 copy (dirname(project.__file__) +sep+ project.__name__+'.py', 47 project.output_build_time_dir + project.__name__+'.py') #copies project.py 48 copy (__file__, project.output_build_time_dir + basename(__file__)) 49 print 'files '+ project.__name__+'.py and '+ basename(__file__)+' copied to '+ project.output_build_time_dir #copies this file 50 #import sys; sys.exit() 51 52 53 #normal screen output is stored in 54 screen_output_name = project.output_build_time_dir + "screen_output.txt" 55 screen_error_name = project.output_build_time_dir + "screen_error.txt" 56 57 #used to catch screen output to file 58 sys.stdout = Screen_Catcher(screen_output_name) 59 sys.stderr = Screen_Catcher(screen_error_name) 44 copy_code_files(project.output_build_time_dir,__file__, 45 dirname(project.__file__)+sep+ project.__name__+'.py' ) 46 47 start_screen_catcher(project.output_build_time_dir) 60 48 61 49 print 'USER: ', project.user … … 125 113 G_off13 = Geospatial_data(file_name = project.offshore_dir_name13 + '.xya') 126 114 G_off14 = Geospatial_data(file_name = project.offshore_dir_name14 + '.xya') 127 115 ''' 116 117 118 print'clip nw', project.clip_poly_nw 119 print'clip e', project.clip_poly_e 120 121 print'reading combined_dir_name' 122 G_offshore_data = Geospatial_data(file_name = project.topographies_dir+'dampier_combined_elevation_final.pts') 123 124 print'reading offshore_dir_name_old' 125 G_offshore_old = Geospatial_data(file_name = project.offshore_dir_name_old + '.pts') 126 127 ''' 128 print 'G_offshore_data_old_nw', 129 G_nw_name = project.topographies_dir + 'nw_old_data' 130 G_offshore_nw = G_offshore_old.clip(project.clip_poly_nw) 131 G_offshore_nw.export_points_file(G_nw_name + '.pts') 132 G_offshore_nw.export_points_file(G_nw_name + '.xya') 133 G_nw = array(G_offshore_nw.get_data_points()) 134 print'shape of arr nw data', G_nw.shape 135 print' max and min of array 0',max(G_nw[:,0]),min(G_nw[:,0]) 136 print' max and min of array 1',max(G_nw[:,1]),min(G_nw[:,1]) 137 138 print 'G_offshore_data_old_e'#, G_offshore_old_nw.get_data_points() 139 G_e_name = project.topographies_dir+'e_old_data' 140 G_offshore_e = G_offshore_old.clip(project.clip_poly_e) 141 G_offshore_e.export_points_file(G_e_name + '.pts') 142 G_offshore_e.export_points_file(G_e_name + '.xya') 143 G_e = array(G_offshore_e.get_data_points()) 144 print'shape of arr e data', G_e.shape 145 print' max and min of array 0',max(G_e[:,0]),min(G_e[:,0]) 146 print' max and min of array 1',max(G_e[:,1]),min(G_e[:,1]) 147 ''' 148 149 print 'G_offshore_data_mid_e'#, G_offshore_old_nw.get_data_points() 150 G_mid_e_name = project.topographies_dir+'mid_e_old_data' 151 G_offshore_mid_e = G_offshore_old.clip(project.clip_poly_mid_e) 152 G_offshore_mid_e.export_points_file(G_mid_e_name + '.pts') 153 G_offshore_mid_e.export_points_file(G_mid_e_name + '.xya') 154 G_mid_e = array(G_offshore_mid_e.get_data_points()) 155 print'shape of arr e data', G_mid_e.shape 156 print' max and min of array 0',max(G_mid_e[:,0]),min(G_mid_e[:,0]) 157 print' max and min of array 1',max(G_mid_e[:,1]),min(G_mid_e[:,1]) 158 159 print 'G_offshore_data_mid_w'#, G_offshore_old_nw.get_data_points() 160 G_mid_w_name = project.topographies_dir+'mid_w_old_data' 161 G_offshore_mid_w = G_offshore_old.clip(project.clip_poly_mid_w) 162 G_offshore_mid_w.export_points_file(G_mid_w_name + '.pts') 163 G_offshore_mid_w.export_points_file(G_mid_w_name + '.xya') 164 G_mid_w = array(G_offshore_mid_w.get_data_points()) 165 print'shape of arr e data', G_mid_w.shape 166 print' max and min of array 0',max(G_mid_w[:,0]),min(G_mid_w[:,0]) 167 print' max and min of array 1',max(G_mid_w[:,1]),min(G_mid_w[:,1]) 168 169 170 print 'G_offshore_data_old_e'#, G_offshore_old_e.get_data_points() 128 171 print'add all geospatial objects' 129 G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4 + G_off5 \ 130 + G_off6 + G_off7 + G_off8 + G_off9 + G_off10 + G_off11 + G_off12 \ 131 + G_off13 + G_off14 132 133 print'clip combined geospatial object by bounding polygon' 134 G.clip(project.bounding_polygon) 172 #G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4 + G_off5 \ 173 # + G_off6 + G_off7 + G_off8 + G_off9 + G_off10 + G_off11 + G_off12 \ 174 # + G_off13 + G_off14 175 176 G = G_offshore_data + G_offshore_mid_w + G_offshore_mid_e 177 print'shape of arr G data', G.get_data_points().shape 178 179 180 #print'clip combined geospatial object by bounding polygon' 181 G_clipped = G.clip(project.bounding_polygon) 135 182 #FIXME: add a clip function to pts 183 print'shape of clipped data', G_clipped.get_data_points().shape 136 184 137 185 print'export combined DEM file' 138 186 if access(project.topographies_time_dir,F_OK) == 0: 139 187 mkdir (project.topographies_time_dir) 140 G.export_points_file(project.combined_time_dir_name + '.pts') 188 G_clipped.export_points_file(project.combined_time_dir_final_name + '.pts') 189 G_clipped.export_points_file(project.combined_time_dir_final_name + '.xya') 141 190 ''' 142 191 #------------------------------------------------------------------------- … … 156 205 if access(project.boundaries_time_dir,F_OK) == 0: 157 206 mkdir (project.boundaries_time_dir) 158 ''' 159 urs2sww(boundaries_in_dir_name,basename_out= project.boundaries_time_dir_name, 160 minlat=project.south_boundary, maxlat=project.north_boundary, 161 minlon= project.west_boundary, maxlon=project.east_boundary, 162 mint=0, maxt= 35000, 163 verbose='true') 164 165 166 urs2sww(boundaries_in_dir_name, basename_out= project.boundaries_time_dir_name, 167 # minlat=project.south, maxlat=project.north, 168 # minlon= project.west, maxlon=project.east, 169 minlat=project.south_boundary, maxlat=project.north_boundary, 170 minlon= project.east_boundary, maxlon=project.west_boundary, 171 mint=0, maxt= 35000, 172 verbose='true') 173 ''' 207 174 208 from caching import cache 175 209 cache(urs2sww, … … 193 227 ) 194 228 # dependencies = source_dir + project.boundary_basename + '.sww') 195 196 197 198 199 200 201 202 203 229 ''' 230 231 232 233 234 235 236 237 -
anuga_work/production/dampier_2006/project.py
r3937 r3940 33 33 34 34 #mesh_name = 'elevation50m' 35 boundaries_name = 'dampier 1'35 boundaries_name = 'dampier' 36 36 boundaries_source = 'mag_9_corrected' 37 37 #boundaries_source = 'test' … … 60 60 offshore_name14 = 'XYDM83' 61 61 62 offshore_old = 'elevation50m' 63 62 64 combined_name ='dampier_combined_elevation' 63 64 65 combined_final_name ='dampier_combined_elevation_final' 66 67 gauge_name = 'dampier_gauges.csv' 65 68 66 69 #Derive subdirectories and filenames … … 83 86 #print 'ctime: ', cctime 84 87 88 output_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep 85 89 output_build_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+build_time+sep 86 90 output_run_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+run_time+sep … … 94 98 95 99 96 #gauge_filename = gaugedir + 'gauge_location_broome.csv' 100 gauges_dir_name = gauges_dir + gauge_name 97 101 98 102 onshore_dir_name = topographies_dir + onshore_name … … 115 119 offshore_dir_name14 = topographies_dir + offshore_name14 116 120 121 offshore_dir_name_old = topographies_dir + offshore_old 122 123 124 117 125 #output dir 118 126 combined_dir_name = topographies_dir + combined_name 119 127 combined_time_dir_name = topographies_time_dir + combined_name 128 combined_time_dir_final_name = topographies_time_dir + combined_final_name 129 120 130 121 131 meshes_dir_name = meshes_dir + scenario_name … … 125 135 boundaries_in_dir_name = boundaries_in_dir + boundaries_name 126 136 boundaries_time_dir_name = boundaries_time_dir + boundaries_name #Used by post processing 127 137 boundaries_dir_name = boundaries_dir + boundaries_name 128 138 129 139 … … 164 174 165 175 bounding_polygon, zone =\ 166 convert_from_latlon_to_utm([p0, p1, p2, p3, p4, p5, p6, p7, p8])176 convert_from_latlon_to_utm([p0, p1, p2, p3, p4, p5, p6, p7, p8]) 167 177 #bounding_polygon, zone =\ 168 178 # convert_from_latlon_to_utm([p1, p2, p3, p4, p5, p6, p7]) … … 192 202 193 203 poly_pipeline = read_polygon(polygons_dir+'pipeline.csv') 204 205 clip_poly_e = read_polygon(polygons_dir+'gap_e.csv') 206 207 clip_poly_nw = read_polygon(polygons_dir+'gap_nw.csv') 208 209 clip_poly_mid_w = read_polygon(polygons_dir+'gap_mid_w.cvs') 210 211 clip_poly_mid_e = read_polygon(polygons_dir+'gap_mid_e.cvs') 194 212 195 213 #Interior regions -
anuga_work/production/dampier_2006/run_dampier.py
r3905 r3940 33 33 34 34 from anuga.geospatial_data.geospatial_data import * 35 from anuga.abstract_2d_finite_volumes.util import Screen_Catcher35 from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files 36 36 from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier 37 37 … … 44 44 #------------------------------------------------------------------------------ 45 45 46 46 47 # filenames 47 48 48 49 #build_time = '20061029_231935_build_tide_24' 49 build_time = '20061030_165746_build_tide_24' 50 #build_time = '20061030_165746_build_tide_24' 51 #build_time = '20061102_215532_build_plus_old_data' 52 #build_time = '20061103_055258_build' 53 build_time = '20061107_063840_build' 50 54 #build_time = '20061025_153643_build_basic' 55 bound_time = '20061102_221245_build' 51 56 52 57 boundaries_name = project.boundaries_name 53 meshes_time_dir_name = project.meshes_time_dir_name+'.msh' 58 #meshes_time_dir_name = project.meshes_time_dir_name+'.msh' 59 meshes_dir_name = project.meshes_dir_name+'.msh' 54 60 #source_dir = project.boundarydir 55 61 #boundaries_time_dir_name = project.boundaries_dir + build_time + sep + boundaries_name 56 62 boundaries_time_dir_name = project.boundaries_time_dir_name 63 boundaries_dir_name = project.boundaries_dir_name 57 64 tide = project.tide 58 65 59 # # creates copy of code in output dir if dir doesn't exist66 # creates copy of code in output dir 60 67 if myid == 0: 61 if access(project.output_run_time_dir,F_OK) == 0: 62 print 'project.output_run_time_dir',dirname(project.output_run_time_dir) 63 mkdir (project.output_run_time_dir,0777) 64 copy (dirname(project.__file__) +sep+ project.__name__+'.py', 65 project.output_run_time_dir + project.__name__+'.py') 66 copy (__file__, project.output_run_time_dir + basename(__file__)) 67 print 'project.output_run_time_dir',project.output_run_time_dir 68 copy_code_files(project.output_run_time_dir,__file__, 69 dirname(project.__file__)+sep+ project.__name__+'.py' ) 68 70 barrier() 69 #normal screen output is stored in 70 screen_output_name = project.output_run_time_dir + "screen_output_%d_%d.txt" %(myid,numprocs) 71 screen_error_name = project.output_run_time_dir + "screen_error_%d_%d.txt" %(myid,numprocs) 72 73 #used to catch screen output to file 74 sys.stdout = Screen_Catcher(screen_output_name) 75 sys.stderr = Screen_Catcher(screen_error_name) 71 start_screen_catcher(project.output_run_time_dir, myid, numprocs) 76 72 77 73 print 'USER: ', project.user 78 74 #sys.exit() 79 75 #-------------------------------------------------------------------------- 80 76 # Create the triangular mesh based on overall clipping polygon with a … … 85 81 86 82 if myid == 0: 87 if access(project.meshes_time_dir,F_OK) == 0: 88 mkdir(project.meshes_time_dir) 83 # if access(project.meshes_time_dir,F_OK) == 0: 84 # mkdir(project.meshes_time_dir) 85 if access(project.meshes_dir,F_OK) == 0: 86 mkdir(project.meshes_dir) 89 87 print 'start create mesh from regions' 90 88 interior_regions = [#[project.karratha_polygon, 25000], 91 89 # [project.cipma_polygon, 1000], 92 [project.poly_pipeline, 1000],93 [project.poly_facility, 10000]]90 [project.poly_pipeline, 5000], 91 [project.poly_facility, 500]] 94 92 # meshes_dir_name = project.meshes_dir_name + '.msh' 95 93 … … 99 97 maximum_triangle_area=100000, 100 98 interior_regions=interior_regions, 101 filename=meshes_time_dir_name, 99 # filename=meshes_time_dir_name, 100 filename=meshes_dir_name, 102 101 use_cache=True, 103 102 verbose=True) … … 110 109 #------------------------------------------------------------------------- 111 110 print 'Setup computational domain' 112 domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True) 111 #domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True) 112 domain = Domain(meshes_dir_name, use_cache=True, verbose=True) 113 113 print domain.statistics() 114 114 … … 128 128 #import sys; sys.exit() 129 129 130 if access(project.boundaries_time_dir,F_OK) == 0:131 mkdir (project.boundaries_time_dir)130 #if access(project.boundaries_time_dir,F_OK) == 0: 131 # mkdir (project.boundaries_time_dir) 132 132 # put above distribute 133 print 'boundary file is: ',boundaries_dir_name 133 134 from caching import cache 134 cache(ferret2sww, 135 (boundaries_in_dir_name, 136 boundaries_time_dir_name), 137 {'verbose': True, 138 'minlat': project.south_boundary, 139 'maxlat': project.north_boundary, 140 'minlon': project.west_boundary, 141 'maxlon': project.east_boundary, 142 # 'minlat': project.south, 143 # 'maxlat': project.north, 144 # 'minlon': project.west, 145 # 'maxlon': project.east, 146 'mint': 0, 'maxt': 35000, 147 'origin': domain.geo_reference.get_origin(), 148 'mean_stage': project.tide, 149 # 'zscale': 1, #Enhance tsunami 150 'fail_on_NaN': False}, 151 verbose = True, 152 ) 153 135 if myid == 0: 136 cache(urs2sww, 137 (boundaries_in_dir_name, 138 # boundaries_time_dir_name), 139 boundaries_dir_name), 140 {'verbose': True, 141 'minlat': project.south_boundary, 142 'maxlat': project.north_boundary, 143 'minlon': project.west_boundary, 144 'maxlon': project.east_boundary, 145 # 'minlat': project.south, 146 # 'maxlat': project.north, 147 # 'minlon': project.west, 148 # 'maxlon': project.east, 149 'mint': 0, 'maxt': 35000, 150 'origin': domain.geo_reference.get_origin(), 151 'mean_stage': project.tide, 152 # 'zscale': 1, #Enhance tsunami 153 'fail_on_NaN': False}, 154 verbose = True, 155 ) 156 barrier() 154 157 155 158 … … 165 168 166 169 domain.set_quantity('elevation', 167 filename = project.topographies_dir + build_time + sep + project.combined_ name + '.pts',170 filename = project.topographies_dir + build_time + sep + project.combined_final_name + '.pts', 168 171 use_cache = True, 169 172 verbose = True, … … 187 190 domain.set_name(project.scenario_name) 188 191 domain.set_datadir(project.output_run_time_dir) 189 domain.set_default_order(2) # associate to the spatial order of the triangle192 domain.set_default_order(2) # Apply second order scheme 190 193 domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm 194 domain.set_store_vertices_uniquely(False) 195 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) 196 domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK) 191 197 192 198 #------------------------------------------------------------------------- … … 198 204 #boundariesname = project.boundaries_dir + '20061101_003322_build'+sep+boundaries_name 199 205 #print'boundariesname',boundariesname 200 Bf = File_boundary(boundaries_time_dir_name + '.sww',206 #Bf = File_boundary(boundaries_time_dir_name + '.sww', 201 207 202 208 #Bf = File_boundary(boundariesname + '.sww', 203 domain, time_thinning=20, use_cache=True, verbose=True) 209 210 print 'domain id', id(domain) 211 Bf = File_boundary(boundaries_dir_name + '.sww', 212 domain, time_thinning=5, use_cache=True, verbose=True) 204 213 205 214 print 'finished reading boundary file' … … 214 223 print'finish set boundary' 215 224 216 217 225 #---------------------------------------------------------------------------- 218 226 # Evolve system through time … … 221 229 t0 = time.time() 222 230 223 for t in domain.evolve(yieldstep = 60, finaltime = 28800): 224 #for t in domain.evolve(yieldstep = 120, finaltime = 28800): 231 #for t in domain.evolve(yieldstep = 60, finaltime = 34000): 232 # domain.write_time() 233 # domain.write_boundary_statistics(tags = 'ocean') 234 235 for t in domain.evolve(yieldstep = 120, finaltime = 9000): 236 #for t in domain.evolve(yieldstep = 120, finaltime = 130): 225 237 domain.write_time() 226 # domain.write_boundary_statistics(tags = 'ocean') 238 domain.write_boundary_statistics(tags = 'ocean') 239 print 'time: ',time.time() 240 241 for t in domain.evolve(yieldstep = 60, finaltime = 28800 242 ,skip_initial_step = True): 243 domain.write_time() 244 domain.write_boundary_statistics(tags = 'ocean') 245 246 for t in domain.evolve(yieldstep = 120, finaltime = 34800 247 ,skip_initial_step = True): 248 domain.write_time() 249 domain.write_boundary_statistics(tags = 'ocean') 227 250 228 251 print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracChangeset
for help on using the changeset viewer.