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