source: anuga_work/production/hobart_2006/run_hobart.py @ 3669

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

reflecting PATH change in project scripts
File size: 10.4 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.py', project.outputtimedir + 'run_hobart.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
71#convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
72
73#creates pts file for onshore DEM - 12.5
74#dem2pts(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,
86        easting_min=project.eastingmin25,
87        easting_max=project.eastingmax25,
88        northing_min=project.northingmin25,
89        northing_max= project.northingmax25,
90        use_cache=True, verbose=True)
91
92#combine_rectangular_points_files(project.onshore_dem_name + '.pts',
93#                                 project.onshore_dem_name_25 + '.pts',
94#                                 project.all_onshore_dem_name + '.pts')
95
96print 'adding data sets'
97G = Geospatial_data(file_name = project.offshore_dem_name_local1 + '.xya')+\
98    Geospatial_data(file_name = project.offshore_dem_name_local2 + '.xya')+\
99    Geospatial_data(file_name = project.offshore_dem_name_local3 + '.xya')+\
100    Geospatial_data(file_name = project.offshore_dem_name_local4 + '.xya')+\
101    Geospatial_data(file_name = project.offshore_dem_name_aho1 + '.xya')+\
102    Geospatial_data(file_name = project.offshore_dem_name_aho2 + '.xya')+\
103    Geospatial_data(file_name = project.offshore_dem_name_aho3 + '.xya')+\
104    Geospatial_data(file_name = project.offshore_dem_name_aho4 + '.xya')+\
105    Geospatial_data(file_name = project.offshore_dem_name_aho5 + '.xya')+\
106    Geospatial_data(file_name = project.offshore_dem_name_aho6 + '.xya')+\
107    Geospatial_data(file_name = project.offshore_dem_name_aho7 + '.xya')+\
108    Geospatial_data(file_name = project.offshore_dem_name_aho8 + '.xya')+\
109    Geospatial_data(file_name = project.offshore_dem_name_aho9 + '.xya')+\
110    Geospatial_data(file_name = project.offshore_dem_name_aho10 + '.xya')+\
111    Geospatial_data(file_name = project.offshore_dem_name_aho11 + '.xya')+\
112    Geospatial_data(file_name = project.offshore_dem_name_aho12 + '.xya')+\
113    Geospatial_data(file_name = project.offshore_dem_name_aho13 + '.xya')+\
114    Geospatial_data(file_name = project.offshore_dem_name_aho14 + '.xya')+\
115    Geospatial_data(file_name = project.offshore_dem_name_aho15 + '.xya')+\
116    Geospatial_data(file_name = project.offshore_dem_name_aho16 + '.xya')+\
117    Geospatial_data(file_name = project.onshore_dem_name_25 + '.pts')
118G.export_points_file(project.combined_dem_name + '.pts')
119#Geospatial_data(file_name = project.all_onshore_dem_name + '.pts')
120#----------------------------------------------------------------------------
121# Create the triangular mesh based on overall clipping polygon with a tagged
122# boundary and interior regions defined in project.py along with
123# resolutions (maximal area of per triangle) for each polygon
124#-------------------------------------------------------------------------------
125
126from anuga.pmesh.mesh_interface import create_mesh_from_regions
127
128# use 75 for onshore components (12.5m DEM)
129island_res = 35000
130hobart_res = 35000
131peninsula_res = 35000
132interior_regions = [[project.poly_hobart, hobart_res],
133                    [project.poly_tasman_peninsula, peninsula_res],
134                    [project.poly_bruny, island_res]]
135
136print 'number of interior regions', len(interior_regions)
137
138from caching import cache
139_ = cache(create_mesh_from_regions,
140          project.polyAll,
141          {'boundary_tags': {'bottom': [0], 'bright': [1],
142                             'topr': [2], 'top': [3], 'left': [4]},
143           'maximum_triangle_area': 250000,
144           'filename': meshname},           
145           #'interior_regions': interior_regions},
146          verbose = True, evaluate=True)
147
148
149#-------------------------------------------------------------------------------                                 
150# Setup computational domain
151#-------------------------------------------------------------------------------                                 
152domain = Domain(meshname, use_cache = False, verbose = True)
153
154print 'Number of triangles = ', len(domain)
155print 'The extent is ', domain.get_extent()
156print domain.statistics()
157
158domain.set_name(project.basename)
159domain.set_datadir(project.outputtimedir)
160domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
161domain.set_minimum_storable_height(0.01)
162domain.set_store_vertices_uniquely(True)  # for writting to sww
163
164#-------------------------------------------------------------------------------                                 
165# Setup initial conditions
166#-------------------------------------------------------------------------------
167
168tide = 0.0
169
170domain.set_quantity('stage', tide)
171domain.set_quantity('friction', 0.0) 
172
173domain.set_quantity('elevation', 
174#                    filename = project.onshore_dem_name + '.pts',
175                    filename = project.combined_dem_name + '.pts',
176#                    filename = project.offshore_dem_name + '.pts',
177                    use_cache = True,
178                    verbose = True,
179                    alpha = 0.1
180                    )
181
182#-------------------------------------------------------------------------------                                 
183# Setup boundary conditions
184#-------------------------------------------------------------------------------
185print 'start ferret2sww'
186from anuga.shallow_water.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
222#Bf = 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
232# for MOST BC
233#domain.set_boundary( {'top': Bd, 'left': Bd,
234#                      'bottom': Bf, 'right': Bf} )
235
236# for testing
237domain.set_boundary( {'topr': Bd, 'left': Bd, 'top': Bd,
238                      'bottom': Bd, 'bright': Bw} )
239
240#-------------------------------------------------------------------------------                                 
241# Evolve system through time
242#-------------------------------------------------------------------------------
243import time
244t0 = time.time()
245
246for t in domain.evolve(yieldstep = 60, finaltime = 120): 
247    domain.write_time()
248    domain.write_boundary_statistics(tags = 'bottom')     
249
250#for t in domain.evolve(yieldstep = 120, finaltime = 12600
251#                       ,skip_initial_step = True):
252#    domain.write_time()
253#    domain.write_boundary_statistics(tags = 'bottom')     
254   
255print 'That took %.2f seconds' %(time.time()-t0)
256
257print 'finished'
Note: See TracBrowser for help on using the repository browser.