Ticket #141: run_dampier.2.py

File run_dampier.2.py, 7.8 KB (added by nick, 17 years ago)

Actually cyclone run file (other one is actually the tornado file)

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=True,
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=True, verbose=True)
103print 'domain id', id(domain) 
104print 'myhash', myhash(domain)     
105       
106print domain.statistics()
107
108boundaries_dir_name=project.boundaries_dir_name
109
110print 'starting to create boundary conditions'
111
112from 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_bathy, 0.)], default = tide)
125    domain.set_quantity('stage', IC)
126    domain.set_quantity('friction', 0.01) 
127    print 'Start Set quantity'
128
129    domain.set_quantity('elevation', 
130#                    filename = project.combined_dir_name + '.pts',
131# MUST USE TXT FILES FOR CACHING TO WORK!
132                    filename = project.combined_dir_name + '.txt',
133#                    filename = project.combined_smallest_dir_name + '.txt',
134                    use_cache = True,
135                    verbose = True,
136                    alpha = 0.1)
137    print 'Finished Set quantity'
138barrier()
139
140
141#------------------------------------------------------
142# Distribute domain to implement parallelism !!!
143#------------------------------------------------------
144
145if numprocs > 1:
146    domain=distribute(domain)
147
148#------------------------------------------------------
149# Set domain parameters
150#------------------------------------------------------
151print 'domain id', id(domain)
152domain.set_name(project.scenario_name)
153domain.set_datadir(project.output_run_time_dir)
154domain.set_default_order(2) # Apply second order scheme
155domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
156domain.set_store_vertices_uniquely(False)
157domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
158domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
159print 'domain id', id(domain)
160#domain.beta_h = 0
161#domain.limit2007 = 1
162
163#-------------------------------------------------------------------------
164# Setup boundary conditions
165#-------------------------------------------------------------------------
166print 'Available boundary tags', domain.get_boundary_tags()
167print 'domain id', id(domain)
168#print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
169
170Bf = Field_boundary(project.boundaries_dir_namea + '.sww',
171                    domain, time_thinning=12, mean_stage=tide, 
172                    use_cache=True, verbose=True)
173
174print 'finished reading boundary file'
175
176Br = Reflective_boundary(domain)
177Bd = Dirichlet_boundary([tide,0,0])
178
179print'set_boundary'
180##domain.set_boundary({'back': Br,
181##                     'side': Bf,
182##                     'ocean': Bf})
183domain.set_boundary({'back': Br,
184                     'side': Bd,
185                     'ocean': Bf}) 
186print'finish set boundary'
187
188#----------------------------------------------------------------------------
189# Evolve system through time
190#----------------------------------------------------------------------------
191
192t0 = time.time()
193
194for t in domain.evolve(yieldstep = 60, finaltime = 29990):
195#    domain.write_time()
196    domain.write_time(track_speeds=True)
197    domain.write_boundary_statistics(tags = 'ocean')
198   
199#    domain.write_time(track_speeds=True)
200
201#for t in domain.evolve(yieldstep = 120, finaltime = 9000):
202#    domain.write_time()
203#    domain.write_boundary_statistics(tags = 'ocean')
204'''     
205for t in domain.evolve(yieldstep = 60, finaltime = 28800
206                       ,skip_initial_step = True):
207    domain.write_time()
208    domain.write_boundary_statistics(tags = 'ocean')   
209
210for t in domain.evolve(yieldstep = 120, finaltime = 34800
211                       ,skip_initial_step = True):
212    domain.write_time()
213    domain.write_boundary_statistics(tags = 'ocean')   
214'''
215x, y = domain.get_maximum_inundation_location()
216q = domain.get_maximum_inundation_elevation()
217
218print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
219
220print 'That took %.2f seconds' %(time.time()-t0)
221