source: production/onslow_2006/run_onslow.py @ 3268

Last change on this file since 3268 was 3264, checked in by sexton, 19 years ago

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