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

Last change on this file since 7078 was 7078, checked in by ole, 15 years ago

Typo

File size: 7.1 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    # Initial run without any event
138    for t in domain.evolve(yieldstep=10, finaltime=60): 
139        print domain.timestepping_statistics()
140        print domain.boundary_statistics(tags='ocean_east')       
141       
142    # Add slide to water surface
143    if allclose(t, 60):
144        domain.add_quantity('stage', tsunami_source)   
145
146    # Continue propagating wave
147    for t in domain.evolve(yieldstep=10, finaltime=5000, 
148                           skip_initial_step = True):
149        print domain.timestepping_statistics()
150        print domain.boundary_statistics(tags='ocean_east')   
151
152if project.scenario == 'fixed_wave':
153    # Save every two mins leading up to wave approaching land
154    for t in domain.evolve(yieldstep=120, finaltime=5000): 
155        print domain.timestepping_statistics()
156        print domain.boundary_statistics(tags='ocean_east')   
157
158    # Save every 30 secs as wave starts inundating ashore
159    for t in domain.evolve(yieldstep=10, finaltime=10000, 
160                           skip_initial_step=True):
161        print domain.timestepping_statistics()
162        print domain.boundary_statistics(tags='ocean_east')
163           
164print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.