source: anuga_work/production/onslow_2006/run_onslow.py @ 4504

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

visualization run for onslow

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