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

Last change on this file since 4003 was 4003, checked in by sexton, 17 years ago

including data and supporting files for Cairns demo + updating scripts

File size: 8.0 KB
Line 
1"""Script for running a tsunami inundation scenario for Broome, WA, 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 project.outputtimedir
6
7The scenario is defined by a triangular mesh created from project.polygon,
8the elevation data and a tsunami wave generated by MOST.
9
10Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
11"""
12
13#-------------------------------------------------------------------------------
14# Import necessary modules
15#-------------------------------------------------------------------------------
16
17# Standard modules
18import os
19import time
20from shutil import copy
21from os import mkdir, access, F_OK
22import sys
23
24# Related major packages
25from anuga.shallow_water import Domain, Reflective_boundary, \
26                            Dirichlet_boundary, Time_boundary, File_boundary
27from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
28from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files
29from anuga.geospatial_data.geospatial_data import *
30
31# Application specific imports
32import project                 # Definition of file names and polygons
33
34#-------------------------------------------------------------------------------
35# Define scenario as either slump or fixed_wave.
36#-------------------------------------------------------------------------------
37scenario = 'slump' # OR 'fixed_wave'
38if access(scenario,F_OK) == 0:
39    mkdir (scenario)
40basename = scenario + 'source'
41
42#-------------------------------------------------------------------------------
43# Preparation of topographic data
44#
45# Convert ASC 2 DEM 2 PTS using source data and store result in source data
46#-------------------------------------------------------------------------------
47
48# filenames
49dem_name = 'cairns' 
50meshname = 'cairns.msh'
51
52# creates DEM from asc data
53convert_dem_from_ascii2netcdf(dem_name, use_cache=True, verbose=True)
54
55#creates pts file for onshore DEM
56dem2pts(dem_name, use_cache=True, verbose=True)
57
58#----------------------------------------------------------------------------
59# Create the triangular mesh based on overall clipping polygon with a tagged
60# boundary and interior regions defined in project.py along with
61# resolutions (maximal area of per triangle) for each polygon
62#-------------------------------------------------------------------------------
63
64from anuga.pmesh.mesh_interface import create_mesh_from_regions
65remainder_res = 10000000
66islands_res = 100000
67cairns_res = 100000
68shallow_res = 500000
69interior_regions = [[project.poly_cairns, cairns_res],
70                    [project.poly_island0, islands_res],
71                    [project.poly_island1, islands_res],
72                    [project.poly_island2, islands_res],
73                    [project.poly_island3, islands_res],
74                    [project.poly_shallow, shallow_res]]
75
76from caching import cache
77_ = cache(create_mesh_from_regions,
78          project.polyAll,
79           {'boundary_tags': {'top': [0], 'ocean_east': [1],
80                              'bottom': [2], 'onshore': [3]},
81           'maximum_triangle_area': remainder_res,
82           'filename': meshname,
83           'interior_regions': interior_regions},
84          verbose = True, evaluate=False)
85print 'created mesh'
86
87#-------------------------------------------------------------------------------                                 
88# Setup computational domain
89#-------------------------------------------------------------------------------                                 
90domain = Domain(meshname, use_cache = True, verbose = True)
91
92print 'Number of triangles = ', len(domain)
93print 'The extent is ', domain.get_extent()
94print domain.statistics()
95
96domain.set_name(basename) 
97domain.set_datadir(scenario)
98domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
99domain.set_minimum_storable_height(0.01)
100
101#-------------------------------------------------------------------------------                                 
102# Setup initial conditions
103#-------------------------------------------------------------------------------
104
105tide = 0.0
106domain.set_quantity('stage', tide)
107domain.set_quantity('friction', 0.0) 
108domain.set_quantity('elevation', 
109                    filename = dem_name + '.pts',
110                    use_cache = True,
111                    verbose = True,
112                    alpha = 0.1
113                    )
114
115#-------------------------------------------------------------------------------                                 
116# Setup information for slump scenario (to be applied 1 min into simulation
117#-------------------------------------------------------------------------------
118if scenario == 'slump':
119    from anuga.shallow_water.smf import slump_tsunami  # Function for submarine slump
120    tsunami_source = slump_tsunami(length=15000.0,
121                                   depth=project.slump_depth,
122                                   slope=6.0,
123                                   thickness=250.0, 
124                                   radius=3330,
125                                   dphi=0.23,
126                                   x0=project.slump_origin[0], 
127                                   y0=project.slump_origin[1], 
128                                   alpha=0.0, 
129                                   domain=domain,
130                                   verbose=True)
131
132
133#-------------------------------------------------------------------------------                                 
134# Setup boundary conditions
135#-------------------------------------------------------------------------------
136print 'Available boundary tags', domain.get_boundary_tags()
137
138Br = Reflective_boundary(domain)
139Bd = Dirichlet_boundary([tide,0,0])
140
141# 7 min square wave starting at 1 min, 10m high
142if scenario == 'fixed_wave':
143    Bw = Time_boundary(domain = domain,
144                       f=lambda t: [(60<t<480)*10, 0, 0])
145    domain.set_boundary( {'ocean_east': Bw, 'bottom': Bd, 'onshore': Bd,
146                           'top': Bd} )
147
148# boundary conditions for slump scenario
149if scenario == 'slump':
150    domain.set_boundary( {'ocean_east': Bd, 'bottom': Bd, 'onshore': Bd,
151                          'top': Bd} )
152
153
154#-------------------------------------------------------------------------------                                 
155# Evolve system through time
156#-------------------------------------------------------------------------------
157import time
158t0 = time.time()
159
160from Numeric import allclose
161from anuga.abstract_2d_finite_volumes.quantity import Quantity
162
163if scenario == 'slump':
164   
165    for t in domain.evolve(yieldstep = 30, finaltime = 60): 
166        domain.write_time()
167        domain.write_boundary_statistics(tags = 'ocean_east')     
168       
169        # add slump
170        thisstagestep = domain.get_quantity('stage')
171        if allclose(t, 60):
172            slump = Quantity(domain)
173            slump.set_values(tsunami_source)
174            domain.set_quantity('stage', slump + thisstagestep)
175
176        # save every two mins leading up to wave approaching land
177        for t in domain.evolve(yieldstep = 120, finaltime = 1000, 
178                               skip_initial_step = True):
179            domain.write_time()
180            domain.write_boundary_statistics(tags = 'ocean_east')
181
182        # save every 30 secs as wave starts inundating ashore
183        for t in domain.evolve(yieldstep = 30, finaltime = 10000, 
184                               skip_initial_step = True):
185            domain.write_time()
186            domain.write_boundary_statistics(tags = 'ocean_east')
187
188if scenario == 'fixed_wave':
189
190    # save every two mins leading up to wave approaching land
191    for t in domain.evolve(yieldstep = 120, finaltime = 1000): 
192        domain.write_time()
193        domain.write_boundary_statistics(tags = 'ocean_east')     
194
195        # save every 30 secs as wave starts inundating ashore
196        for t in domain.evolve(yieldstep = 30, finaltime = 10000, 
197                               skip_initial_step = True):
198            domain.write_time()
199            domain.write_boundary_statistics(tags = 'ocean_east')
200           
201print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.