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

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

Updated evolve in run file so that a tide wave takes over when the sts file finishes

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