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

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

Updates and comments arising from Rajaraman's postings on sourceforge: http://sourceforge.net/forum/forum.php?thread_id=1861165&forum_id=593114

File size: 8.0 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#------------------------------------------------------------------------------
41scenario = 'slide' 
42#scenario = '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#------------------------------------------------------------------------------
97
98domain = Domain(meshname, use_cache=True, verbose=True)
99
100print 'Number of triangles = ', len(domain)
101print 'The extent is ', domain.get_extent()
102print domain.statistics()
103
104domain.set_name(basename) 
105domain.set_datadir(scenario)
106domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
107domain.set_minimum_storable_height(0.01)
108
109#print 'domain.tight_slope_limiters', domain.tight_slope_limiters
110domain.tight_slope_limiters = 0
111print 'domain.tight_slope_limiters', domain.tight_slope_limiters
112
113
114
115#------------------------------------------------------------------------------
116# Setup initial conditions
117#------------------------------------------------------------------------------
118
119tide = 0.0
120domain.set_quantity('stage', tide)
121domain.set_quantity('friction', 0.0) 
122domain.set_quantity('elevation', 
123                    filename=dem_name + '.pts',
124                    use_cache=True,
125                    verbose=True,
126                    alpha=0.1)
127
128
129
130#------------------------------------------------------------------------------
131# Setup information for slide scenario (to be applied 1 min into simulation
132#------------------------------------------------------------------------------
133
134if scenario == 'slide':
135    # Function for submarine slide
136    from anuga.shallow_water.smf import slide_tsunami 
137    tsunami_source = slide_tsunami(length=35000.0,
138                                   depth=project.slide_depth,
139                                   slope=6.0,
140                                   thickness=500.0, 
141                                   x0=project.slide_origin[0], 
142                                   y0=project.slide_origin[1], 
143                                   alpha=0.0, 
144                                   domain=domain,
145                                   verbose=True)
146
147
148#------------------------------------------------------------------------------
149# Setup boundary conditions
150#------------------------------------------------------------------------------
151
152print 'Available boundary tags', domain.get_boundary_tags()
153
154Br = Reflective_boundary(domain)
155Bd = Dirichlet_boundary([tide,0,0])
156
157# 60 min square wave starting at 1 min, 50m high
158if scenario == 'fixed_wave':
159    # FIXME (Ole): Change this to Transmissive Momentum as
160    # suggested by Rajaraman (ticket:208)
161    Bw = Time_boundary(domain = domain,
162                       f=lambda t: [(60<t<3660)*50, 0, 0])
163    domain.set_boundary({'ocean_east': Bw,
164                         'bottom': Bd,
165                         'onshore': Bd,
166                         'top': Bd})
167
168# boundary conditions for slide scenario
169if scenario == 'slide':
170    domain.set_boundary({'ocean_east': Bd,
171                         'bottom': Bd,
172                         'onshore': Bd,
173                         'top': Bd})
174   
175
176#------------------------------------------------------------------------------
177# Evolve system through time
178#------------------------------------------------------------------------------
179
180import time
181t0 = time.time()
182
183from Numeric import allclose
184from anuga.abstract_2d_finite_volumes.quantity import Quantity
185
186if scenario == 'slide':
187   
188    for t in domain.evolve(yieldstep = 10, finaltime = 60): 
189        domain.write_time()
190        domain.write_boundary_statistics(tags = 'ocean_east')     
191       
192    # add slide
193    thisstagestep = domain.get_quantity('stage')
194    if allclose(t, 60):
195        slide = Quantity(domain)
196        slide.set_values(tsunami_source)
197        domain.set_quantity('stage', slide + thisstagestep)
198
199    for t in domain.evolve(yieldstep = 10, finaltime = 5000, 
200                           skip_initial_step = True):
201        domain.write_time()
202        domain.write_boundary_statistics(tags = 'ocean_east')
203
204if scenario == 'fixed_wave':
205
206    # save every two mins leading up to wave approaching land
207    for t in domain.evolve(yieldstep = 120, finaltime = 5000): 
208        domain.write_time()
209        domain.write_boundary_statistics(tags = 'ocean_east')     
210
211    # save every 30 secs as wave starts inundating ashore
212    for t in domain.evolve(yieldstep = 10, finaltime = 10000, 
213                           skip_initial_step = True):
214        domain.write_time()
215        domain.write_boundary_statistics(tags = 'ocean_east')
216           
217print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.