source: anuga_core/documentation/user_manual/demos/cairns/runcairns.py @ 4856

Last change on this file since 4856 was 4856, checked in by sexton, 16 years ago

rename plot_polygons_points to original name of plot_polygons

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