Changeset 5652
- Timestamp:
- Aug 13, 2008, 4:40:45 PM (15 years ago)
- Location:
- anuga_work/production/geraldton
- Files:
-
- 1 added
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/production/geraldton/build_geraldton.py
r5150 r5652 27 27 from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts,start_screen_catcher, copy_code_files 28 28 from anuga.geospatial_data.geospatial_data import * 29 from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files 29 30 30 31 # Application specific imports … … 36 37 #------------------------------------------------------------------------------ 37 38 38 print 'time stamp: ',project.build_time39 40 39 copy_code_files(project.output_build_time_dir,__file__, 41 40 dirname(project.__file__)+sep+ project.__name__+'.py' ) … … 43 42 start_screen_catcher(project.output_build_time_dir) 44 43 45 print 'time stamp: ',project.build_time46 44 print 'USER: ', project.user, project.host 47 48 45 49 46 #------------------------------------------------------------------------------- … … 60 57 island_in_dir_name = project.island_in_dir_name 61 58 offshore_in_dir_name = project.offshore_in_dir_name 59 offshore_in_dir_name1 = project.offshore_in_dir_name1 60 offshore_in_dir_name2 = project.offshore_in_dir_name2 61 offshore_in_dir_name3 = project.offshore_in_dir_name3 62 62 63 63 onshore_dir_name = project.onshore_dir_name … … 65 65 island_dir_name = project.island_dir_name 66 66 offshore_dir_name = project.offshore_dir_name 67 offshore_dir_name1 = project.offshore_dir_name1 68 offshore_dir_name2 = project.offshore_dir_name2 69 offshore_dir_name3 = project.offshore_dir_name3 67 70 68 71 # creates DEM from asc data … … 72 75 73 76 #creates pts file for onshore DEM 74 print "creates pts file for onshore DEM" 75 dem2pts(onshore_dir_name, 76 use_cache=True, 77 verbose=True) 77 dem2pts(onshore_dir_name, use_cache=True, verbose=True) 78 78 79 79 #creates pts file for island DEM … … 82 82 print'create Geospatial data1 objects from topographies' 83 83 G1 = Geospatial_data(file_name = onshore_dir_name + '.pts',verbose=True) 84 print'finished reading in %s.pts' %onshore_dir_name85 84 G2 = Geospatial_data(file_name = coast_in_dir_name + '.txt',verbose=True) 86 85 G3 = Geospatial_data(file_name = island_dir_name + '.pts',verbose=True) 87 G4 = Geospatial_data(file_name = offshore_in_dir_name + '.txt',verbose=True) 86 print'create Geospatial data1 objects from bathymetry' 87 G_off = Geospatial_data(file_name = offshore_in_dir_name + '.txt',verbose=True) 88 G_off1 = Geospatial_data(file_name = offshore_in_dir_name1 + '.txt',verbose=True) 89 G_off2 = Geospatial_data(file_name = offshore_in_dir_name2 + '.txt',verbose=True) 90 G_off3 = Geospatial_data(file_name = offshore_in_dir_name3 + '.txt',verbose=True) 88 91 print'finished reading in files' 89 92 90 #G_onshore_clip = G1.clip(project.poly_all, verbose=True) 91 #print 'finished clip' 92 #G_onshore_clip_clipout=G_onshore_clip.clip_outside(project.poly_cbd, verbose=True) 93 #print 'finished clipout' 94 ## reduce resolution by 10 time (1/10) of the original file. This was determined 95 ## as the resolution for this region is going to be 20000 and the asc grid cellsize is 96 ## 10m there (100). There is 2 20000m triangles in a 200m x 200m grid sqrt of 40000 is about 200, 97 ## a 200m x 200m area would contain 400 points at a 10m grid spacing and we only need 4. So 4/400 is 1/100 or 0.01 98 ## 1/100 reduction is possible but a 1/10 is enough to help with file size and loading 99 #G_onshore_clip_clipout_reduced=G_onshore_clip_clipout.split(0.1, verbose=True) 100 # 101 #G_cbd = G1.clip(projec.poly_cbd,verbose=True) 102 ##G_cbd_reduced = G_cbd.split(0.25, verbose=True) 103 # 104 #G_bathy = G4.clip(project.poly_all) 105 #G_island = G3.clip(project.poly_all) 106 #G_coast = G2.clip(project.poly_all) 107 # 108 #print'add all geospatial objects' 109 #G = G_onshore_clip_clipout_reduced + G_cbd_reduced + G_bathy + G_island + G_coast 110 111 G = G1 + G2 + G3 + G4 112 113 G_clip = G.clip(project.poly_all, verbose=True) 93 print'add all geospatial objects' 94 G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4 114 95 115 96 print'clip combined geospatial object by bounding polygon' 116 #G_clipped = G.clip(project.bounding_polygon) 117 #FIXME: add a clip function to pts 118 #print'shape of clipped data', G_clipped.get_data_points().shape 97 G_clip = G.clip(project.poly_all, verbose=True) 119 98 120 99 print'export combined DEM file' 121 100 if access(project.topographies_dir,F_OK) == 0: 122 101 mkdir (project.topographies_dir) 123 G_clip.export_points_file(project.combined_dir_name + '.txt') 124 #G_clipped.export_points_file(project.combined_dir_name + '.xya') 125 G_small,G_other = G_clip.split(0.01, verbose=True) 126 G_small.export_points_file(project.combined_dir_name_small + '.txt') 102 G_clip.export_points_file(project.combined_dir_name + '.pts') 127 103 128 ''' 129 print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya' 130 G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya') 131 print'split' 132 G_all_1, G_all_2 = G_all.split(.10) 133 print'export 1' 134 G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya') 135 print'export 2' 136 G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya') 137 138 139 #------------------------------------------------------------------------- 140 # Convert URS to SWW file for boundary conditions 141 #------------------------------------------------------------------------- 142 print 'starting to create boundary conditions' 143 boundaries_in_dir_name = project.boundaries_in_dir_name 144 145 from anuga.shallow_water.data_manager import urs2sww 146 147 print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary 148 print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary 149 150 #import sys; sys.exit() 151 152 #if access(project.boundaries_dir,F_OK) == 0: 153 # mkdir (project.boundaries_dir) 154 155 from caching import cache 156 cache(urs2sww, 157 (boundaries_in_dir_name, 158 project.boundaries_dir_name1), 159 {'verbose': True, 160 'minlat': project.south_boundary, 161 'maxlat': project.north_boundary, 162 'minlon': project.west_boundary, 163 'maxlon': project.east_boundary, 164 # 'minlat': project.south, 165 # 'maxlat': project.north, 166 # 'minlon': project.west, 167 # 'maxlon': project.east, 168 'mint': 0, 'maxt': 40000, 169 # 'origin': domain.geo_reference.get_origin(), 170 'mean_stage': project.tide, 171 # 'zscale': 1, #Enhance tsunami 172 'fail_on_NaN': False}, 173 verbose = True, 174 ) 175 # dependencies = source_dir + project.boundary_basename + '.sww') 176 177 ''' 104 import sys 105 sys.exit() 178 106 179 107 … … 182 110 183 111 184 -
anuga_work/production/geraldton/project.py
r5381 r5652 10 10 #from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm 11 11 from anuga.utilities.system_tools import get_user_name, get_host_name 12 from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary 13 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon 12 14 13 15 # file and system info … … 15 17 16 18 home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent dir 19 muxhome = getenv('MUXHOME') 17 20 user = get_user_name() 18 21 host = get_host_name() 22 19 23 # INUNDATIONHOME is the inundation directory, not the data directory. 20 21 #needed when running using mpirun, mpirun doesn't inherit umask from .bashrc22 umask(002)23 24 24 25 #time stuff 25 26 time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir 27 gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir 26 28 build_time = time+'_build' 27 29 run_time = time+'_run' 28 29 tide = 0.75 #??? must check!!! 30 print 'gtime: ', gtime 30 31 31 32 #Making assumptions about the location of scenario data … … 34 35 scenario = 'geraldton_tsunami_scenario' 35 36 36 #Maybe will try to make project a class to allow these parameters to be passed in. 37 38 tide = 0.75 #??? must check!!! 39 37 40 alpha = 0.1 38 41 friction=0.01 39 starttime=10000 40 midtime=21600 41 #finaltime=25000 42 finaltime=10000 43 export_cellsize=50 42 starttime=0 43 finaltime=80000 44 export_cellsize=25 44 45 setup='final' 45 source=' other'46 source='polyline' 46 47 47 48 … … 62 63 yieldstep=60 63 64 64 dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+ str(user)65 dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+ 'alpha' +str(alpha)+'_'+str(user) 65 66 66 67 … … 70 71 71 72 # AHO + DPI data 72 coast_name = 'coastline' 73 offshore_name = 'geraldton_bathy' 73 coast_name = 'XYcoastline_KVP' 74 offshore_name = 'Geraldton_bathymetry' 75 offshore_name1 = 'DPI_Data' 76 offshore_name2 = 'grid250' 77 offshore_name3 = 'Top_bathymetry' 78 74 79 75 80 #final topo name … … 82 87 topographies_dir = home+state+sep+scenario+sep+'anuga'+sep+'topographies'+sep 83 88 84 # input topo file location89 ##input topo file location 85 90 onshore_in_dir_name = topographies_in_dir + onshore_name 86 91 island_in_dir_name = topographies_in_dir + island_name … … 88 93 coast_in_dir_name = topographies_in_dir + coast_name 89 94 offshore_in_dir_name = topographies_in_dir + offshore_name 95 offshore_in_dir_name1 = topographies_in_dir + offshore_name1 96 offshore_in_dir_name2 = topographies_in_dir + offshore_name2 97 offshore_in_dir_name3 = topographies_in_dir + offshore_name3 90 98 91 99 onshore_dir_name = topographies_dir + onshore_name … … 93 101 coast_dir_name = topographies_dir + coast_name 94 102 offshore_dir_name = topographies_dir + offshore_name 103 offshore_dir_name1 = topographies_dir + offshore_name1 104 offshore_dir_name2 = topographies_dir + offshore_name2 105 offshore_dir_name3 = topographies_dir + offshore_name3 95 106 96 107 #final topo files … … 104 115 tide_dir = anuga_dir+'tide_data'+sep 105 116 106 ##if source =='dampier': 107 ## boundaries_name = 'broome_3854_17042007' #Dampier gun 108 ## boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'dampier'+sep+'1_10000'+sep 109 ## 110 ##if source=='onslow': 111 ## boundaries_name = 'broome_3859_16052007' #onslow_hedland_broome gun 112 ## boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'onslow_hedland_broome'+sep+'1_10000'+sep 113 ## 114 ##if source=='exmouth': 115 ## boundaries_name = 'broome_3103_18052007' #exmouth gun 116 ## boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'exmouth'+sep+'1_10000'+sep 117 ## 118 ##if source=='other': 119 ## boundaries_name = 'other' #exmouth gun 120 ## boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'exmouth'+sep+'1_10000'+sep 121 ## 122 ## 123 ###boundaries locations 124 ##boundaries_in_dir_name = boundaries_in_dir + boundaries_name 125 boundaries_in_dir_name = anuga_dir + sep + 'boundaries' + sep + scenario_name + '49' 117 #boundaries_source = '1' 118 119 #boundaries locations 126 120 boundaries_dir = anuga_dir+'boundaries'+sep 127 121 boundaries_dir_name = boundaries_dir + scenario_name 122 boundaries_dir_mux = muxhome 128 123 129 124 #output locations 130 125 output_dir = anuga_dir+'outputs'+sep 131 output_build_time_dir = output_dir+build_time+dir_comment+sep 132 output_run_time_dir = output_dir +run_time+dir_comment+sep 133 #output_run_time_dir = output_dir+sep+run_time+sep 126 output_build_time_dir = anuga_dir+'outputs'+sep+build_time+dir_comment+sep 127 output_run_time_dir = anuga_dir+'outputs'+sep+run_time+dir_comment+sep 134 128 output_run_time_dir_name = output_run_time_dir + scenario_name #Used by post processing 129 130 vertex_filename = output_run_time_dir + 'mesh_vertex.csv' 135 131 136 132 #gauges … … 147 143 148 144 ############################### 149 # Domain definitions145 # Interior region definitions 150 146 ############################### 151 147 152 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon 148 # bounding polygon for study area 149 poly_all = read_polygon(polygons_dir+'poly_all.csv') 150 res_poly_all = 100000*res_factor 153 151 152 poly_neg20_pos20 = read_polygon(polygons_dir+'neg20_pos20.csv') 153 res_neg20_pos20 = 20000*res_factor 154 155 poly_CBD_1km= read_polygon(polygons_dir+'CBD_1km.csv') 156 res_CBD_1km = 800*res_factor 157 158 poly_cbd = read_polygon(polygons_dir+'CBD_500m.csv') 159 res_cbd = 500*res_factor 160 161 poly_wallabi = read_polygon(polygons_dir+'island_wallabi_poly.csv') 162 res_wallabi = 1000*res_factor 163 164 poly_dingiville = read_polygon(polygons_dir+'island_dingiville_poly.csv') 165 res_dingiville = 1000*res_factor 166 167 poly_pelsaert = read_polygon(polygons_dir+'island_pelsaert_poly.csv') 168 res_pelsaert = 1000*res_factor 169 170 171 #plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],figname='boundingpoly2',verbose=False) 172 173 interior_regions = [[poly_neg20_pos20,res_neg20_pos20],[poly_CBD_1km,res_CBD_1km], 174 [poly_cbd,res_cbd],[poly_wallabi,res_wallabi], 175 [poly_dingiville,res_dingiville],[poly_pelsaert, res_pelsaert]] 176 177 trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all) 178 print 'min number triangles', trigs_min 179 180 poly_mainland=read_polygon(topographies_in_dir +'initial_condition.csv') 154 181 155 182 ################################################################### … … 157 184 ################################################################### 158 185 159 # exporting asc grid160 eastingmin = 340000161 eastingmax = 350000162 northingmin = 6273400163 northingmax = 6277700164 186 165 ############################### 166 # Interior region definitions 167 ############################### 187 # Geraldton CBD extract ascii grid 188 xminCBD = 262700 189 xmaxCBD = 267600 190 yminCBD = 6811500 191 ymaxCBD = 6816400 168 192 169 #digitized polygons170 # bounding polygon for study area171 poly_all = read_polygon(polygons_dir+'poly_all_z50.csv')172 res_poly_all = 100000*res_factor173 174 poly_neg20_pos20 = read_polygon(polygons_dir+'neg20_pos20_points.csv')175 res_neg20_pos20 = 20000*res_factor176 177 poly_neg10_pos10= read_polygon(polygons_dir+'neg10_pos10_points.csv')178 res_neg10_pos10 = 2000*res_factor179 180 poly_cbd = read_polygon(polygons_dir+'cbd_500m.csv')181 res_cbd = 500*res_factor182 183 poly_wallabi = read_polygon(polygons_dir+'island_wallabi_poly_points.csv')184 res_wallabi = 20000*res_factor185 186 poly_dingiville = read_polygon(polygons_dir+'island_dingiville_poly_points.csv')187 res_dingiville = 20000*res_factor188 189 poly_pelsaert = read_polygon(polygons_dir+'island_pelsaert_poly_points.csv')190 res_pelsaert = 20000*res_factor191 192 193 #plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],figname='boundingpoly2',verbose=False)194 195 interior_regions = [[poly_neg20_pos20,res_neg20_pos20],[poly_neg10_pos10,res_neg10_pos10],196 [poly_cbd,res_cbd],[poly_wallabi,res_wallabi],197 [poly_dingiville,res_dingiville],[poly_pelsaert, res_pelsaert]]198 199 boundary_tags={'back': [2, 3, 4, 5],200 'side': [1, 6], 'ocean': [0, 7, 8, 9, 10]}201 202 trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)203 204 poly_mainland=read_polygon(polygons_dir+'initial_condition.csv')205 206 print 'min number triangles', trigs_min207 208 -
anuga_work/production/geraldton/run_geraldton.py
r5442 r5652 17 17 # Standard modules 18 18 from os import sep 19 import os 19 20 from os.path import dirname, basename 20 21 from os import mkdir, access, F_OK … … 30 31 from anuga.shallow_water import Field_boundary 31 32 from Numeric import allclose 32 from anuga.shallow_water.data_manager import export_grid 33 from anuga.shallow_water.data_manager import export_grid, create_sts_boundary 33 34 34 35 from anuga.pmesh.mesh_interface import create_mesh_from_regions 35 36 from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters 36 from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier37 #from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier 37 38 from anuga_parallel.parallel_abstraction import get_processor_name 38 39 from anuga.caching import myhash 39 40 from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage 40 41 from anuga.fit_interpolate.benchmark_least_squares import mem_usage 42 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon 43 from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter 44 from Scientific.IO.NetCDF import NetCDFFile 41 45 42 46 # Application specific imports 43 47 import project # Definition of file names and polygons 48 numprocs = 1 49 myid = 0 44 50 45 51 def run_model(**kwargs): … … 62 68 store_parameters(**kwargs) 63 69 64 barrier()70 # barrier() 65 71 66 72 start_screen_catcher(kwargs['output_dir'], myid, numprocs) … … 68 74 print "Processor Name:",get_processor_name() 69 75 70 # filenames 71 # meshes_dir_name = project.meshes_dir_name+'.msh' 72 73 # creates copy of code in output dir 74 print 'min triangles', project.trigs_min, 75 print 'Note: This is generally about 20% less than the final amount' 76 76 77 #----------------------------------------------------------------------- 78 # Domain definitions 79 #----------------------------------------------------------------------- 80 81 # Read in boundary from ordered sts file 82 urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir,project.scenario_name)) 83 84 # Reading the landward defined points, this incorporates the original clipping 85 # polygon minus the 100m contour 86 landward_bounding_polygon = read_polygon(project.polygons_dir+'landward_boundary.txt') 87 88 # Combine sts polyline with landward points 89 bounding_polygon = urs_bounding_polygon + landward_bounding_polygon 90 91 # counting segments 92 N = len(urs_bounding_polygon)-1 93 boundary_tags={'back': [N+2,N+3], 'side': [N,N+1,N+4],'ocean': range(N)} 94 95 77 96 #-------------------------------------------------------------------------- 78 97 # Create the triangular mesh based on overall clipping polygon with a … … 88 107 print 'start create mesh from regions' 89 108 90 create_mesh_from_regions( project.poly_all,91 boundary_tags= project.boundary_tags,109 create_mesh_from_regions(bounding_polygon, 110 boundary_tags=boundary_tags, 92 111 maximum_triangle_area=project.res_poly_all, 93 112 interior_regions=project.interior_regions, 94 113 filename=project.meshes_dir_name+'.msh', 95 use_cache= False,114 use_cache=True, 96 115 verbose=True) 97 barrier() 98 116 # barrier() 117 118 ## covariance_value,alpha = find_optimal_smoothing_parameter (data_file= kwargs['elevation_file'], 119 ## alpha_list=[0.001, 0.01, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5], 120 ## mesh_file = project.meshes_dir_name+'.msh') 121 ## print 'optimal alpha', covariance_value,alpha 99 122 100 123 #------------------------------------------------------------------------- … … 103 126 print 'Setup computational domain' 104 127 105 #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)106 #above don't work107 128 domain = Domain(project.meshes_dir_name+'.msh', use_cache=False, verbose=True) 108 129 print 'memory usage before del domain',mem_usage() … … 123 144 from polygon import Polygon_function 124 145 #following sets the stage/water to be offcoast only 125 IC = Polygon_function( [(project.poly_mainland, -1.0)], default = kwargs['tide'],146 IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'], 126 147 geo_reference = domain.geo_reference) 127 148 domain.set_quantity('stage', IC) 128 # domain.set_quantity('stage', kwargs['tide'])149 #domain.set_quantity('stage',kwargs['tide'] ) 129 150 domain.set_quantity('friction', kwargs['friction']) 130 151 131 print 'Start Set quantity',kwargs[' bathy_file']152 print 'Start Set quantity',kwargs['elevation_file'] 132 153 133 154 domain.set_quantity('elevation', 134 filename = kwargs[' bathy_file'],135 use_cache = True,155 filename = kwargs['elevation_file'], 156 use_cache = False, 136 157 verbose = True, 137 158 alpha = kwargs['alpha']) 138 159 print 'Finished Set quantity' 139 barrier() 140 160 #barrier() 161 162 ## #------------------------------------------------------ 163 ## # Create x,y,z file of mesh vertex!!! 164 ## #------------------------------------------------------ 165 ## coord = domain.get_vertex_coordinates() 166 ## depth = domain.get_quantity('elevation') 167 ## 168 ## # Write vertex coordinates to file 169 ## filename=project.vertex_filename 170 ## fid=open(filename,'w') 171 ## fid.write('x (m), y (m), z(m)\n') 172 ## for i in range(len(coord)): 173 ## pt=coord[i] 174 ## x=pt[0] 175 ## y=pt[1] 176 ## z=depth[i] 177 ## fid.write('%.6f,%.6f,%.6f\n' %(x, y, z)) 178 ## 141 179 142 180 #------------------------------------------------------ … … 157 195 domain.set_store_vertices_uniquely(False) 158 196 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) 197 domain.tight_slope_limiters = 1 159 198 print 'domain id', id(domain) 160 161 199 162 200 #------------------------------------------------------------------------- … … 165 203 print 'Available boundary tags', domain.get_boundary_tags() 166 204 print 'domain id', id(domain) 167 #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww' 168 print'set_boundary' 169 205 206 boundary_urs_out=project.boundaries_dir_name 207 208 print 'Available boundary tags', domain.get_boundary_tags() 209 Bf = Field_boundary(boundary_urs_out+'.sts', # Change from file_boundary 210 domain, mean_stage= project.tide, 211 time_thinning=1, 212 use_cache=True, 213 verbose = True, 214 boundary_polygon=bounding_polygon) 215 216 170 217 Br = Reflective_boundary(domain) 171 218 Bd = Dirichlet_boundary([kwargs['tide'],0,0]) 172 Bo = Dirichlet_boundary([kwargs['tide']+5.0,0,0]) 173 174 if project.source != 'other': 175 Bf = Field_boundary(kwargs['boundary_file'], 176 domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], 177 use_cache=True, verbose=True) 178 print 'finished reading boundary file' 179 domain.set_boundary({'back': Br, 180 'side': Bd, 181 'ocean': Bf}) 182 else: 183 print 'set ocean' 184 domain.set_boundary({'back': Br, 185 'side': Bd, 186 'ocean': Bd}) 187 188 # kwargs['input_start_time']=domain.starttime 189 190 219 220 fid = NetCDFFile(boundary_urs_out+'.sts', 'r') #Open existing file for read 221 sts_time=fid.variables['time'][:]+fid.starttime 222 tmin=min(sts_time) 223 tmax=max(sts_time) 224 fid.close() 225 226 print 'Boundary end time ', tmax-tmin 227 228 ## Bf = Field_boundary(kwargs['boundary_file'], 229 ## domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], 230 ## use_cache=False, verbose=True) 231 232 domain.set_boundary({'back': Br, 233 'side': Bd, 234 'ocean': Bf}) 235 236 kwargs['input_start_time']=domain.starttime 191 237 192 238 print'finish set boundary' … … 197 243 t0 = time.time() 198 244 199 for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']): 245 for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime'] 246 ,skip_initial_step = False): 200 247 domain.write_time() 201 domain.write_boundary_statistics(tags = 'ocean') 202 203 if allclose(t, 240): 204 domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo}) 205 206 if allclose(t, 1440): 207 domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd}) 208 209 # for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = kwargs['midtime'] 210 # ,skip_initial_step = True): 211 # domain.write_time() 212 # domain.write_boundary_statistics(tags = 'ocean') 213 # 214 # for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime'] 215 # ,skip_initial_step = True): 216 # domain.write_time() 217 # domain.write_boundary_statistics(tags = 'ocean') 218 248 domain.write_boundary_statistics(tags = 'ocean') 249 250 if t >= tmax-tmin: 251 print 'changed to tide boundary condition at ocean' 252 domain.set_boundary({'ocean': Bd}) 253 219 254 x, y = domain.get_maximum_inundation_location() 220 255 q = domain.get_maximum_inundation_elevation() … … 229 264 if myid==0: 230 265 store_parameters(**kwargs) 231 barrier266 # barrier 232 267 233 268 print 'memory usage before del domain1',mem_usage() 234 269 235 def export_model(**kwargs):236 #store_parameters(**kwargs)237 238 # print 'memory usage before del domain',mem_usage()239 #del domain240 print 'memory usage after del domain',mem_usage()241 242 swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']243 print'swwfile',swwfile244 245 export_grid(swwfile, extra_name_out = 'town',246 quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation247 #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation248 timestep = None,249 reduction = max,250 cellsize = kwargs['export_cellsize'],251 NODATA_value = -9999,252 easting_min = project.eastingmin,253 easting_max = project.eastingmax,254 northing_min = project.northingmin,255 northing_max = project.northingmax,256 verbose = False,257 origin = None,258 datum = 'GDA94',259 format = 'asc')260 261 inundation_damage(swwfile+'.sww', project.buildings_filename,262 project.buildings_filename_out)263 270 264 271 #------------------------------------------------------------- … … 272 279 kwargs['starttime']=project.starttime 273 280 kwargs['yieldstep']=project.yieldstep 274 kwargs['midtime']=project.midtime275 281 kwargs['finaltime']=project.finaltime 276 282 277 283 kwargs['output_dir']=project.output_run_time_dir 278 # kwargs['bathy_file']=project.combined_dir_name+'.txt' 279 kwargs['bathy_file']=project.combined_dir_name_small + '.txt' 284 kwargs['elevation_file']=project.combined_dir_name+'.pts' 280 285 kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww' 281 286 kwargs['file_name']=project.home+'detail.csv' … … 296 301 if myid==0: 297 302 export_model(**kwargs) 298 barrier303 #barrier
Note: See TracChangeset
for help on using the changeset viewer.