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

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

run scripts for tonga

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