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

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

introduce refined grids around specified sites for Hobart

File size: 8.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 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_clipdata_refine.py', project.outputtimedir + 'run_hobart_clipdata_refine.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]
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_offshore_dem_name_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
115#----------------------------------------------------------------------------
116# Create the triangular mesh based on overall clipping polygon with a tagged
117# boundary and interior regions defined in project.py along with
118# resolutions (maximal area of per triangle) for each polygon
119#-------------------------------------------------------------------------------
120
121from anuga.pmesh.mesh_interface import create_mesh_from_regions
122
123# use 75 for onshore components (12.5m DEM)
124hobart_res = 7500
125bathy_res = 25000
126refine_res = 250
127interior_regions = [[project.poly_site13, refine_res],
128                    [project.poly_kingston, refine_res],
129                    [project.poly_bruny, refine_res],
130                    [project.poly_hobart5, bathy_res]]
131
132print 'number of interior regions', len(interior_regions)
133
134from caching import cache
135_ = cache(create_mesh_from_regions,
136          project.polyAll,
137           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
138                              'e3': [3], 'e4':[4], 'e5': [5],
139                              'e6': [6], 'e7': [7], 'e8': [8],
140                              'e9': [9], 'e10': [10], 'e11': [11],
141                              'e12': [12], 'e13': [13], 'e14': [14]},
142           'maximum_triangle_area': 700000,
143           'filename': meshname,           
144           'interior_regions': interior_regions},
145          verbose = True, evaluate=True)
146
147
148#-------------------------------------------------------------------------------                                 
149# Setup computational domain
150#-------------------------------------------------------------------------------                                 
151domain = Domain(meshname, use_cache = False, verbose = True)
152
153print 'Number of triangles = ', len(domain)
154print 'The extent is ', domain.get_extent()
155print domain.statistics()
156
157domain.set_name(project.basename)
158domain.set_datadir(project.outputtimedir)
159domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
160domain.set_minimum_storable_height(0.01)
161
162#-------------------------------------------------------------------------------                                 
163# Setup initial conditions
164#-------------------------------------------------------------------------------
165
166tide = 0.0
167domain.set_quantity('stage', tide)
168domain.set_quantity('friction', 0.0) 
169domain.set_quantity('elevation', 
170                    filename = project.combined_dem_name + '.pts',
171                    use_cache = True,
172                    verbose = True,
173                    alpha = 0.01
174                    )
175
176#-------------------------------------------------------------------------------                                 
177# Setup boundary conditions
178#-------------------------------------------------------------------------------
179
180print 'Available boundary tags', domain.get_boundary_tags()
181
182Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
183                    domain, verbose = True)
184Br = Reflective_boundary(domain)
185Bd = Dirichlet_boundary([tide,0,0])
186# 7 min square wave starting at 1 min, 6m high
187Bw = Time_boundary(domain = domain,
188                   f=lambda t: [(60<t<480)*10, 0, 0])
189
190#
191domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd,
192                      'e5': Bd,  'e6': Bd, 'e7': Bd, 'e8': Bd, 'e9': Bd,
193                      'e10': Bd, 'e11': Bd, 'e12': Bf, 'e13': Bf, 'e14': Bf} )
194
195#-------------------------------------------------------------------------------                                 
196# Evolve system through time
197#-------------------------------------------------------------------------------
198import time
199t0 = time.time()
200
201for t in domain.evolve(yieldstep = 240, finaltime = 6800): 
202    domain.write_time()
203    domain.write_boundary_statistics(tags = 'e14')     
204
205for t in domain.evolve(yieldstep = 30, finaltime = 15000
206                       ,skip_initial_step = True): 
207    domain.write_time()
208    domain.write_boundary_statistics(tags = 'e14')     
209
210for t in domain.evolve(yieldstep = 240, finaltime = 20000
211                       ,skip_initial_step = True): 
212    domain.write_time()
213    domain.write_boundary_statistics(tags = 'e14') 
214   
215print 'That took %.2f seconds' %(time.time()-t0)
216
217print 'finished'
Note: See TracBrowser for help on using the repository browser.