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