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

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

Cleaned up all of Perth Files

File size: 9.6 KB
Line 
1"""Script for running tsunami inundation scenario for Perth, 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 is combiled into a pts file through build_perth.py
9and a simulated tsunami is generated through an sts file from build_boundary.py
10
11Ole Nielsen and Duncan Gray, GA - 2005, Jane Sexton, Nick Bartzis, GA - 2006
12Ole Nielsen, Jane Sexton and Kristy Van Putten - 2008
13"""
14
15#------------------------------------------------------------------------------
16# Import necessary modules
17#------------------------------------------------------------------------------
18
19# Standard modules
20from os import sep
21import os
22from os.path import dirname, basename
23from os import mkdir, access, F_OK
24from shutil import copy
25import time
26import sys
27
28# Related major packages
29from anuga.shallow_water import Domain
30from anuga.shallow_water import Dirichlet_boundary
31from anuga.shallow_water import File_boundary
32from anuga.shallow_water import Reflective_boundary
33from anuga.shallow_water import Field_boundary
34from Numeric import allclose
35from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
36from anuga.pmesh.mesh_interface import create_mesh_from_regions
37from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
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 polygon import Polygon_function
45   
46# Application specific imports
47import project  # Definition of file names and polygons
48numprocs = 1
49myid = 0
50
51def run_model(**kwargs):
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
61    print 'output_dir',kwargs['output_dir']
62   
63    copy_code_files(kwargs['output_dir'],__file__, 
64             dirname(project.__file__)+sep+ project.__name__+'.py' )
65
66    store_parameters(**kwargs)
67
68    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
69
70    print "Processor Name:",get_processor_name()
71   
72    #-----------------------------------------------------------------------
73    # Domain definitions
74    #-----------------------------------------------------------------------
75
76    # Read in boundary from ordered sts file
77    urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir_event,project.scenario_name))
78
79    # Reading the landward defined points, this incorporates the original clipping
80    # polygon minus the 100m contour
81    landward_bounding_polygon = read_polygon(project.landward_dir)
82
83    # Combine sts polyline with landward points
84    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
85   
86    # counting segments
87    N = len(urs_bounding_polygon)-1
88
89    #boundary tags refer to project.landward 4 points equals 5 segments start at N
90    boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)}
91
92   
93    #--------------------------------------------------------------------------
94    # Create the triangular mesh based on overall clipping polygon with a tagged
95    # boundary and interior regions defined in project.py along with
96    # resolutions (maximal area of per triangle) for each polygon
97    #--------------------------------------------------------------------------
98
99    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
100    # causes problems with the ability to cache set quantity which takes alot of times
101       
102    print 'start create mesh from regions'
103
104    create_mesh_from_regions(bounding_polygon,
105                         boundary_tags=boundary_tags,
106                         maximum_triangle_area=project.res_poly_all,
107                         interior_regions=project.interior_regions,
108                         filename=project.meshes_dir_name,
109                         use_cache=True,
110                         verbose=True)
111   
112    #-------------------------------------------------------------------------
113    # Setup computational domain
114    #-------------------------------------------------------------------------
115    print 'Setup computational domain'
116
117    domain = Domain(project.meshes_dir_name, use_cache=False, verbose=True)
118    print 'memory usage before del domain',mem_usage()
119       
120    print domain.statistics()
121    print 'triangles',len(domain)
122   
123    kwargs['act_num_trigs']=len(domain)
124
125
126    #-------------------------------------------------------------------------
127    # Setup initial conditions
128    #-------------------------------------------------------------------------
129   print 'Setup initial conditions'
130
131    #following sets the stage/water to be offcoast only
132    IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'],
133                             geo_reference = domain.geo_reference)
134    domain.set_quantity('stage', IC)
135    #domain.set_quantity('stage',kwargs['tide'] )
136    domain.set_quantity('friction', kwargs['friction']) 
137   
138    print 'Start Set quantity',kwargs['elevation_file']
139
140    domain.set_quantity('elevation', 
141                        filename = kwargs['elevation_file'],
142                        use_cache = False,
143                        verbose = True,
144                        alpha = kwargs['alpha'])
145    print 'Finished Set quantity'
146
147##   #------------------------------------------------------
148##    # Distribute domain to implement parallelism !!!
149##    #------------------------------------------------------
150##
151##    if numprocs > 1:
152##        domain=distribute(domain)
153
154    #------------------------------------------------------
155    # Set domain parameters
156    #------------------------------------------------------
157    print 'domain id', id(domain)
158    domain.set_name(kwargs['scenario_name'])
159    domain.set_datadir(kwargs['output_dir'])
160    domain.set_default_order(2) # Apply second order scheme
161    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
162    domain.set_store_vertices_uniquely(False)
163    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
164    domain.tight_slope_limiters = 1
165    print 'domain id', id(domain)
166
167    #-------------------------------------------------------------------------
168    # Setup boundary conditions
169    #-------------------------------------------------------------------------
170    print 'Available boundary tags', domain.get_boundary_tags()
171    print 'domain id', id(domain)
172   
173    boundary_urs_out=project.boundaries_dir_event + sep + project.scenario_name
174
175    Br = Reflective_boundary(domain)
176    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
177   
178    print 'Available boundary tags', domain.get_boundary_tags()
179    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
180                   domain, mean_stage= project.tide,
181                   time_thinning=1,
182                   default_boundary=Bd,
183                   use_cache=True,
184                   verbose = True,
185                   boundary_polygon=bounding_polygon)
186
187    domain.set_boundary({'back': Br,
188                         'side': Bd,
189                         'ocean': Bf}) 
190
191    kwargs['input_start_time']=domain.starttime
192
193    print'finish set boundary'
194
195    #----------------------------------------------------------------------------
196    # Evolve system through time
197    #--------------------------------------------------------------------
198    t0 = time.time()
199
200    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
201                       ,skip_initial_step = False): 
202        domain.write_time()
203        domain.write_boundary_statistics(tags = 'ocean')
204
205    x, y = domain.get_maximum_inundation_location()
206    q = domain.get_maximum_inundation_elevation()
207
208    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
209
210    print 'That took %.2f seconds' %(time.time()-t0)
211
212    #kwargs 'completed' must be added to write the final parameters to file
213    kwargs['completed']=str(time.time()-t0)
214     
215    store_parameters(**kwargs)
216     
217    print 'memory usage before del domain1',mem_usage()
218   
219   
220#-------------------------------------------------------------
221if __name__ == "__main__":
222   
223    kwargs={}
224    kwargs['finaltime']=project.finaltime
225    kwargs['output_dir']=project.output_run_time_dir
226    kwargs['elevation_file']=project.combined_dir_name+'.pts'
227    kwargs['scenario_name']=project.scenario_name
228    kwargs['tide']=project.tide
229    kwargs['alpha'] = project.alpha
230    kwargs['friction']=project.friction
231    #kwargs['est_num_trigs']=project.trigs_min
232    #kwargs['num_cpu']=numprocs
233    #kwargs['host']=project.host
234    #kwargs['res_factor']=project.res_factor
235    #kwargs['starttime']=project.starttime
236    #kwargs['yieldstep']=project.yieldstep
237    #kwargs['file_name']=project.home+'detail.csv'
238    #kwargs['ab_time']=project.time
239    #kwargs['res_factor']= project.res_factor
240    #kwargs['user']=project.user
241    #kwargs['time_thinning'] = project.time_thinning
242    #kwargs['dir_comment']=project.dir_comment
243     
244
245    run_model(**kwargs)
246     
247   
Note: See TracBrowser for help on using the repository browser.