source: anuga_work/production/perth/run_perth_250m.py @ 6540

Last change on this file since 6540 was 6032, checked in by kristy, 15 years ago

general maintenance of scripts

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