Changeset 5652 for anuga_work/production/geraldton/run_geraldton.py
- Timestamp:
- Aug 13, 2008, 4:40:45 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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.