source: anuga_work/production/onslow_2006/run_onslow_urs.py @ 4587

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

update model 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_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
34
35from anuga.pmesh.mesh_interface import create_mesh_from_regions
36from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
37from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
38from anuga_parallel.parallel_abstraction import get_processor_name
39from anuga.caching import myhash
40from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
41from anuga.fit_interpolate.benchmark_least_squares import mem_usage
42
43# Application specific imports
44import project_urs                 # Definition of file names and polygons
45
46def run_model(**kwargs):
47   
48#    scenario_name = kwargs['aa_scenario_name']
49   
50    #------------------------------------------------------------------------------
51    # Copy scripts to time stamped output directory and capture screen
52    # output to file
53    #------------------------------------------------------------------------------
54
55    #copy script must be before screen_catcher
56    print 'tide',kwargs['tide']
57    kwargs['est_num_trigs']=project_urs.trigs_min
58    kwargs['num_cpu']=numprocs
59    kwargs['host']=project_urs.host
60    kwargs['res_factor']=project_urs.res_factor
61    kwargs['starttime']=project_urs.starttime
62    kwargs['yieldstep']=project_urs.yieldstep
63    kwargs['finaltime']=project_urs.finaltime
64   
65    kwargs['output_dir']=project_urs.output_run_time_dir
66    kwargs['bathy_file']=project_urs.combined_dir_name + '.pts'
67#    kwargs['bathy_file']=project_urs.combined_small_dir_name + '.txt'
68    kwargs['boundary_file']=project_urs.boundaries_in_dir_name + '.sww'
69    '''
70    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
71       
72    print 'output_dir',kwargs['output_dir']
73    if myid == 0:
74        copy_code_files(kwargs['output_dir'],__file__,
75                 dirname(project_urs.__file__)+sep+ project_urs.__name__+'.py' )
76
77        store_parameters(**kwargs)
78
79    barrier()
80
81    print "Processor Name:",get_processor_name()
82
83    # creates copy of code in output dir
84    print 'min triangles', project_urs.trigs_min,
85    print 'Note: This is generally about 20% less than the final amount'
86
87
88    #--------------------------------------------------------------------------
89    # Create the triangular mesh based on overall clipping polygon with a
90    # tagged
91    # boundary and interior regions defined in project_urs.py along with
92    # resolutions (maximal area of per triangle) for each polygon
93    #--------------------------------------------------------------------------
94    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
95    # causes problems with the ability to cache set quantity which takes alot of times
96
97    if myid == 0:
98   
99        print 'start create mesh from regions'
100
101        create_mesh_from_regions(project_urs.poly_all,
102                             boundary_tags={'back': [4, 5], 'side': [3, 6],
103                                            'ocean': [0, 1, 2]},
104                             maximum_triangle_area=project_urs.res_poly_all,
105                             interior_regions=project_urs.interior_regions,
106                             filename=project_urs.meshes_dir_name+'.msh',
107                             use_cache=False,
108                             verbose=True)
109    barrier()
110   
111
112    #-------------------------------------------------------------------------
113    # Setup computational domain
114    #-------------------------------------------------------------------------
115    print 'Setup computational domain'
116    #from caching import cache
117
118    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
119    #above don't work
120    domain = Domain(project_urs.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_urs.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    #------------------------------------------------------
153    # Distribute domain to implement parallelism !!!
154    #------------------------------------------------------
155
156    if numprocs > 1:
157        domain=distribute(domain)
158
159    #------------------------------------------------------
160    # Set domain parameters
161    #------------------------------------------------------
162    print 'domain id', id(domain)
163    domain.set_name(kwargs['aa_scenario_name'])
164    domain.set_datadir(kwargs['output_dir'])
165    domain.set_default_order(2) # Apply second order scheme
166    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
167    domain.set_store_vertices_uniquely(False)
168    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
169#    domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
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_urs.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']=str(time.time()-t0)
227    '''
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   
236        swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
237        print 'swwfile',swwfile
238   
239        export_grid(swwfile, extra_name_out = 'town',
240                quantities = ['speed'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
241#                quantities = ['elevation','depth','stage','speed'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
242                timestep = None,
243                reduction = max,
244                cellsize = 25,
245                NODATA_value = -9999,
246                easting_min = project_urs.eastingmin,
247                easting_max = project_urs.eastingmax,
248                northing_min = project_urs.northingmin,
249                northing_max = project_urs.northingmax,
250                verbose = True,
251                origin = None,
252                datum = 'WGS84',
253                format = 'asc')
254   
255        buildings_filename = project_urs.buildings_filename
256        buildings_filename_out = project_urs.buildings_filename_out
257               
258        #inundation_damage(swwfile+'.sww', buildings_filename, buildings_filename_out)
259        print '\n Augmented building file written to %s \n' \
260           %buildings_filename_out
261    barrier()
262   
263#-------------------------------------------------------------
264if __name__ == "__main__":
265
266    run_model(file_name=project_urs.home+'detail.csv', aa_scenario_name=project_urs.scenario_name,
267              ab_time=project_urs.time, res_factor= project_urs.res_factor, tide=project_urs.tide, user=project_urs.user,
268              alpha = project_urs.alpha, friction=project_urs.friction,
269              time_thinning = project_urs.time_thinning,
270              dir_comment=project_urs.dir_comment)
271
272
Note: See TracBrowser for help on using the repository browser.