source: anuga_work/production/broome/run_broome_urs.py @ 5001

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

update broome with new structure

File size: 11.5 KB
Line 
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
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
32from anuga.shallow_water.data_manager import export_grid
33
34from anuga.pmesh.mesh_interface import create_mesh_from_regions
35from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
36from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
37from anuga_parallel.parallel_abstraction import get_processor_name
38from anuga.caching import myhash
39from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
40from anuga.fit_interpolate.benchmark_least_squares import mem_usage
41
42# Application specific imports
43import project_urs                 # Definition of file names and polygons
44
45def run_model(**kwargs):
46   
47
48    #------------------------------------------------------------------------------
49    # Copy scripts to time stamped output directory and capture screen
50    # output to file
51    #------------------------------------------------------------------------------
52    print "Processor Name:",get_processor_name()
53
54    #copy script must be before screen_catcher
55    #print kwargs
56    '''
57    print 'tide',kwargs['tide']
58    kwargs['est_num_trigs']=project_urs.trigs_min
59    kwargs['num_cpu']=numprocs
60    kwargs['host']=project_urs.host
61    kwargs['res_factor']=project_urs.res_factor
62    kwargs['starttime']=project_urs.starttime
63    kwargs['yieldstep']=project_urs.yieldstep
64    kwargs['finaltime']=project_urs.finaltime
65   
66    kwargs['output_dir']=project_urs.output_run_time_dir
67    kwargs['bathy_file']=project_urs.combined_dir_name+'.txt'
68#    kwargs['bathy_file']=project_urs.combined_small_dir_name + '.pts'
69    kwargs['boundary_file']=project_urs.boundaries_in_dir_name + '.sww'
70    '''
71    print 'output_dir',kwargs['output_dir']
72    if myid == 0:
73        copy_code_files(kwargs['output_dir'],__file__, 
74                 dirname(project_urs.__file__)+sep+ project_urs.__name__+'.py' )
75
76        store_parameters(**kwargs)
77
78    barrier()
79
80    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
81
82    print "Processor Name:",get_processor_name()
83
84    # filenames
85#    meshes_dir_name = project_urs.meshes_dir_name+'.msh'
86
87    # creates copy of code in output dir
88    print 'min triangles', project_urs.trigs_min,
89    print 'Note: This is generally about 20% less than the final amount'
90
91    #--------------------------------------------------------------------------
92    # Create the triangular mesh based on overall clipping polygon with a
93    # tagged
94    # boundary and interior regions defined in project_urs.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    if myid == 0:
101   
102        print 'start create mesh from regions'
103
104        create_mesh_from_regions(project_urs.poly_all,
105                         boundary_tags=project_urs.boundary_tags,
106                             maximum_triangle_area=project_urs.res_poly_all,
107                             interior_regions=project_urs.interior_regions,
108                             filename=project_urs.meshes_dir_name+'.msh',
109                             use_cache=False,
110                             verbose=True)
111    barrier()
112   
113
114    #-------------------------------------------------------------------------
115    # Setup computational domain
116    #-------------------------------------------------------------------------
117    print 'Setup computational domain'
118
119    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
120    #above don't work
121    domain = Domain(project_urs.meshes_dir_name+'.msh', use_cache=False, verbose=True)
122    print 'memory usage before del domain',mem_usage()
123       
124    print domain.statistics()
125    print 'triangles',len(domain)
126   
127    kwargs['act_num_trigs']=len(domain)
128
129
130    #-------------------------------------------------------------------------
131    # Setup initial conditions
132    #-------------------------------------------------------------------------
133    if myid == 0:
134
135        print 'Setup initial conditions'
136
137        from polygon import Polygon_function
138        #following sets the stage/water to be offcoast only
139#        IC = Polygon_function( [(project_urs.poly_mainland, -1.0)], default = kwargs['tide'],
140#                                 geo_reference = domain.geo_reference)
141#        domain.set_quantity('stage', IC)
142        domain.set_quantity('stage', kwargs['tide'])
143        domain.set_quantity('friction', kwargs['friction']) 
144       
145        print 'Start Set quantity',kwargs['bathy_file']
146
147        domain.set_quantity('elevation', 
148                            filename = kwargs['bathy_file'],
149                            use_cache = True,
150                            verbose = True,
151                            alpha = kwargs['alpha'])
152        print 'Finished Set quantity'
153    barrier()
154
155
156    #------------------------------------------------------
157    # Distribute domain to implement parallelism !!!
158    #------------------------------------------------------
159
160    if numprocs > 1:
161        domain=distribute(domain)
162
163    #------------------------------------------------------
164    # Set domain parameters
165    #------------------------------------------------------
166    print 'domain id', id(domain)
167    domain.set_name(kwargs['aa_scenario_name'])
168    domain.set_datadir(kwargs['output_dir'])
169    domain.set_default_order(2) # Apply second order scheme
170    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
171    domain.set_store_vertices_uniquely(False)
172    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
173    domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
174    print 'domain id', id(domain)
175    domain.beta_h = 0
176    #domain.tight_slope_limiters = 1
177
178    #-------------------------------------------------------------------------
179    # Setup boundary conditions
180    #-------------------------------------------------------------------------
181    print 'Available boundary tags', domain.get_boundary_tags()
182    print 'domain id', id(domain)
183    #print 'Reading Boundary file',project_urs.boundaries_dir_namea + '.sww'
184
185    Bf = Field_boundary(kwargs['boundary_file'],
186                    domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], 
187                    use_cache=True, verbose=True)
188                   
189    kwargs['input_start_time']=domain.starttime
190
191    print 'finished reading boundary file'
192
193    Br = Reflective_boundary(domain)
194    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
195
196    print'set_boundary'
197
198    domain.set_boundary({'back': Br,
199                         'side': Bd,
200                         'ocean': Bf}) 
201    print'finish set boundary'
202
203    #----------------------------------------------------------------------------
204    # Evolve system through time
205    #--------------------------------------------------------------------
206    t0 = time.time()
207
208    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['starttime']): 
209        domain.write_time()
210        domain.write_boundary_statistics(tags = 'ocean')     
211
212    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = 21600
213                       ,skip_initial_step = True): 
214        domain.write_time()
215        domain.write_boundary_statistics(tags = 'ocean')     
216   
217    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']
218                       ,skip_initial_step = True): 
219        domain.write_time()
220        domain.write_boundary_statistics(tags = 'ocean')   
221
222    x, y = domain.get_maximum_inundation_location()
223    q = domain.get_maximum_inundation_elevation()
224
225    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
226
227    print 'That took %.2f seconds' %(time.time()-t0)
228
229    #kwargs 'completed' must be added to write the final parameters to file
230    kwargs['completed']=str(time.time()-t0)
231   
232    if myid==0:
233        store_parameters(**kwargs)
234    barrier
235   
236    print 'memory usage before del domain1',mem_usage()
237   
238def export_model(**kwargs):   
239    #store_parameters(**kwargs)
240   
241    print 'memory usage before del domain',mem_usage()
242    #del domain
243    print 'memory usage after del domain',mem_usage()
244   
245    swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
246    print'swwfile',swwfile
247   
248    export_grid(swwfile, extra_name_out = 'town',
249            quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
250            #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
251            timestep = None,
252            reduction = max,
253            cellsize = 25,
254            NODATA_value = -9999,
255            easting_min = project_urs.eastingmin,
256            easting_max = project_urs.eastingmax,
257            northing_min = project_urs.northingmin,
258            northing_max = project_urs.northingmax,
259            verbose = False,
260            origin = None,
261            datum = 'GDA94',
262            format = 'asc')
263           
264    inundation_damage(swwfile+'.sww', project_urs.buildings_filename, 
265                      project_urs.buildings_filename_out)
266   
267#-------------------------------------------------------------
268if __name__ == "__main__":
269   
270    kwargs={}
271    kwargs['est_num_trigs']=project_urs.trigs_min
272    kwargs['num_cpu']=numprocs
273    kwargs['host']=project_urs.host
274    kwargs['res_factor']=project_urs.res_factor
275    kwargs['starttime']=project_urs.starttime
276    kwargs['yieldstep']=project_urs.yieldstep
277    kwargs['finaltime']=project_urs.finaltime
278   
279    kwargs['output_dir']=project_urs.output_run_time_dir
280    kwargs['bathy_file']=project_urs.combined_dir_name+'.txt'
281#    kwargs['bathy_file']=project_urs.combined_small_dir_name + '.pts'
282    kwargs['boundary_file']=project_urs.boundaries_in_dir_name + '.sww'
283    kwargs['file_name']=project_urs.home+'detail.csv'
284    kwargs['aa_scenario_name']=project_urs.scenario_name
285    kwargs['ab_time']=project_urs.time
286    kwargs['res_factor']= project_urs.res_factor
287    kwargs['tide']=project_urs.tide
288    kwargs['user']=project_urs.user
289    kwargs['alpha'] = project_urs.alpha
290    kwargs['friction']=project_urs.friction
291    kwargs['time_thinning'] = project_urs.time_thinning
292    kwargs['dir_comment']=project_urs.dir_comment
293
294    run_model(**kwargs)
295     
296    if myid==0:
297        export_model(**kwargs)
298    barrier
Note: See TracBrowser for help on using the repository browser.