1 | """Script for running a tsunami inundation scenario for Cairns, QLD 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 directory named after the scenario, i.e |
---|
6 | slide or fixed_wave. |
---|
7 | |
---|
8 | The scenario is defined by a triangular mesh created from project.polygon, |
---|
9 | the elevation data and a tsunami wave generated by a submarine mass failure. |
---|
10 | |
---|
11 | Geoscience Australia, 2004-present |
---|
12 | """ |
---|
13 | |
---|
14 | #------------------------------------------------------------------------------ |
---|
15 | # Import necessary modules |
---|
16 | #------------------------------------------------------------------------------ |
---|
17 | # Standard modules |
---|
18 | import os |
---|
19 | import time |
---|
20 | import sys |
---|
21 | |
---|
22 | # Related major packages |
---|
23 | import anuga |
---|
24 | |
---|
25 | # Application specific imports |
---|
26 | import 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 |
---|
33 | anuga.asc2dem(project.name_stem+'.asc', use_cache=True, verbose=True) |
---|
34 | |
---|
35 | # Create pts file for onshore DEM |
---|
36 | anuga.dem2pts(project.name_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 | #------------------------------------------------------------------------------ |
---|
43 | domain = 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 |
---|
55 | print 'Number of triangles = ', len(domain) |
---|
56 | print 'The extent is ', domain.get_extent() |
---|
57 | print domain.statistics() |
---|
58 | |
---|
59 | #------------------------------------------------------------------------------ |
---|
60 | # Setup parameters of computational domain |
---|
61 | #------------------------------------------------------------------------------ |
---|
62 | domain.set_name('cairns_' + project.scenario) # Name of sww file |
---|
63 | domain.set_datadir('.') # Store sww output here |
---|
64 | domain.set_minimum_storable_height(0.01) # Store only depth > 1cm |
---|
65 | |
---|
66 | #------------------------------------------------------------------------------ |
---|
67 | # Setup initial conditions |
---|
68 | #------------------------------------------------------------------------------ |
---|
69 | tide = 0.0 |
---|
70 | domain.set_quantity('stage', tide) |
---|
71 | domain.set_quantity('friction', 0.0) |
---|
72 | domain.set_quantity('elevation', |
---|
73 | filename=project.name_stem + '.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 | #------------------------------------------------------------------------------ |
---|
81 | if 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 | #------------------------------------------------------------------------------ |
---|
96 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
97 | |
---|
98 | Bd = anuga.Dirichlet_boundary([tide, 0, 0]) # Mean water level |
---|
99 | Bs = anuga.Transmissive_stage_zero_momentum_boundary(domain) # Neutral boundary |
---|
100 | |
---|
101 | if 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 | |
---|
110 | if 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 | #------------------------------------------------------------------------------ |
---|
120 | import time |
---|
121 | t0 = time.time() |
---|
122 | |
---|
123 | from numpy import allclose |
---|
124 | |
---|
125 | if 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 | |
---|
141 | if 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 | |
---|
153 | print 'That took %.2f seconds' %(time.time()-t0) |
---|