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

Last change on this file since 7064 was 7064, checked in by rwilson, 11 years ago

Fiddling with layout of user guide.

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