1 | """Script for running a tsunami inundation scenario for Sydney, NSW, 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.outputdir |
---|
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 | """ |
---|
11 | |
---|
12 | |
---|
13 | #------------------------------------------------------------------------------ |
---|
14 | # Import necessary modules |
---|
15 | #------------------------------------------------------------------------------ |
---|
16 | |
---|
17 | # Standard modules |
---|
18 | import os |
---|
19 | import time |
---|
20 | |
---|
21 | # Related major packages |
---|
22 | from pyvolution.shallow_water import Domain, Dirichlet_boundary |
---|
23 | from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts |
---|
24 | from pyvolution.combine_pts import combine_rectangular_points_files |
---|
25 | from pmesh.mesh_interface import create_mesh_from_regions |
---|
26 | |
---|
27 | # Application specific imports |
---|
28 | import project # Define file names and polygons |
---|
29 | from pyvolution.smf import slump_tsunami # Function for submarine mudslide |
---|
30 | |
---|
31 | |
---|
32 | |
---|
33 | #------------------------------------------------------------------------------ |
---|
34 | # Prepare bathymetric and topographic data |
---|
35 | #------------------------------------------------------------------------------ |
---|
36 | |
---|
37 | # filenames |
---|
38 | coarsedemname = project.coarsedemname + '.pts' |
---|
39 | finedemname = project.finedemname + '.pts' |
---|
40 | combineddemname = project.combineddemname + '.pts' |
---|
41 | meshname = project.meshname+'.msh' |
---|
42 | |
---|
43 | |
---|
44 | # combining the coarse and fine data |
---|
45 | combine_rectangular_points_files(finedemname, |
---|
46 | coarsedemname, |
---|
47 | combineddemname) |
---|
48 | |
---|
49 | |
---|
50 | |
---|
51 | #------------------------------------------------------------------------------ |
---|
52 | # Setup computational domain |
---|
53 | #------------------------------------------------------------------------------ |
---|
54 | |
---|
55 | # Interior regions and mesh resolutions |
---|
56 | interior_res = 5000 |
---|
57 | interior_regions = [[project.northern_polygon, interior_res], |
---|
58 | [project.harbour_polygon, interior_res], |
---|
59 | [project.southern_polygon, interior_res], |
---|
60 | [project.manly_polygon, interior_res], |
---|
61 | [project.top_polygon, interior_res]] |
---|
62 | |
---|
63 | |
---|
64 | create_mesh_from_regions(project.demopoly, |
---|
65 | boundary_tags= {'oceannorth': [0], |
---|
66 | 'onshorenorth1': [1], |
---|
67 | 'onshorenorthwest1': [2], |
---|
68 | 'onshorenorthwest2': [3], |
---|
69 | 'onshorenorth2': [4], |
---|
70 | 'onshorewest': [5], |
---|
71 | 'onshoresouth': [6], |
---|
72 | 'oceansouth': [7], |
---|
73 | 'oceaneast': [8]}, |
---|
74 | maximum_triangle_area=100000, |
---|
75 | filename=meshname, |
---|
76 | interior_regions=interior_regions) |
---|
77 | |
---|
78 | |
---|
79 | #Create shallow water domain |
---|
80 | |
---|
81 | domain = Domain(meshname, |
---|
82 | use_cache = True, |
---|
83 | verbose = True) |
---|
84 | |
---|
85 | |
---|
86 | print 'Number of triangles = ', len(domain) |
---|
87 | print 'The extent is ', domain.get_extent() |
---|
88 | |
---|
89 | domain.set_name(project.basename) |
---|
90 | domain.set_datadir(project.outputdir) |
---|
91 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
92 | |
---|
93 | |
---|
94 | #------------------------------------------------------------------------------ |
---|
95 | # Set up scenario (tsunami_source is a callable object used with set_quantity) |
---|
96 | #------------------------------------------------------------------------------ |
---|
97 | |
---|
98 | tsunami_source = slump_tsunami(length=30000.0, |
---|
99 | depth=400.0, |
---|
100 | slope=6.0, |
---|
101 | thickness=176.0, |
---|
102 | radius=3330, |
---|
103 | dphi=0.23, |
---|
104 | x0=project.slump_origin[0], |
---|
105 | y0=project.slump_origin[1], |
---|
106 | alpha=0.0, |
---|
107 | domain=domain) |
---|
108 | |
---|
109 | |
---|
110 | #------------------------------------------------------------------------------ |
---|
111 | # Setup initial conditions |
---|
112 | #------------------------------------------------------------------------------ |
---|
113 | |
---|
114 | domain.set_quantity('stage', tsunami_source) |
---|
115 | domain.set_quantity('friction', 0.0) |
---|
116 | domain.set_quantity('elevation', |
---|
117 | filename = combineddemname, |
---|
118 | use_cache = True, |
---|
119 | verbose = True) |
---|
120 | |
---|
121 | #------------------------------------------------------------------------------ |
---|
122 | # Setup boundary conditions (all Dirichlet) |
---|
123 | #------------------------------------------------------------------------------ |
---|
124 | |
---|
125 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
126 | |
---|
127 | Bd = Dirichlet_boundary([0.0,0.0,0.0]) |
---|
128 | domain.set_boundary( {'oceannorth': Bd, |
---|
129 | 'onshorenorth1': Bd, |
---|
130 | 'onshorenorthwest1': Bd, |
---|
131 | 'onshorenorthwest2': Bd, |
---|
132 | 'onshorenorth2': Bd, |
---|
133 | 'onshorewest': Bd, |
---|
134 | 'onshoresouth': Bd, |
---|
135 | 'oceansouth': Bd, |
---|
136 | 'oceaneast': Bd} ) |
---|
137 | |
---|
138 | #------------------------------------------------------------------------------ |
---|
139 | # Evolve system through time |
---|
140 | #------------------------------------------------------------------------------ |
---|
141 | |
---|
142 | import time |
---|
143 | t0 = time.time() |
---|
144 | |
---|
145 | for t in domain.evolve(yieldstep = 60, finaltime = 7200): |
---|
146 | print domain.timestepping_statistics() |
---|
147 | print domain.boundary_statistics(tags = 'oceaneast') |
---|
148 | |
---|
149 | print 'That took %.2f seconds' %(time.time()-t0) |
---|