1 | """Example of shallow water wave equation analytical solution |
---|
2 | consists of a flat water surface profile in a parabolic basin |
---|
3 | with linear friction. The analytical solution was derived by Sampson in 2002. |
---|
4 | |
---|
5 | Copyright 2004 |
---|
6 | Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray |
---|
7 | Geoscience Australia |
---|
8 | |
---|
9 | Specific methods pertaining to the 2D shallow water equation |
---|
10 | are imported from shallow_water |
---|
11 | for use with the generic finite volume framework |
---|
12 | |
---|
13 | Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the |
---|
14 | numerical vector named conserved_quantities. |
---|
15 | """ |
---|
16 | |
---|
17 | ###################### |
---|
18 | # Module imports |
---|
19 | # |
---|
20 | |
---|
21 | import sys |
---|
22 | from os import sep |
---|
23 | sys.path.append('..'+sep+'pyvolution') |
---|
24 | |
---|
25 | from shallow_water import Domain, Dirichlet_boundary, gravity, linear_friction |
---|
26 | from math import sqrt, cos, sin, pi, exp |
---|
27 | from mesh_factory import strang_mesh |
---|
28 | from quantity import Quantity |
---|
29 | |
---|
30 | |
---|
31 | ###################### |
---|
32 | # Domain |
---|
33 | # Strang_domain will search through the file and test to see if there are |
---|
34 | # two or three entries. Two entries are for points and three for triangles. |
---|
35 | |
---|
36 | |
---|
37 | points, elements = strang_mesh('yoon_circle.pt') |
---|
38 | domain = Domain(points, elements) |
---|
39 | |
---|
40 | domain.default_order = 2 |
---|
41 | domain.smooth = True |
---|
42 | |
---|
43 | domain.quantities['linear_friction'] = Quantity(domain) |
---|
44 | |
---|
45 | #Reconstruct forcing terms with linear friction instead og manning friction |
---|
46 | domain.forcing_terms = [] |
---|
47 | domain.forcing_terms.append(gravity) |
---|
48 | domain.forcing_terms.append(linear_friction) |
---|
49 | |
---|
50 | print domain.forcing_terms |
---|
51 | |
---|
52 | # Provide file name for storing output |
---|
53 | domain.store = True |
---|
54 | domain.format = 'sww' |
---|
55 | domain.set_name('sampson_strang_second_order') |
---|
56 | |
---|
57 | print 'Number of triangles = ', len(domain) |
---|
58 | |
---|
59 | |
---|
60 | #Reduction operation for get_vertex_values |
---|
61 | from anuga.pyvolution.util import mean |
---|
62 | domain.reduction = mean |
---|
63 | #domain.reduction = min #Looks better near steep slopes |
---|
64 | |
---|
65 | |
---|
66 | ###################### |
---|
67 | #Initial condition |
---|
68 | # |
---|
69 | print 'Initial condition' |
---|
70 | t = 0.0 |
---|
71 | h0 = 10. |
---|
72 | a = 3000. |
---|
73 | g = 9.81 |
---|
74 | tau =0.001 |
---|
75 | B = 5 |
---|
76 | A = 0 |
---|
77 | p = sqrt(8*g*h0/a/a) |
---|
78 | s = sqrt(p*p-tau*tau)/2 |
---|
79 | t = 0. |
---|
80 | |
---|
81 | |
---|
82 | #Set bed-elevation and friction |
---|
83 | def x_slope(x,y): |
---|
84 | n = x.shape[0] |
---|
85 | z = 0*x |
---|
86 | for i in range(n): |
---|
87 | z[i] = h0 - h0*(1.0 -x[i]*x[i]/a/a - y[i]*y[i]/a/a) |
---|
88 | return z |
---|
89 | |
---|
90 | domain.set_quantity('elevation', x_slope) |
---|
91 | domain.set_quantity('linear_friction', tau) |
---|
92 | |
---|
93 | |
---|
94 | #Set the water stage |
---|
95 | def stage(x,y): |
---|
96 | z = x_slope(x,y) |
---|
97 | n = x.shape[0] |
---|
98 | h = 0*x |
---|
99 | for i in range(n): |
---|
100 | h[i] = h0-B*B*exp(-tau*t)/2/g-1/g*(exp(-tau*t/2)*(B*s*cos(s*t) \ |
---|
101 | +tau*B/2*sin(s*t)))*x[i] \ |
---|
102 | -1/g*(exp(-tau*t/2)*(B*s*sin(s*t)-tau*B/2*cos(s*t)))*y[i] |
---|
103 | if h[i] < z[i]: |
---|
104 | h[i] = z[i] |
---|
105 | return h |
---|
106 | |
---|
107 | domain.set_quantity('stage', stage) |
---|
108 | |
---|
109 | |
---|
110 | ############ |
---|
111 | #Boundary |
---|
112 | domain.set_boundary({'exterior': Dirichlet_boundary([0.0, 0.0, 0.0])}) |
---|
113 | |
---|
114 | |
---|
115 | ###################### |
---|
116 | #Evolution |
---|
117 | import time |
---|
118 | t0 = time.time() |
---|
119 | for t in domain.evolve(yieldstep = 10.0, finaltime = 5000): |
---|
120 | domain.write_time() |
---|
121 | |
---|
122 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
123 | |
---|
124 | |
---|