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

Last change on this file since 4058 was 4049, checked in by nick, 18 years ago

updates 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
31from Numeric import allclose
32
33from anuga.pmesh.mesh_interface import create_mesh_from_regions
34
35from anuga.geospatial_data.geospatial_data import *
36from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
37from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
38
39# Application specific imports
40import project                 # Definition of file names and polygons
41
42#------------------------------------------------------------------------------
43# Copy scripts to time stamped output directory and capture screen
44# output to file
45#------------------------------------------------------------------------------
46
47
48# filenames
49
50#build_time = '20061029_231935_build_tide_24'
51#build_time = '20061030_165746_build_tide_24'
52#build_time = '20061102_215532_build_plus_old_data'
53#build_time = '20061103_055258_build'
54build_time = '20061107_063840_build'
55#build_time = '20061025_153643_build_basic'
56bound_time = '20061102_221245_build'
57
58boundaries_name = project.boundaries_name
59#meshes_time_dir_name = project.meshes_time_dir_name+'.msh'
60meshes_dir_name = project.meshes_dir_name+'.msh'
61#source_dir = project.boundarydir
62#boundaries_time_dir_name = project.boundaries_dir + build_time + sep + boundaries_name
63boundaries_time_dir_name = project.boundaries_time_dir_name
64boundaries_dir_name = project.boundaries_dir_name
65tide = project.tide
66
67# creates copy of code in output dir
68if myid == 0:
69    copy_code_files(project.output_run_time_dir,__file__, 
70                 dirname(project.__file__)+sep+ project.__name__+'.py' )
71barrier()
72start_screen_catcher(project.output_run_time_dir, myid, numprocs)
73
74print 'USER: ', project.user
75#sys.exit()
76#--------------------------------------------------------------------------
77# Create the triangular mesh based on overall clipping polygon with a
78# tagged
79# boundary and interior regions defined in project.py along with
80# resolutions (maximal area of per triangle) for each polygon
81#--------------------------------------------------------------------------
82
83if myid == 0:
84   
85    print 'start create mesh from regions'
86    interior_regions = [#[project.karratha_polygon, 25000],
87                    [project.poly_coast, 10000],
88                    [project.poly_pipeline, 2000],
89                    [project.poly_facility, 500]]   
90#                    [project.poly_interior, 1000]]       
91   
92
93    create_mesh_from_regions(project.bounding_polygon,
94                         boundary_tags={'back': [7, 8], 'side': [0, 6],
95                                        'ocean': [1, 2, 3, 4, 5]},
96                         maximum_triangle_area=100000,
97                         interior_regions=interior_regions,
98#                         filename=meshes_time_dir_name,
99                         filename=meshes_dir_name,
100                         use_cache=True,
101                         verbose=True)
102
103# to sync all processors are ready
104barrier()
105
106#-------------------------------------------------------------------------
107# Setup computational domain
108#-------------------------------------------------------------------------
109print 'Setup computational domain'
110#domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True)
111domain = Domain(meshes_dir_name, use_cache=True, verbose=True)
112print domain.statistics()
113
114print 'starting to create boundary conditions'
115boundaries_in_dir_name = project.boundaries_in_dir_name
116
117from anuga.shallow_water.data_manager import urs2sww, ferret2sww
118
119print 'maxlat=project.south_boundary, minlat=project.north_boundary', project.south_boundary,project.north_boundary
120print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
121print ' maxlon=project.east',project.east
122
123print 'origin: domain.geo_reference.get_origin()',domain.geo_reference.get_origin()
124
125#import sys; sys.exit()
126
127#if access(project.boundaries_time_dir,F_OK) == 0:
128#    mkdir (project.boundaries_time_dir)
129# put above distribute
130print 'boundary file is: ',boundaries_dir_name
131from caching import cache
132if myid == 0:
133    cache(ferret2sww,
134          (boundaries_in_dir_name,
135    #       boundaries_time_dir_name),
136           boundaries_dir_name), 
137          {'verbose': True,
138           'minlat': project.south_boundary,
139           'maxlat': project.north_boundary,
140           'minlon': project.west_boundary,
141           'maxlon': project.east_boundary,
142#           'minlat': project.south,
143#           'maxlat': project.north,
144#           'minlon': project.west,
145#           'maxlon': project.east,
146           'mint': 0, 'maxt': 35100,
147           'origin': domain.geo_reference.get_origin(),
148           'mean_stage': project.tide,
149#           'zscale': 1,                 #Enhance tsunami
150           'fail_on_NaN': False},
151           verbose = True,
152           )
153barrier()
154
155
156#-------------------------------------------------------------------------
157# Setup initial conditions
158#-------------------------------------------------------------------------
159if myid == 0:
160
161    print 'Setup initial conditions'
162
163    domain.set_quantity('stage', tide)
164    domain.set_quantity('friction', 0.0) 
165    #combined_time_dir_name = project.topographies_dir+build_time+project.combined_name
166    print 'Start Set quantity'
167
168
169    domain.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)
175    print 'Finished Set quantity'
176barrier()
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    domain.write_time()
237    domain.write_boundary_statistics(tags = 'ocean')
238#    print 'time: ',time.time()
239 
240    if allclose(t, 6000):
241        domain.set_quantity('xmomentum', 0)
242        domain.set_quantity('ymomentum', 0)
243#import sys; sys.exit()
244     
245for t in domain.evolve(yieldstep = 60, finaltime = 28800
246                       ,skip_initial_step = True):
247    domain.write_time()
248    domain.write_boundary_statistics(tags = 'ocean')   
249
250for t in domain.evolve(yieldstep = 120, finaltime = 34800
251                       ,skip_initial_step = True): 
252    domain.write_time()
253    domain.write_boundary_statistics(tags = 'ocean')   
254   
255print 'That took %.2f seconds' %(time.time()-t0)
256
Note: See TracBrowser for help on using the repository browser.