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

Last change on this file since 6889 was 6889, checked in by ole, 11 years ago

Made sure demos in the manual run and did some cleanup

File size: 7.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
11Geoscience Australia, 2004-present
12"""
13
14#------------------------------------------------------------------------------
15# Import necessary modules
16#------------------------------------------------------------------------------
17
18# Standard modules
19import os
20import time
21import sys
22
23# Related major packages
24from anuga.interface import create_domain_from_regions
25from anuga.interface import Reflective_boundary
26from anuga.interface import Dirichlet_boundary
27from anuga.interface import Time_boundary
28from anuga.interface import File_boundary
29from anuga.interface import Transmissive_stage_zero_momentum_boundary
30
31from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
32from anuga.shallow_water.data_manager import dem2pts
33
34from anuga.shallow_water.smf import slide_tsunami 
35
36# Application specific imports
37import project                 # Definition of file names and polygons
38
39
40#------------------------------------------------------------------------------
41# Preparation of topographic data
42# Convert ASC 2 DEM 2 PTS using source data and store result in source data
43#------------------------------------------------------------------------------
44
45# Create DEM from asc data
46convert_dem_from_ascii2netcdf(project.demname, use_cache=True, verbose=True)
47
48# Create pts file for onshore DEM
49dem2pts(project.demname, use_cache=True, verbose=True)
50
51
52#------------------------------------------------------------------------------
53# Create the triangular mesh and domain based on
54# overall clipping polygon with a tagged
55# boundary and interior regions as defined in project.py
56#------------------------------------------------------------------------------
57
58domain = create_domain_from_regions(project.bounding_polygon,
59                                    boundary_tags={'top': [0],
60                                                   'ocean_east': [1],
61                                                   'bottom': [2],
62                                                   'onshore': [3]},
63                                    maximum_triangle_area=project.default_res,
64                                    mesh_filename=project.meshname,
65                                    interior_regions=project.interior_regions,
66                                    use_cache=True,
67                                    verbose=True)
68
69
70# Print some stats about mesh and domain
71print 'Number of triangles = ', len(domain)
72print 'The extent is ', domain.get_extent()
73print domain.statistics()
74                                   
75#------------------------------------------------------------------------------
76# Setup parameters of computational domain
77#------------------------------------------------------------------------------
78
79
80domain.set_name('cairns_' + project.scenario) # Name of sww file
81domain.set_datadir('.')                       # Store sww output here
82domain.set_minimum_storable_height(0.01)      # Store only depth > 1cm
83
84
85#------------------------------------------------------------------------------
86# Setup initial conditions
87#------------------------------------------------------------------------------
88
89tide = 0.0
90domain.set_quantity('stage', tide)
91domain.set_quantity('friction', 0.0) 
92domain.set_quantity('elevation', 
93                    filename=project.demname + '.pts',
94                    use_cache=True,
95                    verbose=True,
96                    alpha=0.1)
97
98
99#------------------------------------------------------------------------------
100# Setup information for slide scenario (to be applied 1 min into simulation
101#------------------------------------------------------------------------------
102
103if project.scenario == 'slide':
104    # Function for submarine slide
105    tsunami_source = slide_tsunami(length=35000.0,
106                                   depth=project.slide_depth,
107                                   slope=6.0,
108                                   thickness=500.0, 
109                                   x0=project.slide_origin[0], 
110                                   y0=project.slide_origin[1], 
111                                   alpha=0.0, 
112                                   domain=domain,
113                                   verbose=True)
114
115
116#------------------------------------------------------------------------------
117# Setup boundary conditions
118#------------------------------------------------------------------------------
119
120print 'Available boundary tags', domain.get_boundary_tags()
121
122
123Bd = Dirichlet_boundary([tide,0,0]) # Mean water level
124Bs = Transmissive_stage_zero_momentum_boundary(domain) # Neutral boundary
125
126
127if project.scenario == 'fixed_wave':
128    # Huge 50m wave starting after 60 seconds and lasting 1 hour.
129    Bw = Time_boundary(domain=domain,
130                       function=lambda t: [(60<t<3660)*50, 0, 0])
131                       
132    domain.set_boundary({'ocean_east': Bw,
133                         'bottom': Bs,
134                         'onshore': Bd,
135                         'top': Bs})
136
137if project.scenario == 'slide':
138    # Boundary conditions for slide scenario
139    domain.set_boundary({'ocean_east': Bd,
140                         'bottom': Bd,
141                         'onshore': Bd,
142                         'top': Bd})
143   
144
145#------------------------------------------------------------------------------
146# Evolve system through time
147#------------------------------------------------------------------------------
148
149import time
150t0 = time.time()
151
152from Numeric import allclose
153from anuga.abstract_2d_finite_volumes.quantity import Quantity
154
155if project.scenario == 'slide':
156   
157    for t in domain.evolve(yieldstep=10, finaltime=60): 
158        print domain.timestepping_statistics()
159        print domain.boundary_statistics(tags='ocean_east')       
160       
161    # Add slide
162    thisstagestep = domain.get_quantity('stage')
163    if allclose(t, 60):
164        slide = Quantity(domain)
165        slide.set_values(tsunami_source)
166        domain.set_quantity('stage', slide + thisstagestep)
167
168    for t in domain.evolve(yieldstep=10, finaltime=5000, 
169                           skip_initial_step = True):
170        print domain.timestepping_statistics()
171        print domain.boundary_statistics(tags='ocean_east')   
172
173       
174if project.scenario == 'fixed_wave':
175
176    # Save every two mins leading up to wave approaching land
177    for t in domain.evolve(yieldstep=120, finaltime=5000): 
178        print domain.timestepping_statistics()
179        print domain.boundary_statistics(tags='ocean_east')   
180
181    # Save every 30 secs as wave starts inundating ashore
182    for t in domain.evolve(yieldstep=10, finaltime=10000, 
183                           skip_initial_step=True):
184        print domain.timestepping_statistics()
185        print domain.boundary_statistics(tags='ocean_east')
186           
187print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.