source: anuga_work/production/perth/run_perth.py @ 5557

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

change script for sts file (polyline method)

File size: 9.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 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
44numprocs = 1
45myid = 0
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.__file__)+sep+ project.__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.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    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
87    # causes problems with the ability to cache set quantity which takes alot of times
88    if myid == 0:
89   
90        print 'start create mesh from regions'
91
92        create_mesh_from_regions(project.poly_all,
93                             boundary_tags=project.boundary_tags,
94                             maximum_triangle_area=project.res_poly_all,
95                             interior_regions=project.interior_regions,
96                             filename=project.meshes_dir_name+'.msh',
97                             use_cache=True,
98                             verbose=True)
99   # barrier()
100   
101
102    #-------------------------------------------------------------------------
103    # Setup computational domain
104    #-------------------------------------------------------------------------
105    print 'Setup computational domain'
106
107    domain = Domain(project.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    #-------------------------------------------------------------------------
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, 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('friction', kwargs['friction']) 
130       
131        print 'Start Set quantity',kwargs['elevation_file']
132
133        domain.set_quantity('elevation', 
134                            filename = kwargs['elevation_file'],
135                            use_cache = False,
136                            verbose = True,
137                            alpha = kwargs['alpha'])
138        print 'Finished Set quantity'
139    #barrier()
140
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
163
164    #-------------------------------------------------------------------------
165    # Setup boundary conditions
166    #-------------------------------------------------------------------------
167    print 'Available boundary tags', domain.get_boundary_tags()
168    print 'domain id', id(domain)
169   
170
171    print 'Available boundary tags', domain.get_boundary_tags()
172    Bf = File_boundary(boundary_urs_out+'.sts', 
173                   domain, time_thinning=1,
174                   use_cache=True,
175                   verbose = True,
176                   boundary_polygon=project.bounding_polygon)
177
178    Br = Reflective_boundary(domain)
179    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
180
181    print 'start reading boundary file'
182
183##    Bf = Field_boundary(kwargs['boundary_file'],
184##                domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'],
185##                use_cache=False, verbose=True)
186
187    domain.set_boundary({'back': Br,
188                         'side': Bd,
189                         'ocean': Bf}) 
190
191    kwargs['input_start_time']=domain.starttime
192
193    print'finish set boundary'
194
195    #----------------------------------------------------------------------------
196    # Evolve system through time
197    #--------------------------------------------------------------------
198    t0 = time.time()
199
200    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
201                       ,skip_initial_step = False): 
202        domain.write_time()
203        domain.write_boundary_statistics(tags = 'ocean')   
204
205    x, y = domain.get_maximum_inundation_location()
206    q = domain.get_maximum_inundation_elevation()
207
208    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
209
210    print 'That took %.2f seconds' %(time.time()-t0)
211
212    #kwargs 'completed' must be added to write the final parameters to file
213    kwargs['completed']=str(time.time()-t0)
214   
215    if myid==0:
216        store_parameters(**kwargs)
217   # barrier
218   
219    print 'memory usage before del domain1',mem_usage()
220   
221   
222#-------------------------------------------------------------
223if __name__ == "__main__":
224   
225    kwargs={}
226    kwargs['est_num_trigs']=project.trigs_min
227    kwargs['num_cpu']=numprocs
228    kwargs['host']=project.host
229    kwargs['res_factor']=project.res_factor
230    kwargs['starttime']=project.starttime
231    kwargs['yieldstep']=project.yieldstep
232    kwargs['finaltime']=project.finaltime
233   
234    kwargs['output_dir']=project.output_run_time_dir
235    kwargs['elevation_file']=project.combined_dir_name+'.pts'
236#    kwargs['elevation_file']=project.combined_small_dir_name + '.pts'
237    kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
238    kwargs['file_name']=project.home+'detail.csv'
239    kwargs['aa_scenario_name']=project.scenario_name
240    kwargs['ab_time']=project.time
241    kwargs['res_factor']= project.res_factor
242    kwargs['tide']=project.tide
243    kwargs['user']=project.user
244    kwargs['alpha'] = project.alpha
245    kwargs['friction']=project.friction
246    kwargs['time_thinning'] = project.time_thinning
247    kwargs['dir_comment']=project.dir_comment
248    kwargs['export_cellsize']=project.export_cellsize
249   
250
251    run_model(**kwargs)
252     
253    if myid==0:
254        export_model(**kwargs)
255    #barrier
Note: See TracBrowser for help on using the repository browser.