source: anuga_work/production/shark_bay_2007/run_shark_bay.py @ 4689

Last change on this file since 4689 was 4689, checked in by ole, 15 years ago

Steep point back log

File size: 9.6 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 Reflective_boundary
29from anuga.shallow_water import Field_boundary
30from Numeric import allclose
31
32from anuga.pmesh.mesh_interface import create_mesh_from_regions
33from anuga.shallow_water.data_manager import start_screen_catcher
34from anuga.shallow_water.data_manager import copy_code_files
35from anuga.shallow_water.data_manager import 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
39
40# Application specific imports
41import project                 # Definition of file names and polygons
42
43def run_model(**kwargs):
44   
45    tide = kwargs['tide']
46    alpha = kwargs['alpha']
47    friction = kwargs['friction']
48    time_thinning = kwargs['time_thinning']
49    scenario_name = kwargs['aa_scenario_name']
50   
51    #------------------------------------------------------------------------------
52    # Copy scripts to time stamped output directory and capture screen
53    # output to file
54    #------------------------------------------------------------------------------
55
56    #copy script must be before screen_catcher
57    print 'tide',tide
58    kwargs['est_num_trigs']=project.trigs_min
59    kwargs['num_cpu']=numprocs
60    kwargs['host']=project.host
61    kwargs['res_factor']=project.res_factor
62    kwargs['starttime']=project.starttime
63    kwargs['yieldstep']=project.yieldstep
64    kwargs['finaltime']=project.finaltime
65   
66    kwargs['output_dir']=project.output_run_time_dir
67    kwargs['bathy_file']=project.combined_dir_name + '.pts'
68#    kwargs['bathy_file']=project.combined_small_dir_name + '.pts'
69    kwargs['boundary_file']=project.boundary_name + '.sww'
70#    kwargs['Completed']=''
71
72    print 'output_dir',kwargs['output_dir']
73    if myid == 0:
74        copy_code_files(kwargs['output_dir'],__file__, 
75                        dirname(project.__file__)+sep+project.__name__+'.py' )
76
77        store_parameters(**kwargs)
78
79    barrier()
80
81    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
82
83    print 'Processor Name:', get_processor_name()
84
85    #--------------------------------------------------------------------------
86    # Create the triangular mesh based on overall clipping polygon with a
87    # tagged
88    # boundary and interior regions defined in project.py along with
89    # resolutions (maximal area of per triangle) for each polygon
90    #--------------------------------------------------------------------------
91
92    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
93    # causes problems with the ability to cache set quantity which takes alot of times
94
95    # FIXME (Ole): Move to build_shark_bay.py
96    if myid == 0:
97   
98        print 'start create mesh from regions'
99       
100        create_mesh_from_regions(project.bounding_polygon,
101                                 boundary_tags=project.boundary_tags,
102                                 maximum_triangle_area=project.res_bounding_polygon,
103                                 interior_regions=project.interior_regions,
104                                 filename=project.mesh_name+'.msh',
105                                 use_cache=False,
106                                 verbose=True)
107
108    #-------------------------------------------------------------------------
109    # Setup computational domain
110    #-------------------------------------------------------------------------
111    print 'Setup computational domain'
112
113    #domain = cache(Domain, (mesh_name), {'use_cache':True, 'verbose':True}, verbose=True)
114    #above don't work
115    domain = Domain(project.mesh_name+'.msh',
116                    use_cache=False, verbose=True)
117     
118    print domain.statistics()
119    print 'triangles',len(domain)
120   
121    kwargs['act_num_trigs']=len(domain)
122
123    #-------------------------------------------------------------------------
124    # Setup initial conditions
125    #-------------------------------------------------------------------------
126    if myid == 0:
127
128        IC = project.tide
129        #if project.tide == 0.0:
130        #    print 'Tide is MSL'
131        #    IC = 0.0
132        #else:   
133        #    print 'Create onshore polygon'
134        #    from polygon import Polygon_function
135        #    #following sets the stage/water to be offcoast only
136        #    IC = Polygon_function([(project.onshore_polygon, -1.0)],
137        #                          default = tide,
138        #                          geo_reference = domain.geo_reference)
139           
140
141        print 'Set initial conditions'
142        domain.set_quantity('stage', IC)
143        domain.set_quantity('friction', friction)
144        print 'Set elevation'
145        domain.set_quantity('elevation', 
146                            filename = kwargs['bathy_file'],
147                            use_cache = True,
148                            verbose = True,
149                            alpha = alpha)
150       
151
152    #------------------------------------------------------
153    # Distribute domain to implement parallelism !!!
154    #------------------------------------------------------
155    barrier()
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(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    print 'domain id', id(domain)
171    domain.beta_h = 0
172    #domain.tight_slope_limiters = 1
173   
174
175    #-------------------------------------------------------------------------
176    # Setup boundary conditions
177    #-------------------------------------------------------------------------
178    print 'Available boundary tags', domain.get_boundary_tags()
179    print 'domain id', id(domain)
180    #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
181
182    Bf = Field_boundary(kwargs['boundary_file'],
183                        domain,
184                        time_thinning=time_thinning,
185                        mean_stage=tide,
186                        momentum_scale=project.momentum_scale,                       
187                        use_cache=False,
188                        verbose=True)
189                   
190    kwargs['input_start_time']=domain.starttime
191
192    print 'finished reading boundary file'
193
194    Br = Reflective_boundary(domain)
195    Bd = Dirichlet_boundary([tide,0,0])
196
197    print'set_boundary'
198
199    #domain.set_boundary({'back': Br,
200    #                     'side': Bd,
201    #                     'ocean': Bf})
202    domain.set_boundary({'tide': Bd,
203    #                     'ocean': Dirichlet_boundary([3,0,0])})
204                         'ocean': Bf})     
205    print'finish set boundary'
206
207    #----------------------------------------------------------------------------
208    # Evolve system through time
209    #----------------------------------------------------------------------------
210
211    t0 = time.time()
212
213    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['starttime']): 
214        domain.write_time()
215        domain.write_boundary_statistics(tags = 'ocean')     
216
217
218    # Reset stage in study area
219    print 'Reset initial condition in study area'
220    domain.set_quantity('stage', project.tide,
221                        polygon=project.poly_tsunami_approach_area)
222   
223
224    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = 18000
225                       ,skip_initial_step = True): 
226        domain.write_time()
227        domain.write_boundary_statistics(tags = 'ocean')     
228   
229    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']
230                       ,skip_initial_step = True): 
231        domain.write_time()
232        domain.write_boundary_statistics(tags = 'ocean')   
233
234    x, y = domain.get_maximum_inundation_location()
235    q = domain.get_maximum_inundation_elevation()
236
237    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
238
239    print 'That took %.2f seconds' %(time.time()-t0)
240   
241    kwargs['completed']=str(time.time()-t0)
242   
243    store_parameters(**kwargs)
244
245#-------------------------------------------------------------
246if __name__ == "__main__":
247
248    run_model(file_name=project.home+'detail.csv', aa_scenario_name=project.scenario_name,
249              ab_time=project.time, res_factor= project.res_factor, tide=project.tide, user=project.user,
250              alpha = project.alpha, friction=project.friction,
251              time_thinning = project.time_thinning,
252              dir_comment=project.dir_comment)
253
254
Note: See TracBrowser for help on using the repository browser.