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 |
---|
7 | ANU |
---|
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 | from anuga.pyvolution.shallow_water import Domain, Transmissive_boundary, Reflective_boundary,\ |
---|
25 | Dirichlet_boundary, gravity, linear_friction |
---|
26 | from math import sqrt, cos, sin, pi, exp |
---|
27 | from anuga.pyvolution.mesh_factory import rectangular_cross |
---|
28 | from anuga.pyvolution.quantity import Quantity |
---|
29 | from anuga.utilities.polygon import inside_polygon |
---|
30 | from Numeric import asarray |
---|
31 | from least_squares import Interpolation |
---|
32 | |
---|
33 | #------------------------------- |
---|
34 | # Domain |
---|
35 | n = 100 |
---|
36 | m = 100 |
---|
37 | lenx = 10000.0 |
---|
38 | leny = 10000.0 |
---|
39 | origin = (-5000.0, -5000.0) |
---|
40 | |
---|
41 | points, elements, boundary = rectangular_cross(m, n, lenx, leny, origin) |
---|
42 | domain = Domain(points, elements, boundary) |
---|
43 | |
---|
44 | #---------------- |
---|
45 | # Order of scheme |
---|
46 | domain.default_order = 1 |
---|
47 | domain.smooth = True |
---|
48 | |
---|
49 | domain.quantities['linear_friction'] = Quantity(domain) |
---|
50 | |
---|
51 | #--------------------------------------------------------------------------- |
---|
52 | # Reconstruct forcing terms with linear friction instead of manning friction |
---|
53 | domain.forcing_terms = [] |
---|
54 | domain.forcing_terms.append(gravity) |
---|
55 | domain.forcing_terms.append(linear_friction) |
---|
56 | |
---|
57 | print domain.forcing_terms |
---|
58 | |
---|
59 | #------------------------------------- |
---|
60 | # Provide file name for storing output |
---|
61 | domain.store = True |
---|
62 | domain.format = 'sww' |
---|
63 | domain.filename = 'sampson_first_order_cross' |
---|
64 | |
---|
65 | print 'Number of triangles = ', len(domain) |
---|
66 | |
---|
67 | #----------------------------------------- |
---|
68 | #Reduction operation for get_vertex_values |
---|
69 | from anuga.utilities.numerical_tools import mean |
---|
70 | domain.reduction = mean |
---|
71 | #domain.reduction = min #Looks better near steep slopes |
---|
72 | |
---|
73 | |
---|
74 | #----------------- |
---|
75 | #Initial condition |
---|
76 | # |
---|
77 | print 'Initial condition' |
---|
78 | t = 0.0 |
---|
79 | h0 = 10. |
---|
80 | a = 3000. |
---|
81 | g = 9.81 |
---|
82 | tau =0.001 |
---|
83 | B = 5 |
---|
84 | A = 0 |
---|
85 | p = sqrt(8*g*h0/a/a) |
---|
86 | s = sqrt(p*p-tau*tau)/2 |
---|
87 | t = 0. |
---|
88 | |
---|
89 | #------------------------------ |
---|
90 | #Set bed-elevation and friction |
---|
91 | def x_slope(x,y): |
---|
92 | n = x.shape[0] |
---|
93 | z = 0*x |
---|
94 | for i in range(n): |
---|
95 | z[i] = h0 - h0*(1.0 -x[i]*x[i]/a/a - y[i]*y[i]/a/a) |
---|
96 | return z |
---|
97 | |
---|
98 | domain.set_quantity('elevation', x_slope) |
---|
99 | domain.set_quantity('linear_friction', tau) |
---|
100 | |
---|
101 | #------------------- |
---|
102 | #Set the water stage |
---|
103 | def stage(x,y): |
---|
104 | z = x_slope(x,y) |
---|
105 | n = x.shape[0] |
---|
106 | h = 0*x |
---|
107 | for i in range(n): |
---|
108 | h[i] = h0-B*B*exp(-tau*t)/2/g-1/g*(exp(-tau*t/2)*(B*s*cos(s*t) \ |
---|
109 | +tau*B/2*sin(s*t)))*x[i] \ |
---|
110 | -1/g*(exp(-tau*t/2)*(B*s*sin(s*t)-tau*B/2*cos(s*t)))*y[i] |
---|
111 | if h[i] < z[i]: |
---|
112 | h[i] = z[i] |
---|
113 | return h |
---|
114 | |
---|
115 | domain.set_quantity('stage', stage) |
---|
116 | |
---|
117 | |
---|
118 | #--------- |
---|
119 | # Boundary |
---|
120 | print 'Boundary conditions' |
---|
121 | R = Reflective_boundary(domain) |
---|
122 | T = Transmissive_boundary(domain) |
---|
123 | D = Dirichlet_boundary([0.0, 0.0, 0.0]) |
---|
124 | domain.set_boundary({'left': D, 'right': D, 'top': D, 'bottom': D}) |
---|
125 | |
---|
126 | #--------------------------------------------- |
---|
127 | # Find triangle that contains the point points |
---|
128 | # and print to file |
---|
129 | points = [0.,0.] |
---|
130 | for n in range(len(domain.triangles)): |
---|
131 | t = domain.triangles[n] |
---|
132 | tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]] |
---|
133 | |
---|
134 | if inside_polygon(points,tri): |
---|
135 | print 'Point is within triangle with vertices '+'%s'%tri |
---|
136 | n_point = n |
---|
137 | |
---|
138 | print 'n_point = ',n_point |
---|
139 | t = domain.triangles[n_point] |
---|
140 | tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]] |
---|
141 | |
---|
142 | filename=domain.filename |
---|
143 | file = open(filename,'w') |
---|
144 | |
---|
145 | #---------- |
---|
146 | # Evolution |
---|
147 | import time |
---|
148 | t0 = time.time() |
---|
149 | |
---|
150 | Stage = domain.quantities['stage'] |
---|
151 | Xmomentum = domain.quantities['xmomentum'] |
---|
152 | Ymomentum = domain.quantities['ymomentum'] |
---|
153 | |
---|
154 | for t in domain.evolve(yieldstep = 20.0, finaltime = 10000.0 ): |
---|
155 | domain.write_time() |
---|
156 | |
---|
157 | tri_array = asarray(tri) |
---|
158 | t_array = asarray([[0,1,2]]) |
---|
159 | interp = Interpolation(tri_array,t_array,[points]) |
---|
160 | |
---|
161 | |
---|
162 | stage = Stage.get_values(location='centroids',indices=[n_point])[0] |
---|
163 | xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])[0] |
---|
164 | ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])[0] |
---|
165 | file.write( '%10.6f %10.6f %10.6f %10.6f\n'%(t,stage,xmomentum,ymomentum) ) |
---|
166 | |
---|
167 | file.close() |
---|
168 | print 'That took %.2f seconds' %(time.time()-t0) |
---|