1 | """Script for running a tsunami inundation scenario for Broome, WA, Australia. |
---|
2 | |
---|
3 | Source data such as elevation and boundary data is assumed to be available in |
---|
4 | directories specified by project.py |
---|
5 | The output sww file is stored in project.outputtimedir |
---|
6 | |
---|
7 | The scenario is defined by a triangular mesh created from project.polygon, |
---|
8 | the elevation data and a tsunami wave generated by MOST. |
---|
9 | |
---|
10 | Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006 |
---|
11 | """ |
---|
12 | |
---|
13 | #------------------------------------------------------------------------------- |
---|
14 | # Import necessary modules |
---|
15 | #------------------------------------------------------------------------------- |
---|
16 | |
---|
17 | # Standard modules |
---|
18 | import os |
---|
19 | import time |
---|
20 | from shutil import copy |
---|
21 | from os import mkdir, access, F_OK |
---|
22 | import sys |
---|
23 | |
---|
24 | # Related major packages |
---|
25 | from anuga.shallow_water import Domain, Reflective_boundary, \ |
---|
26 | Dirichlet_boundary, Time_boundary, File_boundary |
---|
27 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts |
---|
28 | from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files |
---|
29 | from anuga.geospatial_data.geospatial_data import * |
---|
30 | |
---|
31 | # Application specific imports |
---|
32 | import project # Definition of file names and polygons |
---|
33 | |
---|
34 | #------------------------------------------------------------------------------- |
---|
35 | # Define scenario as either slump or fixed_wave. |
---|
36 | #------------------------------------------------------------------------------- |
---|
37 | scenario = 'slump' # OR 'fixed_wave' |
---|
38 | if access(scenario,F_OK) == 0: |
---|
39 | mkdir (scenario) |
---|
40 | basename = 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 |
---|
49 | dem_name = 'cairns' |
---|
50 | meshname = 'cairns.msh' |
---|
51 | |
---|
52 | # creates DEM from asc data |
---|
53 | convert_dem_from_ascii2netcdf(dem_name, use_cache=True, verbose=True) |
---|
54 | |
---|
55 | #creates pts file for onshore DEM |
---|
56 | dem2pts(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 | |
---|
64 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
65 | remainder_res = 10000000 |
---|
66 | islands_res = 100000 |
---|
67 | cairns_res = 100000 |
---|
68 | shallow_res = 500000 |
---|
69 | interior_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 | |
---|
76 | from 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) |
---|
85 | print 'created mesh' |
---|
86 | |
---|
87 | #------------------------------------------------------------------------------- |
---|
88 | # Setup computational domain |
---|
89 | #------------------------------------------------------------------------------- |
---|
90 | domain = Domain(meshname, use_cache = True, verbose = True) |
---|
91 | |
---|
92 | print 'Number of triangles = ', len(domain) |
---|
93 | print 'The extent is ', domain.get_extent() |
---|
94 | print domain.statistics() |
---|
95 | |
---|
96 | domain.set_name(basename) |
---|
97 | domain.set_datadir(scenario) |
---|
98 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
99 | domain.set_minimum_storable_height(0.01) |
---|
100 | |
---|
101 | #------------------------------------------------------------------------------- |
---|
102 | # Setup initial conditions |
---|
103 | #------------------------------------------------------------------------------- |
---|
104 | |
---|
105 | tide = 0.0 |
---|
106 | domain.set_quantity('stage', tide) |
---|
107 | domain.set_quantity('friction', 0.0) |
---|
108 | domain.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 | #------------------------------------------------------------------------------- |
---|
118 | if 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 | #------------------------------------------------------------------------------- |
---|
136 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
137 | |
---|
138 | Br = Reflective_boundary(domain) |
---|
139 | Bd = Dirichlet_boundary([tide,0,0]) |
---|
140 | |
---|
141 | # 7 min square wave starting at 1 min, 10m high |
---|
142 | if 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 |
---|
149 | if 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 | #------------------------------------------------------------------------------- |
---|
157 | import time |
---|
158 | t0 = time.time() |
---|
159 | |
---|
160 | from Numeric import allclose |
---|
161 | from anuga.abstract_2d_finite_volumes.quantity import Quantity |
---|
162 | |
---|
163 | if 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 | |
---|
188 | if 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 | |
---|
201 | print 'That took %.2f seconds' %(time.time()-t0) |
---|