source: production/onslow_2006/run_onslow.py @ 2556

Last change on this file since 2556 was 2556, checked in by nick, 18 years ago

new onslow with metadata control

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