1 | """Script for running a tsunami inundation scenario for Onslow, 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 simulated submarine landslide. |
---|
9 | |
---|
10 | Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006 |
---|
11 | """ |
---|
12 | #-------------------------------------------------------------------------------# Import necessary modules |
---|
13 | #------------------------------------------------------------------------------- |
---|
14 | |
---|
15 | # Standard modules |
---|
16 | from os import sep |
---|
17 | from os.path import dirname, basename |
---|
18 | import time |
---|
19 | |
---|
20 | # Related major packages |
---|
21 | from anuga.shallow_water import Domain, Reflective_boundary, \ |
---|
22 | Dirichlet_boundary, Time_boundary, File_boundary |
---|
23 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, \ |
---|
24 | dem2pts, ferret2sww |
---|
25 | from shutil import copy |
---|
26 | from os import mkdir, access, F_OK |
---|
27 | from anuga.geospatial_data.geospatial_data import * |
---|
28 | import sys |
---|
29 | from anuga.abstract_2d_finite_volumes.util import Screen_Catcher |
---|
30 | |
---|
31 | from anuga.fit_interpolate.fit import fit_to_mesh_file |
---|
32 | |
---|
33 | # Application specific imports |
---|
34 | import project # Definition of file names and polygons |
---|
35 | |
---|
36 | #------------------------------------------------------------------------------- |
---|
37 | # Copy scripts to time stamped output directory and capture screen |
---|
38 | # output to file |
---|
39 | #------------------------------------------------------------------------------- |
---|
40 | |
---|
41 | # creates copy of code in output dir if dir doesn't exist |
---|
42 | if access(project.outputtimedir,F_OK) == 0 : |
---|
43 | mkdir (project.outputtimedir) |
---|
44 | copy (dirname(project.__file__) +sep+ project.__name__+'.py', project.outputtimedir + project.__name__+'.py') |
---|
45 | copy (__file__, project.outputtimedir + basename(__file__)) |
---|
46 | print 'project.outputtimedir',project.outputtimedir |
---|
47 | |
---|
48 | # normal screen output is stored in |
---|
49 | #screen_output_name = project.outputtimedir + "screen_output.txt" |
---|
50 | #screen_error_name = project.outputtimedir + "screen_error.txt" |
---|
51 | |
---|
52 | # used to catch screen output to file |
---|
53 | #sys.stdout = Screen_Catcher(screen_output_name) |
---|
54 | #sys.stderr = Screen_Catcher(screen_error_name) |
---|
55 | print 'USER: ', project.user |
---|
56 | |
---|
57 | #------------------------------------------------------------------------------- |
---|
58 | # Preparation of topographic data |
---|
59 | # |
---|
60 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
61 | # Do for coarse and fine data |
---|
62 | # Fine pts file to be clipped to area of interest |
---|
63 | #------------------------------------------------------------------------------- |
---|
64 | |
---|
65 | # filenames |
---|
66 | meshname = project.meshname + '.tsh' |
---|
67 | mesh_elevname = project.mesh_elevname + '.tsh' |
---|
68 | source_dir = project.boundarydir |
---|
69 | |
---|
70 | |
---|
71 | #------------------------------------------------------------------------------- |
---|
72 | # Create the triangular mesh based on overall clipping polygon with a tagged |
---|
73 | # boundary and interior regions defined in project.py along with |
---|
74 | # resolutions (maximal area of per triangle) for each polygon |
---|
75 | #------------------------------------------------------------------------------- |
---|
76 | |
---|
77 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
78 | |
---|
79 | region_res = 5000000 |
---|
80 | coast_res = 500000 |
---|
81 | pt_hedland_res = 500000 |
---|
82 | interior_regions = [[project.poly_pt_hedland, pt_hedland_res], |
---|
83 | [project.poly_region, region_res]] |
---|
84 | |
---|
85 | print 'number of interior regions', len(interior_regions) |
---|
86 | |
---|
87 | print 'start create mesh from regions' |
---|
88 | from caching import cache |
---|
89 | _ = cache(create_mesh_from_regions, |
---|
90 | project.polyAll, |
---|
91 | {'boundary_tags': {'topright': [0], 'topleft': [1], |
---|
92 | 'left': [2], 'bottom0': [3], |
---|
93 | 'bottom1': [4], 'bottom2': [5], |
---|
94 | 'bottom3': [6], 'right': [7]}, |
---|
95 | 'maximum_triangle_area': 5000000, |
---|
96 | 'filename': meshname, |
---|
97 | 'interior_regions': interior_regions}, |
---|
98 | verbose = True, evaluate=True) |
---|
99 | |
---|
100 | cache(fit_to_mesh_file,(meshname, |
---|
101 | 'pt_hedland_combined_elevation_31204' + '.pts', |
---|
102 | mesh_elevname), |
---|
103 | {'verbose': True} |
---|
104 | #,evaluate = True |
---|
105 | ,verbose = False |
---|
106 | ) |
---|
107 | #------------------------------------------------------------------------------- |
---|
108 | # Setup computational domain |
---|
109 | #------------------------------------------------------------------------------- |
---|
110 | domain = Domain(mesh_elevname, use_cache = False, verbose = True) |
---|
111 | |
---|
112 | print domain.statistics() |
---|
113 | print 'Number of triangles = ', len(domain) |
---|
114 | print 'The extent is ', domain.get_extent() |
---|
115 | print domain.statistics() |
---|
116 | |
---|
117 | domain.set_name(project.basename) |
---|
118 | domain.set_datadir(project.outputtimedir) |
---|
119 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
120 | |
---|
121 | #------------------------------------------------------------------------------- |
---|
122 | # Setup initial conditions |
---|
123 | #------------------------------------------------------------------------------- |
---|
124 | |
---|
125 | tide = 0. |
---|
126 | #high |
---|
127 | #tide = 3.6 |
---|
128 | #low |
---|
129 | #tide = -3.9 |
---|
130 | |
---|
131 | domain.set_quantity('stage', tide) |
---|
132 | domain.set_quantity('friction', 0.0) |
---|
133 | print 'hi and file',project.combined_dem_name + '.pts' |
---|
134 | |
---|
135 | #domain.set_quantity('elevation', |
---|
136 | # filename = project.combined_dem_name + '.pts', |
---|
137 | # use_cache = True, |
---|
138 | # verbose = True, |
---|
139 | # alpha = 0.1 |
---|
140 | # ) |
---|
141 | |
---|
142 | #------------------------------------------------------------------------------- |
---|
143 | # Setup boundary conditions (all reflective) |
---|
144 | #------------------------------------------------------------------------------- |
---|
145 | print 'start ferret2sww' |
---|
146 | # skipped as results in file SU-AU_clipped is correct for all WA |
---|
147 | |
---|
148 | |
---|
149 | south = project.south |
---|
150 | north = project.north |
---|
151 | west = project.west |
---|
152 | east = project.east |
---|
153 | |
---|
154 | #note only need to do when an SWW file for the MOST boundary doesn't exist |
---|
155 | cache(ferret2sww, |
---|
156 | (project.boundary_basename, |
---|
157 | project.boundary_basename+'_'+project.basename), |
---|
158 | {'verbose': True, |
---|
159 | 'minlat': south, |
---|
160 | 'maxlat': north, |
---|
161 | 'minlon': west, |
---|
162 | 'maxlon': east, |
---|
163 | # 'origin': project.mesh_origin, |
---|
164 | 'origin': domain.geo_reference.get_origin(), |
---|
165 | 'mean_stage': tide, |
---|
166 | 'zscale': 10, #Enhance tsunami |
---|
167 | 'fail_on_NaN': False, |
---|
168 | 'inverted_bathymetry': True}, |
---|
169 | evaluate = True, |
---|
170 | verbose = True, |
---|
171 | dependencies = source_dir + project.boundary_basename + '.sww') |
---|
172 | |
---|
173 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
174 | |
---|
175 | Bf = File_boundary(project.boundary_basename+'_'+project.basename + '.sww', |
---|
176 | domain, verbose = True) |
---|
177 | Br = Reflective_boundary(domain) |
---|
178 | Bd = Dirichlet_boundary([tide,0,0]) |
---|
179 | domain.set_boundary( {'topright': Bf,'topleft': Bf, 'left': Bd, 'bottom0': Bd, |
---|
180 | 'bottom1': Bd, 'bottom2': Bd, 'bottom3': Bd, |
---|
181 | 'right': Bd}) |
---|
182 | |
---|
183 | #------------------------------------------------------------------------------- |
---|
184 | # Evolve system through time |
---|
185 | #------------------------------------------------------------------------------- |
---|
186 | import time |
---|
187 | t0 = time.time() |
---|
188 | |
---|
189 | for t in domain.evolve(yieldstep = 240, finaltime = 10800): |
---|
190 | domain.write_time() |
---|
191 | domain.write_boundary_statistics(tags = 'topright') |
---|
192 | |
---|
193 | for t in domain.evolve(yieldstep = 120, finaltime = 16200 |
---|
194 | ,skip_initial_step = True): |
---|
195 | domain.write_time() |
---|
196 | domain.write_boundary_statistics(tags = 'topright') |
---|
197 | |
---|
198 | for t in domain.evolve(yieldstep = 60, finaltime = 21600 |
---|
199 | ,skip_initial_step = True): |
---|
200 | domain.write_time() |
---|
201 | domain.write_boundary_statistics(tags = 'topright') |
---|
202 | |
---|
203 | for t in domain.evolve(yieldstep = 120, finaltime = 27000 |
---|
204 | ,skip_initial_step = True): |
---|
205 | domain.write_time() |
---|
206 | domain.write_boundary_statistics(tags = 'topright') |
---|
207 | |
---|
208 | for t in domain.evolve(yieldstep = 240, finaltime = 36000 |
---|
209 | ,skip_initial_step = True): |
---|
210 | domain.write_time() |
---|
211 | domain.write_boundary_statistics(tags = 'topright') |
---|
212 | |
---|
213 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
214 | |
---|
215 | print 'finished' |
---|