source: anuga_work/production/onslow_2008/run_onslow.py @ 5635

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

To re-run Onslow with new boundary input

File size: 11.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
43from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
44
45# Application specific imports
46import project                 # Definition of file names and polygons
47numprocs = 1
48myid = 0
49
50def run_model(**kwargs):
51   
52
53    #------------------------------------------------------------------------------
54    # Copy scripts to time stamped output directory and capture screen
55    # output to file
56    #------------------------------------------------------------------------------
57    print "Processor Name:",get_processor_name()
58
59    #copy script must be before screen_catcher
60    #print kwargs
61
62    print 'output_dir',kwargs['output_dir']
63    if myid == 0:
64        copy_code_files(kwargs['output_dir'],__file__, 
65                 dirname(project.__file__)+sep+ project.__name__+'.py' )
66
67        store_parameters(**kwargs)
68
69   # barrier()
70
71    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
72
73    print "Processor Name:",get_processor_name()
74
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.polygons_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], 'side': [N,N+4],'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=True,
114                             verbose=True)
115   # barrier()
116
117##        covariance_value,alpha = find_optimal_smoothing_parameter (data_file= kwargs['elevation_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['elevation_file']
152
153        domain.set_quantity('elevation', 
154                            filename = kwargs['elevation_file'],
155                            use_cache = False,
156                            verbose = True,
157                            alpha = kwargs['alpha'])
158        print 'Finished Set quantity'
159    #barrier()
160
161##    #------------------------------------------------------
162##    # Create x,y,z file of mesh vertex!!!
163##    #------------------------------------------------------
164##        coord = domain.get_vertex_coordinates()
165##        depth = domain.get_quantity('elevation')
166##       
167##        # Write vertex coordinates to file
168##        filename=project.vertex_filename
169##        fid=open(filename,'w')
170##        fid.write('x (m), y (m), z(m)\n')
171##        for i in range(len(coord)):
172##            pt=coord[i]
173##            x=pt[0]
174##            y=pt[1]
175##            z=depth[i]
176##            fid.write('%.6f,%.6f,%.6f\n' %(x, y, z))
177##
178
179    #------------------------------------------------------
180    # Distribute domain to implement parallelism !!!
181    #------------------------------------------------------
182
183    if numprocs > 1:
184        domain=distribute(domain)
185
186    #------------------------------------------------------
187    # Set domain parameters
188    #------------------------------------------------------
189    print 'domain id', id(domain)
190    domain.set_name(kwargs['aa_scenario_name'])
191    domain.set_datadir(kwargs['output_dir'])
192    domain.set_default_order(2) # Apply second order scheme
193    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
194    domain.set_store_vertices_uniquely(False)
195    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
196    domain.tight_slope_limiters = 1
197    print 'domain id', id(domain)
198
199    #-------------------------------------------------------------------------
200    # Setup boundary conditions
201    #-------------------------------------------------------------------------
202    print 'Available boundary tags', domain.get_boundary_tags()
203    print 'domain id', id(domain)
204   
205    boundary_urs_out=project.boundaries_dir_name
206   
207    print 'Available boundary tags', domain.get_boundary_tags()
208    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
209                   domain, mean_stage= project.tide,
210                   time_thinning=1,
211                   use_cache=True,
212                   verbose = True,
213                   boundary_polygon=bounding_polygon)
214
215    Br = Reflective_boundary(domain)
216    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
217
218    print dir(Bf)
219    print 'start reading boundary file'
220   
221
222##    Bf = Field_boundary(kwargs['boundary_file'],
223##                domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'],
224##                use_cache=False, verbose=True)
225
226    domain.set_boundary({'back': Br,
227                         'side': Bd,
228                         'ocean': Bf}) 
229
230    kwargs['input_start_time']=domain.starttime
231
232    print'finish set boundary'
233
234    #----------------------------------------------------------------------------
235    # Evolve system through time
236    #--------------------------------------------------------------------
237    t0 = time.time()
238
239    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
240                       ,skip_initial_step = False): 
241        domain.write_time()
242        domain.write_boundary_statistics(tags = 'ocean')
243           
244    x, y = domain.get_maximum_inundation_location()
245    q = domain.get_maximum_inundation_elevation()
246
247    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
248
249    print 'That took %.2f seconds' %(time.time()-t0)
250
251    #kwargs 'completed' must be added to write the final parameters to file
252    kwargs['completed']=str(time.time()-t0)
253   
254    if myid==0:
255        store_parameters(**kwargs)
256   # barrier
257   
258    print 'memory usage before del domain1',mem_usage()
259   
260   
261#-------------------------------------------------------------
262if __name__ == "__main__":
263   
264    kwargs={}
265    kwargs['est_num_trigs']=project.trigs_min
266    kwargs['num_cpu']=numprocs
267    kwargs['host']=project.host
268    kwargs['res_factor']=project.res_factor
269    kwargs['starttime']=project.starttime
270    kwargs['yieldstep']=project.yieldstep
271    kwargs['finaltime']=project.finaltime
272   
273    kwargs['output_dir']=project.output_run_time_dir
274    kwargs['elevation_file']=project.combined_dir_name+'.pts'
275    kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
276    kwargs['file_name']=project.home+'detail.csv'
277    kwargs['aa_scenario_name']=project.scenario_name
278    kwargs['ab_time']=project.time
279    kwargs['res_factor']= project.res_factor
280    kwargs['tide']=project.tide
281    kwargs['user']=project.user
282    kwargs['alpha'] = project.alpha
283    kwargs['friction']=project.friction
284    kwargs['time_thinning'] = project.time_thinning
285    kwargs['dir_comment']=project.dir_comment
286    kwargs['export_cellsize']=project.export_cellsize
287   
288
289    run_model(**kwargs)
290     
291    if myid==0:
292        export_model(**kwargs)
293    #barrier
Note: See TracBrowser for help on using the repository browser.