1 | """Script for running tsunami inundation scenario for Dampier, 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.output_run_time_dir |
---|
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 Jane Sexton, Nick Bartzis, GA - 2006 |
---|
11 | """ |
---|
12 | |
---|
13 | #------------------------------------------------------------------------------ |
---|
14 | # Import necessary modules |
---|
15 | #------------------------------------------------------------------------------ |
---|
16 | |
---|
17 | # Standard modules |
---|
18 | from os import sep |
---|
19 | from os.path import dirname, basename |
---|
20 | from os import mkdir, access, F_OK |
---|
21 | from shutil import copy |
---|
22 | import time |
---|
23 | import sys |
---|
24 | |
---|
25 | |
---|
26 | # Related major packages |
---|
27 | from anuga.shallow_water import Domain |
---|
28 | from anuga.shallow_water import Dirichlet_boundary |
---|
29 | from anuga.shallow_water import File_boundary |
---|
30 | from anuga.shallow_water import Reflective_boundary |
---|
31 | from Numeric import allclose |
---|
32 | |
---|
33 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
34 | |
---|
35 | from anuga.geospatial_data.geospatial_data import * |
---|
36 | from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files |
---|
37 | from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier |
---|
38 | |
---|
39 | # Application specific imports |
---|
40 | import project # Definition of file names and polygons |
---|
41 | |
---|
42 | #------------------------------------------------------------------------------ |
---|
43 | # Copy scripts to time stamped output directory and capture screen |
---|
44 | # output to file |
---|
45 | #------------------------------------------------------------------------------ |
---|
46 | |
---|
47 | |
---|
48 | # filenames |
---|
49 | |
---|
50 | #build_time = '20061029_231935_build_tide_24' |
---|
51 | #build_time = '20061030_165746_build_tide_24' |
---|
52 | #build_time = '20061102_215532_build_plus_old_data' |
---|
53 | #build_time = '20061103_055258_build' |
---|
54 | build_time = '20061107_063840_build' |
---|
55 | #build_time = '20061025_153643_build_basic' |
---|
56 | bound_time = '20061102_221245_build' |
---|
57 | |
---|
58 | boundaries_name = project.boundaries_name |
---|
59 | #meshes_time_dir_name = project.meshes_time_dir_name+'.msh' |
---|
60 | meshes_dir_name = project.meshes_dir_name+'.msh' |
---|
61 | #source_dir = project.boundarydir |
---|
62 | #boundaries_time_dir_name = project.boundaries_dir + build_time + sep + boundaries_name |
---|
63 | boundaries_time_dir_name = project.boundaries_time_dir_name |
---|
64 | boundaries_dir_name = project.boundaries_dir_name |
---|
65 | tide = project.tide |
---|
66 | |
---|
67 | # creates copy of code in output dir |
---|
68 | if myid == 0: |
---|
69 | copy_code_files(project.output_run_time_dir,__file__, |
---|
70 | dirname(project.__file__)+sep+ project.__name__+'.py' ) |
---|
71 | barrier() |
---|
72 | start_screen_catcher(project.output_run_time_dir, myid, numprocs) |
---|
73 | |
---|
74 | print 'USER: ', project.user |
---|
75 | #sys.exit() |
---|
76 | #-------------------------------------------------------------------------- |
---|
77 | # Create the triangular mesh based on overall clipping polygon with a |
---|
78 | # tagged |
---|
79 | # boundary and interior regions defined in project.py along with |
---|
80 | # resolutions (maximal area of per triangle) for each polygon |
---|
81 | #-------------------------------------------------------------------------- |
---|
82 | |
---|
83 | if myid == 0: |
---|
84 | |
---|
85 | print 'start create mesh from regions' |
---|
86 | interior_regions = [#[project.karratha_polygon, 25000], |
---|
87 | [project.poly_coast, 10000], |
---|
88 | [project.poly_pipeline, 2000], |
---|
89 | [project.poly_facility, 500]] |
---|
90 | # [project.poly_interior, 1000]] |
---|
91 | |
---|
92 | |
---|
93 | create_mesh_from_regions(project.bounding_polygon, |
---|
94 | boundary_tags={'back': [7, 8], 'side': [0, 6], |
---|
95 | 'ocean': [1, 2, 3, 4, 5]}, |
---|
96 | maximum_triangle_area=100000, |
---|
97 | interior_regions=interior_regions, |
---|
98 | # filename=meshes_time_dir_name, |
---|
99 | filename=meshes_dir_name, |
---|
100 | use_cache=True, |
---|
101 | verbose=True) |
---|
102 | |
---|
103 | # to sync all processors are ready |
---|
104 | barrier() |
---|
105 | |
---|
106 | #------------------------------------------------------------------------- |
---|
107 | # Setup computational domain |
---|
108 | #------------------------------------------------------------------------- |
---|
109 | print 'Setup computational domain' |
---|
110 | #domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True) |
---|
111 | domain = Domain(meshes_dir_name, use_cache=True, verbose=True) |
---|
112 | print domain.statistics() |
---|
113 | |
---|
114 | print 'starting to create boundary conditions' |
---|
115 | boundaries_in_dir_name = project.boundaries_in_dir_name |
---|
116 | |
---|
117 | from anuga.shallow_water.data_manager import urs2sww, ferret2sww |
---|
118 | |
---|
119 | print 'maxlat=project.south_boundary, minlat=project.north_boundary', project.south_boundary,project.north_boundary |
---|
120 | print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary |
---|
121 | print ' maxlon=project.east',project.east |
---|
122 | |
---|
123 | print 'origin: domain.geo_reference.get_origin()',domain.geo_reference.get_origin() |
---|
124 | |
---|
125 | #import sys; sys.exit() |
---|
126 | |
---|
127 | #if access(project.boundaries_time_dir,F_OK) == 0: |
---|
128 | # mkdir (project.boundaries_time_dir) |
---|
129 | # put above distribute |
---|
130 | print 'boundary file is: ',boundaries_dir_name |
---|
131 | from caching import cache |
---|
132 | if myid == 0: |
---|
133 | cache(ferret2sww, |
---|
134 | (boundaries_in_dir_name, |
---|
135 | # boundaries_time_dir_name), |
---|
136 | boundaries_dir_name), |
---|
137 | {'verbose': True, |
---|
138 | 'minlat': project.south_boundary, |
---|
139 | 'maxlat': project.north_boundary, |
---|
140 | 'minlon': project.west_boundary, |
---|
141 | 'maxlon': project.east_boundary, |
---|
142 | # 'minlat': project.south, |
---|
143 | # 'maxlat': project.north, |
---|
144 | # 'minlon': project.west, |
---|
145 | # 'maxlon': project.east, |
---|
146 | 'mint': 0, 'maxt': 35100, |
---|
147 | 'origin': domain.geo_reference.get_origin(), |
---|
148 | 'mean_stage': project.tide, |
---|
149 | # 'zscale': 1, #Enhance tsunami |
---|
150 | 'fail_on_NaN': False}, |
---|
151 | verbose = True, |
---|
152 | ) |
---|
153 | barrier() |
---|
154 | |
---|
155 | |
---|
156 | #------------------------------------------------------------------------- |
---|
157 | # Setup initial conditions |
---|
158 | #------------------------------------------------------------------------- |
---|
159 | if myid == 0: |
---|
160 | |
---|
161 | print 'Setup initial conditions' |
---|
162 | |
---|
163 | domain.set_quantity('stage', tide) |
---|
164 | domain.set_quantity('friction', 0.0) |
---|
165 | #combined_time_dir_name = project.topographies_dir+build_time+project.combined_name |
---|
166 | print 'Start Set quantity' |
---|
167 | |
---|
168 | |
---|
169 | domain.set_quantity('elevation', |
---|
170 | filename = project.topographies_dir + build_time + sep + project.combined_final_name + '.pts', |
---|
171 | use_cache = True, |
---|
172 | verbose = True, |
---|
173 | alpha = 0.1) |
---|
174 | #domain.set_quantity('elevation', -50) |
---|
175 | print 'Finished Set quantity' |
---|
176 | barrier() |
---|
177 | |
---|
178 | |
---|
179 | #------------------------------------------------------ |
---|
180 | # Distribute domain to implement parallelism !!! |
---|
181 | #------------------------------------------------------ |
---|
182 | |
---|
183 | if numprocs > 1: |
---|
184 | domain=distribute(domain) |
---|
185 | |
---|
186 | #------------------------------------------------------ |
---|
187 | # Set domain parameters |
---|
188 | #------------------------------------------------------ |
---|
189 | |
---|
190 | domain.set_name(project.scenario_name) |
---|
191 | domain.set_datadir(project.output_run_time_dir) |
---|
192 | domain.set_default_order(2) # Apply second order scheme |
---|
193 | domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm |
---|
194 | domain.set_store_vertices_uniquely(False) |
---|
195 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
196 | domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK) |
---|
197 | |
---|
198 | #------------------------------------------------------------------------- |
---|
199 | # Setup boundary conditions |
---|
200 | #------------------------------------------------------------------------- |
---|
201 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
202 | |
---|
203 | print 'Reading Boundary file' |
---|
204 | #boundariesname = project.boundaries_dir + '20061101_003322_build'+sep+boundaries_name |
---|
205 | #print'boundariesname',boundariesname |
---|
206 | #Bf = File_boundary(boundaries_time_dir_name + '.sww', |
---|
207 | |
---|
208 | #Bf = File_boundary(boundariesname + '.sww', |
---|
209 | |
---|
210 | print 'domain id', id(domain) |
---|
211 | Bf = File_boundary(boundaries_dir_name + '.sww', |
---|
212 | domain, time_thinning=5, use_cache=True, verbose=True) |
---|
213 | |
---|
214 | print 'finished reading boundary file' |
---|
215 | |
---|
216 | Br = Reflective_boundary(domain) |
---|
217 | Bd = Dirichlet_boundary([tide,0,0]) |
---|
218 | |
---|
219 | print'set_boundary' |
---|
220 | domain.set_boundary({'back': Br, |
---|
221 | 'side': Bd, |
---|
222 | 'ocean': Bf}) |
---|
223 | print'finish set boundary' |
---|
224 | |
---|
225 | #---------------------------------------------------------------------------- |
---|
226 | # Evolve system through time |
---|
227 | #---------------------------------------------------------------------------- |
---|
228 | |
---|
229 | t0 = time.time() |
---|
230 | |
---|
231 | #for t in domain.evolve(yieldstep = 60, finaltime = 34000): |
---|
232 | # domain.write_time() |
---|
233 | # domain.write_boundary_statistics(tags = 'ocean') |
---|
234 | |
---|
235 | for t in domain.evolve(yieldstep = 120, finaltime = 9000): |
---|
236 | domain.write_time() |
---|
237 | domain.write_boundary_statistics(tags = 'ocean') |
---|
238 | # print 'time: ',time.time() |
---|
239 | |
---|
240 | if allclose(t, 6000): |
---|
241 | domain.set_quantity('xmomentum', 0) |
---|
242 | domain.set_quantity('ymomentum', 0) |
---|
243 | #import sys; sys.exit() |
---|
244 | |
---|
245 | for t in domain.evolve(yieldstep = 60, finaltime = 28800 |
---|
246 | ,skip_initial_step = True): |
---|
247 | domain.write_time() |
---|
248 | domain.write_boundary_statistics(tags = 'ocean') |
---|
249 | |
---|
250 | for t in domain.evolve(yieldstep = 120, finaltime = 34800 |
---|
251 | ,skip_initial_step = True): |
---|
252 | domain.write_time() |
---|
253 | domain.write_boundary_statistics(tags = 'ocean') |
---|
254 | |
---|
255 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
256 | |
---|