source: production/onslow_2006/run_onslow.py @ 2554

Last change on this file since 2554 was 2554, checked in by ole, 18 years ago

cosmetic

File size: 7.6 KB
Line 
1"""Script for running a tsunami inundation scenario for Onslow, 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.outputdir
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 Nick Bartzis, GA - 2006
11"""
12
13
14#------------------------------------------------------------------------------
15# Import necessary modules
16#------------------------------------------------------------------------------
17
18# Standard modules
19import os
20import time
21
22# Related major packages
23from pyvolution.shallow_water import Domain, Reflective_boundary, \
24                            Dirichlet_boundary, Time_boundary, File_boundary
25from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
26from pyvolution.combine_pts import combine_rectangular_points_files
27from pyvolution.pmesh2domain import pmesh_to_domain_instance
28
29# Application specific imports
30import project                 # Definition of file names and polygons
31from smf import slump_tsunami  # Function for submarine mudslide
32
33
34#------------------------------------------------------------------------------
35# Preparation of topographic data
36#
37# Convert ASC 2 DEM 2 PTS using source data and store result in source data
38# Do for coarse and fine data
39# Fine pts file to be clipped to area of interest
40#------------------------------------------------------------------------------
41
42# filenames
43coarsedemname = project.coarsedemname
44
45onshore_dem_name = project.onshore_dem_name
46
47meshname = project.meshname+'.msh'
48
49source_dir = project.boundarydir
50
51# coarse data
52convert_dem_from_ascii2netcdf(coarsedemname, use_cache=True, verbose=True)
53dem2pts(coarsedemname, use_cache=True, verbose=True)
54
55# fine data (clipping the points file to smaller area)
56convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
57dem2pts(onshore_dem_name,
58        easting_min=project.eastingmin,
59        easting_max=project.eastingmax,
60        northing_min=project.northingmin,
61        northing_max= project.northingmax,
62        use_cache=True, 
63        verbose=True)
64
65
66# combining the coarse and fine data
67combine_rectangular_points_files(project.coarsedemname + '.pts',
68                                 project.onshore_dem_name + '.pts',
69                                 project.combineddemname + '.pts')
70
71
72#------------------------------------------------------------------------------
73# Create the triangular mesh based on overall clipping polygon with a tagged
74# boundary and interior regions defined in project.py along with
75# resolutions (maximal area of per triangle) for each polygon
76#------------------------------------------------------------------------------
77from pmesh.mesh_interface import create_mesh_from_regions
78
79# original
80interior_res = 50000
81interior_regions = [[project.poly_onslow, interior_res],
82                    [project.poly_thevenard, interior_res],
83                    [project.poly_direction, interior_res]]
84                    #[project.testpoly, interior_res]]
85print 'number of interior regions', len(interior_regions)
86
87from caching import cache
88_ = cache(create_mesh_from_regions,
89          project.polyAll,
90          {'boundary_tags': {'top': [0], 'topleft': [1],
91                             'left': [2], 'bottom': [3],
92                             'bottomright': [4], 'topright': [5]},
93           'maximum_triangle_area': 1000000,
94           'filename': meshname,           
95           'interior_regions': interior_regions},
96          verbose = True)
97
98
99#------------------------------------------------------------------------------
100# Setup computational domain
101#------------------------------------------------------------------------------
102
103domain = pmesh_to_domain_instance(meshname, Domain,
104                                  use_cache = True,
105                                  verbose = True)
106
107print 'Number of triangles = ', len(domain)
108print 'The extent is ', domain.get_extent()
109print domain.statistics()
110
111domain.set_name(project.basename)
112domain.set_datadir(project.outputdir)
113domain.set_quantities_to_be_stored(['stage'])
114
115
116#------------------------------------------------------------------------------
117# Set up scenario (tsunami_source is a callable object used with set_quantity)
118#------------------------------------------------------------------------------
119'''
120tsunami_source = slump_tsunami(length=30000.0,
121                               depth=400.0,
122                               slope=6.0,
123                               thickness=176.0,
124                               radius=3330,
125                               dphi=0.23,
126                               x0=project.slump_origin[0],
127                               y0=project.slump_origin[1],
128                               alpha=0.0,
129                               domain=domain)
130
131'''
132#------------------------------------------------------------------------------
133# Setup initial conditions
134#------------------------------------------------------------------------------
135
136tide = 0.
137
138domain.set_quantity('stage', tide)
139domain.set_quantity('friction', 0.0) 
140
141domain.set_quantity('elevation', 
142#                    0.
143#                    filename = project.onshore_dem_name + '.pts',
144                    filename = project.combineddemname + '.pts',
145#                    filename = project.coarsedemname + '.pts',
146                    use_cache = True,
147                    verbose = True
148                    )
149
150
151#------------------------------------------------------------------------------
152# Setup boundary conditions (all reflective)
153#------------------------------------------------------------------------------
154
155from pyvolution.data_manager import ferret2sww
156
157south = project.south
158north = project.north
159west = project.west
160east = project.east
161
162cache(ferret2sww,
163      (source_dir + project.boundary_basename,
164       source_dir + project.boundary_basename), 
165      {'verbose': True,
166#       'minlat': south - 1,
167#       'maxlat': north + 1,
168#       'minlon': west - 1,
169#       'maxlon': east + 1,
170       'minlat': south,
171       'maxlat': north,
172       'minlon': west,
173       'maxlon': east,
174#       'origin': project.mesh_origin,
175       'origin': domain.geo_reference.get_origin(),
176       'mean_stage': tide,
177       'zscale': 1,                 #Enhance tsunami
178       'fail_on_NaN': False,
179       'inverted_bathymetry': True},
180      #evaluate = True,
181       verbose = True)
182
183
184print 'Available boundary tags', domain.get_boundary_tags()
185
186Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
187                    domain, verbose = True)
188Br = Reflective_boundary(domain)
189Bd = Dirichlet_boundary([tide,0,0])
190
191
192# 7 min square wave starting at 1 min, 6m high
193Bw = Time_boundary(domain = domain,
194                   f=lambda t: [(60<t<480)*6, 0, 0])
195
196domain.set_boundary( {'top': Bf, 'topleft': Bf,
197                             'left': Br, 'bottom': Br,
198                             'bottomright': Br, 'topright': Br} )
199
200
201#------------------------------------------------------------------------------
202# Evolve system through time
203#------------------------------------------------------------------------------
204
205import time
206t0 = time.time()
207
208for t in domain.evolve(yieldstep = 50, finaltime = 50): 
209    domain.write_time()
210    domain.write_boundary_statistics(tags = 'top')     
211   
212print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.