source: production/onslow_2006/run_onslow.py @ 2916

Last change on this file since 2916 was 2902, checked in by nick, 19 years ago

update onslow and pt hedland
File size: 10.2 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.outputtimedir
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 *
37import sys
38from pyvolution.util import Screen_Catcher
39
40#-------------------------------------------------------------------------------
41# Preparation of topographic data
42#
43# Convert ASC 2 DEM 2 PTS using source data and store result in source data
44# Do for coarse and fine data
45# Fine pts file to be clipped to area of interest
46#-------------------------------------------------------------------------------
47
48# filenames
49#coarsedemname = project.coarsedemname
50
51onshore_dem_name = project.onshore_dem_name
52
53offshore_points = project.offshore_dem_name
54
55meshname = project.meshname+'.msh'
56
57source_dir = project.boundarydir
58
59# creates copy of code in output dir if dir doesn't exist
60if access(project.outputtimedir,F_OK) == 0 :
61    mkdir (project.outputtimedir)
62copy (project.codedirname, project.outputtimedir + project.codename)
63copy (project.codedir + 'run_onslow.py', project.outputtimedir + 'run_onslow.py')
64print'output dir', project.outputtimedir
65
66#normal screen output is stored in
67screen_output_name = project.outputtimedir + "screen_output.txt"
68screen_error_name = project.outputtimedir + "screen_error.txt"
69
70#used to catch screen output to file
71sys.stdout = Screen_Catcher(screen_output_name)
72#sys.stderr = Screen_Catcher(screen_output_name)
73sys.stderr = Screen_Catcher(screen_error_name)
74
75'''
76copied_files = False
77
78# files to be used
79files_used = [onshore_dem_name, offshore_points,]
80
81if sys.platform != 'win32':   
82    copied_files = True
83    for name in file_list:
84        copy(name, )
85'''   
86#print' most file', project.MOST_dir + project.boundary_basename+'_ha.nc'
87#if access(project.MOST_dir + project.boundary_basename+'_ha.nc',F_OK) == 1 :
88#    print' most file', project.MOST_dir + project.boundary_basename
89
90
91# fine data (clipping the points file to smaller area)
92# creates DEM from asc data
93convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
94
95#creates pts file from DEM
96dem2pts(onshore_dem_name,
97        easting_min=project.eastingmin,
98        easting_max=project.eastingmax,
99        northing_min=project.northingmin,
100        northing_max= project.northingmax,
101        use_cache=True, 
102        verbose=True)
103
104print'create G1'
105G1 = Geospatial_data(file_name = project.offshore_dem_name + '.xya')
106
107print'create G2'
108G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
109
110print'add G1+G2'
111G = G1 + G2
112
113print'export G'
114G.export_points_file(project.combined_dem_name + '.pts')
115
116
117#-------------------------------------------------------------------------------                                 
118# Create the triangular mesh based on overall clipping polygon with a tagged
119# boundary and interior regions defined in project.py along with
120# resolutions (maximal area of per triangle) for each polygon
121#-------------------------------------------------------------------------------
122
123from pmesh.mesh_interface import create_mesh_from_regions
124'''
125# original
126interior_res = 5000
127high_res = 1500
128interior_regions = [[project.poly_onslow, high_res],
129                    [project.poly_thevenard, interior_res],
130                    [project.poly_coast, interior_res]]
131'''
132#new
133region_res = 50000
134coast_res = 2500
135onslow_res = 500
136interior_regions = [[project.poly_onslow, onslow_res],
137                    [project.poly_coast, coast_res],
138                    [project.poly_region, region_res]]
139
140print 'number of interior regions', len(interior_regions)
141
142from caching import cache
143_ = cache(create_mesh_from_regions,
144          project.polyAll,
145          {'boundary_tags': {'top': [0], 'topleft': [1],
146                             'topleft1': [2], 'bottomleft': [3],
147                             'bottom': [4], 'bottomright': [5],
148                             'topright':[6]},
149           'maximum_triangle_area': 200000,
150           'filename': meshname,           
151           'interior_regions': interior_regions},
152          verbose = True)
153
154
155#-------------------------------------------------------------------------------                                 
156# Setup computational domain
157#-------------------------------------------------------------------------------                                 
158
159domain = pmesh_to_domain_instance(meshname, Domain,
160                                  use_cache = False,
161                                  verbose = True)
162
163print 'Number of triangles = ', len(domain)
164print 'The extent is ', domain.get_extent()
165print domain.statistics()
166
167domain.set_name(project.basename)
168domain.set_datadir(project.outputtimedir)
169domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
170
171#-------------------------------------------------------------------------------
172# Set up scenario (tsunami_source is a callable object used with set_quantity)
173#-------------------------------------------------------------------------------
174'''
175tsunami_source = slump_tsunami(length=30000.0,
176                               depth=400.0,
177                               slope=6.0,
178                               thickness=176.0,
179                               radius=3330,
180                               dphi=0.23,
181                               x0=project.slump_origin[0],
182                               y0=project.slump_origin[1],
183                               alpha=0.0,
184                               domain=domain)
185
186'''
187#-------------------------------------------------------------------------------                                 
188# Setup initial conditions
189#-------------------------------------------------------------------------------
190
191tide = 0.0
192
193domain.set_quantity('stage', tide)
194domain.set_quantity('friction', 0.0) 
195print 'hi and file',project.combined_dem_name + '.pts'
196
197domain.set_quantity('elevation', 
198#                    0.
199#                    filename = project.onshore_dem_name + '.pts',
200                    filename = project.combined_dem_name + '.pts',
201#                    filename = project.offshore_dem_name + '.pts',
202                    use_cache = True,
203                    verbose = True,
204                    alpha = 0.1
205                    )
206
207print 'hi1'
208
209#-------------------------------------------------------------------------------                                 
210# Setup boundary conditions (all reflective)
211#-------------------------------------------------------------------------------
212print 'start ferret2sww'
213from pyvolution.data_manager import ferret2sww
214
215south = project.south
216north = project.north
217west = project.west
218east = project.east
219
220#note only need to do when an SWW file for the MOST boundary doesn't exist
221cache(ferret2sww,
222      (source_dir + project.boundary_basename,
223       source_dir + project.boundary_basename), 
224#      (project.MOST_dir + project.boundary_basename,
225#       source_dir + project.boundary_basename),
226      {'verbose': True,
227# note didn't work with the below
228#       'minlat': south - 1,
229#       'maxlat': north + 1,
230#       'minlon': west - 1,
231#       'maxlon': east + 1,
232       'minlat': south,
233       'maxlat': north,
234       'minlon': west,
235       'maxlon': east,
236#       'origin': project.mesh_origin,
237       'origin': domain.geo_reference.get_origin(),
238       'mean_stage': tide,
239       'zscale': 1,                 #Enhance tsunami
240       'fail_on_NaN': False,
241       'inverted_bathymetry': True},
242      #evaluate = True,
243       verbose = True)
244
245
246print 'Available boundary tags', domain.get_boundary_tags()
247
248Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
249                    domain, verbose = True)
250Br = Reflective_boundary(domain)
251Bd = Dirichlet_boundary([tide,0,0])
252
253
254# 7 min square wave starting at 1 min, 6m high
255Bw = Time_boundary(domain = domain,
256                   f=lambda t: [(60<t<480)*6, 0, 0])
257
258domain.set_boundary( {'top': Bf, 'topleft': Bf,
259                             'topleft1': Bf, 'bottomleft': Bd,
260                             'bottom': Br, 'bottomright': Br, 'topright': Bd} )
261
262#-------------------------------------------------------------------------------                                 
263# Evolve system through time
264#-------------------------------------------------------------------------------
265import time
266t0 = time.time()
267
268for t in domain.evolve(yieldstep = 240, finaltime = 7200): 
269    domain.write_time()
270    domain.write_boundary_statistics(tags = 'top')     
271
272for t in domain.evolve(yieldstep = 120, finaltime = 12600
273                       ,skip_initial_step = True): 
274    domain.write_time()
275    domain.write_boundary_statistics(tags = 'top')     
276
277for t in domain.evolve(yieldstep = 60, finaltime = 19800
278                       ,skip_initial_step = True): 
279    domain.write_time()
280    domain.write_boundary_statistics(tags = 'top')     
281   
282for t in domain.evolve(yieldstep = 120, finaltime = 25200
283                       ,skip_initial_step = True): 
284    domain.write_time()
285    domain.write_boundary_statistics(tags = 'top')     
286
287for t in domain.evolve(yieldstep = 240, finaltime = 36000
288                       ,skip_initial_step = True): 
289    domain.write_time()
290    domain.write_boundary_statistics(tags = 'top')     
291 
292print 'That took %.2f seconds' %(time.time()-t0)
293
294print 'finished'
Note: See TracBrowser for help on using the repository browser.