source: production/onslow_2006/run_onslow.py @ 2629

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

small onslow changes

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