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

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

used to make unit tests to check boundary data to evolved data

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