source: anuga_work/production/hobart_2006/run_hobart_usepts.py @ 3852

Last change on this file since 3852 was 3795, checked in by sexton, 17 years ago

incorporating new interior region for hobart

File size: 10.1 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_usepts.py', project.outputtimedir + 'run_hobart_usepts.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
63#onshore_dem_name = project.onshore_dem_name
64onshore_offshore_dem_name = project.onshore_offshore_dem_name
65onshore_offshore_dem_name_25 = project.onshore_offshore_dem_name_25
66meshname = project.meshname+'.msh'
67source_dir = project.boundarydir
68
69copied_files = False
70
71# create DEM from 50m asc data
72convert_dem_from_ascii2netcdf(onshore_offshore_dem_name, use_cache=True, verbose=True)
73
74# creates pts file for combined 50m DEM and make a Geospatial data object
75dem2pts(onshore_offshore_dem_name, use_cache=True, verbose=True)
76G = Geospatial_data(file_name = project.onshore_offshore_dem_name + '.pts')
77
78# clip 50m pts based on interior regions - want 50m data OUTSIDE of these polygons
79clip_regions = [project.poly_hobart1, project.poly_hobart2, \
80                project.poly_hobart3, project.poly_hobart4]
81
82# set up initial value as Geospatial data object
83U = Geospatial_data(clip_regions[0])
84allG = G.clip_outside(U)
85for j in clip_regions[1:]:
86    U = Geospatial_data(j)
87    allG += G.clip_outside(U)
88
89print 'created outside clipping'
90# clip 50m to be inside bounding polygon
91allG = allG.clip(Geospatial_data(project.polyAll))
92print 'finished 50m clipping'
93
94# create DEM from asc data - 25m data
95convert_dem_from_ascii2netcdf(onshore_offshore_dem_name_25, use_cache=True, verbose=True)
96
97# creates pts file for onshore DEM - 25 and make a Geospatial data object
98dem2pts(onshore_dem_name_offshore_25,
99        easting_min=project.eastingmin25_3,
100        easting_max=project.eastingmax25_3,
101        northing_min=project.northingmin25_3,
102        northing_max= project.northingmax25_3,
103        use_cache=True, verbose=True)
104G = Geospatial_data(file_name = project.onshore_offshore_dem_name_25 + '.pts')
105
106# clip 25m pts based on interior regions - want 25m data INSIDE these polygons
107for j in clip_regions[1:]:
108    allG += G.clip(j)
109
110print 'finished 25m clipping'
111
112allG.export_points_file(project.combined_dem_name_3 + '.pts')
113print 'exported points'
114'''
115print 'adding offshore data sets'
116G = Geospatial_data(file_name = project.offshore_dem_name_local1 + '.xya')+\
117    Geospatial_data(file_name = project.offshore_dem_name_local2 + '.xya')+\
118    Geospatial_data(file_name = project.offshore_dem_name_local3 + '.xya')+\
119    Geospatial_data(file_name = project.offshore_dem_name_local4 + '.xya')+\
120    Geospatial_data(file_name = project.offshore_dem_name_aho1 + '.xya')+\
121    Geospatial_data(file_name = project.offshore_dem_name_aho2 + '.xya')+\
122    Geospatial_data(file_name = project.offshore_dem_name_aho3 + '.xya')+\
123    Geospatial_data(file_name = project.offshore_dem_name_aho4 + '.xya')+\
124    Geospatial_data(file_name = project.offshore_dem_name_aho5 + '.xya')+\
125    Geospatial_data(file_name = project.offshore_dem_name_aho6 + '.xya')+\
126    Geospatial_data(file_name = project.offshore_dem_name_aho7 + '.xya')+\
127    Geospatial_data(file_name = project.offshore_dem_name_aho8 + '.xya')+\
128    Geospatial_data(file_name = project.offshore_dem_name_aho9 + '.xya')+\
129    Geospatial_data(file_name = project.offshore_dem_name_aho10 + '.xya')+\
130    Geospatial_data(file_name = project.offshore_dem_name_aho11 + '.xya')+\
131    Geospatial_data(file_name = project.offshore_dem_name_aho12 + '.xya')+\
132    Geospatial_data(file_name = project.offshore_dem_name_aho13 + '.xya')+\
133    Geospatial_data(file_name = project.offshore_dem_name_aho14 + '.xya')+\
134    Geospatial_data(file_name = project.offshore_dem_name_aho15 + '.xya')+\
135    Geospatial_data(file_name = project.offshore_dem_name_aho16 + '.xya')#+\
136    #eospatial_data(file_name = project.onshore_dem_name_25 + '.pts')
137G.export_points_file(project.combined_dem_name_3 + '.pts')
138'''
139#----------------------------------------------------------------------------
140# Create the triangular mesh based on overall clipping polygon with a tagged
141# boundary and interior regions defined in project.py along with
142# resolutions (maximal area of per triangle) for each polygon
143#-------------------------------------------------------------------------------
144
145from anuga.pmesh.mesh_interface import create_mesh_from_regions
146
147# use 75 for onshore components (12.5m DEM)
148hobart_res = 7500
149bathy_res = 50000
150interior_regions = [[project.poly_hobart1, hobart_res],
151                    [project.poly_hobart2, hobart_res],
152                    [project.poly_hobart3, hobart_res],
153                    [project.poly_hobart5, bathy_res]]
154
155print 'number of interior regions', len(interior_regions)
156
157from caching import cache
158_ = cache(create_mesh_from_regions,
159          project.polyAll,
160           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
161                              'e3': [3], 'e4':[4], 'e5': [5],
162                              'e6': [6], 'e7': [7], 'e8': [8],
163                              'e9': [9], 'e10': [10], 'e11': [11],
164                              'e12': [12], 'e13': [13], 'e14': [14],
165                              'e15': [15]},
166           'maximum_triangle_area': 500000,
167           'filename': meshname,           
168           'interior_regions': interior_regions},
169          verbose = True, evaluate=True)
170
171
172#-------------------------------------------------------------------------------                                 
173# Setup computational domain
174#-------------------------------------------------------------------------------                                 
175domain = Domain(meshname, use_cache = False, verbose = True)
176
177print 'Number of triangles = ', len(domain)
178print 'The extent is ', domain.get_extent()
179print domain.statistics()
180
181domain.set_name(project.basename)
182domain.set_datadir(project.outputtimedir)
183domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
184domain.set_minimum_storable_height(0.01)
185
186#-------------------------------------------------------------------------------                                 
187# Setup initial conditions
188#-------------------------------------------------------------------------------
189
190tide = 0.0
191domain.set_quantity('stage', tide)
192domain.set_quantity('friction', 0.0) 
193domain.set_quantity('elevation', 
194                    filename = project.combined_dem_name + '.pts',
195                    use_cache = True,
196                    verbose = True,
197                    alpha = 0.01
198                    )
199
200#-------------------------------------------------------------------------------                                 
201# Setup boundary conditions
202#-------------------------------------------------------------------------------
203
204print 'Available boundary tags', domain.get_boundary_tags()
205
206Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
207                    domain, verbose = True)
208Br = Reflective_boundary(domain)
209Bd = Dirichlet_boundary([tide,0,0])
210# 7 min square wave starting at 1 min, 6m high
211Bw = Time_boundary(domain = domain,
212                   f=lambda t: [(60<t<480)*10, 0, 0])
213
214#
215domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd,
216                      'e5': Bd,  'e6': Bd, 'e7': Bd, 'e8': Bd, 'e9': Bd,
217                      'e10': Bd, 'e11': Bd, 'e12': Bf, 'e13': Bf, 'e14': Bf,
218                      'e15': Bf} )
219
220#-------------------------------------------------------------------------------                                 
221# Evolve system through time
222#-------------------------------------------------------------------------------
223import time
224t0 = time.time()
225
226for t in domain.evolve(yieldstep = 240, finaltime = 6800): 
227    domain.write_time()
228    domain.write_boundary_statistics(tags = 'e14')     
229
230for t in domain.evolve(yieldstep = 30, finaltime = 9000
231                       ,skip_initial_step = True): 
232    domain.write_time()
233    domain.write_boundary_statistics(tags = 'e14')     
234
235for t in domain.evolve(yieldstep = 240, finaltime = 15000
236                       ,skip_initial_step = True): 
237    domain.write_time()
238    domain.write_boundary_statistics(tags = 'e14') 
239   
240print 'That took %.2f seconds' %(time.time()-t0)
241
242print 'finished'
Note: See TracBrowser for help on using the repository browser.