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