source: anuga_work/production/busselton/2006/run_busselton.py @ 4995

Last change on this file since 4995 was 4995, checked in by nick, 16 years ago

updates to busselton

File size: 8.0 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# 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.geospatial_data.geospatial_data import *
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.utilities.polygon import plot_polygons, polygon_area
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
46start_screen_catcher(project.output_run_time_dir, myid, numprocs)
47
48print 'time stamp: ',project.gtime
49print 'USER: ', project.user
50print 'min triangles', project.trigs_min
51
52# filenames
53boundaries_name = project.scenario
54meshes_dir_name = project.meshes_dir_name+'.msh'
55#boundaries_time_dir_name = project.boundaries_time_dir_name
56boundaries_dir_name = project.boundaries_dir_name
57
58tide = project.tide
59
60# creates copy of code in output dir
61if myid == 0:
62    copy_code_files(project.output_run_time_dir,__file__, 
63                 dirname(project.__file__)+sep+ project.__name__+'.py' )
64barrier()
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
73if myid == 0:
74
75    print 'start create mesh from regions'
76#    interior_regions = [[project.poly_busselton1, project.res_bus1],
77#                    [project.poly_busselton2, project.res_bus2],
78#                    [project.poly_busselton3, project.res_bus3]]   
79
80#    from anuga.utilities.polygon import plot_polygons
81#    figname = project.output_run_time_dir + 'poly_pic'
82#    plot_polygons([project.poly_coast,project.poly_pipeline,project.poly_facility,
83#               project.bounding_polygon],
84#               figname,
85#               verbose = True)
86#    import sys; sys.exit()
87
88    print 'start create mesh from regions'
89    create_mesh_from_regions(project.poly_all,
90                         boundary_tags={'back': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 14],
91                                        'side': [10,13], 'ocean': [11, 12]},
92                         maximum_triangle_area=project.res_poly_all,
93                         interior_regions=project.interior_regions,
94#                         filename=meshes_time_dir_name,
95                         filename=meshes_dir_name,
96                         use_cache=True,
97                         verbose=True)
98
99# to sync all processors are ready
100barrier()
101
102#-------------------------------------------------------------------------
103# Setup computational domain
104#-------------------------------------------------------------------------
105print 'Setup computational domain'
106#domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True)
107domain = Domain(meshes_dir_name, use_cache=True, verbose=True)
108print domain.statistics()
109'''
110print 'starting to create boundary conditions'
111boundaries_in_dir_name = project.boundaries_in_dir_name
112
113from anuga.shallow_water.data_manager import urs2sww, ferret2sww
114
115# put above distribute
116print 'boundary file is: ',boundaries_dir_name
117from caching import cache
118if myid == 0:
119    cache(ferret2sww,
120          (boundaries_in_dir_name,
121    #       boundaries_time_dir_name),
122           boundaries_dir_name),
123          {'verbose': True,
124           'minlat': project.south_boundary,
125           'maxlat': project.north_boundary,
126           'minlon': project.west_boundary,
127           'maxlon': project.east_boundary,
128#           'minlat': project.south,
129#           'maxlat': project.north,
130#           'minlon': project.west,
131#           'maxlon': project.east,
132           'mint': 0, 'maxt': 35100,
133           'origin': domain.geo_reference.get_origin(),
134           'mean_stage': project.tide,
135#           'zscale': 1,                 #Enhance tsunami
136           'fail_on_NaN': False},
137           verbose = True,
138           )
139barrier()
140'''
141
142#-------------------------------------------------------------------------
143# Setup initial conditions
144#-------------------------------------------------------------------------
145if myid == 0:
146
147    print 'Setup initial conditions'
148
149    domain.set_quantity('stage', tide)
150    domain.set_quantity('friction', 0.01) 
151    #combined_time_dir_name = project.topographies_dir+build_time+project.combined_name
152    print 'Start Set quantity'
153
154    domain.set_quantity('elevation', 
155#                    filename = project.combined_smaller_dir_name+'.xya',
156                    filename = project.combined_dir_name+'.pts',
157                    use_cache = True,
158                    verbose = True,
159                    alpha = 0.1)
160    print 'Finished Set quantity'
161barrier()
162
163#------------------------------------------------------
164# Distribute domain to implement parallelism !!!
165#------------------------------------------------------
166
167if numprocs > 1:
168    domain=distribute(domain)
169
170#------------------------------------------------------
171# Set domain parameters
172#------------------------------------------------------
173
174domain.set_name(project.scenario_name)
175domain.set_datadir(project.output_run_time_dir)
176domain.set_default_order(2) # Apply second order scheme
177domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
178domain.set_store_vertices_uniquely(False)
179domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
180domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
181
182#-------------------------------------------------------------------------
183# Setup boundary conditions
184#-------------------------------------------------------------------------
185print 'Available boundary tags', domain.get_boundary_tags()
186
187print 'Reading Boundary file'
188print 'domain id', id(domain)
189#Bf = File_boundary(boundaries_dir_name + '.sww',
190#                  domain, time_thinning=5, use_cache=True, verbose=True)
191
192print 'finished reading boundary file'
193Br = Reflective_boundary(domain)
194Bd = Dirichlet_boundary([tide,0,0])
195Bo = Dirichlet_boundary([tide+5.0,0,0])
196
197print'set_boundary'
198domain.set_boundary({'back': Br,
199                     'side': Bd,
200                     'ocean': Bd}) 
201print'finish set boundary'
202
203#----------------------------------------------------------------------------
204# Evolve system through time
205#----------------------------------------------------------------------------
206
207t0 = time.time()
208
209for t in domain.evolve(yieldstep = 60, finaltime = 10000):
210    domain.write_time()
211    domain.write_boundary_statistics(tags = 'ocean')
212    if allclose(t, 120):
213        domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
214
215    if allclose(t, 1020):
216        domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
217
218   
219print 'That took %.2f seconds' %(time.time()-t0)
220
Note: See TracBrowser for help on using the repository browser.