source: anuga_work/production/exmouth_2006/run_exmouth.py @ 4546

Last change on this file since 4546 was 4546, checked in by nick, 17 years ago

update scripts

File size: 10.8 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.py
5The output sww file is stored in project.output_run_time_dir
6
7The scenario is defined by a triangular mesh created from project.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                 # Definition of file names and polygons
44
45def run_model(**kwargs):
46   
47    #scenario_name = kwargs['aa_scenario_name']
48   
49    #------------------------------------------------------------------------------
50    # Copy scripts to time stamped output directory and capture screen
51    # output to file
52    #------------------------------------------------------------------------------
53
54    #copy script must be before screen_catcher
55    print 'tide',kwargs['tide']
56    kwargs['est_num_trigs']=project.trigs_min
57    kwargs['num_cpu']=numprocs
58    kwargs['host']=project.host
59    kwargs['res_factor']=project.res_factor
60    kwargs['starttime']=project.starttime
61    kwargs['yieldstep']=project.yieldstep
62    kwargs['finaltime']=project.finaltime
63   
64    kwargs['output_dir']=project.output_run_time_dir
65    kwargs['bathy_file']=project.combined_dir_name+'.txt'
66#    kwargs['bathy_file']=project.combined_small_dir_name + '.txt'
67    kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
68   
69    print 'output_dir',kwargs['output_dir']
70    if myid == 0:
71        copy_code_files(kwargs['output_dir'],__file__, 
72                 dirname(project.__file__)+sep+ project.__name__+'.py' )
73
74        store_parameters(**kwargs)
75
76    barrier()
77
78    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
79
80    print "Processor Name:",get_processor_name()
81
82    # filenames
83#    meshes_dir_name = project.meshes_dir_name+'.msh'
84
85    # creates copy of code in output dir
86    print 'min triangles', project.trigs_min,
87    print 'Note: This is generally about 20% less than the final amount'
88
89    #--------------------------------------------------------------------------
90    # Create the triangular mesh based on overall clipping polygon with a
91    # tagged
92    # boundary and interior regions defined in project.py along with
93    # resolutions (maximal area of per triangle) for each polygon
94    #--------------------------------------------------------------------------
95
96    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
97    # causes problems with the ability to cache set quantity which takes alot of times
98    if myid == 0:
99   
100        print 'start create mesh from regions'
101
102        create_mesh_from_regions(project.poly_all,
103                         boundary_tags={'back': [3,4,5], 'side': [2,6],
104                                        'ocean': [0,1,7]},
105                             maximum_triangle_area=project.res_poly_all,
106                             interior_regions=project.interior_regions,
107                             filename=project.meshes_dir_name+'.msh',
108                             use_cache=False,
109                             verbose=True)
110        copy(project.meshes_dir_name+'.msh',project.output_dir+project.scenario_name+'.msh')
111    barrier()
112
113    #-------------------------------------------------------------------------
114    # Setup computational domain
115    #-------------------------------------------------------------------------
116    print 'Setup computational domain'
117
118    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
119    #above don't work
120    domain = Domain(project.meshes_dir_name+'.msh', use_cache=False, verbose=True)
121       
122    print domain.statistics()
123    print 'triangles',len(domain)
124   
125    kwargs['act_num_trigs']=len(domain)
126
127    #-------------------------------------------------------------------------
128    # Setup initial conditions
129    #-------------------------------------------------------------------------
130    if myid == 0:
131
132        print 'Setup initial conditions'
133
134        from polygon import Polygon_function
135        #following sets the stage/water to be offcoast only
136        IC = Polygon_function( [(project.poly_mainland, -1.0)], default = kwargs['tide'],
137                                 geo_reference = domain.geo_reference)
138        domain.set_quantity('stage', IC)
139        domain.set_quantity('friction', kwargs['friction']) 
140       
141        print 'Start Set quantity'
142
143        domain.set_quantity('elevation', 
144                            filename = kwargs['bathy_file'],
145                            use_cache = True,
146                            verbose = True,
147                            alpha = kwargs['alpha'])
148        print 'Finished Set quantity'
149    barrier()
150
151    #------------------------------------------------------
152    # Distribute domain to implement parallelism !!!
153    #------------------------------------------------------
154
155    if numprocs > 1:
156        domain=distribute(domain)
157
158    #------------------------------------------------------
159    # Set domain parameters
160    #------------------------------------------------------
161    print 'domain id', id(domain)
162    domain.set_name(kwargs['aa_scenario_name'])
163    domain.set_datadir(kwargs['output_dir'])
164    domain.set_default_order(2) # Apply second order scheme
165    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
166    domain.set_store_vertices_uniquely(False)
167    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
168    #domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
169    #print 'domain id', id(domain)
170    domain.beta_h = 0 #sets the surface of the triangle to follow the bathy
171    #domain.H0=0.01 #controls the flux limiter (limiter2007)
172    #domain.limit2007 = 1 #minimises creep
173
174    #-------------------------------------------------------------------------
175    # Setup boundary conditions
176    #-------------------------------------------------------------------------
177    print 'Available boundary tags', domain.get_boundary_tags()
178    print 'domain id', id(domain)
179    #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
180
181    Bf = Field_boundary(kwargs['boundary_file'],
182                    domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], 
183                    use_cache=True, verbose=True)
184                   
185    kwargs['input_start_time']=domain.starttime
186
187    print 'finished reading boundary file'
188
189    Br = Reflective_boundary(domain)
190    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
191
192    print'set_boundary'
193
194    domain.set_boundary({'back': Br,
195                         'side': Bd,
196                         'ocean': Bf}) 
197    print'finish set boundary'
198
199    #----------------------------------------------------------------------------
200    # Evolve system through time
201    #----------------------------------------------------------------------------
202
203    t0 = time.time()
204
205    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['starttime']): 
206        domain.write_time()
207        domain.write_boundary_statistics(tags = 'ocean')     
208
209    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = 21600
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
219    x, y = domain.get_maximum_inundation_location()
220    q = domain.get_maximum_inundation_elevation()
221
222    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
223
224    print 'That took %.2f seconds' %(time.time()-t0)
225
226    #kwargs 'completed' must be added to write the final parameters to file
227    kwargs['completed']=str(time.time()-t0)
228    if myid == 0:
229        store_parameters(**kwargs)
230       
231        print 'memory usage before del domain',mem_usage()
232        del domain
233        print 'memory usage after del domain',mem_usage()
234       
235        swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
236        export_grid(swwfile, extra_name_out = 'town',
237#                quantities = ['depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
238                quantities = ['elevation','depth','stage','speed'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
239                timestep = None,
240                reduction = max,
241                cellsize = 25,
242                NODATA_value = -9999,
243                easting_min = project.eastingmin,
244                easting_max = project.eastingmax,
245                northing_min = project.northingmin,
246                northing_max = project.northingmax,
247                verbose = True,
248                origin = None,
249                datum = 'WGS84',
250                format = 'asc')
251   
252        buildings_filename = project.buildings_filename
253        buildings_filename_out = project.buildings_filename_out
254               
255        inundation_damage(swwfile+'.sww', buildings_filename, buildings_filename_out)
256        print '\n Augmented building file written to %s \n' %buildings_filename_out
257   
258    barrier()
259   
260#-------------------------------------------------------------
261if __name__ == "__main__":
262
263    run_model(file_name=project.home+'detail.csv', aa_scenario_name=project.scenario_name,
264              ab_time=project.time, res_factor= project.res_factor, tide=project.tide, user=project.user,
265              alpha = project.alpha, friction=project.friction,
266              time_thinning = project.time_thinning,
267              dir_comment=project.dir_comment)
268
269
Note: See TracBrowser for help on using the repository browser.