source: anuga_work/production/busselton/run_busselton.py @ 5480

Last change on this file since 5480 was 5480, checked in by kristy, 14 years ago

Took out a comment

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.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
36#from 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
45numprocs = 1
46myid = 0
47
48def run_model(**kwargs):
49   
50
51    #------------------------------------------------------------------------------
52    # Copy scripts to time stamped output directory and capture screen
53    # output to file
54    #------------------------------------------------------------------------------
55    print "Processor Name:",get_processor_name()
56
57    #copy script must be before screen_catcher
58    #print kwargs
59
60    print 'output_dir',kwargs['output_dir']
61    if myid == 0:
62        copy_code_files(kwargs['output_dir'],__file__, 
63                 dirname(project.__file__)+sep+ project.__name__+'.py' )
64
65        store_parameters(**kwargs)
66
67    #barrier()
68
69    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
70
71    print "Processor Name:",get_processor_name()
72
73    # filenames
74#    meshes_dir_name = project.meshes_dir_name+'.msh'
75
76    # creates copy of code in output dir
77    print 'min triangles', project.trigs_min,
78    print 'Note: This is generally about 20% less than the final amount'
79
80    #--------------------------------------------------------------------------
81    # Create the triangular mesh based on overall clipping polygon with a
82    # tagged
83    # boundary and interior regions defined in project.py along with
84    # resolutions (maximal area of per triangle) for each polygon
85    #--------------------------------------------------------------------------
86
87    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
88    # causes problems with the ability to cache set quantity which takes alot of times
89    if myid == 0:
90   
91        print 'start create mesh from regions'
92       
93        create_mesh_from_regions(project.poly_all,
94                             boundary_tags=project.boundary_tags,
95                             maximum_triangle_area=project.res_poly_all,
96                             interior_regions=project.interior_regions,
97                             filename=project.meshes_dir_name+'.msh',
98                             use_cache=False,
99                             verbose=False)
100    #barrier()
101   
102
103    #-------------------------------------------------------------------------
104    # Setup computational domain
105    #-------------------------------------------------------------------------
106    print 'Setup computational domain'
107
108    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
109    #above don't work
110    domain = Domain(project.meshes_dir_name+'.msh', use_cache=False, verbose=True)
111    print 'memory usage before del domain',mem_usage()
112       
113    print domain.statistics()
114    print 'triangles',len(domain)
115   
116    kwargs['act_num_trigs']=len(domain)
117
118
119    #-------------------------------------------------------------------------
120    # Setup initial conditions
121    #-------------------------------------------------------------------------
122    if myid == 0:
123
124        print 'Setup initial conditions'
125
126        from polygon import Polygon_function
127        #following sets the stage/water to be offcoast only
128        IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'],
129                                 geo_reference = domain.geo_reference)
130        domain.set_quantity('stage', IC)
131#        domain.set_quantity('stage', kwargs['tide'])
132        domain.set_quantity('friction', kwargs['friction']) 
133       
134        print 'Start Set quantity',kwargs['bathy_file']
135
136        domain.set_quantity('elevation', 
137                            filename = kwargs['bathy_file'],
138                            use_cache = True,
139                            verbose = True,
140                            alpha = kwargs['alpha'])
141        print 'Finished Set quantity'
142    #barrier()
143
144
145    #------------------------------------------------------
146    # Distribute domain to implement parallelism !!!
147    #------------------------------------------------------
148
149    if numprocs > 1:
150        domain=distribute(domain)
151
152    #------------------------------------------------------
153    # Set domain parameters
154    #------------------------------------------------------
155    print 'domain id', id(domain)
156    domain.set_name(kwargs['aa_scenario_name'])
157    domain.set_datadir(kwargs['output_dir'])
158    domain.set_default_order(2) # Apply second order scheme
159    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
160    domain.set_store_vertices_uniquely(False)
161    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
162    domain.tight_slope_limiters = 1
163    print 'domain id', id(domain)
164
165    #-------------------------------------------------------------------------
166    # Setup boundary conditions
167    #-------------------------------------------------------------------------
168    print 'Available boundary tags', domain.get_boundary_tags()
169    print 'domain id', id(domain)
170    #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
171
172    print'set_boundary'
173
174    Br = Reflective_boundary(domain)
175    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
176    Bw = Dirichlet_boundary([kwargs['tide']+10.0,0,0]) # To be deleted for FESA runs
177
178    Bf = Field_boundary(kwargs['boundary_file'],
179                    domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], 
180                    use_cache=True, verbose=True)
181    print 'finished reading boundary file'
182    domain.set_boundary({'back': Bd,
183                             'side': Bd,
184                             'ocean': Bd}) #changed from Bf to Bd for large wave
185
186    kwargs['input_start_time']=domain.starttime
187
188    print'finish set boundary'
189
190    #----------------------------------------------------------------------------
191    # Evolve system through time
192    #--------------------------------------------------------------------
193    t0 = time.time()
194
195    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
196                            ,skip_initial_step = False ): 
197        domain.write_time()
198        domain.write_boundary_statistics(tags = 'ocean')
199
200        # To be deleted for FESA runs
201        if allclose(t, 240):
202            domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bw})
203
204        if allclose(t, 1440):
205            domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
206
207
208
209    x, y = domain.get_maximum_inundation_location()
210    q = domain.get_maximum_inundation_elevation()
211
212    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
213
214    print 'That took %.2f seconds' %(time.time()-t0)
215
216    #kwargs 'completed' must be added to write the final parameters to file
217    kwargs['completed']=str(time.time()-t0)
218   
219    if myid==0:
220        store_parameters(**kwargs)
221    #barrier
222   
223    print 'memory usage before del domain1',mem_usage()
224
225 #-------------------------------------------------------------
226if __name__ == "__main__":
227   
228    kwargs={}
229    kwargs['est_num_trigs']=project.trigs_min
230    kwargs['num_cpu']=numprocs
231    kwargs['host']=project.host
232    kwargs['res_factor']=project.res_factor
233    kwargs['starttime']=project.starttime
234    kwargs['yieldstep']=project.yieldstep
235    kwargs['finaltime']=project.finaltime
236   
237    kwargs['output_dir']=project.output_run_time_dir
238    kwargs['bathy_file']=project.combined_dir_name+'.txt'
239#    kwargs['bathy_file']=project.combined_small_dir_name + '.pts'
240    kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
241    kwargs['file_name']=project.home+'detail.csv'
242    kwargs['aa_scenario_name']=project.scenario_name
243    kwargs['ab_time']=project.time
244    kwargs['res_factor']= project.res_factor
245    kwargs['tide']=project.tide
246    kwargs['user']=project.user
247    kwargs['alpha'] = project.alpha
248    kwargs['friction']=project.friction
249    kwargs['time_thinning'] = project.time_thinning
250    kwargs['dir_comment']=project.dir_comment
251    kwargs['export_cellsize']=project.export_cellsize
252   
253
254    run_model(**kwargs)
255     
256    if myid==0:
257        export_model(**kwargs)
258    #barrier
259   
Note: See TracBrowser for help on using the repository browser.