Changeset 6178 for anuga_work/production/patong/run_patong.py
- Timestamp:
- Jan 15, 2009, 5:45:46 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/production/patong/run_patong.py
r6175 r6178 23 23 from os import sep 24 24 import os 25 from os.path import dirname, basename 26 from os import mkdir, access, F_OK 27 from shutil import copy 25 from os.path import dirname 28 26 import time 29 import sys30 27 31 28 # Related major packages … … 35 32 from anuga.shallow_water import Reflective_boundary 36 33 from anuga.shallow_water import Field_boundary 37 from Numeric import allclose38 34 from anuga.shallow_water.data_manager import export_grid, create_sts_boundary, csv2building_polygons 39 from anuga. pmesh.mesh_interface import create_mesh_from_regions35 from anuga.shallow_water import create_domain_from_regions 40 36 from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters 41 37 from anuga.caching import myhash … … 51 47 myid = 0 52 48 53 def run_model(**kwargs):54 55 #------------------------------------------------------------------------------56 # Copy scripts to time stamped output directory and capture screen57 # output to file58 #------------------------------------------------------------------------------59 60 #copy script must be before screen_catcher61 62 print 'output_dir',kwargs['output_dir']63 64 copy_code_files(kwargs['output_dir'],__file__,65 dirname(project.__file__)+sep+ project.__name__+'.py' )66 67 store_parameters(**kwargs)68 69 start_screen_catcher(kwargs['output_dir'], myid, numprocs)70 49 71 50 72 #----------------------------------------------------------------------- 73 # Domain definitions 74 #----------------------------------------------------------------------- 51 #------------------------------------------------------------------------------ 52 # Copy scripts to time stamped output directory and capture screen 53 # output to file 54 #------------------------------------------------------------------------------ 75 55 76 # Read in boundary from ordered sts file 77 urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir,project.scenario_name)) 56 #copy script must be before screen_catcher 78 57 79 # Reading the landward defined points, this incorporates the original clipping 80 # polygon minus the 100m contour 81 landward_bounding_polygon = read_polygon(project.landward_dir) 58 output_dir = project.output_run_time_dir 59 print 'output_dir', output_dir 82 60 83 # Combine sts polyline with landward points 84 bounding_polygon = urs_bounding_polygon + landward_bounding_polygon 85 86 # counting segments 87 N = len(urs_bounding_polygon)-1 61 #copy_code_files(output_dir,__file__, 62 # dirname(project.__file__)+sep+ project.__name__+'.py' ) 88 63 89 # boundary tags refer to project.landward 4 points equals 5 segments start at N 90 boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)} 91 92 #-------------------------------------------------------------------------- 93 # Create the triangular mesh based on overall clipping polygon with a tagged 94 # boundary and interior regions defined in project.py along with 95 # resolutions (maximal area of per triangle) for each polygon 96 #-------------------------------------------------------------------------- 97 98 # IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it 99 # causes problems with the ability to cache set quantity which takes alot of times 100 101 print 'start create mesh from regions' 102 103 create_mesh_from_regions(bounding_polygon, 104 boundary_tags=boundary_tags, 105 maximum_triangle_area=project.res_poly_all, 106 interior_regions=project.interior_regions, 107 filename=project.meshes_dir_name, 108 use_cache=True, 109 verbose=True) 110 111 #------------------------------------------------------------------------- 112 # Setup computational domain 113 #------------------------------------------------------------------------- 114 print 'Setup computational domain' 115 116 domain = Domain(project.meshes_dir_name, use_cache=False, verbose=True) 117 print 'memory usage before del domain',mem_usage() 118 119 print domain.statistics() 120 print 'triangles',len(domain) 121 122 kwargs['act_num_trigs']=len(domain) 64 #start_screen_catcher(output_dir, myid, numprocs) 123 65 124 66 125 #------------------------------------------------------------------------- 126 # Setup initial conditions 127 #------------------------------------------------------------------------- 128 print 'Setup initial conditions' 67 #----------------------------------------------------------------------- 68 # Domain definitions 69 #----------------------------------------------------------------------- 129 70 130 # sets the initial stage in the offcoast region only 131 IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'], 132 geo_reference = domain.geo_reference) 133 domain.set_quantity('stage', IC) 134 #domain.set_quantity('stage',kwargs['tide'] ) 135 domain.set_quantity('friction', kwargs['friction']) 136 137 print 'Start Set quantity',kwargs['elevation_file'] 71 # Read in boundary from ordered sts file 72 urs_boundary_name = os.path.join(project.boundaries_dir, 73 project.scenario_name) 74 urs_bounding_polygon = create_sts_boundary(urs_boundary_name) 138 75 139 domain.set_quantity('elevation', 140 filename = kwargs['elevation_file'], 141 use_cache = False, 142 verbose = True, 143 alpha = kwargs['alpha']) 76 # Reading the landward defined points, this incorporates the original clipping 77 # polygon minus the 100m contour 78 landward_bounding_polygon = read_polygon(project.landward_dir) 144 79 145 # Add buildings from file 146 building_polygons, building_heights = csv2building_polygons(project.building_polygon_file) 80 # Combine sts polyline with landward points 81 bounding_polygon = urs_bounding_polygon + landward_bounding_polygon 147 82 148 L = [] 149 for key in building_polygons: 150 poly = building_polygons[key] 151 elev = building_heights[key] 152 L.append((poly, elev)) 153 domain.add_quantity('elevation', Polygon_function(L, default=0.0)) 154 155 print 'Finished Set quantity' 83 # Counting segments 84 N = len(urs_bounding_polygon)-1 156 85 157 ## #------------------------------------------------------ 158 ## # Distribute domain to implement parallelism !!! 159 ## #------------------------------------------------------ 160 ## 161 ## if numprocs > 1: 162 ## domain=distribute(domain) 86 # boundary tags refer to project.landward 4 points equals 5 segments start at N 87 boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)} 163 88 164 #------------------------------------------------------ 165 # Set domain parameters 166 #------------------------------------------------------ 167 print 'domain id', id(domain) 168 domain.set_name(kwargs['scenario_name']) 169 domain.set_datadir(kwargs['output_dir']) 170 domain.set_default_order(2) # Apply second order scheme 171 domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm 172 domain.set_store_vertices_uniquely(False) 173 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) 174 domain.tight_slope_limiters = 1 175 print 'domain id', id(domain) 89 #-------------------------------------------------------------------------- 90 # Create the computational domain based on overall clipping polygon with 91 # a tagged boundary and interior regions defined in project.py along with 92 # resolutions (maximal area of per triangle) for each polygon 93 #-------------------------------------------------------------------------- 176 94 177 #------------------------------------------------------------------------- 178 # Setup boundary conditions179 #-------------------------------------------------------------------------180 print 'Available boundary tags', domain.get_boundary_tags()181 print 'domain id', id(domain)182 183 boundary_urs_out=project.boundaries_dir + sep + project.scenario_name95 domain = create_domain_from_regions(bounding_polygon, 96 boundary_tags=boundary_tags, 97 maximum_triangle_area=project.res_poly_all, 98 interior_regions=project.interior_regions, 99 mesh_filename=project.meshes_dir_name, 100 use_cache=True, 101 verbose=True) 184 102 185 Br = Reflective_boundary(domain) 186 Bd = Dirichlet_boundary([kwargs['tide'],0,0]) 187 188 print 'Available boundary tags', domain.get_boundary_tags() 189 Bf = Field_boundary(boundary_urs_out+'.sts', # Change from file_boundary 190 domain, mean_stage= project.tide, 191 time_thinning=1, 192 time_limit=project.finaltime, 193 default_boundary=Bd, 194 use_cache=True, 195 verbose=True, 196 boundary_polygon=bounding_polygon) 103 print 'memory usage after creation of domain', mem_usage() 104 print domain.statistics() 197 105 198 domain.set_boundary({'back': Br,199 'side': Bd,200 'ocean': Bf})201 106 202 kwargs['input_start_time']=domain.starttime 107 #------------------------------------------------------------------------- 108 # Setup initial conditions 109 #------------------------------------------------------------------------- 110 print 'Setup initial conditions' 203 111 204 print'finish set boundary' 112 # sets the initial stage in the offcoast region only 113 IC = Polygon_function([(project.poly_mainland, 0)], 114 default=project.tide, 115 geo_reference=domain.geo_reference) 116 domain.set_quantity('stage', IC) 117 domain.set_quantity('friction', project.friction) 118 domain.set_quantity('elevation', 119 filename=project.combined_dir_name+'.pts', 120 use_cache=True, 121 verbose=True, 122 alpha=project.alpha) 205 123 206 #---------------------------------------------------------------------------- 207 # Evolve system through time 208 #-------------------------------------------------------------------- 209 t0 = time.time() 124 # Add buildings from file 125 print 'Creating building polygons' 126 building_polygons, building_heights = csv2building_polygons(project.building_polygon_file) 210 127 211 for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime'] 212 ,skip_initial_step = False): 213 domain.write_time() 214 domain.write_boundary_statistics(tags = 'ocean') 128 L = [] 129 for key in building_polygons: 130 poly = building_polygons[key] 131 elev = building_heights[key] 132 L.append((poly, elev)) 215 133 216 # these outputs should be checked with the resultant inundation map 217 x, y = domain.get_maximum_inundation_location() 218 q = domain.get_maximum_inundation_elevation() 219 print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q) 134 buildings = Polygon_function(L, default=0.0) 135 domain.add_quantity('elevation', buildings) 220 136 221 print 'Simulation took %.2f seconds' %(time.time()-t0)222 137 223 #kwargs 'completed' must be added to write the final parameters to file 224 kwargs['completed']=str(time.time()-t0) 225 226 store_parameters(**kwargs) 227 228 229 230 #------------------------------------------------------------- 231 if __name__ == "__main__": 232 233 kwargs={} 234 kwargs['file_name']=project.dir_comment 235 kwargs['finaltime']=project.finaltime 236 kwargs['output_dir']=project.output_run_time_dir 237 kwargs['elevation_file']=project.combined_dir_name+'.pts' 238 kwargs['scenario_name']=project.scenario_name 239 kwargs['tide']=project.tide 240 kwargs['alpha'] = project.alpha 241 kwargs['friction']=project.friction 242 #kwargs['num_cpu']=numprocs 243 #kwargs['host']=project.host 244 #kwargs['starttime']=project.starttime 245 #kwargs['yieldstep']=project.yieldstep 246 #kwargs['user']=project.user 247 #kwargs['time_thinning'] = project.time_thinning 248 249 run_model(**kwargs) 138 #------------------------------------------------------ 139 # Distribute domain to implement parallelism !!! 140 #------------------------------------------------------ 141 142 # if numprocs > 1: 143 # domain=distribute(domain) 144 145 #------------------------------------------------------ 146 # Set domain parameters 147 #------------------------------------------------------ 148 domain.set_name(project.scenario_name) 149 domain.set_datadir(output_dir) 150 domain.set_default_order(2) # Apply second order scheme 151 domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm 152 153 #------------------------------------------------------------------------- 154 # Setup boundary conditions 155 #------------------------------------------------------------------------- 156 print 'Set boundary conditions' 157 158 Br = Reflective_boundary(domain) 159 Bd = Dirichlet_boundary([project.tide,0,0]) 160 Bf = Field_boundary(urs_boundary_name+'.sts', 161 domain, 162 mean_stage= project.tide, 163 time_thinning=project.time_thinning, 164 default_boundary=Bd, 165 use_cache=True, 166 verbose=True, 167 boundary_polygon=bounding_polygon) 168 169 domain.set_boundary({'back': Br, 170 'side': Bd, 171 'ocean': Bf}) 172 173 174 #---------------------------------------------------------------------------- 175 # Evolve system through time 176 #-------------------------------------------------------------------- 177 t0 = time.time() 178 179 for t in domain.evolve(yieldstep=project.yieldstep, 180 finaltime=project.finaltime) 181 print domain.timestepping_statistics() 182 print domain.boundary_statistics(tags='ocean') 183 184 print 'Simulation took %.2f seconds' %(time.time()-t0) 250 185 251 186
Note: See TracChangeset
for help on using the changeset viewer.