source: anuga_work/production/dampier_2006/run_dampier.py @ 4430

Last change on this file since 4430 was 4430, checked in by nick, 17 years ago

update dampier

File size: 7.8 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
19from os.path import dirname, basename
20from os import mkdir, access, F_OK
21from shutil import copy
22import time
23import sys
24
25# Related major packages
26from anuga.shallow_water import Domain
27from anuga.shallow_water import Dirichlet_boundary
28from anuga.shallow_water import File_boundary
29from anuga.shallow_water import Reflective_boundary
30from anuga.shallow_water import Field_boundary
31from Numeric import allclose
32
33from anuga.pmesh.mesh_interface import create_mesh_from_regions
34from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
35from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
36from anuga_parallel.parallel_abstraction import get_processor_name
37from anuga.caching import myhash
38# Application specific imports
39import project                 # Definition of file names and polygons
40
41#------------------------------------------------------------------------------
42# Copy scripts to time stamped output directory and capture screen
43# output to file
44#------------------------------------------------------------------------------
45
46#copy script must be before screen_catcher
47if myid == 0:
48    copy_code_files(project.output_run_time_dir,__file__, 
49                 dirname(project.__file__)+sep+ project.__name__+'.py' )
50barrier()
51
52start_screen_catcher(project.output_run_time_dir, myid, numprocs)
53
54print "Processor Name:",get_processor_name()
55barrier()
56
57# filenames
58#boundaries_name = project.boundaries_name
59meshes_dir_name = project.meshes_dir_name+'.msh'
60#boundaries_dir_name = project.boundaries_dir_name
61
62tide = project.tide
63
64# creates copy of code in output dir
65
66
67print 'USER: ', project.user
68print 'min triangles', project.trigs_min,
69print 'Note: This is generally about 20% less than the final amount'
70
71#--------------------------------------------------------------------------
72# Create the triangular mesh based on overall clipping polygon with a
73# tagged
74# boundary and interior regions defined in project.py along with
75# resolutions (maximal area of per triangle) for each polygon
76#--------------------------------------------------------------------------
77
78if myid == 0:
79   
80    print 'start create mesh from regions'
81
82    create_mesh_from_regions(project.poly_all,
83                             boundary_tags={'back': [2,3], 'side': [0, 1, 4],
84                                            'ocean': [5]},
85                             maximum_triangle_area=project.res_poly_all,
86                             interior_regions=project.interior_regions,
87                             filename=meshes_dir_name,
88                             use_cache=False,
89                             verbose=True)
90
91# to sync all processors are ready
92barrier()
93
94#-------------------------------------------------------------------------
95# Setup computational domain
96#-------------------------------------------------------------------------
97print 'Setup computational domain'
98#from caching import cache
99
100#domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
101#above don't work
102domain = Domain(meshes_dir_name, use_cache=False, verbose=True)
103print 'domain id', id(domain) 
104print 'myhash', myhash(domain)     
105       
106print domain.statistics()
107
108#boundaries_dir_name=project.boundaries_dir_name
109
110#print 'starting to create boundary conditions'
111
112#from anuga.shallow_water.data_manager import urs2sww
113
114
115#-------------------------------------------------------------------------
116# Setup initial conditions
117#-------------------------------------------------------------------------
118if myid == 0:
119
120    print 'Setup initial conditions'
121
122    from polygon import *
123    #following sets the stage/water to be offcoast only
124    IC = Polygon_function( [(project.poly_mainland, 0.)], default = tide)
125#    IC = Polygon_function( [(project.poly_mainland, 0.)], default = 200)
126    domain.set_quantity('stage', IC)
127    domain.set_quantity('friction', 0.01) 
128    print 'Start Set quantity'
129
130    domain.set_quantity('elevation', 
131#                    filename = project.combined_dir_name + '.pts',
132# MUST USE TXT FILES FOR CACHING TO WORK!
133                    filename = project.combined_dir_name + '.txt',
134#                    filename = project.combined_smallest_dir_name + '.txt',
135                    use_cache = True,
136                    verbose = True,
137                    alpha = 0.1)
138    print 'Finished Set quantity'
139barrier()
140
141
142#------------------------------------------------------
143# Distribute domain to implement parallelism !!!
144#------------------------------------------------------
145
146if numprocs > 1:
147    domain=distribute(domain)
148
149#------------------------------------------------------
150# Set domain parameters
151#------------------------------------------------------
152print 'domain id', id(domain)
153domain.set_name(project.scenario_name)
154domain.set_datadir(project.output_run_time_dir)
155domain.set_default_order(2) # Apply second order scheme
156domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
157domain.set_store_vertices_uniquely(False)
158domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
159domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
160print 'domain id', id(domain)
161domain.beta_h = 0
162#domain.limit2007 = 1
163
164#-------------------------------------------------------------------------
165# Setup boundary conditions
166#-------------------------------------------------------------------------
167print 'Available boundary tags', domain.get_boundary_tags()
168print 'domain id', id(domain)
169#print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
170
171Bf = Field_boundary(project.boundaries_dir_namea + '.sww',
172                    domain, time_thinning=12, mean_stage=tide, 
173                    use_cache=True, verbose=True)
174
175print 'finished reading boundary file'
176
177Br = Reflective_boundary(domain)
178Bd = Dirichlet_boundary([tide,0,0])
179
180print'set_boundary'
181##domain.set_boundary({'back': Br,
182##                     'side': Bf,
183##                     'ocean': Bf})
184domain.set_boundary({'back': Br,
185                     'side': Bd,
186                     'ocean': Bf}) 
187print'finish set boundary'
188
189#----------------------------------------------------------------------------
190# Evolve system through time
191#----------------------------------------------------------------------------
192
193t0 = time.time()
194
195for t in domain.evolve(yieldstep = 60, finaltime = 29950):
196    domain.write_time()
197#    domain.write_time(track_speeds=True)
198    domain.write_boundary_statistics(tags = 'ocean')
199   
200#    domain.write_time(track_speeds=True)
201
202#for t in domain.evolve(yieldstep = 120, finaltime = 9000):
203#    domain.write_time()
204#    domain.write_boundary_statistics(tags = 'ocean')
205'''     
206for t in domain.evolve(yieldstep = 60, finaltime = 28800
207                       ,skip_initial_step = True):
208    domain.write_time()
209    domain.write_boundary_statistics(tags = 'ocean')   
210
211for t in domain.evolve(yieldstep = 120, finaltime = 34800
212                       ,skip_initial_step = True):
213    domain.write_time()
214    domain.write_boundary_statistics(tags = 'ocean')   
215'''
216x, y = domain.get_maximum_inundation_location()
217q = domain.get_maximum_inundation_elevation()
218
219print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
220
221print 'That took %.2f seconds' %(time.time()-t0)
222
Note: See TracBrowser for help on using the repository browser.