source: production/hobart_2006/run_hobart.py @ 3559

Last change on this file since 3559 was 3559, checked in by sexton, 18 years ago

starting on Hobart script

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