source: anuga_work/production/australia_ph2/dampier_ph2/run_dampier_ph2_250m.py @ 6200

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

New model for phase 2 work

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