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 | Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton and |
---|
12 | Nick Bartzis, GA - 2006 |
---|
13 | """ |
---|
14 | |
---|
15 | #------------------------------------------------------------------------------ |
---|
16 | # Import necessary modules |
---|
17 | #------------------------------------------------------------------------------ |
---|
18 | |
---|
19 | # Standard modules |
---|
20 | import os |
---|
21 | import time |
---|
22 | import sys |
---|
23 | |
---|
24 | # Related major packages |
---|
25 | from anuga.shallow_water import Domain |
---|
26 | from anuga.shallow_water import Reflective_boundary |
---|
27 | from anuga.shallow_water import Dirichlet_boundary |
---|
28 | from anuga.shallow_water import Time_boundary |
---|
29 | from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary |
---|
30 | from anuga.abstract_2d_finite_volumes.util import file_function |
---|
31 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
32 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf |
---|
33 | from anuga.shallow_water.data_manager import dem2pts |
---|
34 | from anuga.geospatial_data.geospatial_data import Geospatial_data |
---|
35 | from anuga.shallow_water.shallow_water_domain import Inflow |
---|
36 | |
---|
37 | # Application specific imports |
---|
38 | import project # Definition of file names and polygons |
---|
39 | |
---|
40 | |
---|
41 | #------------------------------------------------------------------------------ |
---|
42 | # Define scenario as either slide or fixed_wave. |
---|
43 | #------------------------------------------------------------------------------ |
---|
44 | scenario = 'slide' |
---|
45 | #scenario = 'fixed_wave' |
---|
46 | |
---|
47 | if os.access(scenario, os.F_OK) == 0: |
---|
48 | os.mkdir(scenario) |
---|
49 | basename = scenario + 'source' |
---|
50 | |
---|
51 | |
---|
52 | #------------------------------------------------------------------------------ |
---|
53 | # Preparation of topographic data |
---|
54 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
55 | #------------------------------------------------------------------------------ |
---|
56 | |
---|
57 | # Filenames |
---|
58 | basename = 'blackbutt' |
---|
59 | meshname = basename+'.msh' |
---|
60 | |
---|
61 | elevation_file = 'ALL_Grd_Bui_Comb_Smth_Rds_XYZPts_only' |
---|
62 | # Create pts file |
---|
63 | #G = Geospatial_data(elevation_file+'.csv') |
---|
64 | #G.export_points_file(basename+'.pts') |
---|
65 | |
---|
66 | #------------------------------------------------------------------------------ |
---|
67 | # Create the triangular mesh based on overall clipping polygon with a tagged |
---|
68 | # boundary and interior regions defined in project.py along with |
---|
69 | # resolutions (maximal area of per triangle) for each polygon |
---|
70 | #------------------------------------------------------------------------------ |
---|
71 | |
---|
72 | |
---|
73 | W = 302500 |
---|
74 | S = 6171740 |
---|
75 | N = 6172270 |
---|
76 | E = 304000 |
---|
77 | |
---|
78 | bounding_polygon = [[W, S], [E, S], [E, N], [W, N]] |
---|
79 | |
---|
80 | channel_polygon = [[W, S+100], [E, S+100], [E, N-100], [W, N-100]] |
---|
81 | |
---|
82 | interior_regions = [ [channel_polygon, 100] ] # 100 m^2 |
---|
83 | create_mesh_from_regions(bounding_polygon, |
---|
84 | boundary_tags={'south': [0], |
---|
85 | 'east': [1], |
---|
86 | 'north': [2], |
---|
87 | 'west': [3]}, |
---|
88 | maximum_triangle_area=400, |
---|
89 | interior_regions=interior_regions, |
---|
90 | filename=meshname, |
---|
91 | use_cache=False, |
---|
92 | verbose=True) |
---|
93 | |
---|
94 | |
---|
95 | #------------------------------------------------------------------------------ |
---|
96 | # Setup computational domain |
---|
97 | #------------------------------------------------------------------------------ |
---|
98 | |
---|
99 | domain = Domain(meshname, use_cache=True, verbose=True) |
---|
100 | |
---|
101 | print domain.statistics() |
---|
102 | domain.set_name(basename) |
---|
103 | |
---|
104 | |
---|
105 | #------------------------------------------------------------------------------ |
---|
106 | # Setup initial conditions |
---|
107 | #------------------------------------------------------------------------------ |
---|
108 | |
---|
109 | tide = -100 |
---|
110 | domain.set_quantity('stage', tide) |
---|
111 | domain.set_quantity('friction', 0.0, location = 'centroids') |
---|
112 | domain.set_quantity('elevation', |
---|
113 | filename=basename+'.pts', |
---|
114 | use_cache=True, |
---|
115 | verbose=True, |
---|
116 | alpha=0.01) |
---|
117 | |
---|
118 | |
---|
119 | |
---|
120 | #------------------------------------------------------------------------------ |
---|
121 | # Setup specialised forcing terms |
---|
122 | #------------------------------------------------------------------------------ |
---|
123 | |
---|
124 | hydrograph = Inflow(center=(300, 300), radius=50, |
---|
125 | flow=file_function('Q/QPMF_Rot_Sub13.tms')) |
---|
126 | |
---|
127 | domain.forcing_terms.append(hydrograph) |
---|
128 | |
---|
129 | |
---|
130 | #------------------------------------------------------------------------------ |
---|
131 | # Setup boundary conditions |
---|
132 | #------------------------------------------------------------------------------ |
---|
133 | |
---|
134 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
135 | |
---|
136 | Br = Reflective_boundary(domain) |
---|
137 | Bd = Dirichlet_boundary([tide,0,0]) |
---|
138 | |
---|
139 | |
---|
140 | # boundary conditions for slide scenario |
---|
141 | domain.set_boundary({'west': Bd, |
---|
142 | 'south': Bd, |
---|
143 | 'east': Transmissive_Momentum_Set_Stage_boundary(domain, lambda t: tide), |
---|
144 | 'north': Bd}) |
---|
145 | |
---|
146 | |
---|
147 | #------------------------------------------------------------------------------ |
---|
148 | # Evolve system through time |
---|
149 | #------------------------------------------------------------------------------ |
---|
150 | |
---|
151 | import time |
---|
152 | t0 = time.time() |
---|
153 | |
---|
154 | for t in domain.evolve(yieldstep = 1, finaltime = 3000): |
---|
155 | domain.write_time() |
---|
156 | #domain.write_boundary_statistics(tags='east', quantities='stage') |
---|
157 | |
---|
158 | if t > 2538: |
---|
159 | hydrograph.flow = 0.0 |
---|