source: anuga_work/development/cairns_demo/runcairns.py @ 4862

Last change on this file since 4862 was 4862, checked in by duncan, 16 years ago

ticket:#224

File size: 8.5 KB
Line 
1
2"""Script for running a tsunami inundation scenario for Cairns, QLD Australia.
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
12Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton and
13Nick Bartzis, GA - 2006
14"""
15
16#------------------------------------------------------------------------------
17# Import necessary modules
18#------------------------------------------------------------------------------
19
20# Standard modules
21import os
22import time
23import sys
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
33from anuga.shallow_water.data_manager import dem2pts
34from anuga.fit_interpolate.search_functions import search_times, \
35     reset_search_times
36from anuga.fit_interpolate.general_fit_interpolate import \
37     get_build_quadtree_time
38
39# Application specific imports
40import project                 # Definition of file names and polygons
41
42
43#------------------------------------------------------------------------------
44# Define scenario as either slide or fixed_wave.
45#------------------------------------------------------------------------------
46scenario = 'slide' 
47#scenario = 'fixed_wave'
48
49if os.access(scenario, os.F_OK) == 0:
50    os.mkdir(scenario)
51basename = scenario + 'source'
52
53
54#------------------------------------------------------------------------------
55# Preparation of topographic data
56# Convert ASC 2 DEM 2 PTS using source data and store result in source data
57#------------------------------------------------------------------------------
58
59# Filenames
60dem_name = 'cairns' 
61meshname = 'cairns.msh'
62
63# Create DEM from asc data
64convert_dem_from_ascii2netcdf(dem_name, use_cache=True, verbose=True)
65
66# Create pts file for onshore DEM
67dem2pts(dem_name, use_cache=True, verbose=True)
68
69
70#------------------------------------------------------------------------------
71# Create the triangular mesh based on overall clipping polygon with a tagged
72# boundary and interior regions defined in project.py along with
73# resolutions (maximal area of per triangle) for each polygon
74#------------------------------------------------------------------------------
75
76remainder_res = 10000000
77islands_res = 100000
78cairns_res = 100000
79shallow_res = 500000
80interior_regions = [[project.poly_cairns, cairns_res],
81                    [project.poly_island0, islands_res],
82                    [project.poly_island1, islands_res],
83                    [project.poly_island2, islands_res],
84                    [project.poly_island3, islands_res],
85                    [project.poly_shallow, shallow_res]]
86
87create_mesh_from_regions(project.bounding_polygon,
88                         boundary_tags={'top': [0],
89                                        'ocean_east': [1],
90                                        'bottom': [2],
91                                        'onshore': [3]},
92                         maximum_triangle_area=remainder_res,
93                         filename=meshname,
94                         interior_regions=interior_regions,
95                         use_cache=True,
96                         verbose=True)
97
98
99#------------------------------------------------------------------------------
100# Setup computational domain
101#------------------------------------------------------------------------------
102
103from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
104from anuga.caching import cache
105
106##domain = cache(Domain(meshname, use_cache=True, verbose=True)
107
108domain = cache(pmesh_to_domain_instance,
109               (meshname, Domain),
110               dependencies = [meshname])
111
112print 'Number of triangles = ', len(domain)
113print 'The extent is ', domain.get_extent()
114print domain.statistics()
115
116domain.set_name(basename) 
117domain.set_datadir(scenario)
118domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
119domain.set_minimum_storable_height(0.01)
120
121#print 'domain.tight_slope_limiters', domain.tight_slope_limiters
122domain.tight_slope_limiters = 0
123print 'domain.tight_slope_limiters', domain.tight_slope_limiters
124
125
126
127domain.points_file_block_line_size = 50000
128
129#------------------------------------------------------------------------------
130# Setup initial conditions
131#------------------------------------------------------------------------------
132
133tide = 0.0
134domain.set_quantity('stage', tide)
135domain.set_quantity('friction', 0.0)
136
137
138t0 = time.time()
139domain.set_quantity('elevation', 
140                    filename=dem_name + '.pts',
141                    use_cache=False,
142                    verbose=True,
143                    alpha=0.1)
144
145print 'Fitting the elevation data took %.2f seconds' %(time.time()-t0)
146
147search_one_cell_time, search_more_cells_time = search_times()
148reset_search_times()
149print "search_one_cell_time",search_one_cell_time
150print "search_more_cells_time", search_more_cells_time
151print "build_quadtree_time", get_build_quadtree_time()
152import sys; sys.exit() 
153#------------------------------------------------------------------------------
154# Setup information for slide scenario (to be applied 1 min into simulation
155#------------------------------------------------------------------------------
156
157if scenario == 'slide':
158    # Function for submarine slide
159    from anuga.shallow_water.smf import slide_tsunami 
160    tsunami_source = slide_tsunami(length=35000.0,
161                                   depth=project.slide_depth,
162                                   slope=6.0,
163                                   thickness=500.0, 
164                                   x0=project.slide_origin[0], 
165                                   y0=project.slide_origin[1], 
166                                   alpha=0.0, 
167                                   domain=domain,
168                                   verbose=True)
169
170
171#------------------------------------------------------------------------------
172# Setup boundary conditions
173#------------------------------------------------------------------------------
174
175print 'Available boundary tags', domain.get_boundary_tags()
176
177Br = Reflective_boundary(domain)
178Bd = Dirichlet_boundary([tide,0,0])
179
180# 60 min square wave starting at 1 min, 50m high
181if scenario == 'fixed_wave':
182    Bw = Transmissive_Momentum_Set_Stage_boundary(domain = domain,
183                          f=lambda t: [(60<t<3660)*50, 0, 0])
184    domain.set_boundary({'ocean_east': Bw,
185                         'bottom': Bd,
186                         'onshore': Bd,
187                         'top': Bd})
188
189# boundary conditions for slide scenario
190if scenario == 'slide':
191    domain.set_boundary({'ocean_east': Bd,
192                         'bottom': Bd,
193                         'onshore': Bd,
194                         'top': Bd})
195   
196
197#------------------------------------------------------------------------------
198# Evolve system through time
199#------------------------------------------------------------------------------
200
201t0 = time.time()
202
203from Numeric import allclose
204from anuga.abstract_2d_finite_volumes.quantity import Quantity
205
206if scenario == 'slide':
207   
208    for t in domain.evolve(yieldstep = 10, finaltime = 60): 
209        domain.write_time()
210        domain.write_boundary_statistics(tags = 'ocean_east')     
211       
212    # add slide
213    thisstagestep = domain.get_quantity('stage')
214    if allclose(t, 60):
215        slide = Quantity(domain)
216        slide.set_values(tsunami_source)
217        domain.set_quantity('stage', slide + thisstagestep)
218
219    for t in domain.evolve(yieldstep = 10, finaltime = 5000, 
220                           skip_initial_step = True):
221        domain.write_time()
222        domain.write_boundary_statistics(tags = 'ocean_east')
223
224if scenario == 'fixed_wave':
225
226    # save every two mins leading up to wave approaching land
227    for t in domain.evolve(yieldstep = 120, finaltime = 5000): 
228        domain.write_time()
229        domain.write_boundary_statistics(tags = 'ocean_east')     
230
231    # save every 30 secs as wave starts inundating ashore
232    for t in domain.evolve(yieldstep = 10, finaltime = 10000, 
233                           skip_initial_step = True):
234        domain.write_time()
235        domain.write_boundary_statistics(tags = 'ocean_east')
236           
237print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.