1 | """Example of shallow water wave equation analytical solution of the |
---|
2 | one-dimensional Thacker and Greenspan wave run-up treated as a two-dimensional solution. |
---|
3 | |
---|
4 | Copyright 2004 |
---|
5 | Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray |
---|
6 | Geoscience Australia |
---|
7 | |
---|
8 | Specific methods pertaining to the 2D shallow water equation |
---|
9 | are imported from shallow_water |
---|
10 | for use with the generic finite volume framework |
---|
11 | |
---|
12 | Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the |
---|
13 | numerical vector named conserved_quantities. |
---|
14 | """ |
---|
15 | |
---|
16 | #------------------- |
---|
17 | # Module imports |
---|
18 | import sys |
---|
19 | from math import sqrt, cos, sin, pi |
---|
20 | from numpy import asarray |
---|
21 | from time import localtime, strftime, gmtime |
---|
22 | |
---|
23 | |
---|
24 | #------------------- |
---|
25 | # Anuga Imports |
---|
26 | #------------------- |
---|
27 | from anuga.shallow_water_balanced.swb_domain import Domain, Transmissive_boundary, Reflective_boundary,\ |
---|
28 | Dirichlet_boundary, Time_boundary |
---|
29 | |
---|
30 | #from anuga.interface import Domain, Transmissive_boundary, Reflective_boundary,\ |
---|
31 | # Dirichlet_boundary |
---|
32 | |
---|
33 | from anuga.utilities.polygon import inside_polygon, is_inside_triangle |
---|
34 | from anuga.abstract_2d_finite_volumes.mesh_factory import strang_mesh |
---|
35 | |
---|
36 | |
---|
37 | #------------------- |
---|
38 | #Convenience functions |
---|
39 | #------------------- |
---|
40 | def imag(a): |
---|
41 | return a.imag |
---|
42 | |
---|
43 | def real(a): |
---|
44 | return a.real |
---|
45 | |
---|
46 | |
---|
47 | #-------------------- |
---|
48 | # Domain |
---|
49 | # Strang_domain will search through the file and test to see if there are |
---|
50 | # two or three entries. Two entries are for points and three for triangles. |
---|
51 | |
---|
52 | #points, elements = strang_mesh('Run-up.pt') |
---|
53 | points, elements = strang_mesh('strang_7389.pt') |
---|
54 | domain = Domain(points, elements) |
---|
55 | |
---|
56 | print 'extent ',domain.get_extent() |
---|
57 | |
---|
58 | #---------------- |
---|
59 | # Order of scheme |
---|
60 | # Good compromise between |
---|
61 | # limiting and CFL |
---|
62 | #--------------- |
---|
63 | domain.set_default_order(2) |
---|
64 | domain.set_timestepping_method(2) |
---|
65 | domain.set_beta(0.7) |
---|
66 | domain.set_CFL(0.6) |
---|
67 | |
---|
68 | |
---|
69 | #------------------- |
---|
70 | #Set a default tagging |
---|
71 | #------------------- |
---|
72 | epsilon = 1.0e-12 |
---|
73 | |
---|
74 | for id, face in domain.boundary: |
---|
75 | domain.boundary[(id,face)] = 'external' |
---|
76 | if domain.get_vertex_coordinate(id,(face+1)%3)[0] < -200.0+1.0e-10: |
---|
77 | domain.boundary[(id,face)] = 'left' |
---|
78 | |
---|
79 | #--------------------- |
---|
80 | # Provide file name for storing output |
---|
81 | #--------------------- |
---|
82 | domain.store = True |
---|
83 | domain.format = 'sww' |
---|
84 | |
---|
85 | time = strftime('%Y%m%d_%H%M%S',localtime()) |
---|
86 | output_file= 'carrier_wave_runup_'+time |
---|
87 | |
---|
88 | domain.set_name(output_file) |
---|
89 | |
---|
90 | print "Number of triangles = ", len(domain) |
---|
91 | |
---|
92 | #----------------------- |
---|
93 | #Define the boundary condition |
---|
94 | #----------------------- |
---|
95 | def stage_setup(x,t): |
---|
96 | vh = 0 |
---|
97 | alpha = 0.1 |
---|
98 | eta = 0.1 |
---|
99 | a = 1.5*sqrt(1.+0.9*eta) |
---|
100 | l_0 = 200. |
---|
101 | ii = complex(0,1) |
---|
102 | g = 9.81 |
---|
103 | v_0 = sqrt(g*l_0*alpha) |
---|
104 | v1 = 0. |
---|
105 | |
---|
106 | sigma_max = 100. |
---|
107 | sigma_min = -100. |
---|
108 | for j in range (1,50): |
---|
109 | sigma0 = (sigma_max+sigma_min)/2. |
---|
110 | lambda_prime = 2./a*(t/sqrt(l_0/alpha/g)+v1) |
---|
111 | sigma_prime = sigma0/a |
---|
112 | const = (1.-ii*lambda_prime)**2+sigma_prime**2 |
---|
113 | |
---|
114 | v1 = 8.*eta/a*imag(1./const**(3./2.) \ |
---|
115 | -3./4.*(1.-ii*lambda_prime)/const**(5./2.)) |
---|
116 | |
---|
117 | x1 = -v1**2/2.-a**2*sigma_prime**2/16.+eta \ |
---|
118 | *real(1.-2.*(5./4.-ii*lambda_prime) \ |
---|
119 | /const**(3./2.)+3./2.*(1.-ii*lambda_prime)**2 \ |
---|
120 | /const**(5./2.)) |
---|
121 | |
---|
122 | neta1 = x1 + a*a*sigma_prime**2/16. |
---|
123 | |
---|
124 | v_star1 = v1*v_0 |
---|
125 | x_star1 = x1*l_0 |
---|
126 | neta_star1 = neta1*alpha*l_0 |
---|
127 | stage = neta_star1 |
---|
128 | z = stage - x_star1*alpha |
---|
129 | uh = z*v_star1 |
---|
130 | |
---|
131 | if x_star1-x > 0: |
---|
132 | sigma_max = sigma0 |
---|
133 | else: |
---|
134 | sigma_min = sigma0 |
---|
135 | |
---|
136 | if abs(abs(sigma0)-100.) < 10: |
---|
137 | |
---|
138 | # solution does not converge because bed is dry |
---|
139 | stage = 0. |
---|
140 | uh = 0. |
---|
141 | z = 0. |
---|
142 | |
---|
143 | return [stage, uh, vh] |
---|
144 | |
---|
145 | def boundary_stage(t): |
---|
146 | x = -200 |
---|
147 | return stage_setup(x,t) |
---|
148 | |
---|
149 | |
---|
150 | #--------------------- |
---|
151 | #Initial condition |
---|
152 | #--------------------- |
---|
153 | print 'Initial condition' |
---|
154 | t_star1 = 0.0 |
---|
155 | slope = -0.1 |
---|
156 | |
---|
157 | #Set bed-elevation and friction(None) |
---|
158 | def x_slope(x,y): |
---|
159 | n = x.shape[0] |
---|
160 | z = 0*x |
---|
161 | for i in range(n): |
---|
162 | z[i] = -slope*x[i] |
---|
163 | return z |
---|
164 | |
---|
165 | domain.set_quantity('elevation', x_slope) |
---|
166 | |
---|
167 | #Set the water depth |
---|
168 | print 'Initial water depth' |
---|
169 | def stage(x,y): |
---|
170 | z = x_slope(x,y) |
---|
171 | n = x.shape[0] |
---|
172 | w = 0*x |
---|
173 | for i in range(n): |
---|
174 | w[i], uh, vh = stage_setup(x[i],t_star1) |
---|
175 | h = w[i] - z[i] |
---|
176 | if h < 0: |
---|
177 | h = 0 |
---|
178 | w[i] = z[i] |
---|
179 | return w |
---|
180 | |
---|
181 | domain.set_quantity('stage', stage) |
---|
182 | |
---|
183 | |
---|
184 | ##################### |
---|
185 | #Set up boundary conditions |
---|
186 | Br = Reflective_boundary(domain) |
---|
187 | Bw = Time_boundary(domain, boundary_stage) |
---|
188 | domain.set_boundary({'left': Bw, 'external': Br}) |
---|
189 | |
---|
190 | import time |
---|
191 | visualize = True |
---|
192 | if visualize: |
---|
193 | from anuga.visualiser import RealtimeVisualiser |
---|
194 | vis = RealtimeVisualiser(domain) |
---|
195 | vis.render_quantity_height("elevation", zScale=3.0, offset = 0.01, dynamic=False) |
---|
196 | vis.render_quantity_height("stage", zScale = 3.0, dynamic=True, opacity = 0.6, wireframe=False) |
---|
197 | #vis.colour_height_quantity('stage', (lambda q:q['stage'], 1.0, 2.0)) |
---|
198 | vis.colour_height_quantity('stage', (0.4, 0.6, 0.4)) |
---|
199 | vis.start() |
---|
200 | time.sleep(2.0) |
---|
201 | |
---|
202 | #domain.visualise = True |
---|
203 | |
---|
204 | ###################### |
---|
205 | #Evolution |
---|
206 | |
---|
207 | t0 = time.time() |
---|
208 | |
---|
209 | |
---|
210 | for t in domain.evolve(yieldstep = 1., finaltime = 100): |
---|
211 | domain.write_time() |
---|
212 | print boundary_stage(domain.time) |
---|
213 | if visualize: vis.update() |
---|
214 | |
---|
215 | if visualize: vis.evolveFinished() |
---|
216 | |
---|
217 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
218 | |
---|