source: anuga_work/production/onslow_2006/run_onslow_grad.py @ 5442

Last change on this file since 5442 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: 9.5 KB
Line 
1"""Script for running tsunami inundation scenario for Onslow, WA, Australia.
2
3Source data such as elevation and boundary data is assumed to be available in
4directories specified by project_grad_grad.py
5The output sww file is stored in project_grad_grad.output_run_time_dir
6
7The scenario is defined by a triangular mesh created from project_grad_grad.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
32
33from anuga.pmesh.mesh_interface import create_mesh_from_regions
34from anuga.shallow_water.data_manager import start_screen_catcher
35from anuga.shallow_water.data_manager import copy_code_files
36from anuga.shallow_water.data_manager import 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
40
41from anuga.shallow_water.data_manager import export_grid
42from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
43
44# Application specific imports
45import project_grad                 # Definition of file names and polygons
46
47def run_model(**kwargs):
48   
49    tide = kwargs['tide']
50    alpha = kwargs['alpha']
51    friction = kwargs['friction']
52    time_thinning = kwargs['time_thinning']
53    scenario_name = kwargs['aa_scenario_name']
54   
55    #------------------------------------------------------------------------------
56    # Copy scripts to time stamped output directory and capture screen
57    # output to file
58    #------------------------------------------------------------------------------
59
60    #copy script must be before screen_catcher
61    print 'tide',tide
62    kwargs['est_num_trigs']=project_grad.trigs_min
63    kwargs['num_cpu']=numprocs
64    kwargs['host']=project_grad.host
65    kwargs['res_factor']=project_grad.res_factor
66    kwargs['starttime']=project_grad.starttime
67    kwargs['yieldstep']=project_grad.yieldstep
68    kwargs['finaltime']=project_grad.finaltime
69   
70    kwargs['output_dir']=project_grad.output_run_time_dir
71    kwargs['bathy_file']=project_grad.combined_dir_name + '.pts'
72    kwargs['boundary_file']=project_grad.boundaries_in_dir+sep+project_grad.boundaries_name+'_40000.sww'
73#    kwargs['Completed']=''
74
75    print 'output_dir',kwargs['output_dir']
76    if myid == 0:
77        copy_code_files(kwargs['output_dir'],__file__, 
78                        dirname(project_grad.__file__)+sep+project_grad.__name__+'.py' )
79
80        store_parameters(**kwargs)
81
82    barrier()
83
84    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
85
86    print 'Processor Name:', get_processor_name()
87
88    #--------------------------------------------------------------------------
89    # Create the triangular mesh based on overall clipping polygon with a
90    # tagged
91    # boundary and interior regions defined in project_grad.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
98    # FIXME (Ole): Move to build_shark_bay.py
99    if myid == 0:
100   
101        print 'start create mesh from regions'
102       
103        create_mesh_from_regions(project_grad.bounding_polygon,
104                                 boundary_tags=project_grad.boundary_tags,
105                                 maximum_triangle_area=project_grad.res_bounding_polygon,
106                                 interior_regions=project_grad.interior_regions,
107                                 filename=project_grad.mesh_name+'.msh',
108                                 use_cache=True,
109                                 verbose=True)
110
111    #-------------------------------------------------------------------------
112    # Setup computational domain
113    #-------------------------------------------------------------------------
114    print 'Setup computational domain'
115
116    #domain = cache(Domain, (mesh_name), {'use_cache':True, 'verbose':True}, verbose=True)
117    #above don't work
118    domain = Domain(project_grad.mesh_name+'.msh',
119                    use_cache=True, 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        IC = project_grad.tide
132        #if project_grad.tide == 0.0:
133        #    print 'Tide is MSL'
134        #    IC = 0.0
135        #else:   
136        #    print 'Create onshore polygon'
137        #    from polygon import Polygon_function
138        #    #following sets the stage/water to be offcoast only
139        #    IC = Polygon_function([(project_grad.onshore_polygon, -1.0)],
140        #                          default = tide,
141        #                          geo_reference = domain.geo_reference)
142           
143
144        print 'Set initial conditions'
145        domain.set_quantity('stage', IC)
146        domain.set_quantity('friction', friction)
147        print 'Set elevation'
148        domain.set_quantity('elevation', 
149                            filename = kwargs['bathy_file'],
150                            use_cache = True,
151                            verbose = True,
152                            alpha = project_grad.alpha)
153       
154
155    #------------------------------------------------------
156    # Distribute domain to implement parallelism !!!
157    #------------------------------------------------------
158    barrier()
159    if numprocs > 1:
160        domain=distribute(domain)
161
162    #------------------------------------------------------
163    # Set domain parameters
164    #------------------------------------------------------
165    print 'domain id', id(domain)
166    domain.set_name(scenario_name)
167    domain.set_datadir(kwargs['output_dir'])
168    domain.set_default_order(2) # Apply second order scheme
169    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
170    domain.set_store_vertices_uniquely(False)
171    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
172    domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
173    print 'domain id', id(domain)
174
175    domain.tight_slope_limiters = 0
176   
177
178    #-------------------------------------------------------------------------
179    # Setup boundary conditions
180    #-------------------------------------------------------------------------
181    print 'Available boundary tags', domain.get_boundary_tags()
182    print 'domain id', id(domain)
183    #print 'Reading Boundary file',project_grad.boundaries_dir_namea + '.sww'
184
185    Bf = Field_boundary(kwargs['boundary_file'],
186                        domain,
187                        time_thinning=time_thinning,
188                        mean_stage=tide,
189                        use_cache=True,
190                        verbose=True)
191                   
192    kwargs['input_start_time']=domain.starttime
193
194    print 'finished reading boundary file'
195
196    Br = Reflective_boundary(domain)
197    Bd = Dirichlet_boundary([tide,0,0])
198
199    print'set_boundary'
200
201    domain.set_boundary({'back': Br,
202                         'side': Bd,
203                         'ocean': Bf})
204
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    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = 21600
218                       ,skip_initial_step = True): 
219        domain.write_time()
220        domain.write_boundary_statistics(tags = 'ocean')     
221   
222    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']
223                       ,skip_initial_step = True): 
224        domain.write_time()
225        domain.write_boundary_statistics(tags = 'ocean')   
226
227    x, y = domain.get_maximum_inundation_location()
228    q = domain.get_maximum_inundation_elevation()
229
230    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
231
232    print 'That took %.2f seconds' %(time.time()-t0)
233   
234    kwargs['completed']=str(time.time()-t0)
235   
236    store_parameters(**kwargs)
237
238#-------------------------------------------------------------
239if __name__ == "__main__":
240
241    run_model(file_name=project_grad.home+'detail.csv', aa_scenario_name=project_grad.scenario_name,
242              ab_time=project_grad.time, res_factor= project_grad.res_factor, tide=project_grad.tide, user=project_grad.user,
243              alpha = project_grad.alpha, friction=project_grad.friction,
244              time_thinning = project_grad.time_thinning,
245              dir_comment=project_grad.dir_comment)
246
247
Note: See TracBrowser for help on using the repository browser.