source: anuga_work/development/boxingday08/run_boxingday_most.py @ 5545

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

updates of scripts used for boxingday validation to include boundary ordering

File size: 8.4 KB
Line 
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# Define scenario as either slide or fixed_wave.
45#------------------------------------------------------------------------------
46
47#scenario = 'poor_simulation'
48scenario = 'good_simulation'
49
50if os.access(scenario, os.F_OK) == 0:
51    os.mkdir(scenario)
52
53timestamp = strftime('%Y%m%d_%H%M%S',localtime())
54basename = scenario[:4] + '_most'
55
56
57#------------------------------------------------------------------------------
58# Preparation of topographic data
59# Convert ASC 2 DEM 2 PTS using source data and store result in source data
60#------------------------------------------------------------------------------
61
62if scenario == 'good_simualation':
63     msg = 'Must use combine_good_data to create bathymetry .pts file'
64     assert os.path.exists(project.good_combined_dir_name+'.pts'), msg
65
66#------------------------------------------------------------------------------
67# Create the triangular mesh based on overall clipping polygon with a tagged
68# boundary and interior regions defined in project.py along with
69# resolutions (maximal area of per triangle) for each polygon
70#------------------------------------------------------------------------------
71extent_res     = 1000000.0
72contour20m_res = 50000.0
73island_res     = 5000.0
74bay_res        = 2000.0
75patong_res     = 400.0
76
77
78
79# make polygon that contains land that does not affect result.
80interior_regions = [[project.patong, patong_res],
81                    [project.bay, bay_res],
82                    [project.contour20m, contour20m_res]]#,
83                    #[project.island_north, island_res],
84                    #[project.island_south, island_res],
85                    #[project.island_south2, island_res]]
86
87#for coarse run to test gauges
88#extent_res = 1000000.0 #causes inaccurate stage
89#interior_regions=None
90
91from Numeric import arange,allclose
92boundary_tags={'ocean': arange(0,41).tolist(), 'otherocean': [41,44], 'land': [42,43]}
93
94
95#trigs_min = number_mesh_triangles(interior_regions, project.bounding_polygon,extent_res)
96#print 'Minimum number of traingles ', trigs_min
97
98# filenames
99meshname = project.meshname + '_most.tsh'
100mesh_elevname = project.mesh_elevname  + '_most.tsh'
101
102print 'start create mesh from regions'
103
104_ = cache(create_mesh_from_regions,
105          project.bounding_polygon,
106          {'boundary_tags': boundary_tags,
107           'maximum_triangle_area': extent_res,
108           'filename': meshname,           
109           'interior_regions': interior_regions},
110          verbose = True#,
111          #dependencies = ['run_boxingday_most.py']
112          #, evaluate=True
113          )
114
115#------------------------------------------------------------------------------
116# Setup computational domain
117#------------------------------------------------------------------------------
118print 'Converting mesh to domain'
119
120#domain = Domain(meshname, use_cache=False, verbose=True)
121
122domain = cache(pmesh_to_domain_instance,
123               (meshname, Domain),
124               dependencies = [meshname])
125
126print 'Number of triangles = ', len(domain)
127print 'The extent is ', domain.get_extent()
128print domain.statistics()
129
130domain.set_name(basename+timestamp) 
131domain.set_datadir(scenario)
132domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
133domain.set_minimum_storable_height(0.01)
134
135domain.tight_slope_limiters = 1
136domain.set_default_order(2)
137print 'domain.tight_slope_limiters', domain.tight_slope_limiters
138
139domain.points_file_block_line_size = 50000
140
141#------------------------------------------------------------------------------
142# Setup initial conditions
143#------------------------------------------------------------------------------
144     
145tide = 0.35
146domain.set_quantity('stage', tide)
147domain.set_quantity('friction', 0.01) 
148
149if scenario == 'poor_simulation':
150     domain.set_quantity('elevation', 
151                         filename=project.poor_combined_dir_name + '.pts',
152                         use_cache=True,
153                         verbose=True,
154                         alpha=0.1)
155
156if scenario == 'good_simulation':
157     domain.set_quantity('elevation', 
158                         filename=project.good_combined_dir_name + '.pts',
159                         use_cache=True,
160                         verbose=True,
161                         alpha=0.1)
162#------------------------------------------------------------------------------
163# Setup boundary conditions
164#------------------------------------------------------------------------------
165
166print 'start ferret2sww'
167# Note only need to do when an SWW file for the MOST boundary doesn't exist or the bounding polygon has changed
168
169if os.path.exists(project.boundary_most_out+'.sww'):
170    print 'sww boundary file already created ensure boundary polygon has not changed'
171else:
172    south = project.south
173    north = project.north
174    west = project.west
175    east = project.east
176
177    #note only need to do when an SWW file for the MOST boundary doesn't exist
178    cache(ferret2sww,
179          (project.boundary_most_in,
180           project.boundary_most_out), 
181          {'verbose': True,
182           'minlat': south,
183           'maxlat': north,
184           'minlon': west,
185           'maxlon': east,
186           #'origin': project.mesh_origin,
187           'origin': domain.geo_reference.get_origin(),
188           'mean_stage': tide,
189           'zscale': 1.0,                 #Enhance tsunami
190           'fail_on_NaN': False,
191           'inverted_bathymetry': True},
192           evaluate = True,
193           verbose = True,
194          dependencies = 'run_boxingday_most.py')
195
196print 'Available boundary tags', domain.get_boundary_tags()
197Bf = File_boundary(project.boundary_most_out+'.sww', 
198                    domain, time_thinning=1, use_cache=True,verbose = True)
199Br = Reflective_boundary(domain)
200Bd = Dirichlet_boundary([tide,0.0,0.0])
201
202domain.set_boundary({'ocean': Bf,
203                     'otherocean': Bd,
204                     'land': Br,
205                     'both': Bd})
206
207#------------------------------------------------------------------------------
208# Evolve system through time
209#------------------------------------------------------------------------------
210import time
211t0 = time.time()
212
213#from Numeric import allclose
214#from anuga.abstract_2d_finite_volumes.quantity import Quantity
215
216# Add new loop that uses larger yieldstep until wave first reaches a point of
217# the ANUGA boundary. Or find a way to clip MOST sww boundary file to only
218# start when boundary stage first becomes non-zero.
219
220for t in domain.evolve(yieldstep = 20.0, finaltime = 18000.0, 
221                       skip_initial_step = False):
222    if allclose(t,10800.0):
223        print 'Changing urs file boundary to dirichlet. Urs gauges only have 3 hours of data'
224        domain.set_boundary({'ocean': Bd,
225                     'otherocean': Bd,
226                     'land': Br,
227                     'both': Bd})
228   
229    domain.write_time()
230print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.