1 | """Example of shallow water wave equation analytical solution |
---|
2 | consists of a parabolic profile in a parabolic basin. Analytical |
---|
3 | solutiuon to this problem was derived by Carrier and Greenspan |
---|
4 | and used by Yoon and Chou. |
---|
5 | |
---|
6 | Copyright 2005 |
---|
7 | Christopher Zoppou, Stephen Roberts, ANU, Geoscience Australia |
---|
8 | |
---|
9 | """ |
---|
10 | |
---|
11 | #--------------- |
---|
12 | # Module imports |
---|
13 | #import sys |
---|
14 | #from os import sep |
---|
15 | #sys.path.append('..'+sep+'pyvolution') |
---|
16 | from pyvolution.shallow_water import Domain, Transmissive_boundary, Reflective_boundary,\ |
---|
17 | Dirichlet_boundary |
---|
18 | from math import sqrt, cos, sin, pi |
---|
19 | from pyvolution.mesh_factory import rectangular_cross |
---|
20 | from utilities.polygon import inside_polygon |
---|
21 | from Numeric import asarray |
---|
22 | from pyvolution.least_squares import Interpolation |
---|
23 | |
---|
24 | #------------------------------- |
---|
25 | # Domain |
---|
26 | n = 100 |
---|
27 | m = 100 |
---|
28 | delta_x = 80.0 |
---|
29 | delta_y = 80.0 |
---|
30 | lenx = delta_x*n |
---|
31 | leny = delta_y*m |
---|
32 | origin = (-4000.0, -4000.0) |
---|
33 | |
---|
34 | points, elements, boundary = rectangular_cross(m, n, lenx, leny, origin) |
---|
35 | domain = Domain(points, elements, boundary) |
---|
36 | |
---|
37 | #---------------- |
---|
38 | # Order of scheme |
---|
39 | domain.default_order = 1 |
---|
40 | |
---|
41 | domain.smooth = True |
---|
42 | |
---|
43 | #------------------------------------- |
---|
44 | # Provide file name for storing output |
---|
45 | domain.store = False |
---|
46 | domain.format = 'sww' |
---|
47 | domain.filename = 'yoon_mesh_second_order_cross' |
---|
48 | print 'Number of triangles = ', len(domain) |
---|
49 | |
---|
50 | #---------------------------------------------------------- |
---|
51 | # Decide which quantities are to be stored at each timestep |
---|
52 | domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum'] |
---|
53 | |
---|
54 | #------------------------------------------ |
---|
55 | # Reduction operation for get_vertex_values |
---|
56 | from util import mean |
---|
57 | domain.reduction = mean #domain.reduction = min #Looks better near steep slopes |
---|
58 | |
---|
59 | #------------------ |
---|
60 | # Initial condition |
---|
61 | print 'Initial condition' |
---|
62 | t = 0.0 |
---|
63 | D0 = 1. |
---|
64 | L = 2500. |
---|
65 | R0 = 2000. |
---|
66 | g = 9.81 |
---|
67 | |
---|
68 | A = (L**4 - R0**4)/(L**4 + R0**4) |
---|
69 | omega = 2./L*sqrt(2.*g*D0) |
---|
70 | T = pi/omega |
---|
71 | |
---|
72 | #------------------ |
---|
73 | # Set bed elevation |
---|
74 | def x_slope(x,y): |
---|
75 | n = x.shape[0] |
---|
76 | z = 0*x |
---|
77 | for i in range(n): |
---|
78 | r = sqrt(x[i]*x[i] + y[i]*y[i]) |
---|
79 | z[i] = -D0*(1.-r*r/L/L) |
---|
80 | return z |
---|
81 | domain.set_quantity('elevation', x_slope) |
---|
82 | |
---|
83 | #---------------------------- |
---|
84 | # Set the initial water level |
---|
85 | def level(x,y): |
---|
86 | z = x_slope(x,y) |
---|
87 | n = x.shape[0] |
---|
88 | h = 0*x |
---|
89 | for i in range(n): |
---|
90 | r = sqrt(x[i]*x[i] + y[i]*y[i]) |
---|
91 | h[i] = D0*((sqrt(1-A*A))/(1.-A*cos(omega*t)) |
---|
92 | -1.-r*r/L/L*((1.-A*A)/((1.-A*cos(omega*t))**2)-1.)) |
---|
93 | if h[i] < z[i]: |
---|
94 | h[i] = z[i] |
---|
95 | return h |
---|
96 | domain.set_quantity('stage', level) |
---|
97 | |
---|
98 | #--------- |
---|
99 | # Boundary |
---|
100 | print 'Boundary conditions' |
---|
101 | R = Reflective_boundary(domain) |
---|
102 | T = Transmissive_boundary(domain) |
---|
103 | D = Dirichlet_boundary([0.0, 0.0, 0.0]) |
---|
104 | domain.set_boundary({'left': D, 'right': D, 'top': D, 'bottom': D}) |
---|
105 | |
---|
106 | #--------------------------------------------- |
---|
107 | # Find triangle that contains the point points |
---|
108 | # and print to file |
---|
109 | points = [0.,0.] |
---|
110 | for n in range(len(domain.triangles)): |
---|
111 | t = domain.triangles[n] |
---|
112 | tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]] |
---|
113 | |
---|
114 | if inside_polygon(points,tri): |
---|
115 | print 'Point is within triangle with vertices '+'%s'%tri |
---|
116 | n_point = n |
---|
117 | |
---|
118 | print 'n_point = ',n_point |
---|
119 | t = domain.triangles[n_point] |
---|
120 | tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]] |
---|
121 | |
---|
122 | filename=domain.filename |
---|
123 | file = open(filename,'w') |
---|
124 | |
---|
125 | #---------- |
---|
126 | # Evolution |
---|
127 | import time |
---|
128 | t0 = time.time() |
---|
129 | |
---|
130 | time_array = [] |
---|
131 | stage_array = [] |
---|
132 | Stage = domain.quantities['stage'] |
---|
133 | Xmomentum = domain.quantities['xmomentum'] |
---|
134 | Ymomentum = domain.quantities['ymomentum'] |
---|
135 | |
---|
136 | for t in domain.evolve(yieldstep = 20.0, finaltime = 17700.0 ): |
---|
137 | domain.write_time() |
---|
138 | |
---|
139 | tri_array = asarray(tri) |
---|
140 | t_array = asarray([[0,1,2]]) |
---|
141 | interp = Interpolation(tri_array,t_array,[points]) |
---|
142 | |
---|
143 | |
---|
144 | stage = Stage.get_values(location='centroids',indices=[n_point])[0] |
---|
145 | xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])[0] |
---|
146 | ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])[0] |
---|
147 | file.write( '%10.6f %10.6f %10.6f %10.6f\n'%(t,stage,xmomentum,ymomentum) ) |
---|
148 | |
---|
149 | time_array.append(t) |
---|
150 | stage_array.append(stage) |
---|
151 | |
---|
152 | file.close() |
---|
153 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
154 | |
---|
155 | |
---|
156 | from pylab import * |
---|
157 | ion() |
---|
158 | hold(False) |
---|
159 | plot(time_array, stage_array, 'r.-') |
---|
160 | #title('Gauge %s' %name) |
---|
161 | xlabel('time(s)') |
---|
162 | ylabel('stage (m)') |
---|
163 | legend(('Observed', 'Modelled'), shadow=True, loc='upper left') |
---|
164 | #savefig(name, dpi = 300) |
---|
165 | |
---|
166 | #raw_input('Next') |
---|
167 | show() |
---|
168 | |
---|
169 | |
---|