source: anuga_work/production/hobart_2006/run_hobart_test.py @ 5212

Last change on this file since 5212 was 4063, checked in by sexton, 17 years ago

updates for nsw slide modelling and cairns demo

File size: 12.6 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 tsunami wave generated by MOST.
9
10Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
11"""
12#-------------------------------------------------------------------------------# Import necessary modules
13#-------------------------------------------------------------------------------
14
15# Standard modules
16import os
17import time
18from shutil import copy
19from os import mkdir, access, F_OK
20import sys
21
22# Related major packages
23from anuga.shallow_water import Domain, Reflective_boundary, \
24                            Dirichlet_boundary, Time_boundary, File_boundary
25from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
26from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files
27from anuga.geospatial_data.geospatial_data import *
28from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
29
30# Application specific imports
31import project                 # Definition of file names and polygons
32
33#-------------------------------------------------------------------------------
34# Copy scripts to time stamped output directory and capture screen
35# output to file
36#-------------------------------------------------------------------------------
37
38# creates copy of code in output dir if dir doesn't exist
39if access(project.outputtimedir,F_OK) == 0 :
40    mkdir (project.outputtimedir)
41copy (project.codedirname, project.outputtimedir + project.codename)
42copy (project.codedir + 'run_hobart_test.py', project.outputtimedir + 'run_hobart_test.py')
43print'output dir', project.outputtimedir
44
45#normal screen output is stored in
46screen_output_name = project.outputtimedir + "screen_output.txt"
47screen_error_name = project.outputtimedir + "screen_error.txt"
48
49#used to catch screen output to file
50sys.stdout = Screen_Catcher(screen_output_name)
51#sys.stderr = Screen_Catcher(screen_output_name)
52sys.stderr = Screen_Catcher(screen_error_name)
53
54print 'USER:    ', project.user
55
56#-------------------------------------------------------------------------------
57# Preparation of topographic data
58#
59# Convert ASC 2 DEM 2 PTS using source data and store result in source data
60#-------------------------------------------------------------------------------
61
62# filenames
63onshore_dem_name = project.onshore_dem_name
64onshore_dem_name_25 = project.onshore_dem_name_25
65meshname = project.meshname+'.msh'
66source_dir = project.boundarydir
67
68copied_files = False
69
70# creates DEM from asc data - 12.5m
71convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
72
73#creates pts file for onshore DEM - 12.5
74dem2pts(onshore_dem_name,
75        easting_min=project.eastingmin,
76        easting_max=project.eastingmax,
77        northing_min=project.northingmin,
78        northing_max= project.northingmax,
79        use_cache=True, verbose=True)
80
81# create DEM from asc data - 25m data
82convert_dem_from_ascii2netcdf(onshore_dem_name_25, use_cache=True, verbose=True)
83
84#creates pts file for onshore DEM - 25
85dem2pts(onshore_dem_name_25, use_cache=True, verbose=True)
86
87combine_rectangular_points_files(project.onshore_dem_name + '.pts',
88                                 project.onshore_dem_name_25 + '.pts',
89                                 project.all_onshore_dem_name + '.pts')
90#print 'local offshore data sets'
91#G1 = Geospatial_data(file_name = project.offshore_dem_name_local1 + '.xya')
92#G2 = Geospatial_data(file_name = project.offshore_dem_name_local2 + '.xya')
93#G3 = Geospatial_data(file_name = project.offshore_dem_name_local3 + '.xya')
94#G4 = Geospatial_data(file_name = project.offshore_dem_name_local4 + '.xya')
95#G11 = G1 + G2 + G3 + G4
96#print 'aho offshore data sets'
97#G5 = Geospatial_data(file_name = project.offshore_dem_name_aho1 + '.xya')
98#G6 = Geospatial_data(file_name = project.offshore_dem_name_aho2 + '.xya')
99#G7 = Geospatial_data(file_name = project.offshore_dem_name_aho3 + '.xya')
100#G8 = Geospatial_data(file_name = project.offshore_dem_name_aho4 + '.xya')
101#G9 = Geospatial_data(file_name = project.offshore_dem_name_aho5 + '.xya')
102#G10 = Geospatial_data(file_name = project.offshore_dem_name_aho6 + '.xya')
103#G11 = Geospatial_data(file_name = project.offshore_dem_name_aho7 + '.xya')
104#G12 = Geospatial_data(file_name = project.offshore_dem_name_aho8 + '.xya')
105#G13 = Geospatial_data(file_name = project.offshore_dem_name_aho9 + '.xya')
106#G14 = Geospatial_data(file_name = project.offshore_dem_name_aho10 + '.xya')
107#G15 = Geospatial_data(file_name = project.offshore_dem_name_aho11 + '.xya')
108#G16 = Geospatial_data(file_name = project.offshore_dem_name_aho12 + '.xya')
109#G17 = Geospatial_data(file_name = project.offshore_dem_name_aho13 + '.xya')
110#G18 = Geospatial_data(file_name = project.offshore_dem_name_aho14 + '.xya')
111#G19 = Geospatial_data(file_name = project.offshore_dem_name_aho15 + '.xya')
112#G20 = Geospatial_data(file_name = project.offshore_dem_name_aho16 + '.xya')
113#G12 = G5+G6+G7+G8+G9+G10+G11+G12+G13+G14+G15+G16+G17+G18+G19+G20
114#print 'local onshore dem'
115#G21 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
116print 'coast'
117#G22 = Geospatial_data(file_name = project.coast_dem_name + '.xya')
118#print 'adding data sets'
119#G = G11 + G12 + G21 + G22
120#G = G1+G2+G3+G4+G5+G6+G7+G8+G9+G10+G11+G12+G13+G14+G15+G16+G17+G18+G19+G20+G21+G22
121#G = G5+G6+G7+G8+G9+G10+G11+G12+G13+G14+G15+G16+G17+G18+G19+G20+G21+G22
122#G = G21+G22+G5+G6+G7+G8+G9+G10+G11+G12+G13+G14+G15+G16+G17+G18+G19+G20
123#G = G21+G22+G5
124#G = Geospatial_data(file_name = project.onshore_dem_name + '.pts') #+ \
125#    Geospatial_data(file_name = project.coast_dem_name + '.xya') #+ \
126#    Geospatial_data(file_name = project.offshore_dem_name_aho1 + '.xya')
127G = Geospatial_data(file_name = project.offshore_dem_name_local1 + '.xya')+\
128    Geospatial_data(file_name = project.offshore_dem_name_local2 + '.xya')+\
129    Geospatial_data(file_name = project.offshore_dem_name_local3 + '.xya')+\
130    Geospatial_data(file_name = project.offshore_dem_name_local4 + '.xya')+\
131    Geospatial_data(file_name = project.offshore_dem_name_aho1 + '.xya')+\
132    Geospatial_data(file_name = project.offshore_dem_name_aho2 + '.xya')+\
133    Geospatial_data(file_name = project.offshore_dem_name_aho3 + '.xya')+\
134    Geospatial_data(file_name = project.offshore_dem_name_aho4 + '.xya')+\
135    Geospatial_data(file_name = project.offshore_dem_name_aho5 + '.xya')+\
136    Geospatial_data(file_name = project.offshore_dem_name_aho6 + '.xya')+\
137    Geospatial_data(file_name = project.offshore_dem_name_aho7 + '.xya')+\
138    Geospatial_data(file_name = project.offshore_dem_name_aho8 + '.xya')+\
139    Geospatial_data(file_name = project.offshore_dem_name_aho9 + '.xya')+\
140    Geospatial_data(file_name = project.offshore_dem_name_aho10 + '.xya')+\
141    Geospatial_data(file_name = project.offshore_dem_name_aho11 + '.xya')+\
142    Geospatial_data(file_name = project.offshore_dem_name_aho12 + '.xya')+\
143    Geospatial_data(file_name = project.offshore_dem_name_aho13 + '.xya')+\
144    Geospatial_data(file_name = project.offshore_dem_name_aho14 + '.xya')+\
145    Geospatial_data(file_name = project.offshore_dem_name_aho15 + '.xya')+\
146    Geospatial_data(file_name = project.offshore_dem_name_aho16 + '.xya')+\
147    Geospatial_data(file_name = project.all_onshore_dem_name + '.pts')
148G.export_points_file(project.combined_dem_name_4 + '.pts')
149
150#-------------------------------------------------------------------------------                                 
151# Create the triangular mesh based on overall clipping polygon with a tagged
152# boundary and interior regions defined in project.py along with
153# resolutions (maximal area of per triangle) for each polygon
154#-------------------------------------------------------------------------------
155
156from anuga.pmesh.mesh_interface import create_mesh_from_regions
157
158# use 75 for onshore components (12.5m DEM)
159hobart_res = 2500
160interior_regions = [[project.poly_hobart1, hobart_res],
161                    [project.poly_hobart2, hobart_res],
162                    [project.poly_hobart3, hobart_res]]
163
164print 'number of interior regions', len(interior_regions)
165
166from caching import cache
167_ = cache(create_mesh_from_regions,
168          project.polyAll,
169           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
170                              'e3': [3], 'e4':[4], 'e5': [5],
171                              'e6': [6], 'e7': [7], 'e8': [8],
172                              'e9': [9], 'e10': [10], 'e11': [11],
173                              'e12': [12], 'e13': [13], 'e14': [14]},
174           'maximum_triangle_area': 750000,
175           'filename': meshname,           
176           'interior_regions': interior_regions},
177          verbose = True, evaluate=True)
178
179
180#-------------------------------------------------------------------------------                                 
181# Setup computational domain
182#-------------------------------------------------------------------------------                                 
183domain = Domain(meshname, use_cache = False, verbose = True)
184
185print 'Number of triangles = ', len(domain)
186print 'The extent is ', domain.get_extent()
187print domain.statistics()
188
189domain.set_name(project.basename)
190domain.set_datadir(project.outputtimedir)
191domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
192domain.set_minimum_storable_height(0.01)
193domain.set_store_vertices_uniquely(True)  # for writting to sww
194
195#-------------------------------------------------------------------------------                                 
196# Setup initial conditions
197#-------------------------------------------------------------------------------
198
199tide = 0.0
200
201domain.set_quantity('stage', tide)
202domain.set_quantity('friction', 0.0) 
203
204domain.set_quantity('elevation', 
205#                    filename = project.onshore_dem_name + '.pts',
206                    filename = project.combined_dem_name_4 + '.pts',
207#                    filename = project.offshore_dem_name + '.pts',
208                    use_cache = True,
209                    verbose = True,
210                    alpha = 0.1
211                    )
212
213#-------------------------------------------------------------------------------                                 
214# Setup boundary conditions (all reflective)
215#-------------------------------------------------------------------------------
216'''
217print 'start ferret2sww'
218from anuga.shallow_water.data_manager import ferret2sww
219south = project.south
220north = project.north
221west = project.west
222east = project.east
223
224#note only need to do when an SWW file for the MOST boundary doesn't exist
225cache(ferret2sww,
226      (source_dir + project.boundary_basename,
227       source_dir + project.boundary_basename),
228#      (project.MOST_dir + project.boundary_basename,
229#       source_dir + project.boundary_basename),
230      {'verbose': True,
231# note didn't work with the below
232#       'minlat': south - 1,
233#       'maxlat': north + 1,
234#       'minlon': west - 1,
235#       'maxlon': east + 1,
236       'minlat': south,
237       'maxlat': north,
238       'minlon': west,
239       'maxlon': east,
240#       'origin': project.mesh_origin,
241       'origin': domain.geo_reference.get_origin(),
242       'mean_stage': tide,
243       'zscale': 1,                 #Enhance tsunami
244       'fail_on_NaN': False,
245       'inverted_bathymetry': True},
246      #evaluate = True,
247       verbose = True,
248       dependencies = source_dir + project.boundary_basename + '.sww')
249
250'''
251print 'Available boundary tags', domain.get_boundary_tags()
252
253Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
254                    domain, verbose = True)
255Br = Reflective_boundary(domain)
256Bd = Dirichlet_boundary([tide,0,0])
257# 7 min square wave starting at 1 min, 6m high
258Bw = Time_boundary(domain = domain,
259                   f=lambda t: [(60<t<480)*6, 0, 0])
260
261
262domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd,
263                      'e5': Bd,  'e6': Bd, 'e7': Bd, 'e8': Bd, 'e9': Bd,
264                      'e10': Bd, 'e11': Bd, 'e12': Bf, 'e13': Bf, 'e14': Bf} )
265
266#-------------------------------------------------------------------------------                                 
267# Evolve system through time
268#-------------------------------------------------------------------------------
269import time
270t0 = time.time()
271
272for t in domain.evolve(yieldstep = 240, finaltime = 7200): 
273    domain.write_time()
274    domain.write_boundary_statistics(tags = 'e13')     
275
276for t in domain.evolve(yieldstep = 30, finaltime = 16000
277                       ,skip_initial_step = True): 
278    domain.write_time()
279    domain.write_boundary_statistics(tags = 'e13')     
280
281for t in domain.evolve(yieldstep = 240, finaltime = 25000
282                       ,skip_initial_step = True):
283    domain.write_time()
284   
285print 'That took %.2f seconds' %(time.time()-t0)
286
287print 'finished'
Note: See TracBrowser for help on using the repository browser.