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

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

update to dampier

File size: 9.2 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 submarine landslide.
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
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
31
32from anuga.pmesh.mesh_interface import create_mesh_from_regions
33
34from anuga.geospatial_data.geospatial_data import *
35from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
36from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
37
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
47# filenames
48
49#build_time = '20061029_231935_build_tide_24'
50#build_time = '20061030_165746_build_tide_24'
51#build_time = '20061102_215532_build_plus_old_data'
52#build_time = '20061103_055258_build'
53build_time = '20061107_063840_build'
54#build_time = '20061025_153643_build_basic'
55bound_time = '20061102_221245_build'
56
57boundaries_name = project.boundaries_name
58#meshes_time_dir_name = project.meshes_time_dir_name+'.msh'
59meshes_dir_name = project.meshes_dir_name+'.msh'
60#source_dir = project.boundarydir
61#boundaries_time_dir_name = project.boundaries_dir + build_time + sep + boundaries_name
62boundaries_time_dir_name = project.boundaries_time_dir_name
63boundaries_dir_name = project.boundaries_dir_name
64tide = project.tide
65
66# creates copy of code in output dir
67if myid == 0:
68    copy_code_files(project.output_run_time_dir,__file__, 
69                 dirname(project.__file__)+sep+ project.__name__+'.py' )
70barrier()
71start_screen_catcher(project.output_run_time_dir, myid, numprocs)
72
73print 'USER: ', project.user
74#sys.exit()
75#--------------------------------------------------------------------------
76# Create the triangular mesh based on overall clipping polygon with a
77# tagged
78# boundary and interior regions defined in project.py along with
79# resolutions (maximal area of per triangle) for each polygon
80#--------------------------------------------------------------------------
81
82if myid == 0:
83#    if access(project.meshes_time_dir,F_OK) == 0:
84#        mkdir(project.meshes_time_dir)
85    if access(project.meshes_dir,F_OK) == 0:
86        mkdir(project.meshes_dir)
87    print 'start create mesh from regions'
88    interior_regions = [#[project.karratha_polygon, 25000],
89#                    [project.cipma_polygon, 1000],
90                    [project.poly_pipeline, 5000],
91                    [project.poly_facility, 500]]   
92#    meshes_dir_name = project.meshes_dir_name + '.msh'
93
94    create_mesh_from_regions(project.bounding_polygon,
95                         boundary_tags={'back': [7, 8], 'side': [0, 6],
96                                        'ocean': [1, 2, 3, 4, 5]},
97                         maximum_triangle_area=100000,
98                         interior_regions=interior_regions,
99#                         filename=meshes_time_dir_name,
100                         filename=meshes_dir_name,
101                         use_cache=True,
102                         verbose=True)
103
104# to sync all processors are ready
105barrier()
106
107#-------------------------------------------------------------------------
108# Setup computational domain
109#-------------------------------------------------------------------------
110print 'Setup computational domain'
111#domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True)
112domain = Domain(meshes_dir_name, use_cache=True, verbose=True)
113print domain.statistics()
114
115
116
117print 'starting to create boundary conditions'
118boundaries_in_dir_name = project.boundaries_in_dir_name
119
120from anuga.shallow_water.data_manager import urs2sww, ferret2sww
121
122print 'maxlat=project.south_boundary, minlat=project.north_boundary', project.south_boundary,project.north_boundary
123print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
124print ' maxlon=project.east',project.east
125
126print 'origin: domain.geo_reference.get_origin()',domain.geo_reference.get_origin()
127
128#import sys; sys.exit()
129
130#if access(project.boundaries_time_dir,F_OK) == 0:
131#    mkdir (project.boundaries_time_dir)
132# put above distribute
133print 'boundary file is: ',boundaries_dir_name
134from caching import cache
135if myid == 0:
136    cache(urs2sww,
137          (boundaries_in_dir_name,
138    #       boundaries_time_dir_name),
139           boundaries_dir_name), 
140          {'verbose': True,
141           'minlat': project.south_boundary,
142           'maxlat': project.north_boundary,
143           'minlon': project.west_boundary,
144           'maxlon': project.east_boundary,
145#           'minlat': project.south,
146#           'maxlat': project.north,
147#           'minlon': project.west,
148#           'maxlon': project.east,
149           'mint': 0, 'maxt': 35000,
150           'origin': domain.geo_reference.get_origin(),
151           'mean_stage': project.tide,
152#           'zscale': 1,                 #Enhance tsunami
153           'fail_on_NaN': False},
154           verbose = True,
155           )
156barrier()
157
158
159#-------------------------------------------------------------------------
160# Setup initial conditions
161#-------------------------------------------------------------------------
162print 'Setup initial conditions'
163
164domain.set_quantity('stage', tide)
165domain.set_quantity('friction', 0.0) 
166#combined_time_dir_name = project.topographies_dir+build_time+project.combined_name
167print 'Start Set quantity'
168
169domain.set_quantity('elevation', 
170                    filename = project.topographies_dir + build_time + sep + project.combined_final_name + '.pts',
171                    use_cache = True,
172                    verbose = True,
173                    alpha = 0.1)
174#domain.set_quantity('elevation', -50)
175print 'Finished Set quantity'
176
177
178
179#------------------------------------------------------
180# Distribute domain to implement parallelism !!!
181#------------------------------------------------------
182
183if numprocs > 1:
184    domain=distribute(domain)
185
186#------------------------------------------------------
187# Set domain parameters
188#------------------------------------------------------
189
190domain.set_name(project.scenario_name)
191domain.set_datadir(project.output_run_time_dir)
192domain.set_default_order(2) # Apply second order scheme
193domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
194domain.set_store_vertices_uniquely(False)
195domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
196domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
197
198#-------------------------------------------------------------------------
199# Setup boundary conditions
200#-------------------------------------------------------------------------
201print 'Available boundary tags', domain.get_boundary_tags()
202
203print 'Reading Boundary file'
204#boundariesname = project.boundaries_dir + '20061101_003322_build'+sep+boundaries_name
205#print'boundariesname',boundariesname
206#Bf = File_boundary(boundaries_time_dir_name + '.sww',
207
208#Bf = File_boundary(boundariesname + '.sww',
209
210print 'domain id', id(domain)
211Bf = File_boundary(boundaries_dir_name + '.sww',
212                  domain, time_thinning=5, use_cache=True, verbose=True)
213
214print 'finished reading boundary file'
215
216Br = Reflective_boundary(domain)
217Bd = Dirichlet_boundary([tide,0,0])
218
219print'set_boundary'
220domain.set_boundary({'back': Br,
221                     'side': Bd,
222                     'ocean': Bf}) 
223print'finish set boundary'
224
225#----------------------------------------------------------------------------
226# Evolve system through time
227#----------------------------------------------------------------------------
228
229t0 = time.time()
230
231#for t in domain.evolve(yieldstep = 60, finaltime = 34000):
232#    domain.write_time()
233#    domain.write_boundary_statistics(tags = 'ocean')
234
235for t in domain.evolve(yieldstep = 120, finaltime = 9000):
236#for t in domain.evolve(yieldstep = 120, finaltime = 130):
237    domain.write_time()
238    domain.write_boundary_statistics(tags = 'ocean')
239    print 'time: ',time.time()
240     
241for t in domain.evolve(yieldstep = 60, finaltime = 28800
242                       ,skip_initial_step = True):
243    domain.write_time()
244    domain.write_boundary_statistics(tags = 'ocean')   
245
246for t in domain.evolve(yieldstep = 120, finaltime = 34800
247                       ,skip_initial_step = True): 
248    domain.write_time()
249    domain.write_boundary_statistics(tags = 'ocean')   
250   
251print 'That took %.2f seconds' %(time.time()-t0)
252
Note: See TracBrowser for help on using the repository browser.