source: anuga_work/development/boxingday08/run_boxingday.py @ 5247

Last change on this file since 5247 was 5247, checked in by jakeman, 16 years ago

First commit of 2008 Boxing day validation files

File size: 8.6 KB
RevLine 
[5247]1"""Script for running the 2004 boxing Day tsunami inundation scenario for
2Phuket, Thailand.
3
4Source data such as elevation and boundary data is assumed to be available in
5directories specified by project.py
6The output sww file is stored in directory named after the scenario, i.e
7slide or fixed_wave.
8
9The scenario is defined by a triangular mesh created from project.polygon,
10the elevation data and a tsunami wave generated by a submarine mass failure.
11
12Author: John Jakeman, The Australian National University (2008)
13"""
14
15#------------------------------------------------------------------------------
16# Import necessary modules
17#------------------------------------------------------------------------------
18
19# Standard modules
20import os
21import time
22import sys
23from time import localtime, strftime
24
25# Related major packages
26from anuga.shallow_water import Domain
27from anuga.shallow_water import Reflective_boundary
28from anuga.shallow_water import Dirichlet_boundary
29from anuga.shallow_water import Time_boundary
30from anuga.shallow_water import File_boundary
31from anuga.pmesh.mesh_interface import create_mesh_from_regions
32from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, ferret2sww
33from anuga.shallow_water.data_manager import dem2pts
34from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
35from anuga.utilities.polygon import number_mesh_triangles
36from anuga.fit_interpolate.fit import fit_to_mesh_file
37from anuga.caching import cache
38from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
39
40# Application specific imports
41import project                 # Definition of file names and polygons
42
43
44#------------------------------------------------------------------------------
45# Define scenario as either slide or fixed_wave.
46#------------------------------------------------------------------------------
47
48#scenario = 'poor_simulation'
49scenario = 'good_simulation'
50
51if os.access(scenario, os.F_OK) == 0:
52    os.mkdir(scenario)
53
54timestamp = strftime('%Y%m%d_%H%M%S',localtime())
55basename = scenario + 'source'
56
57
58#------------------------------------------------------------------------------
59# Preparation of topographic data
60# Convert ASC 2 DEM 2 PTS using source data and store result in source data
61#------------------------------------------------------------------------------
62
63if scenario == 'good_simualation':
64     msg = 'Must use combine_good_data to create bathymetry .pts file'
65     assert os.path.exists(project.good_combined_dir_name+'.pts'), msg
66
67#------------------------------------------------------------------------------
68# Create the triangular mesh based on overall clipping polygon with a tagged
69# boundary and interior regions defined in project.py along with
70# resolutions (maximal area of per triangle) for each polygon
71#------------------------------------------------------------------------------
72#extent_res     = 10000000.0
73#contour20m_res = 50000.0
74#island_res     = 5000.0
75#bay_res        = 2000.0
76#patong_res     = 100.0
77
78extent_res     = 10000000.0
79contour20m_res = 1000000.0
80island_res     = 1000000.0
81bay_res        = 1000000.0
82patong_res     = 1000000.0
83
84
85# make polygon that contains land that does not affect result.
86
87interior_regions=None
88#interior_regions = [[project.patong, patong_res],
89#                    [project.bay, bay_res],
90#                    [project.contour20m, contour20m_res],
91#                    [project.island_north, island_res],
92#                    [project.island_south, island_res],
93#                    [project.island_south2, island_res]]
94
95# need to include runup_res regions
96
97# There are 9 boundaries. FIX ME: CHECK TAGS ARE CORRECT. I.E. THAT I HAVE
98#CHOSEN CORECT ORDER AS SEEEN BY ANUGA. IS THERE A WAY TO DO This
99print 'must determine which boundaries are which'
100
101#boundary_tags={'ocean': [1,2], 'otherocean': [0,3, 7, 8], 'land': [5, 6], 'both': [4]}
102
103boundary_tags =  {'most':[10,11,12,13,14,15],
104                  'land':[0,1,2,3,4,5,6,7,8],
105                  'shallow':[9,16]}
106
107#trigs_min = number_mesh_triangles(interior_regions, project.bounding_polygon,extent_res)
108#print 'Minimum number of traingles ', trigs_min
109
110# filenames
111meshname = project.meshname + '.tsh'
112mesh_elevname = project.mesh_elevname  + '.tsh'
113
114print 'start create mesh from regions'
115
116_ = cache(create_mesh_from_regions,
117          project.bounding_polygon,
118          {'boundary_tags': boundary_tags,
119           'maximum_triangle_area': extent_res,
120           'filename': meshname,           
121           'interior_regions': interior_regions},
122          verbose = True,
123          dependencies = ['run_boxingday.py']
124          #, evaluate=True
125          )
126
127
128#------------------------------------------------------------------------------
129# Setup computational domain
130#------------------------------------------------------------------------------
131print 'Converting mesh to domain'
132
133#domain = Domain(meshname, use_cache=False, verbose=True)
134
135domain = cache(pmesh_to_domain_instance,
136               (meshname, Domain),
137               dependencies = [meshname])
138
139print 'Number of triangles = ', len(domain)
140print 'The extent is ', domain.get_extent()
141print domain.statistics()
142
143domain.set_name(basename+'friction'+timestamp) 
144domain.set_datadir(scenario)
145domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
146domain.set_minimum_storable_height(0.01)
147
148domain.beta_h = 0
149domain.tight_slope_limiters = 1
150domain.set_default_order(2)
151print 'domain.tight_slope_limiters', domain.tight_slope_limiters
152
153domain.points_file_block_line_size = 50000
154
155#------------------------------------------------------------------------------
156# Setup initial conditions
157#------------------------------------------------------------------------------
158     
159tide = 0.0
160domain.set_quantity('stage', tide)
161domain.set_quantity('friction', 0.01) 
162
163if scenario == 'poor_simulation':
164     domain.set_quantity('elevation', 
165                         filename=project.poor_combined_dir_name + '.pts',
166                         use_cache=True,
167                         verbose=True,
168                         alpha=0.1)
169
170if scenario == 'good_simulation':
171     domain.set_quantity('elevation', 
172                         2.0,#filename=project.good_combined_dir_name + '.pts',
173                         use_cache=True,
174                         verbose=True,
175                         alpha=0.1)
176#------------------------------------------------------------------------------
177# Setup boundary conditions
178#------------------------------------------------------------------------------
179
180print 'start ferret2sww'
181south = project.south
182north = project.north
183west = project.west
184east = project.east
185
186# Note only need to do when an SWW file for the MOST boundary doesn't exist or the bounding polygon has changed
187
188if os.path.exists(project.boundary_most_out+'.sww'):
189    print 'sww boundary file already created ensure boundary polygon has not changed'
190else:
191     cache(ferret2sww,
192           (project.boundary_most_in,
193            project.boundary_most_out), 
194           {'verbose': True,
195            'minlat': south,
196            'maxlat': north,
197            'minlon': west,
198            'maxlon': east,
199            #'origin': project.mesh_origin,
200            'origin': domain.geo_reference.get_origin(),
201            'mean_stage': tide,
202            'zscale': 1.0,                 #Enhance tsunami
203            'fail_on_NaN': False,
204            'inverted_bathymetry': True},
205           #evaluate = True,
206           verbose = True,
207           dependencies = 'run_boxingday.py')
208
209print 'Available boundary tags', domain.get_boundary_tags()
210Bf = File_boundary(project.boundary_most_out+'.sww', 
211                    domain, verbose = True)
212Br = Reflective_boundary(domain)
213Bd = Dirichlet_boundary([tide,0,0])
214
215#domain.set_boundary({'ocean': Bf,
216#                     'otherocean': Bd,
217#                     'land': Br,
218#                     'both': Bd})
219
220domain.set_boundary({'most': Bf,
221                     'land': Br,
222                     'shallow': Bd})
223
224#------------------------------------------------------------------------------
225# Evolve system through time
226#------------------------------------------------------------------------------
227import time
228t0 = time.time()
229
230#from Numeric import allclose
231#from anuga.abstract_2d_finite_volumes.quantity import Quantity
232
233# Add new loop that uses larger yieldstep until wave first reaches a point of
234# the ANUGA boundary. Or find a way to clip MOST sww boundary file to only
235# start when boundary stage first becomes non-zero.
236
237for t in domain.evolve(yieldstep = 60, finaltime = 18000, 
238                       skip_initial_step = False):
239     domain.write_time()
240     
241# Add new loop that uses smaller yieldstep when wave approaches shore?
242
243           
244print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.