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

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

updating Cairns demo for the usermanual

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