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

Last change on this file since 5587 was 5442, checked in by ole, 16 years ago

Retired h-limiter and beta_h as per ticket:194.
All unit tests and validation tests pass.
File size: 10.7 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    barrier()
111
112    #-------------------------------------------------------------------------
113    # Setup computational domain
114    #-------------------------------------------------------------------------
115    print 'Setup computational domain'
116
117    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
118    #above don't work
119    domain = Domain(project.meshes_dir_name+'.msh', use_cache=False, verbose=True)
120       
121    print domain.statistics()
122    print 'triangles',len(domain)
123   
124    kwargs['act_num_trigs']=len(domain)
125
126    #-------------------------------------------------------------------------
127    # Setup initial conditions
128    #-------------------------------------------------------------------------
129    if myid == 0:
130
131        print 'Setup initial conditions'
132
133        from polygon import Polygon_function
134        #following sets the stage/water to be offcoast only
135        IC = Polygon_function( [(project.poly_mainland, -1.0)], default = kwargs['tide'],
136                                 geo_reference = domain.geo_reference)
137        domain.set_quantity('stage', IC)
138        domain.set_quantity('friction', kwargs['friction']) 
139       
140        print 'Start Set quantity'
141
142        domain.set_quantity('elevation', 
143                            filename = kwargs['bathy_file'],
144                            use_cache = True,
145                            verbose = True,
146                            alpha = kwargs['alpha'])
147        print 'Finished Set quantity'
148    barrier()
149
150    #------------------------------------------------------
151    # Distribute domain to implement parallelism !!!
152    #------------------------------------------------------
153
154    if numprocs > 1:
155        domain=distribute(domain)
156
157    #------------------------------------------------------
158    # Set domain parameters
159    #------------------------------------------------------
160    print 'domain id', id(domain)
161    domain.set_name(kwargs['aa_scenario_name'])
162    domain.set_datadir(kwargs['output_dir'])
163    domain.set_default_order(2) # Apply second order scheme
164    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
165    domain.set_store_vertices_uniquely(False)
166    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
167    #domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
168    #print 'domain id', id(domain)
169
170    #domain.H0=0.01 #controls the flux limiter (limiter2007)
171    domain.tight_slope_limiters = 0 #minimises creep
172
173    #-------------------------------------------------------------------------
174    # Setup boundary conditions
175    #-------------------------------------------------------------------------
176    print 'Available boundary tags', domain.get_boundary_tags()
177    print 'domain id', id(domain)
178    #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
179
180    Bf = Field_boundary(kwargs['boundary_file'],
181                    domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], 
182                    use_cache=True, verbose=True)
183                   
184    kwargs['input_start_time']=domain.starttime
185
186    print 'finished reading boundary file'
187
188    Br = Reflective_boundary(domain)
189    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
190
191    print'set_boundary'
192
193    domain.set_boundary({'back': Br,
194                         'side': Bd,
195                         'ocean': Bf}) 
196    print'finish set boundary'
197
198    #----------------------------------------------------------------------------
199    # Evolve system through time
200    #----------------------------------------------------------------------------
201
202    t0 = time.time()
203
204    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['starttime']): 
205        domain.write_time()
206        domain.write_boundary_statistics(tags = 'ocean')     
207
208    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = 21600
209                       ,skip_initial_step = True): 
210        domain.write_time()
211        domain.write_boundary_statistics(tags = 'ocean')     
212   
213    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']
214                       ,skip_initial_step = True): 
215        domain.write_time()
216        domain.write_boundary_statistics(tags = 'ocean')   
217
218    x, y = domain.get_maximum_inundation_location()
219    q = domain.get_maximum_inundation_elevation()
220
221    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
222
223    print 'That took %.2f seconds' %(time.time()-t0)
224
225    #kwargs 'completed' must be added to write the final parameters to file
226    kwargs['completed']=str(time.time()-t0)
227   
228   
229    if myid == 0:
230        store_parameters(**kwargs)
231       
232        print 'memory usage before del domain',mem_usage()
233        #del domain
234        print 'memory usage after del domain',mem_usage()
235       
236        swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
237        export_grid(swwfile, extra_name_out = 'town',
238                quantities = ['speed'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
239#                quantities = ['elevation','depth','speed','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
240                timestep = None,
241                reduction = max,
242                cellsize = 25,
243                NODATA_value = -9999,
244                easting_min = project.eastingmin,
245                easting_max = project.eastingmax,
246                northing_min = project.northingmin,
247                northing_max = project.northingmax,
248                verbose = True,
249                origin = None,
250                datum = 'WGS84',
251                format = 'asc')
252   
253        buildings_filename = project.buildings_filename
254        buildings_filename_out = project.buildings_filename_out
255               
256        #inundation_damage(swwfile+'.sww', buildings_filename, buildings_filename_out)
257        print '\n Augmented building file written to %s \n' %buildings_filename_out
258   
259    barrier()
260   
261#-------------------------------------------------------------
262if __name__ == "__main__":
263
264    run_model(file_name=project.home+'detail.csv', aa_scenario_name=project.scenario_name,
265              ab_time=project.time, res_factor= project.res_factor, tide=project.tide, user=project.user,
266              alpha = project.alpha, friction=project.friction,
267              time_thinning = project.time_thinning,
268              dir_comment=project.dir_comment)
269
270
Note: See TracBrowser for help on using the repository browser.