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

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

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