1 | """Script for running a tsunami inundation scenario for Broome, 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 | #------------------------------------------------------------------------------ |
---|
13 | # Import necessary modules |
---|
14 | #------------------------------------------------------------------------------ |
---|
15 | |
---|
16 | # Standard modules |
---|
17 | from os import sep |
---|
18 | from os.path import dirname, basename |
---|
19 | import time |
---|
20 | |
---|
21 | # Related major packages |
---|
22 | from anuga.shallow_water import Domain |
---|
23 | from anuga.shallow_water import Dirichlet_boundary |
---|
24 | from anuga.shallow_water import File_boundary |
---|
25 | from anuga.shallow_water import Reflective_boundary |
---|
26 | |
---|
27 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
28 | |
---|
29 | from shutil import copy |
---|
30 | from os import mkdir, access, F_OK |
---|
31 | from anuga.geospatial_data.geospatial_data import * |
---|
32 | import sys |
---|
33 | from anuga.abstract_2d_finite_volumes.util import Screen_Catcher |
---|
34 | |
---|
35 | # Application specific imports |
---|
36 | import project # Definition of file names and polygons |
---|
37 | |
---|
38 | from anuga_parallel.parallel_api import distribute, myid |
---|
39 | |
---|
40 | #------------------------------------------------------------------------------ |
---|
41 | # Copy scripts to time stamped output directory and capture screen |
---|
42 | # output to file |
---|
43 | #------------------------------------------------------------------------------ |
---|
44 | |
---|
45 | # filenames |
---|
46 | meshname = project.meshname+'.msh' |
---|
47 | source_dir = project.boundarydir |
---|
48 | |
---|
49 | |
---|
50 | |
---|
51 | if myid == 0: |
---|
52 | # creates copy of code in output dir if dir doesn't exist |
---|
53 | if access(project.outputtimedir,F_OK) == 0: |
---|
54 | mkdir (project.outputtimedir) |
---|
55 | copy (dirname(project.__file__) +sep+ project.__name__+'.py', |
---|
56 | project.outputtimedir + project.__name__+'.py') |
---|
57 | copy (__file__, project.outputtimedir + basename(__file__)) |
---|
58 | print 'project.outputtimedir',project.outputtimedir |
---|
59 | |
---|
60 | # normal screen output is stored in |
---|
61 | screen_output_name = project.outputtimedir + 'screen_output.txt' |
---|
62 | screen_error_name = project.outputtimedir + 'screen_error.txt' |
---|
63 | |
---|
64 | # used to catch screen output to file |
---|
65 | #sys.stdout = Screen_Catcher(screen_output_name) |
---|
66 | #sys.stderr = Screen_Catcher(screen_error_name) |
---|
67 | print 'USER: ', project.user |
---|
68 | |
---|
69 | |
---|
70 | #-------------------------------------------------------------------------- |
---|
71 | # Create the triangular mesh based on overall clipping polygon with a tagged |
---|
72 | # boundary and interior regions defined in project.py along with |
---|
73 | # resolutions (maximal area of per triangle) for each polygon |
---|
74 | #-------------------------------------------------------------------------- |
---|
75 | |
---|
76 | |
---|
77 | resolution = 4000 |
---|
78 | interior_regions = [[project.neil1_polygon, resolution], |
---|
79 | [project.neil2_polygon, 64000]] |
---|
80 | |
---|
81 | print 'number of interior regions', len(interior_regions) |
---|
82 | |
---|
83 | |
---|
84 | print 'start create mesh from regions' |
---|
85 | from caching import cache |
---|
86 | |
---|
87 | meshname = project.meshname + '_%d.msh' %myid |
---|
88 | _ = cache(create_mesh_from_regions, |
---|
89 | project.bounding_polygon, |
---|
90 | {'boundary_tags': {'back': [7, 8], 'side': [0, 6], |
---|
91 | 'ocean': [1, 2, 3, 4, 5]}, |
---|
92 | 'maximum_triangle_area': 100000, |
---|
93 | 'filename': meshname}, |
---|
94 | #'interior_regions': interior_regions}, |
---|
95 | verbose = True, |
---|
96 | evaluate = False) |
---|
97 | |
---|
98 | #------------------------------------------------------------------------- |
---|
99 | # Setup computational domain |
---|
100 | #------------------------------------------------------------------------- |
---|
101 | |
---|
102 | domain = Domain(meshname, use_cache = True, verbose = True) |
---|
103 | print domain.statistics() |
---|
104 | |
---|
105 | |
---|
106 | domain.set_name(project.basename) |
---|
107 | domain.set_datadir(project.outputtimedir) |
---|
108 | |
---|
109 | |
---|
110 | #------------------------------------------------------------------------- |
---|
111 | # Setup initial conditions |
---|
112 | #------------------------------------------------------------------------- |
---|
113 | |
---|
114 | tide = 0. |
---|
115 | |
---|
116 | domain.set_quantity('stage', tide) |
---|
117 | domain.set_quantity('friction', 0.0) |
---|
118 | domain.set_quantity('elevation', |
---|
119 | filename = project.datadir + project.basename + '.pts', |
---|
120 | use_cache = False, |
---|
121 | verbose = True, |
---|
122 | alpha = 0.1 |
---|
123 | ) |
---|
124 | |
---|
125 | #------------------------------------------------------------------------- |
---|
126 | # Setup boundary conditions |
---|
127 | #------------------------------------------------------------------------- |
---|
128 | |
---|
129 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
130 | |
---|
131 | Bf = File_boundary(source_dir + project.boundary_basename + '.sww', |
---|
132 | domain, verbose = True) |
---|
133 | Br = Reflective_boundary(domain) |
---|
134 | Bd = Dirichlet_boundary([tide,0,0]) |
---|
135 | domain.set_boundary({'back': Br, |
---|
136 | 'side': Bd, |
---|
137 | #'ocean': None}) # Bind this one later |
---|
138 | 'ocean': Bf}) # Bind this one later |
---|
139 | else: |
---|
140 | domain = None |
---|
141 | |
---|
142 | |
---|
143 | #-------------------------------------------------------------------------- |
---|
144 | # Create the parallel domain |
---|
145 | #-------------------------------------------------------------------------- |
---|
146 | domain = distribute(domain, verbose=True) |
---|
147 | |
---|
148 | |
---|
149 | # Set those boundaries that can't be communicated automatically |
---|
150 | Bf = File_boundary(source_dir + project.boundary_basename + '.sww', |
---|
151 | domain, verbose = True) |
---|
152 | boundary_map = domain.boundary_map |
---|
153 | boundary_map['ocean'] = Bf |
---|
154 | |
---|
155 | print boundary_map |
---|
156 | domain.set_boundary(boundary_map) |
---|
157 | |
---|
158 | |
---|
159 | |
---|
160 | #---------------------------------------------------------------------------- |
---|
161 | # Evolve system through time |
---|
162 | #---------------------------------------------------------------------------- |
---|
163 | import time |
---|
164 | t0 = time.time() |
---|
165 | |
---|
166 | for t in domain.evolve(yieldstep = 60, finaltime = 28600): |
---|
167 | domain.write_time() |
---|
168 | domain.write_boundary_statistics(tags = 'ocean') |
---|
169 | |
---|
170 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
171 | |
---|