source: anuga_work/production/dampier_2006/run_dampier_simple.py @ 5442

Last change on this file since 5442 was 5442, checked in by ole, 16 years ago

Retired h-limiter and beta_h as per ticket:194.
All unit tests and validation tests pass.

File size: 10.8 KB
RevLine 
[4307]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
[4357]19from os.path import dirname, basename,
[4307]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 Numeric import allclose
31
32from anuga.pmesh.mesh_interface import create_mesh_from_regions
33from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
34from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
35from anuga_parallel.parallel_abstraction import get_processor_name
36from anuga.caching import myhash
37# Application specific imports
[4311]38import project                 # Definition of file names and polygons
[4307]39
40#------------------------------------------------------------------------------
41# Copy scripts to time stamped output directory and capture screen
42# output to file
43#------------------------------------------------------------------------------
44
[4311]45start_screen_catcher(project.output_run_time_dir, myid, numprocs)
[4307]46print "Processor Name:",get_processor_name()
47
48# filenames
49#boundaries_name = project.boundaries_name
50meshes_dir_name = 'small.tsh' # this will be local
51#meshes_dir_name = project.meshes_dir_name+'.msh'
52#boundaries_dir_name = project.boundaries_dir_name
53
[4311]54tide = project.tide
[4307]55
56# creates copy of code in output dir
57if myid == 0:
58    copy_code_files(project.output_run_time_dir,__file__, 
59                 dirname(project.__file__)+sep+ project.__name__+'.py' )
60barrier()
[4311]61
62print 'USER: ', project.user
63print 'min triangles', project.trigs_min,
[4307]64print 'Note: This is generally about 20% less than the final amount'
65
66#--------------------------------------------------------------------------
67# Create the triangular mesh based on overall clipping polygon with a
68# tagged
69# boundary and interior regions defined in project.py along with
70# resolutions (maximal area of per triangle) for each polygon
71#--------------------------------------------------------------------------
72'''
73poly = [[0,0],[0,100],[100,100],[100,0]]
74
75create_mesh_from_regions(poly,
76                             boundary_tags={'back': [0], 'side': [1,3],
77                                            'ocean': [2]},
78                             maximum_triangle_area=1,
79                             interior_regions=None,
80                             filename=meshes_dir_name,
81                             use_cache=True,
82                             verbose=True)
83
84sys.exit()
85
86if myid == 0:
87   
88    print 'start create mesh from regions'
89#    print 'project.res_poly_all',project.res_poly_all
90#    print 'project.interior_regions_test',project.interior_regions_test
91#    print 'meshes_dir_name',meshes_dir_name
92    create_mesh_from_regions(project.poly_all,
93                             boundary_tags={'back': [2,3], 'side': [0, 1, 4],
94                                            'ocean': [5]},
95                             maximum_triangle_area=project.res_poly_all,
96                             interior_regions=project.interior_regions_test,
97                             filename=meshes_dir_name,
98                             use_cache=True,
99                             verbose=True)
100
101# to sync all processors are ready
102barrier()
103
104#all= [[517303.61,7760858.88],[517477.06,7733020.49],[425462.92,7733020.49],[425462.92,7760338.54]]
105#all= [[500000,7760800],[500000,7733000],[450000,7733000],[450000,7760800]]
106#all= [[469000,7760800],[490000,7750000],[459000,7750000]]
107all= [[479000,7760800],[469000,7750000],[459000,7760800]] #swapped edge of triangle
108#all= [[469000,7760800],[490000,7750000],[459000,7750000]]
109create_mesh_from_regions(all,
110                             boundary_tags={'ocean': [ 2],'side': [0, 1]},
111#                             boundary_tags={'back': [2,3], 'side': [0, 1, 4],'ocean': [5]},
112                             maximum_triangle_area=50000,
113#                             interior_regions=project.interior_regions_test,
114                             filename=meshes_dir_name,
115                             use_cache=False,
116                             verbose=True)
117'''                           
[4311]118#all = [[469000,7760000],[470000,7758000],[468000,7758000]]
119all = [[470000,7760000],[469000,7758000],[468000,7760000]]
[4307]120create_mesh_from_regions(all,
[4311]121                             boundary_tags={'ocean': [0,1,2]},
122#                             boundary_tags={'ocean': [ 2],'side': [0, 1]},
123                             maximum_triangle_area=50000,
[4307]124                             filename=meshes_dir_name,
[4311]125                             use_cache=False,
[4307]126                             verbose=True)
[4311]127
[4307]128#-------------------------------------------------------------------------
129# Setup computational domain
130#-------------------------------------------------------------------------
131print 'Setup computational domain'
132domain = Domain( meshes_dir_name, verbose=True)
133'''
134from caching import cache
135domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
136print 'domain id', id(domain)
137print 'myhash', myhash(domain)     
138'''       
139print domain.statistics()
140
141boundaries_dir_name=project.boundaries_dir_name
142
143print 'starting to create boundary conditions'
144
145from anuga.shallow_water.data_manager import urs2sww
146'''
147# put above distribute
148print 'boundary file is: ',project.boundaries_dir_name
149from caching import cache
150if myid == 0:
151    cache(urs2sww,
152          (project.boundaries_in_dir_name,
153           project.boundaries_dir_name),
154          {'verbose': True,
155           'minlat': project.south_boundary,
156           'maxlat': project.north_boundary,
157           'minlon': project.west_boundary,
158           'maxlon': project.east_boundary,
159           'mint': 0, 'maxt': 35100,
160           'origin': domain.geo_reference.get_origin(),
161           'mean_stage': project.tide,
162#           'zscale': 1,                 #Enhance tsunami
163           'fail_on_NaN': False},
164           verbose = True,
165           )
166barrier()
167'''
168
169#-------------------------------------------------------------------------
170# Setup initial conditions
171#-------------------------------------------------------------------------
172'''
173
174if myid == 0:
175
176    print 'Setup initial conditions'
177
178    from polygon import *
179    #following sets the stage/water to be offcoast only
180    IC = Polygon_function( [(project.poly_bathy, 0.)], default = tide)
181    domain.set_quantity('stage', IC)
182    domain.set_quantity('friction', 0.01)
183    print 'Start Set quantity'
184
185    domain.set_quantity('elevation',
186#                    filename = project.combined_dir_name + '.pts',
187# MUST USE TXT FILES FOR CACHING TO WORK!
188                    filename = project.combined_smallest_dir_name + '.txt',
189                    use_cache = True,
190                    verbose = True,
191                    alpha = 0.1)
192    print 'Finished Set quantity'
193barrier()
194'''
195print 'Start Set quantity'
196
[4311]197#domain.set_quantity('elevation',
[4307]198#                    filename = project.combined_dir_name + '.pts',
199# MUST USE TXT FILES FOR CACHING TO WORK!
[4311]200#                    filename = project.combined_smallest_dir_name + '.txt',
201#                    use_cache = True,
202#                    verbose = True,
203#                    alpha = 0.1)
204domain.set_quantity('elevation', -42.3)
205
[4307]206print 'Finished Set quantity'
207
208#------------------------------------------------------
209# Distribute domain to implement parallelism !!!
210#------------------------------------------------------
211'''
212if numprocs > 1:
213    domain=distribute(domain)
214'''
215#------------------------------------------------------
216# Set domain parameters
217#------------------------------------------------------
218
219domain.set_quantity('stage', 2.4)
220domain.set_quantity('friction', 0.01)
221domain.set_name(project.scenario_name)
222domain.set_datadir(project.output_run_time_dir)
223
224print 'domain id', id(domain)
225'''
226domain.set_name(project.scenario_name)
227domain.set_datadir(project.output_run_time_dir)
228domain.set_default_order(2) # Apply second order scheme
229domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
230domain.set_store_vertices_uniquely(False)
231domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
232domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
233print 'domain id', id(domain)
234'''
235#-------------------------------------------------------------------------
236# Setup boundary conditions
237#-------------------------------------------------------------------------
238print 'Available boundary tags', domain.get_boundary_tags()
239print 'domain id', id(domain)
240print 'Reading Boundary file',project.boundaries_dir_name4 + '.sww'
241#Bf = File_boundary(boundaries_dir_name + '.sww',
242Bf = File_boundary(project.boundaries_dir_name4 + '.sww',
243                  domain, time_thinning=24, use_cache=False, verbose=True)
244
245print 'finished reading boundary file'
246
247Br = Reflective_boundary(domain)
248Bd = Dirichlet_boundary([tide,0,0])
249
250print'set_boundary'
251##domain.set_boundary({'back': Br,
252##                     'side': Bf,
253##                     'ocean': Bf})
254domain.set_boundary({#'back': Br,
[4311]255                     'side': Bd,
256                    'ocean': Bf
[4307]257                    })
258print'finish set boundary'
259
260#----------------------------------------------------------------------------
261# Evolve system through time
262#----------------------------------------------------------------------------
263
264t0 = time.time()
265
[4311]266for t in domain.evolve(yieldstep = 60, finaltime = 3400):
[4307]267    domain.write_time()
268    domain.write_boundary_statistics(tags = 'ocean')
269
270#for t in domain.evolve(yieldstep = 120, finaltime = 9000):
271#    domain.write_time()
272#    domain.write_boundary_statistics(tags = 'ocean')
273'''     
274for t in domain.evolve(yieldstep = 60, finaltime = 28800
275                       ,skip_initial_step = True):
276    domain.write_time()
277    domain.write_boundary_statistics(tags = 'ocean')   
278
279for t in domain.evolve(yieldstep = 120, finaltime = 34800
280                       ,skip_initial_step = True):
281    domain.write_time()
282    domain.write_boundary_statistics(tags = 'ocean')   
283'''
284x, y = domain.get_maximum_inundation_location()
285q = domain.get_maximum_inundation_elevation()
286
287print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
288
289print 'That took %.2f seconds' %(time.time()-t0)
290
291
292
293
Note: See TracBrowser for help on using the repository browser.