source: anuga_work/production/busselton/run_busselton.py @ 5575

Last change on this file since 5575 was 5575, checked in by kristy, 16 years ago

Addition of Polyline information with sts files

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