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

Last change on this file since 7816 was 7816, checked in by hudson, 12 years ago

buildings demo added

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