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 tsunami generated by a subduction zone earthquake. |
---|
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 |
---|
25 | from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files |
---|
26 | from shutil import copy |
---|
27 | from os import mkdir, access, F_OK |
---|
28 | from anuga.geospatial_data.geospatial_data import * |
---|
29 | import sys |
---|
30 | from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files |
---|
31 | |
---|
32 | # Application specific imports |
---|
33 | import project # Definition of file names and polygons |
---|
34 | |
---|
35 | #------------------------------------------------------------------------------- |
---|
36 | # Copy scripts to time stamped output directory and capture screen |
---|
37 | # output to file |
---|
38 | #------------------------------------------------------------------------------- |
---|
39 | |
---|
40 | copy_code_files(project.outputtimedir,__file__,dirname(project.__file__)+sep+ project.__name__+'.py' ) |
---|
41 | myid = 0 |
---|
42 | numprocs = 1 |
---|
43 | start_screen_catcher(project.outputtimedir, myid, numprocs) |
---|
44 | |
---|
45 | print 'USER: ', project.user |
---|
46 | |
---|
47 | #------------------------------------------------------------------------------- |
---|
48 | # Preparation of topographic data |
---|
49 | # |
---|
50 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
51 | # Do for coarse and fine data |
---|
52 | # Fine pts file to be clipped to area of interest |
---|
53 | #------------------------------------------------------------------------------- |
---|
54 | |
---|
55 | # filenames |
---|
56 | meshname = project.meshname+'.msh' |
---|
57 | source_dir = project.boundarydir |
---|
58 | |
---|
59 | # fine data (clipping the points file to smaller area) |
---|
60 | # creates DEM from asc data |
---|
61 | convert_dem_from_ascii2netcdf(project.onshore_dem_name, use_cache=True, verbose=True) |
---|
62 | |
---|
63 | #creates pts file from DEM |
---|
64 | dem2pts(project.onshore_dem_name, |
---|
65 | easting_min=project.eastingmin, |
---|
66 | easting_max=project.eastingmax, |
---|
67 | northing_min=project.northingmin, |
---|
68 | northing_max= project.northingmax, |
---|
69 | use_cache=True, |
---|
70 | verbose=True) |
---|
71 | |
---|
72 | print 'create offshore' |
---|
73 | G11= Geospatial_data(file_name = project.offshore_dem_name0 + '.xya')+\ |
---|
74 | Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')+\ |
---|
75 | Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')+\ |
---|
76 | Geospatial_data(file_name = project.offshore_dem_name3 + '.xya')+\ |
---|
77 | Geospatial_data(file_name = project.offshore_dem_name4 + '.xya')+\ |
---|
78 | Geospatial_data(file_name = project.offshore_dem_name5 + '.xya')+\ |
---|
79 | Geospatial_data(file_name = project.offshore_dem_name6 + '.xya')+\ |
---|
80 | Geospatial_data(file_name = project.offshore_dem_name7 + '.xya')+\ |
---|
81 | Geospatial_data(file_name = project.offshore_dem_name8 + '.xya')+\ |
---|
82 | Geospatial_data(file_name = project.offshore_dem_name9 + '.xya')+\ |
---|
83 | Geospatial_data(file_name = project.offshore_dem_name10 + '.xya') |
---|
84 | G12= Geospatial_data(file_name = project.offshore_dem_name11 + '.xya')+\ |
---|
85 | Geospatial_data(file_name = project.offshore_dem_name12 + '.xya')+\ |
---|
86 | Geospatial_data(file_name = project.offshore_dem_name13 + '.xya')+\ |
---|
87 | Geospatial_data(file_name = project.offshore_dem_name14 + '.xya')+\ |
---|
88 | Geospatial_data(file_name = project.offshore_dem_name15 + '.xya')+\ |
---|
89 | Geospatial_data(file_name = project.offshore_dem_name16 + '.xya')+\ |
---|
90 | Geospatial_data(file_name = project.offshore_dem_name17 + '.xya')+\ |
---|
91 | Geospatial_data(file_name = project.offshore_dem_name18 + '.xya')+\ |
---|
92 | Geospatial_data(file_name = project.offshore_dem_name19 + '.xya')+\ |
---|
93 | Geospatial_data(file_name = project.offshore_dem_name20 + '.xya') |
---|
94 | G13= Geospatial_data(file_name = project.offshore_dem_name21 + '.xya')+\ |
---|
95 | Geospatial_data(file_name = project.offshore_dem_name22 + '.xya')+\ |
---|
96 | Geospatial_data(file_name = project.offshore_dem_name23 + '.xya')+\ |
---|
97 | Geospatial_data(file_name = project.offshore_dem_name24 + '.xya')+\ |
---|
98 | Geospatial_data(file_name = project.offshore_dem_name25 + '.xya')+\ |
---|
99 | Geospatial_data(file_name = project.offshore_dem_name26 + '.xya')+\ |
---|
100 | Geospatial_data(file_name = project.offshore_dem_name27 + '.xya')+\ |
---|
101 | Geospatial_data(file_name = project.offshore_dem_name28 + '.xya')+\ |
---|
102 | Geospatial_data(file_name = project.offshore_dem_name29 + '.xya') |
---|
103 | G14= Geospatial_data(file_name = project.offshore_dem_name30 + '.xya')+\ |
---|
104 | Geospatial_data(file_name = project.offshore_dem_name31 + '.xya')+\ |
---|
105 | Geospatial_data(file_name = project.offshore_dem_name32 + '.xya')+\ |
---|
106 | Geospatial_data(file_name = project.offshore_dem_name33 + '.xya')+\ |
---|
107 | Geospatial_data(file_name = project.offshore_dem_name34 + '.xya')+\ |
---|
108 | Geospatial_data(file_name = project.offshore_dem_name35 + '.xya')+\ |
---|
109 | Geospatial_data(file_name = project.offshore_dem_name36 + '.xya')+\ |
---|
110 | Geospatial_data(file_name = project.offshore_dem_name37 + '.xya')+\ |
---|
111 | Geospatial_data(file_name = project.offshore_dem_name38 + '.xya')+\ |
---|
112 | Geospatial_data(file_name = project.offshore_dem_name39 + '.xya') |
---|
113 | G15= Geospatial_data(file_name = project.offshore_dem_name40 + '.xya')+\ |
---|
114 | Geospatial_data(file_name = project.offshore_dem_name41 + '.xya')+\ |
---|
115 | Geospatial_data(file_name = project.offshore_dem_name42 + '.xya')+\ |
---|
116 | Geospatial_data(file_name = project.offshore_dem_name43 + '.xya')+\ |
---|
117 | Geospatial_data(file_name = project.offshore_dem_name44 + '.xya')+\ |
---|
118 | Geospatial_data(file_name = project.offshore_dem_name45 + '.xya')+\ |
---|
119 | Geospatial_data(file_name = project.offshore_dem_name46 + '.xya')+\ |
---|
120 | Geospatial_data(file_name = project.offshore_dem_name47 + '.xya')+\ |
---|
121 | Geospatial_data(file_name = project.offshore_dem_name48 + '.xya')+\ |
---|
122 | Geospatial_data(file_name = project.offshore_dem_name49 + '.xya')+\ |
---|
123 | Geospatial_data(file_name = project.offshore_dem_name50 + '.xya') |
---|
124 | |
---|
125 | print 'create onshore' |
---|
126 | G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts') |
---|
127 | print 'create coast' |
---|
128 | G3 = Geospatial_data(file_name = project.coast_dem_name + '.xya') |
---|
129 | print 'add' |
---|
130 | G = G11 + G12 + G13 + G14 + G15 + G2 + G3 |
---|
131 | print 'export G' |
---|
132 | G.export_points_file(project.combined_dem_name + '.pts') |
---|
133 | |
---|
134 | #------------------------------------------------------------------------------- |
---|
135 | # Create the triangular mesh based on overall clipping polygon with a tagged |
---|
136 | # boundary and interior regions defined in project.py along with |
---|
137 | # resolutions (maximal area of per triangle) for each polygon |
---|
138 | #------------------------------------------------------------------------------- |
---|
139 | |
---|
140 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
141 | |
---|
142 | region_res = 500000 |
---|
143 | coast_res = 500 |
---|
144 | pt_hedland_res = 5000 |
---|
145 | interior_regions = [[project.poly_pt_hedland, pt_hedland_res], |
---|
146 | [project.poly_region, region_res]] |
---|
147 | |
---|
148 | print 'number of interior regions', len(interior_regions) |
---|
149 | |
---|
150 | from anuga.utilities.polygon import plot_polygons |
---|
151 | if sys.platform == 'win32': |
---|
152 | #figname = project.outputtimedir + 'pt_hedland_polys' |
---|
153 | figname = 'pt_hedland_polys_test' |
---|
154 | plot_polygons([project.polyAll,project.poly_pt_hedland,project.poly_region], |
---|
155 | figname, |
---|
156 | verbose = True) |
---|
157 | |
---|
158 | print 'start create mesh from regions' |
---|
159 | from caching import cache |
---|
160 | _ = cache(create_mesh_from_regions, |
---|
161 | project.polyAll, |
---|
162 | {'boundary_tags': {'topright': [0], 'topleft': [1], |
---|
163 | 'left': [2], 'bottom0': [3], |
---|
164 | 'bottom1': [4], 'bottom2': [5], |
---|
165 | 'bottom3': [6], 'right': [7]}, |
---|
166 | 'maximum_triangle_area': 250000, |
---|
167 | 'filename': meshname, |
---|
168 | 'interior_regions': interior_regions}, |
---|
169 | verbose = True, evaluate=True) |
---|
170 | |
---|
171 | #------------------------------------------------------------------------------- |
---|
172 | # Setup computational domain |
---|
173 | #------------------------------------------------------------------------------- |
---|
174 | domain = Domain(meshname, use_cache = False, verbose = True) |
---|
175 | |
---|
176 | print domain.statistics() |
---|
177 | print 'Number of triangles = ', len(domain) |
---|
178 | print 'The extent is ', domain.get_extent() |
---|
179 | print domain.statistics() |
---|
180 | |
---|
181 | domain.set_name(project.basename) |
---|
182 | domain.set_datadir(project.outputtimedir) |
---|
183 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
184 | |
---|
185 | #------------------------------------------------------------------------------- |
---|
186 | # Setup initial conditions |
---|
187 | #------------------------------------------------------------------------------- |
---|
188 | |
---|
189 | tide = 0. |
---|
190 | #high |
---|
191 | #tide = 3.6 |
---|
192 | #low |
---|
193 | #tide = -3.9 |
---|
194 | |
---|
195 | domain.set_quantity('stage', tide) |
---|
196 | domain.set_quantity('friction', 0.0) |
---|
197 | print 'hi and file',project.combined_dem_name + '.pts' |
---|
198 | |
---|
199 | domain.set_quantity('elevation', |
---|
200 | filename = project.combined_dem_name + '.pts', |
---|
201 | use_cache = True, |
---|
202 | verbose = True, |
---|
203 | alpha = 0.1 |
---|
204 | ) |
---|
205 | |
---|
206 | #------------------------------------------------------------------------------- |
---|
207 | # Setup boundary conditions (all reflective) |
---|
208 | #------------------------------------------------------------------------------- |
---|
209 | print 'start ferret2sww' |
---|
210 | # skipped as results in file SU-AU_clipped is correct for all WA |
---|
211 | |
---|
212 | from anuga.pyvolution.data_manager import ferret2sww |
---|
213 | |
---|
214 | south = project.south |
---|
215 | north = project.north |
---|
216 | west = project.west |
---|
217 | east = project.east |
---|
218 | |
---|
219 | #note only need to do when an SWW file for the MOST boundary doesn't exist |
---|
220 | cache(ferret2sww, |
---|
221 | (source_dir + project.boundary_basename, |
---|
222 | source_dir + project.boundary_basename+'_'+project.basename), |
---|
223 | {'verbose': True, |
---|
224 | 'minlat': south, |
---|
225 | 'maxlat': north, |
---|
226 | 'minlon': west, |
---|
227 | 'maxlon': east, |
---|
228 | # 'origin': project.mesh_origin, |
---|
229 | 'origin': domain.geo_reference.get_origin(), |
---|
230 | 'mean_stage': tide, |
---|
231 | 'zscale': 1, #Enhance tsunami |
---|
232 | 'fail_on_NaN': False, |
---|
233 | 'inverted_bathymetry': True}, |
---|
234 | evaluate = True, |
---|
235 | verbose = True, |
---|
236 | dependencies = source_dir + project.boundary_basename + '.sww') |
---|
237 | |
---|
238 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
239 | |
---|
240 | Bf = File_boundary(source_dir + project.boundary_basename + '.sww', |
---|
241 | domain, verbose = True) |
---|
242 | Br = Reflective_boundary(domain) |
---|
243 | Bd = Dirichlet_boundary([tide,0,0]) |
---|
244 | domain.set_boundary( {'topright': Bf,'topleft': Bf, 'left': Bd, 'bottom0': Bd, |
---|
245 | 'bottom1': Bd, 'bottom2': Bd, 'bottom3': Bd, |
---|
246 | 'right': Bd}) |
---|
247 | |
---|
248 | #------------------------------------------------------------------------------- |
---|
249 | # Evolve system through time |
---|
250 | #------------------------------------------------------------------------------- |
---|
251 | import time |
---|
252 | t0 = time.time() |
---|
253 | |
---|
254 | for t in domain.evolve(yieldstep = 240, finaltime = 10800): |
---|
255 | domain.write_time() |
---|
256 | domain.write_boundary_statistics(tags = 'topright') |
---|
257 | |
---|
258 | for t in domain.evolve(yieldstep = 120, finaltime = 16200 |
---|
259 | ,skip_initial_step = True): |
---|
260 | domain.write_time() |
---|
261 | domain.write_boundary_statistics(tags = 'topright') |
---|
262 | |
---|
263 | for t in domain.evolve(yieldstep = 60, finaltime = 21600 |
---|
264 | ,skip_initial_step = True): |
---|
265 | domain.write_time() |
---|
266 | domain.write_boundary_statistics(tags = 'topright') |
---|
267 | |
---|
268 | for t in domain.evolve(yieldstep = 120, finaltime = 27000 |
---|
269 | ,skip_initial_step = True): |
---|
270 | domain.write_time() |
---|
271 | domain.write_boundary_statistics(tags = 'topright') |
---|
272 | |
---|
273 | for t in domain.evolve(yieldstep = 240, finaltime = 36000 |
---|
274 | ,skip_initial_step = True): |
---|
275 | domain.write_time() |
---|
276 | domain.write_boundary_statistics(tags = 'topright') |
---|
277 | |
---|
278 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
279 | |
---|
280 | print 'finished' |
---|