source: anuga_work/production/onslow_2008/run_onslow.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: 11.4 KB
RevLine 
[5635]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
[5645]44from Scientific.IO.NetCDF import NetCDFFile
[5635]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
[5645]85    landward_bounding_polygon = read_polygon(project.boundaries_dir+'landward_bounding_polygon.txt')
[5635]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
[5645]218    fid = NetCDFFile(boundary_urs_out+'.sts', 'r')    #Open existing file for read
219    sts_time=fid.variables['time'][:]+fid.starttime
220    tmin=min(sts_time)
221    tmax=max(sts_time)
222    fid.close()
[5635]223
[5645]224    print 'Boundary end time ', tmax-tmin
225     
[5635]226##    Bf = Field_boundary(kwargs['boundary_file'],
227##                domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'],
228##                use_cache=False, verbose=True)
229
230    domain.set_boundary({'back': Br,
231                         'side': Bd,
232                         'ocean': Bf}) 
233
234    kwargs['input_start_time']=domain.starttime
235
236    print'finish set boundary'
237
238    #----------------------------------------------------------------------------
239    # Evolve system through time
240    #--------------------------------------------------------------------
241    t0 = time.time()
242
243    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
244                       ,skip_initial_step = False): 
245        domain.write_time()
246        domain.write_boundary_statistics(tags = 'ocean')
[5645]247
248        if t >= tmax-tmin:
249            print 'changed to tide boundary condition at ocean'
250            domain.set_boundary({'ocean': Bd})
[5635]251           
252    x, y = domain.get_maximum_inundation_location()
253    q = domain.get_maximum_inundation_elevation()
254
255    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
256
257    print 'That took %.2f seconds' %(time.time()-t0)
258
259    #kwargs 'completed' must be added to write the final parameters to file
260    kwargs['completed']=str(time.time()-t0)
261   
262    if myid==0:
263        store_parameters(**kwargs)
264   # barrier
265   
266    print 'memory usage before del domain1',mem_usage()
267   
268   
269#-------------------------------------------------------------
270if __name__ == "__main__":
271   
272    kwargs={}
273    kwargs['est_num_trigs']=project.trigs_min
274    kwargs['num_cpu']=numprocs
275    kwargs['host']=project.host
276    kwargs['res_factor']=project.res_factor
277    kwargs['starttime']=project.starttime
278    kwargs['yieldstep']=project.yieldstep
279    kwargs['finaltime']=project.finaltime
280   
281    kwargs['output_dir']=project.output_run_time_dir
282    kwargs['elevation_file']=project.combined_dir_name+'.pts'
283    kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
284    kwargs['file_name']=project.home+'detail.csv'
285    kwargs['aa_scenario_name']=project.scenario_name
286    kwargs['ab_time']=project.time
287    kwargs['res_factor']= project.res_factor
288    kwargs['tide']=project.tide
289    kwargs['user']=project.user
290    kwargs['alpha'] = project.alpha
291    kwargs['friction']=project.friction
292    kwargs['time_thinning'] = project.time_thinning
293    kwargs['dir_comment']=project.dir_comment
294    kwargs['export_cellsize']=project.export_cellsize
295   
296
297    run_model(**kwargs)
298     
299    if myid==0:
300        export_model(**kwargs)
301    #barrier
Note: See TracBrowser for help on using the repository browser.