source: anuga_work/production/pt_hedland_2006/run_pt_hedland_urs.py @ 4943

Last change on this file since 4943 was 4943, checked in by nick, 16 years ago

update run_pt_hedland

File size: 11.8 KB
RevLine 
[4434]1"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
2
3Source data such as elevation and boundary data is assumed to be available in
4directories specified by project_urs.py
5The output sww file is stored in project_urs.output_run_time_dir
6
7The scenario is defined by a triangular mesh created from project_urs.polygon,
8the elevation data and a simulated tsunami generated with URS code.
9
10Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
11"""
12
13#------------------------------------------------------------------------------
14# Import necessary modules
15#------------------------------------------------------------------------------
16
17# Standard modules
18from os import sep,umask
19from os.path import dirname, basename
20from os import mkdir, access, F_OK
21from shutil import copy
22import time
23import sys
24
25# Related major packages
26from anuga.shallow_water import Domain
27from anuga.shallow_water import Dirichlet_boundary
28from anuga.shallow_water import File_boundary
29from anuga.shallow_water import Reflective_boundary
30from anuga.shallow_water import Field_boundary
31from Numeric import allclose
[4542]32from anuga.shallow_water.data_manager import export_grid
[4434]33
34from anuga.pmesh.mesh_interface import create_mesh_from_regions
[4542]35from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
[4434]36from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
37from anuga_parallel.parallel_abstraction import get_processor_name
38from anuga.caching import myhash
[4542]39from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
[4587]40from anuga.fit_interpolate.benchmark_least_squares import mem_usage
[4434]41# Application specific imports
42import project_urs                 # Definition of file names and polygons
43
44def run_model(**kwargs):
45   
[4544]46    print 'tide',kwargs['tide']
[4942]47#    kwargs['est_num_trigs']=project_urs.trigs_min
48#    kwargs['num_cpu']=numprocs
49#    kwargs['host']=project_urs.host
50#    kwargs['res_factor']=project_urs.res_factor
51#    kwargs['starttime']=project_urs.starttime
52#    kwargs['yieldstep']=project_urs.yieldstep
53#    kwargs['finaltime']=project_urs.finaltime
54#   
55#    kwargs['output_dir']=project_urs.output_run_time_dir
56#    kwargs['bathy_file']=project_urs.combined_dir_name + '.pts'
57#    kwargs['boundary_file']=project_urs.boundaries_in_dir_name + '.sww'
[4434]58
[4544]59    #------------------------------------------------------------------------------
60    # Copy scripts to time stamped output directory and capture screen
61    # output to file
62    #------------------------------------------------------------------------------
63
[4434]64    print 'output_dir',kwargs['output_dir']
65    if myid == 0:
66        copy_code_files(kwargs['output_dir'],__file__, 
67                 dirname(project_urs.__file__)+sep+ project_urs.__name__+'.py' )
68
69        store_parameters(**kwargs)
70
71    barrier()
[4587]72   
73    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
74   
[4434]75
76    print "Processor Name:",get_processor_name()
77
78    # creates copy of code in output dir
79    print 'min triangles', project_urs.trigs_min,
80    print 'Note: This is generally about 20% less than the final amount'
81
82    #--------------------------------------------------------------------------
83    # Create the triangular mesh based on overall clipping polygon with a
84    # tagged
85    # boundary and interior regions defined in project_urs.py along with
86    # resolutions (maximal area of per triangle) for each polygon
87    #--------------------------------------------------------------------------
88
89    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
90    # causes problems with the ability to cache set quantity which takes alot of times
91    if myid == 0:
92   
93        print 'start create mesh from regions'
94
95        create_mesh_from_regions(project_urs.poly_all,
96                             boundary_tags={'back': [3, 4, 5, 6], 'side': [2, 7],
97                                            'ocean': [0, 1]},
98                             maximum_triangle_area=project_urs.res_poly_all,
99                             interior_regions=project_urs.interior_regions,
100                             filename=project_urs.meshes_dir_name+'.msh',
101                             use_cache=False,
102                             verbose=True)
103    barrier()
104
105    #-------------------------------------------------------------------------
106    # Setup computational domain
107    #-------------------------------------------------------------------------
108    print 'Setup computational domain'
109
110    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
111    #above don't work
112    domain = Domain(project_urs.meshes_dir_name+'.msh', use_cache=False, verbose=True)
113       
114    print domain.statistics()
115    print 'triangles',len(domain)
116   
117    kwargs['act_num_trigs']=len(domain)
118
119    #-------------------------------------------------------------------------
120    # Setup initial conditions
121    #-------------------------------------------------------------------------
122    if myid == 0:
123
124        print 'Setup initial conditions'
125
126        from polygon import Polygon_function
127        #following sets the stage/water to be offcoast only
[4587]128        #IC = Polygon_function( [(project_urs.poly_mainland, -50.0)], default = kwargs['tide'],
129        #                         geo_reference = domain.geo_reference)
130        #domain.set_quantity('stage', IC)
131        domain.set_quantity('stage', kwargs['tide'])
[4544]132        domain.set_quantity('friction', kwargs['friction']) 
[4434]133       
[4587]134        print 'Start Set quantity', kwargs['bathy_file'],kwargs['alpha']
135       
[4434]136        domain.set_quantity('elevation', 
137                            filename = kwargs['bathy_file'],
138                            use_cache = True,
139                            verbose = True,
[4544]140                            alpha = kwargs['alpha'])
[4434]141        print 'Finished Set quantity'
142    barrier()
143
144    #------------------------------------------------------
145    # Distribute domain to implement parallelism !!!
146    #------------------------------------------------------
147
148    if numprocs > 1:
149        domain=distribute(domain)
150
151    #------------------------------------------------------
152    # Set domain parameters
153    #------------------------------------------------------
154    print 'domain id', id(domain)
[4544]155    domain.set_name(kwargs['aa_scenario_name'])
[4434]156    domain.set_datadir(kwargs['output_dir'])
157    domain.set_default_order(2) # Apply second order scheme
158    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
159    domain.set_store_vertices_uniquely(False)
160    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
[4943]161    domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
[4544]162    domain.beta_h = 0 #sets the surface of the triangle to follow the bathy
163    #domain.H0=0.01 #controls the flux limiter (limiter2007)
[4943]164    domain.tight_slope_limiters = 1 #minimises creep
[4434]165
166    #-------------------------------------------------------------------------
167    # Setup boundary conditions
168    #-------------------------------------------------------------------------
169    print 'Available boundary tags', domain.get_boundary_tags()
170    print 'domain id', id(domain)
171    #print 'Reading Boundary file',project_urs.boundaries_dir_namea + '.sww'
172
173    Bf = Field_boundary(kwargs['boundary_file'],
[4544]174                    domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], 
[4434]175                    use_cache=True, verbose=True)
176                   
177    kwargs['input_start_time']=domain.starttime
178
179    print 'finished reading boundary file'
180
181    Br = Reflective_boundary(domain)
[4544]182    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
[4434]183
184    print'set_boundary'
185
186    domain.set_boundary({'back': Br,
187                         'side': Bd,
188                         'ocean': Bf}) 
189    print'finish set boundary'
190
191    #----------------------------------------------------------------------------
192    # Evolve system through time
193    #----------------------------------------------------------------------------
194
195    t0 = time.time()
196
197    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['starttime']): 
198        domain.write_time()
199        domain.write_boundary_statistics(tags = 'ocean')     
200
201    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = 21600
202                       ,skip_initial_step = True): 
203        domain.write_time()
204        domain.write_boundary_statistics(tags = 'ocean')     
205   
206    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']
207                       ,skip_initial_step = True): 
208        domain.write_time()
209        domain.write_boundary_statistics(tags = 'ocean')   
210
211    x, y = domain.get_maximum_inundation_location()
212    q = domain.get_maximum_inundation_elevation()
213
214    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
215
216    print 'That took %.2f seconds' %(time.time()-t0)
217   
218    kwargs['completed']=str(time.time()-t0)
219   
[4467]220    if myid == 0:
221        store_parameters(**kwargs)
[4587]222       
223        print 'memory usage before del domain',mem_usage()
224        del domain
225        print 'memory usage after del domain',mem_usage()
226   
227        swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
228        print'swwfile',swwfile
[4542]229
[4587]230        export_grid(swwfile, extra_name_out = 'town',
231                    quantities = ['speed','depth','stage','elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
232                    timestep = None,
233                    reduction = max,
234                    cellsize = 25,
235                    NODATA_value = -9999,
236                    easting_min = project_urs.eastingmin,
237                    easting_max = project_urs.eastingmax,
238                    northing_min = project_urs.northingmin,
239                    northing_max = project_urs.northingmax,
240                    verbose = True,
241                    origin = None,
242                    datum = 'WGS84',
243                    format = 'asc')
244       
245        buildings_filename = project_urs.buildings_filename
246        buildings_filename_out = project_urs.buildings_filename_out
247                   
248        inundation_damage(swwfile+'.sww', buildings_filename, buildings_filename_out)
249        print '\n Augmented building file written to %s \n' %buildings_filename_out
[4542]250   
[4467]251    barrier()
[4434]252#-------------------------------------------------------------
253if __name__ == "__main__":
254
[4942]255#    run_model(file_name=project_urs.home+'detail.csv', aa_scenario_name=project_urs.scenario_name,
256#              ab_time=project_urs.time, res_factor= project_urs.res_factor, tide=project_urs.tide, user=project_urs.user,
257#              alpha = project_urs.alpha, friction=project_urs.friction,
258#              time_thinning = project_urs.time_thinning,
259#              dir_comment=project_urs.dir_comment)
260    kwargs={}
261    kwargs['est_num_trigs']=project_urs.trigs_min
262    kwargs['num_cpu']=numprocs
263    kwargs['host']=project_urs.host
264    kwargs['res_factor']=project_urs.res_factor
265    kwargs['starttime']=project_urs.starttime
266    kwargs['yieldstep']=project_urs.yieldstep
267    kwargs['finaltime']=project_urs.finaltime
268   
269    kwargs['output_dir']=project_urs.output_run_time_dir
270#    kwargs['bathy_file']=project_urs.combined_dir_name+'.txt'
271    kwargs['bathy_file']=project_urs.combined_dir_name + '.pts'
272    kwargs['boundary_file']=project_urs.boundaries_in_dir_name + '.sww'
273    kwargs['file_name']=project_urs.home+'detail.csv'
274    kwargs['aa_scenario_name']=project_urs.scenario_name
275    kwargs['ab_time']=project_urs.time
276    kwargs['res_factor']= project_urs.res_factor
277    kwargs['tide']=project_urs.tide
278    kwargs['user']=project_urs.user
279    kwargs['alpha'] = project_urs.alpha
280    kwargs['friction']=project_urs.friction
281    kwargs['time_thinning'] = project_urs.time_thinning
282    kwargs['dir_comment']=project_urs.dir_comment
[4434]283
[4942]284    run_model(**kwargs)
285     
286    if myid==0:
287        export_model(**kwargs)
288    barrier
[4434]289
[4942]290
Note: See TracBrowser for help on using the repository browser.