source: anuga_work/production/tonga/run_fangauta.py @ 5212

Last change on this file since 5212 was 5212, checked in by herve, 16 years ago

update_fanguata

File size: 14.5 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
24from math import radians
25
26# Related major packages
27from anuga.shallow_water import Time_boundary
28from anuga.shallow_water import Domain
29from anuga.shallow_water import Dirichlet_boundary
30from anuga.shallow_water import Transmissive_boundary
31from anuga.shallow_water import File_boundary
32from anuga.shallow_water import Reflective_boundary
33from anuga.shallow_water import Field_boundary
34from Numeric import allclose
35from anuga.shallow_water.data_manager import export_grid
36from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
37from anuga.pmesh.mesh_interface import create_mesh_from_regions
38from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
39from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
40from anuga_parallel.parallel_abstraction import get_processor_name
41from anuga.caching import myhash
42from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
43from anuga.fit_interpolate.benchmark_least_squares import mem_usage
44
45# Application specific imports
46import project_fangauta                 # Definition of file names and polygons
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_fangauta.__file__)+sep+ project_fangauta.__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_fangauta.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    if myid == 0:
88   
89        print 'start create mesh from regions'
90
91        create_mesh_from_regions(project_fangauta.poly_all,
92                             boundary_tags=project_fangauta.boundary_tags,
93                             maximum_triangle_area=project_fangauta.res_poly_all,
94                             filename=project_fangauta.meshes_dir_name+'.msh',
95                             use_cache=False,
96                             verbose=True)
97    barrier()
98
99    scenario='fixed_wave'
100   
101    #-------------------------------------------------------------------------
102    # Setup computational domain
103    #-------------------------------------------------------------------------
104    print 'Setup computational domain'
105
106    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
107    #above don't work
108    domain = Domain(project_fangauta.meshes_dir_name+'.msh', use_cache=False, verbose=True)
109    print 'memory usage before del domain',mem_usage()
110       
111    print domain.statistics()
112    print 'triangles',len(domain)
113   
114    kwargs['act_num_trigs']=len(domain)
115   
116    #-------------------------------------------------------------------------
117    # Setup initial conditions
118    #-------------------------------------------------------------------------
119    if myid == 0:
120
121        print 'Setup initial conditions'
122
123        from polygon import Polygon_function
124        #following sets the stage/water to be offcoast only
125#        IC = Polygon_function( [(project.poly_mainland, -1.0)], default = kwargs['tide'],
126#                                 geo_reference = domain.geo_reference)
127#        domain.set_quantity('stage', IC)
128        domain.set_quantity('stage',kwargs['tide'] )
129#        domain.set_quantity('stage', kwargs['tide'])
130        domain.set_quantity('friction', kwargs['friction']) 
131       
132        print 'Start Set quantity',kwargs['bathy_file']
133
134        domain.set_quantity('elevation', 
135                            filename = kwargs['bathy_file'],
136                            use_cache = False,
137                            verbose = True,
138                            alpha = kwargs['alpha'])
139        print 'Finished Set quantity'
140    barrier()
141   
142    #------------------------------------------------------
143    # Distribute domain to implement parallelism !!!
144    #------------------------------------------------------
145
146    if numprocs > 1:
147        domain=distribute(domain)
148
149    #------------------------------------------------------
150    # Set domain parameters
151    #------------------------------------------------------
152    print 'domain id', id(domain)
153    domain.set_name(kwargs['aa_scenario_name'])
154    domain.set_datadir(kwargs['output_dir'])
155    domain.set_default_order(2) # Apply second order scheme
156    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
157    domain.set_store_vertices_uniquely(False)
158    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
159    domain.tight_slope_limiters = 1
160    #domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
161    print 'domain id', id(domain)
162    domain.beta_h = 0
163
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    Bt = Transmissive_boundary(domain)
172    Br = Reflective_boundary(domain)
173    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
174    Bo = Dirichlet_boundary([kwargs['tide']+5.0,0,0])
175   
176    from math import sin, pi
177    Bw = Time_boundary(domain = domain,
178                                f=lambda t: [0.5*sin(t*radians(2*pi/60)), 0, 0])
179       
180    domain.set_boundary({'land':Bt,'mouth':Bw})
181
182    #----------------------------------------------------------------------------
183    # Evolve system through time
184    #--------------------------------------------------------------------
185    import time
186    t0 = time.time()
187   
188    for t in domain.evolve(yieldstep=10, finaltime = 10800): 
189            domain.write_time()
190            domain.write_boundary_statistics(tags = 'mouth')
191                           
192
193    for t in domain.evolve(10,finaltime=21600,skip_initial_step=True):
194                           
195        Bw = Time_boundary(domain = domain,
196                                        f=lambda t: [0.5*sin(t*radians(2*pi/180)), 0, 0])
197        domain.set_boundary({'land':Bt,'mouth':Bw})
198        domain.write_time()
199        domain.write_boundary_statistics(tags = 'mouth')
200                           
201
202    for t in domain.evolve(10,finaltime=32400,skip_initial_step=True):
203        Bw = Time_boundary(domain = domain,
204                                        f=lambda t: [0.5*sin(t*radians(2*pi/300)), 0, 0])
205        domain.set_boundary({'land':Bt,'mouth':Bw})
206        domain.write_time()
207        domain.write_boundary_statistics(tags = 'mouth')
208                           
209
210    for t in domain.evolve(10,finaltime=43200,skip_initial_step=True):
211        Bw = Time_boundary(domain = domain,
212                                        f=lambda t: [0.5*sin(t*radians(2*pi/420)), 0, 0])
213        domain.set_boundary({'land':Bt,'mouth':Bw})
214        domain.write_time()
215        domain.write_boundary_statistics(tags = 'mouth')
216                           
217     
218    for t in domain.evolve(10,finaltime=54001,skip_initial_step=True):
219        Bw = Time_boundary(domain = domain,
220                                        f=lambda t: [0.5*sin(t*radians(2*pi/540)), 0, 0])
221        domain.set_boundary({'land':Bt,'mouth':Bw})
222        domain.write_time()
223        domain.write_boundary_statistics(tags = 'mouth')
224                           
225                         
226    for t in domain.evolve(10,finaltime=64800,skip_initial_step=True):
227        Bw = Time_boundary(domain = domain,
228                                        f=lambda t: [0.5*sin(t*radians(2*pi/660)), 0, 0])
229        domain.set_boundary({'land':Bt,'mouth':Bw})
230        domain.write_time()
231        domain.write_boundary_statistics(tags = 'mouth')
232                           
233
234    for t in domain.evolve(10,finaltime=75600,skip_initial_step=True):
235        Bw = Time_boundary(domain = domain,
236                                        f=lambda t: [0.5*sin(t*radians(2*pi/780)), 0, 0])
237        domain.set_boundary({'land':Bt,'mouth':Bw})
238        domain.write_time()
239        domain.write_boundary_statistics(tags = 'mouth')
240                           
241
242    for t in domain.evolve(10,finaltime=86400,skip_initial_step=True):
243        Bw = Time_boundary(domain = domain,
244                                        f=lambda t: [0.5*sin(t*radians(2*pi/900)), 0, 0])
245        domain.set_boundary({'land':Bt,'mouth':Bw})
246        domain.write_time()
247        domain.write_boundary_statistics(tags = 'mouth')
248                           
249
250    for t in domain.evolve(10,finaltime=97200,skip_initial_step=True):
251        Bw = Time_boundary(domain = domain,
252                                        f=lambda t: [0.5*sin(t*radians(2*pi/1020)), 0, 0])
253        domain.set_boundary({'land':Bt,'mouth':Bw})
254        domain.write_time()
255        domain.write_boundary_statistics(tags = 'mouth')
256                           
257
258    for t in domain.evolve(10,finaltime=108000,skip_initial_step=True):
259        Bw = Time_boundary(domain = domain,
260                                        f=lambda t: [0.5*sin(t*radians(2*pi/1140)), 0, 0])
261        domain.set_boundary({'land':Bt,'mouth':Bw})
262        domain.write_time()
263        domain.write_boundary_statistics(tags = 'mouth')
264                           
265    for t in domain.evolve(10,finaltime=11800,skip_initial_step=True):
266        Bw = Time_boundary(domain = domain,
267                                        f=lambda t: [0.5*sin(t*radians(2*pi/1260)), 0, 0])
268        domain.set_boundary({'land':Bt,'mouth':Bw})
269        domain.write_time()
270        domain.write_boundary_statistics(tags = 'mouth')
271       
272    x, y = domain.get_maximum_inundation_location()
273    q = domain.get_maximum_inundation_elevation()
274
275    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
276
277    print 'That took %.2f seconds' %(time.time()-t0)
278
279    #kwargs 'completed' must be added to write the final parameters to file
280    kwargs['completed']=str(time.time()-t0)
281   
282    if myid==0:
283        store_parameters(**kwargs)
284    barrier
285   
286    print 'memory usage before del domain1',mem_usage()
287
288def export_model(**kwargs):   
289    #store_parameters(**kwargs)
290   
291#    print 'memory usage before del domain',mem_usage()
292    #del domain
293    print 'memory usage after del domain',mem_usage()
294   
295    swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
296    print'swwfile',swwfile
297   
298    export_grid(swwfile, extra_name_out = 'town',
299            quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
300            #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
301            timestep = None,
302            reduction = max,
303            cellsize = kwargs['export_cellsize'],
304            NODATA_value = -1E-030,
305            easting_min = project_fangauta.eastingmin,
306            easting_max = project_fangauta.eastingmax,
307            northing_min = project_fangauta.northingmin,
308            northing_max = project_fangauta.northingmax,
309            verbose = False,
310            origin = None,
311            datum = 'WGS84',
312            format = 'asc')
313   
314#-------------------------------------------------------------
315if __name__ == "__main__":
316   
317    kwargs={}
318    kwargs['est_num_trigs']=project_fangauta.trigs_min
319    kwargs['num_cpu']=numprocs
320    kwargs['host']=project_fangauta.host
321    kwargs['res_factor']=project_fangauta.res_factor
322    kwargs['starttime']=project_fangauta.starttime
323    kwargs['yieldstep']=project_fangauta.yieldstep
324    kwargs['midtime']=project_fangauta.midtime
325    kwargs['finaltime']=project_fangauta.finaltime
326    kwargs['output_dir']=project_fangauta.output_run_time_dir
327    kwargs['bathy_file']=project_fangauta.combined_dir_name+'.txt'
328#    kwargs['bathy_file']=project.combined_small_dir_name + '.pts'
329    kwargs['boundary_file']=project_fangauta.boundaries_in_dir_name + '.sww'
330    kwargs['file_name']=project_fangauta.home+'detail.csv'
331    kwargs['aa_scenario_name']=project_fangauta.scenario_name
332    kwargs['ab_time']=project_fangauta.time
333    kwargs['res_factor']= project_fangauta.res_factor
334    kwargs['tide']=project_fangauta.tide
335    kwargs['user']=project_fangauta.user
336    kwargs['alpha'] = project_fangauta.alpha
337    kwargs['friction']=project_fangauta.friction
338    kwargs['time_thinning'] = project_fangauta.time_thinning
339    kwargs['dir_comment']=project_fangauta.dir_comment
340    kwargs['export_cellsize']=project_fangauta.export_cellsize
341
342    run_model(**kwargs)
343     
344    if myid==0:
345        export_model(**kwargs)
346    barrier
Note: See TracBrowser for help on using the repository browser.