source: trunk/anuga_work/development/boxingday08/run_boxingday_urs.py @ 7884

Last change on this file since 7884 was 5442, checked in by ole, 16 years ago

Retired h-limiter and beta_h as per ticket:194.
All unit tests and validation tests pass.

File size: 8.0 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] + '_urs'
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# make polygon that contains land that does not affect result.
78
79interior_regions = [[project.patong, patong_res],
80                    [project.bay, bay_res],
81                    [project.contour20m, contour20m_res]]#,
82                    #[project.island_north, island_res],
83                    #[project.island_south, island_res],
84                    #[project.island_south2, island_res]]
85
86#for coarse run to test gauges
87#extent_res = 10000000.0
88#interior_regions=None
89
90from Numeric import arange,allclose
91boundary_tags={'ocean': arange(0,41).tolist(), 'otherocean': [41,44], 'land': [42,43]}
92
93
94#trigs_min = number_mesh_triangles(interior_regions, project.bounding_polygon,extent_res)
95#print 'Minimum number of traingles ', trigs_min
96
97# filenames
98meshname = project.meshname + '_cloud.tsh'
99mesh_elevname = project.mesh_elevname  + '_cloud.tsh'
100
101print 'start create mesh from regions'
102
103_ = cache(create_mesh_from_regions,
104          project.bounding_polygon,
105          {'boundary_tags': boundary_tags,
106           'maximum_triangle_area': extent_res,
107           'filename': meshname,           
108           'interior_regions': interior_regions},
109          verbose = True,
110          #dependencies = ['project.py']
111          dependencies = ['run_boxingday_urs.py']
112          #, evaluate=True
113          )
114
115
116#------------------------------------------------------------------------------
117# Setup computational domain
118#------------------------------------------------------------------------------
119print 'Converting mesh to domain'
120
121#domain = Domain(meshname, use_cache=False, verbose=True)
122
123domain = cache(pmesh_to_domain_instance,
124               (meshname, Domain),
125               dependencies = [meshname])
126
127print 'Number of triangles = ', len(domain)
128print 'The extent is ', domain.get_extent()
129print domain.statistics()
130
131domain.set_name(basename+timestamp)
132#domain.set_name(basename)
133domain.set_datadir(scenario)
134domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
135domain.set_minimum_storable_height(0.01)
136
137domain.tight_slope_limiters = 1
138domain.set_default_order(2)
139print 'domain.tight_slope_limiters', domain.tight_slope_limiters
140
141domain.points_file_block_line_size = 50000
142
143#------------------------------------------------------------------------------
144# Setup initial conditions
145#------------------------------------------------------------------------------
146     
147tide = 0.0
148domain.set_quantity('stage', tide)
149domain.set_quantity('friction', 0.01) 
150
151if scenario == 'poor_simulation':
152     domain.set_quantity('elevation', 
153                         filename=project.poor_combined_dir_name + '.pts',
154                         use_cache=True,
155                         verbose=True,
156                         alpha=0.1)
157
158if scenario == 'good_simulation':
159     domain.set_quantity('elevation', 
160                         filename=project.good_combined_dir_name + '.pts',
161                         use_cache=True,
162                         verbose=True,
163                         alpha=0.1)
164#------------------------------------------------------------------------------
165# Setup boundary conditions
166#------------------------------------------------------------------------------
167boundary_urs_in='data/boxing'
168boundary_urs_out='data/urs_boundary_condition'
169print 'start urs2sww'
170# Note only need to do when an SWW file for the MOST boundary doesn't exist or the bounding polygon has changed
171
172if os.path.exists(boundary_urs_out+'.sww'):
173    print 'sww boundary file already created ensure boundary polygon has not changed'
174else:
175    from anuga.shallow_water.data_manager import urs_ungridded2sww
176    urs_ungridded2sww(basename_in=boundary_urs_in, basename_out=boundary_urs_out,
177                      verbose=True,
178                      #mint=0, maxt=10800,
179                      zscale=1)
180
181print 'Available boundary tags', domain.get_boundary_tags()
182Bf = File_boundary(boundary_urs_out+'.sww', 
183                   domain, time_thinning=1,
184                   use_cache=True,verbose = True)
185
186Br = Reflective_boundary(domain)
187Bd = Dirichlet_boundary([tide,0,0])
188
189domain.set_boundary({'ocean': Bf,
190                     'otherocean': Bd,
191                     'land': Br,
192                     'both': Bd})
193
194#------------------------------------------------------------------------------
195# Evolve system through time
196#------------------------------------------------------------------------------
197import time
198t0 = time.time()
199
200#from Numeric import allclose
201#from anuga.abstract_2d_finite_volumes.quantity import Quantity
202
203# Add new loop that uses larger yieldstep until wave first reaches a point of
204# the ANUGA boundary. Or find a way to clip MOST sww boundary file to only
205# start when boundary stage first becomes non-zero.
206
207for t in domain.evolve(yieldstep = 20.0, finaltime = 18000.0, 
208                       skip_initial_step = False):
209    if allclose(t,10800.0):
210        print 'Changing urs file boundary to dirichlet. Urs gauges only have 3 hours of data'
211        domain.set_boundary({'ocean': Bd,
212                     'otherocean': Bd,
213                     'land': Br,
214                     'both': Bd})
215   
216    domain.write_time()
217           
218print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.