source: anuga_work/production/dampier_2006/run_dampier_urs.py @ 4631

Last change on this file since 4631 was 4631, checked in by ole, 17 years ago

Refactored limit2007 to tight_slope_limiters (using eclipse)

File size: 10.9 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,umask
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# Application specific imports
42import project_urs                 # Definition of file names and polygons
43
44def run_model(**kwargs):
45   
46    tide = kwargs['tide']
47    alpha = kwargs['alpha']
48    friction = kwargs['friction']
49    time_thinning = kwargs['time_thinning']
50    scenario_name = kwargs['aa_scenario_name']
51   
52    #------------------------------------------------------------------------------
53    # Copy scripts to time stamped output directory and capture screen
54    # output to file
55    #------------------------------------------------------------------------------
56
57    #copy script must be before screen_catcher
58    print 'tide',tide
59    kwargs['est_num_trigs']=project_urs.trigs_min
60    kwargs['num_cpu']=numprocs
61    kwargs['host']=project_urs.host
62    kwargs['res_factor']=project_urs.res_factor
63    kwargs['starttime']=project_urs.starttime
64    kwargs['yieldstep']=project_urs.yieldstep
65    kwargs['finaltime']=project_urs.finaltime
66   
67    kwargs['output_dir']=project_urs.output_run_time_dir
68    kwargs['bathy_file']=project_urs.combined_dir_name + '.txt'
69#    kwargs['bathy_file']=project_urs.combined_small_dir_name + '.pts'
70    kwargs['boundary_file']=project_urs.boundaries_in_dir_name + '.sww'
71#    kwargs['Completed']=''
72   
73    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
74
75    print 'output_dir',kwargs['output_dir']
76    if myid == 0:
77        copy_code_files(kwargs['output_dir'],__file__, 
78                 dirname(project_urs.__file__)+sep+ project_urs.__name__+'.py' )
79
80        store_parameters(**kwargs)
81
82    barrier()
83
84    print "Processor Name:",get_processor_name()
85
86    # creates copy of code in output dir
87    print 'min triangles', project_urs.trigs_min,
88    print 'Note: This is generally about 20% less than the final amount'
89
90    #--------------------------------------------------------------------------
91    # Create the triangular mesh based on overall clipping polygon with a
92    # tagged
93    # boundary and interior regions defined in project_urs.py along with
94    # resolutions (maximal area of per triangle) for each polygon
95    #--------------------------------------------------------------------------
96
97    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
98    # causes problems with the ability to cache set quantity which takes alot of times
99    if myid == 0:
100   
101        print 'start create mesh from regions'
102
103        create_mesh_from_regions(project_urs.poly_all,
104                             boundary_tags={'back': [2,3], 'side': [0, 1, 4],
105                                            'ocean': [5]},
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    # 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_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 = tide,
137                                 geo_reference = domain.geo_reference)
138        domain.set_quantity('stage', IC)
139        domain.set_quantity('friction', 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 = 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(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
171    #domain.tight_slope_limiters = 1
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_urs.boundaries_dir_namea + '.sww'
179
180    Bf = Field_boundary(kwargs['boundary_file'],
181                    domain, time_thinning=time_thinning, mean_stage=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([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']=str(time.time()-t0)
226     
227    if myid == 0:
228 #       store_parameters(**kwargs)
229       
230        print 'memory usage before del domain',mem_usage()
231#        del domain
232        print 'memory usage after del domain',mem_usage()
233
234        swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
235        export_grid(swwfile, extra_name_out = 'town',
236                quantities = ['elevation','depth','stage','speed'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
237                timestep = None,
238                reduction = max,
239                cellsize = 25,
240                NODATA_value = -9999,
241                easting_min = project_urs.eastingmin,
242                easting_max = project_urs.eastingmax,
243                northing_min = project_urs.northingmin,
244                northing_max = project_urs.northingmax,
245                verbose = True,
246                origin = None,
247                datum = 'WGS84',
248                format = 'asc')
249   
250        buildings_filename = project_urs.buildings_filename
251        buildings_filename_out = project_urs.buildings_filename_out
252               
253        inundation_damage(swwfile+'.sww', buildings_filename, buildings_filename_out)
254        print '\n Augmented building file written to %s \n' %buildings_filename_out
255   
256    barrier()
257#-------------------------------------------------------------
258if __name__ == "__main__":
259
260    run_model(file_name=project_urs.home+'detail.csv', aa_scenario_name=project_urs.scenario_name,
261              ab_time=project_urs.time, res_factor= project_urs.res_factor, tide=project_urs.tide, user=project_urs.user,
262              alpha = project_urs.alpha, friction=project_urs.friction,
263              time_thinning = project_urs.time_thinning,
264              dir_comment=project_urs.dir_comment)
265
266
Note: See TracBrowser for help on using the repository browser.