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

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

updated sts boundary

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