source: anuga_work/production/FSM/run_kolonia.py @ 5442

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

Retired h-limiter and beta_h as per ticket:194.
All unit tests and validation tests pass.

File size: 11.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 Time_boundary
27from anuga.shallow_water import Domain
28from anuga.shallow_water import Dirichlet_boundary, Transmissive_boundary
29from anuga.shallow_water import File_boundary
30from anuga.shallow_water import Reflective_boundary
31from anuga.shallow_water import Field_boundary
32from Numeric import allclose
33from anuga.shallow_water.data_manager import export_grid
34from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
35from anuga.pmesh.mesh_interface import create_mesh_from_regions
36from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,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
40from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
41from anuga.fit_interpolate.benchmark_least_squares import mem_usage
42
43# Application specific imports
44import project                 # Definition of file names and polygons
45
46def run_model(**kwargs):
47   
48
49    #------------------------------------------------------------------------------
50    # Copy scripts to time stamped output directory and capture screen
51    # output to file
52    #------------------------------------------------------------------------------
53    print "Processor Name:",get_processor_name()
54
55    #copy script must be before screen_catcher
56    #print kwargs
57
58    print 'output_dir',kwargs['output_dir']
59    if myid == 0:
60        copy_code_files(kwargs['output_dir'],__file__, 
61                 dirname(project.__file__)+sep+ project.__name__+'.py' )
62
63        store_parameters(**kwargs)
64
65    barrier()
66
67    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
68
69    print "Processor Name:",get_processor_name()
70
71    # filenames
72#    meshes_dir_name = project.meshes_dir_name+'.msh'
73
74    # creates copy of code in output dir
75    print 'min triangles', project.trigs_min,
76    print 'Note: This is generally about 20% less than the final amount'
77
78    #--------------------------------------------------------------------------
79    # Create the triangular mesh based on overall clipping polygon with a
80    # tagged
81    # boundary and interior regions defined in project.py along with
82    # resolutions (maximal area of per triangle) for each polygon
83    #--------------------------------------------------------------------------
84   
85    if myid == 0:
86   
87        print 'start create mesh from regions'
88
89        create_mesh_from_regions(project.poly_all,
90                             boundary_tags=project.boundary_tags,
91                             maximum_triangle_area=project.res_poly_all,
92                             interior_regions=project.interior_regions,
93                             filename=project.meshes_dir_name+'.msh',
94                             use_cache=True,
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.meshes_dir_name+'.msh', use_cache=True, 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
162    #------------------------------------------------------------------------------
163    # Setup information for slide scenario (to be applied 1 min into simulation
164    #------------------------------------------------------------------------------
165
166    if scenario == 'slide':
167        # Function for submarine slide
168        from anuga.shallow_water.smf import slide_tsunami 
169        tsunami_source = slide_tsunami(length=35000.0,
170                                       depth=project.slide_depth,
171                                       slope=10.0,
172                                       thickness=800.0, 
173                                       x0=project.slide_origin[0], 
174                                       y0=project.slide_origin[1], 
175                                       alpha=0.0, 
176                                       domain=domain,
177                                       verbose=True)
178   
179    #-------------------------------------------------------------------------
180    # Setup boundary conditions
181    #-------------------------------------------------------------------------
182    print 'Available boundary tags', domain.get_boundary_tags()
183    print 'domain id', id(domain)
184    #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
185
186    Bt = Transmissive_boundary(domain)
187    Br = Reflective_boundary(domain)
188    Bd = Dirichlet_boundary([0,0,0])
189
190    if scenario == 'fixed_wave':
191    # FIXME (Ole): Change this to Transmissive Momentum as
192    # suggested by Rajaraman (ticket:208)
193        from math import sin, pi
194       
195        Bw = Time_boundary(domain = domain,
196                           f=lambda t: [(120<t<220)*3.0*sin(t*pi*2/10), 0, 0])
197        domain.set_boundary({'lagoon_west':Bt,
198                             'Land':Bt,
199                             'Reef_edge':Bw,
200                             'East&land':Bt})
201
202    # boundary conditions for slide scenario
203    if scenario == 'slide':
204        domain.set_boundary({'ocean':Bd,'land':Bd})
205
206    #----------------------------------------------------------------------------
207    # Evolve system through time
208    #--------------------------------------------------------------------
209    t0 = time.time()
210
211    for t in domain.evolve(kwargs['yieldstep'], finaltime = kwargs['finaltime']): 
212        domain.write_time()
213        domain.write_boundary_statistics(tags = 'Reef_edge')     
214
215       
216    x, y = domain.get_maximum_inundation_location()
217    q = domain.get_maximum_inundation_elevation()
218
219    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
220
221    print 'That took %.2f seconds' %(time.time()-t0)
222
223    #kwargs 'completed' must be added to write the final parameters to file
224    kwargs['completed']=str(time.time()-t0)
225   
226    if myid==0:
227        store_parameters(**kwargs)
228    barrier
229   
230    print 'memory usage before del domain1',mem_usage()
231
232def export_model(**kwargs):   
233    #store_parameters(**kwargs)
234   
235#    print 'memory usage before del domain',mem_usage()
236    #del domain
237    print 'memory usage after del domain',mem_usage()
238   
239    swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
240    print'swwfile',swwfile
241   
242    export_grid(swwfile, extra_name_out = 'FSM',
243            quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
244            #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
245            timestep = None,
246            reduction = max,
247            cellsize = 25,
248            NODATA_value = -1E-030,
249            easting_min = project.eastingmin,
250            easting_max = project.eastingmax,
251            northing_min = project.northingmin,
252            northing_max = project.northingmax,
253            verbose = False,
254            origin = None,
255            datum = 'WGS84',
256            format = 'txt')
257   
258#-------------------------------------------------------------
259if __name__ == "__main__":
260   
261    kwargs={}
262    kwargs['est_num_trigs']=project.trigs_min
263    kwargs['num_cpu']=numprocs
264    kwargs['host']=project.host
265    kwargs['res_factor']=project.res_factor
266    kwargs['starttime']=project.starttime
267    kwargs['yieldstep']=project.yieldstep
268    kwargs['midtime']=project.midtime
269    kwargs['finaltime']=project.finaltime
270    kwargs['output_dir']=project.output_run_time_dir
271    kwargs['bathy_file']=project.combined_dir_name+'.txt'
272#    kwargs['bathy_file']=project.combined_small_dir_name + '.pts'
273    kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
274    kwargs['file_name']=project.home+'detail.csv'
275    kwargs['aa_scenario_name']=project.scenario_name
276    kwargs['ab_time']=project.time
277    kwargs['res_factor']= project.res_factor
278    kwargs['tide']=project.tide
279    kwargs['user']=project.user
280    kwargs['alpha'] = project.alpha
281    kwargs['friction']=project.friction
282    kwargs['time_thinning'] = project.time_thinning
283    kwargs['dir_comment']=project.dir_comment
284    kwargs['export_cellsize']=project.export_cellsize
285   
286
287    run_model(**kwargs)
288     
289    if myid==0:
290        export_model(**kwargs)
291    barrier
Note: See TracBrowser for help on using the repository browser.