1 | from scipy import sin, cos, sqrt, linspace, pi, dot |
---|
2 | from Numeric import zeros, Float, array |
---|
3 | from gaussPivot import * |
---|
4 | from analytical_prescription import * |
---|
5 | from parameter import * |
---|
6 | import os, time, csv, pprint |
---|
7 | from domain_cg import * |
---|
8 | from config import g, epsilon |
---|
9 | from rootsearch import * |
---|
10 | from bisect import * |
---|
11 | |
---|
12 | #Analytical computations################################################################# |
---|
13 | def root_g(a,b,t): |
---|
14 | dx = 0.01 |
---|
15 | def g(u): |
---|
16 | return u + 2.0*A*pi/T*sin(2.0*pi/T*(t+u)) |
---|
17 | while 1: |
---|
18 | x1,x2 = rootsearch(g,a,b,dx) |
---|
19 | if x1 != None: |
---|
20 | a = x2 |
---|
21 | root = bisect(g,x1,x2,1) |
---|
22 | else: |
---|
23 | break |
---|
24 | return root |
---|
25 | def shore(t): |
---|
26 | a = -1.0#-0.2#-1.0 |
---|
27 | b = 1.0#0.2#1.0 |
---|
28 | #dx = 0.01 |
---|
29 | u = root_g(a,b,t) |
---|
30 | xi = -0.5*u*u + A*cos(2.0*pi/T*(t+u)) |
---|
31 | position = 1.0 + xi |
---|
32 | return position, u |
---|
33 | |
---|
34 | |
---|
35 | |
---|
36 | |
---|
37 | #Numerical computations################################################################### |
---|
38 | def newtonRaphson2(f,q,tol=1.0e-15): ##1.0e-9 may be too large |
---|
39 | for i in range(30): |
---|
40 | h = 1.0e-4 ##1.0e-4 may be too large. |
---|
41 | n = len(q) |
---|
42 | jac = zeros((n,n),Float) |
---|
43 | if 1.0+q[0]-x<0.0: |
---|
44 | temp1 = 1.0+q[0]-x |
---|
45 | q[0] = q[0]-temp1 |
---|
46 | q[1] = v |
---|
47 | return q |
---|
48 | f0 = f(q) |
---|
49 | for i in range(n): |
---|
50 | temp = q[i] |
---|
51 | q[i] = temp + h |
---|
52 | f1 = f(q) |
---|
53 | q[i] = temp |
---|
54 | jac[:,i] = (f1 - f0)/h |
---|
55 | if sqrt(dot(f0,f0)/len(q)) < tol: return q |
---|
56 | dq = gaussPivot(jac,-f0) |
---|
57 | q = q + dq |
---|
58 | if sqrt(dot(dq,dq)) < tol*max(max(abs(q)),1.0): return q |
---|
59 | print 'Too many iterations' |
---|
60 | |
---|
61 | def elevation(X): |
---|
62 | N = len(X) |
---|
63 | z = zeros(N,Float) |
---|
64 | for i in range(N): |
---|
65 | z[i] = (h_0/L)*X[i] - h_0 |
---|
66 | return z |
---|
67 | |
---|
68 | def height(X): |
---|
69 | N = len(X) |
---|
70 | z = zeros(N,Float) |
---|
71 | for i in range(N): |
---|
72 | z[i] = h_0 - (h_0/L)*X[i] |
---|
73 | return z |
---|
74 | |
---|
75 | def velocity(X): |
---|
76 | N = len(X) |
---|
77 | return zeros(N,Float) |
---|
78 | |
---|
79 | boundary = { (0,0): 'left',(N-1,1): 'right'} |
---|
80 | |
---|
81 | domain = Domain(points,boundary) |
---|
82 | domain.order = 2 |
---|
83 | domain.set_timestepping_method('rk2') |
---|
84 | domain.set_CFL(1.0) |
---|
85 | domain.beta = 1.0 |
---|
86 | domain.set_limiter("minmod") |
---|
87 | |
---|
88 | |
---|
89 | def f_CG(t): |
---|
90 | timing = t*sqrt(g*h_0)/L |
---|
91 | w, u = prescribe(0.0,timing) |
---|
92 | wO = w*h_0 |
---|
93 | uO = u*sqrt(g*h_0) |
---|
94 | zO = -h_0 |
---|
95 | hO = wO - zO |
---|
96 | pO = uO * hO |
---|
97 | #[ 'stage', 'xmomentum', 'elevation', 'height', 'velocity'] |
---|
98 | return [wO, pO, zO, hO, uO] |
---|
99 | def f_JOHNS(t): |
---|
100 | timing = t*sqrt(g*h_0)/L |
---|
101 | w, u = prescribe_at_O_JOHNS(timing) |
---|
102 | wO = w*h_0 |
---|
103 | uO = u*sqrt(g*h_0) |
---|
104 | zO = -h_0 |
---|
105 | hO = wO - zO |
---|
106 | pO = uO * hO |
---|
107 | #[ 'stage', 'xmomentum', 'elevation', 'height', 'velocity'] |
---|
108 | return [wO, pO, zO, hO, uO] |
---|
109 | |
---|
110 | T1 = Time_boundary(domain,f_CG) |
---|
111 | D2 = Dirichlet_boundary([50.5, 0.0, 50.5, 0.0, 0.0]) |
---|
112 | domain.set_boundary({'left':T1,'right':D2}) |
---|
113 | |
---|
114 | domain.set_quantity('height',height) |
---|
115 | domain.set_quantity('elevation',elevation) |
---|
116 | domain.set_quantity('velocity',velocity) |
---|
117 | |
---|
118 | |
---|
119 | Ver = domain.vertices |
---|
120 | n_V = len(Ver) |
---|
121 | AnalitW_V = zeros((n_V,2), Float) |
---|
122 | AnalitP_V = zeros((n_V,2), Float) |
---|
123 | AnalitZ_V = zeros((n_V,2), Float) |
---|
124 | AnalitH_V = zeros((n_V,2), Float) |
---|
125 | AnalitU_V = zeros((n_V,2), Float) |
---|
126 | |
---|
127 | Cen = domain.centroids |
---|
128 | n_C = len(Cen) |
---|
129 | AnalitW_C = zeros(n_C, Float) |
---|
130 | AnalitP_C = zeros(n_C, Float) |
---|
131 | AnalitZ_C = zeros(n_C, Float) |
---|
132 | AnalitH_C = zeros(n_C, Float) |
---|
133 | AnalitU_C = zeros(n_C, Float) |
---|
134 | |
---|
135 | waktu = 10.0 #3.0*60.0 |
---|
136 | WAKTU = 12690.0 #20.0*Tp #Note: Tp=15.0*60.0 |
---|
137 | yieldstep = finaltime = waktu |
---|
138 | t0 = time.time() |
---|
139 | counter=1 |
---|
140 | |
---|
141 | shorelines_numerical_cg = zeros(int(WAKTU/waktu), Float) |
---|
142 | shorelines_analytical = zeros(int(WAKTU/waktu), Float) |
---|
143 | time_instants = zeros(int(WAKTU/waktu), Float) |
---|
144 | |
---|
145 | |
---|
146 | while finaltime < WAKTU+0.1: |
---|
147 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
148 | domain.write_time() |
---|
149 | time_instants[counter-1] = domain.time |
---|
150 | |
---|
151 | Stage = domain.quantities['stage'] |
---|
152 | Momentum = domain.quantities['xmomentum'] |
---|
153 | Elevation = domain.quantities['elevation'] |
---|
154 | Height = domain.quantities['height'] |
---|
155 | Velocity = domain.quantities['velocity'] |
---|
156 | |
---|
157 | StageV = Stage.vertex_values |
---|
158 | MomV = Momentum.vertex_values |
---|
159 | ElevationV = Elevation.vertex_values |
---|
160 | HeightV = Height.vertex_values |
---|
161 | VelV = Velocity.vertex_values |
---|
162 | |
---|
163 | StageC = Stage.centroid_values |
---|
164 | MomC = Momentum.centroid_values |
---|
165 | ElevationC = Elevation.centroid_values |
---|
166 | HeightC = Height.centroid_values |
---|
167 | VelC = Velocity.centroid_values |
---|
168 | |
---|
169 | table = zeros((len(Ver.flat),6),Float) |
---|
170 | for r in range(len(Ver.flat)): |
---|
171 | for c in range(6): |
---|
172 | if c==0: |
---|
173 | table[r][c] = Ver.flat[r] |
---|
174 | elif c==1: |
---|
175 | table[r][c] = StageV.flat[r] |
---|
176 | elif c==2: |
---|
177 | table[r][c] = MomV.flat[r] |
---|
178 | elif c==3: |
---|
179 | table[r][c] = ElevationV.flat[r] |
---|
180 | elif c==4: |
---|
181 | table[r][c] = HeightV.flat[r] |
---|
182 | else: |
---|
183 | table[r][c] = VelV.flat[r] |
---|
184 | |
---|
185 | outname = "%s%04i%s%f%s" %("numerical_cg_", counter, "_", domain.time, ".csv") |
---|
186 | outfile = open(outname, 'w') |
---|
187 | writer = csv.writer(outfile) |
---|
188 | for row in table: |
---|
189 | writer.writerow(row) |
---|
190 | outfile.close() |
---|
191 | |
---|
192 | |
---|
193 | for s in range(2*n_V): |
---|
194 | heiL = HeightV.flat[s] |
---|
195 | momR = MomV.flat[s+1] |
---|
196 | if heiL >= 1e-6: |
---|
197 | if abs(momR)==0.0: #<1e-15: |
---|
198 | break |
---|
199 | #print "s+1=",s+1 |
---|
200 | shorelines_numerical_cg[counter-1] = Ver.flat[s+1] |
---|
201 | #print "shorelines_numerical_cg=",shorelines_numerical_cg[counter-1] |
---|
202 | |
---|
203 | |
---|
204 | #print "Now the ANALYTIC" |
---|
205 | pos_shore, vel_shore = shore(domain.time*sqrt(g*h_0)/L) #dimensionless |
---|
206 | pos_shore = pos_shore*L |
---|
207 | vel_shore = vel_shore*sqrt(g*h_0) |
---|
208 | |
---|
209 | shorelines_analytical[counter-1] = pos_shore |
---|
210 | |
---|
211 | #The following is for calculating the error at centroids |
---|
212 | for i in range(n_C): |
---|
213 | x = Cen[i] |
---|
214 | if x < pos_shore: |
---|
215 | sta, vel = prescribe(x/L,domain.time*sqrt(g*h_0)/L) #dimensionless |
---|
216 | AnalitW_C[i] = sta*h_0 |
---|
217 | AnalitU_C[i] = vel #It needs dimensionalisation |
---|
218 | AnalitZ_C[i] = (h_0/L)*x - h_0 |
---|
219 | AnalitH_C[i] = AnalitW_C[i] - AnalitZ_C[i] |
---|
220 | AnalitP_C[i] = AnalitH_C[i]*AnalitU_C[i] #It needs dimensionalisation |
---|
221 | else: |
---|
222 | AnalitW_C[i] = (h_0/L)*x - h_0 |
---|
223 | AnalitU_C[i] = 0.0 |
---|
224 | AnalitZ_C[i] = (h_0/L)*x - h_0 |
---|
225 | AnalitH_C[i] = 0.0 |
---|
226 | AnalitP_C[i] = 0.0 |
---|
227 | AnalitU_C = AnalitU_C*sqrt(g*h_0) #This is the dimensionalisation |
---|
228 | AnalitP_C = AnalitP_C*sqrt(g*h_0) #This is the dimensionalisation |
---|
229 | |
---|
230 | error_W = (1.0/n_C)*sum(abs(StageC-AnalitW_C)) |
---|
231 | error_P = (1.0/n_C)*sum(abs(MomC-AnalitP_C)) |
---|
232 | error_U = (1.0/n_C)*sum(abs(VelC-AnalitU_C)) |
---|
233 | table = zeros((1,4),Float) |
---|
234 | for r in range(1): |
---|
235 | for c in range(4): |
---|
236 | if c==0: |
---|
237 | table[r][c] = domain.time |
---|
238 | elif c==1: |
---|
239 | table[r][c] = error_W |
---|
240 | elif c==2: |
---|
241 | table[r][c] = error_P |
---|
242 | else: |
---|
243 | table[r][c] = error_U |
---|
244 | outname = "%s%04i%s%f%s" %("error_cg_", counter, "_", domain.time, ".csv") |
---|
245 | outfile = open(outname, 'w') |
---|
246 | writer = csv.writer(outfile) |
---|
247 | for row in table: |
---|
248 | writer.writerow(row) |
---|
249 | outfile.close() |
---|
250 | |
---|
251 | |
---|
252 | #The following is for ploting the quantities at vertex values |
---|
253 | for i in range(n_V): |
---|
254 | vector_x = Ver[i] |
---|
255 | for k in range(2): |
---|
256 | x = vector_x[k] |
---|
257 | if x < pos_shore: |
---|
258 | sta, vel = prescribe(x/L,domain.time*sqrt(g*h_0)/L) #dimensionless |
---|
259 | AnalitW_V[i,k] = sta*h_0 |
---|
260 | AnalitU_V[i,k] = vel #It needs dimensionalisation |
---|
261 | AnalitZ_V[i,k] = (h_0/L)*x - h_0 |
---|
262 | AnalitH_V[i,k] = AnalitW_V[i,k] - AnalitZ_V[i,k] |
---|
263 | AnalitP_V[i,k] = AnalitH_V[i,k]*AnalitU_V[i,k] #It needs dimensionalisation |
---|
264 | else: |
---|
265 | AnalitW_V[i,k] = (h_0/L)*x - h_0 |
---|
266 | AnalitU_V[i,k] = 0.0 |
---|
267 | AnalitZ_V[i,k] = (h_0/L)*x - h_0 |
---|
268 | AnalitH_V[i,k] = 0.0 |
---|
269 | AnalitP_V[i,k] = 0.0 |
---|
270 | AnalitU_V = AnalitU_V*sqrt(g*h_0) #This is the dimensionalisation |
---|
271 | AnalitP_V = AnalitP_V*sqrt(g*h_0) #This is the dimensionalisation |
---|
272 | |
---|
273 | |
---|
274 | table = zeros((len(Ver.flat),6),Float) |
---|
275 | for r in range(len(Ver.flat)): |
---|
276 | for c in range(6): |
---|
277 | if c==0: |
---|
278 | table[r][c] = Ver.flat[r] |
---|
279 | elif c==1: |
---|
280 | table[r][c] = AnalitW_V.flat[r] |
---|
281 | elif c==2: |
---|
282 | table[r][c] = AnalitP_V.flat[r] |
---|
283 | elif c==3: |
---|
284 | table[r][c] = AnalitZ_V.flat[r] |
---|
285 | elif c==4: |
---|
286 | table[r][c] = AnalitH_V.flat[r] |
---|
287 | else: |
---|
288 | table[r][c] = AnalitU_V.flat[r] |
---|
289 | |
---|
290 | outname = "%s%04i%s%f%s" %("analytical_", counter, "_", domain.time, ".csv") |
---|
291 | outfile = open(outname, 'w') |
---|
292 | writer = csv.writer(outfile) |
---|
293 | for row in table: |
---|
294 | writer.writerow(row) |
---|
295 | outfile.close() |
---|
296 | |
---|
297 | |
---|
298 | from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion |
---|
299 | hold(False) |
---|
300 | clf() |
---|
301 | |
---|
302 | plot1 = subplot(311) |
---|
303 | plot(Ver/1e+4,AnalitW_V,'b-', Ver/1e+4,StageV,'g-', Ver/1e+4,ElevationV,'k-') |
---|
304 | #plot(Ver,StageV, Ver,ElevationV) |
---|
305 | #plot(Ver/L,StageV/h_0, Ver/L,ElevationV/h_0) |
---|
306 | #xlabel('Position') |
---|
307 | ylabel('Stage') |
---|
308 | #plot1.set_xlim([0.0,1.2]) |
---|
309 | plot1.set_ylim([-6.0,6.0])#([-9.0e-3,9.0e-3]) |
---|
310 | #legend(('Analytical Solution', 'Numerical Solution', 'Discretized Bed'), |
---|
311 | # 'upper right', shadow=False) |
---|
312 | |
---|
313 | plot2 = subplot(312) |
---|
314 | plot(Ver/1e+4,AnalitP_V,'b-', Ver/1e+4,MomV,'g-') |
---|
315 | #plot(Ver/L, VelV/sqrt(g*h_0)) |
---|
316 | #xlabel('Position') |
---|
317 | ylabel('Momentum') |
---|
318 | #plot2.set_xlim([0.0,1.2]) |
---|
319 | #legend(('Analytical Solution','Numerical Solution'), |
---|
320 | # 'upper right', shadow=False) |
---|
321 | |
---|
322 | plot3 = subplot(313) |
---|
323 | plot(Ver/1e+4,AnalitU_V,'b-', Ver/1e+4,VelV,'g-') |
---|
324 | #plot(Ver/L, VelV/sqrt(g*h_0)) |
---|
325 | xlabel('Position / 10,000') |
---|
326 | ylabel('Velocity') |
---|
327 | #plot2.set_xlim([0.0,1.2]) |
---|
328 | legend(('Analytical Solution','Numerical Solution'), |
---|
329 | 'upper center', shadow=False) |
---|
330 | |
---|
331 | #filename = "%s%04i%s%f%s" %("numerical_cg_", counter,"_", domain.time, ".png") |
---|
332 | #savefig(filename) |
---|
333 | #show() |
---|
334 | #raw_input("Press ENTER to continue") |
---|
335 | |
---|
336 | counter = counter+1 |
---|
337 | finaltime = finaltime + waktu |
---|
338 | |
---|
339 | |
---|
340 | table = zeros((int(WAKTU/waktu), 2),Float) |
---|
341 | for r in range(int(WAKTU/waktu)): |
---|
342 | for c in range(2): |
---|
343 | if c==0: |
---|
344 | table[r][c] = time_instants[r] |
---|
345 | else: |
---|
346 | table[r][c] = shorelines_numerical_cg[r] |
---|
347 | |
---|
348 | outname = "%s" %("shore_numerical_cg.csv") |
---|
349 | outfile = open(outname, 'w') |
---|
350 | writer = csv.writer(outfile) |
---|
351 | for row in table: |
---|
352 | writer.writerow(row) |
---|
353 | outfile.close() |
---|
354 | |
---|
355 | |
---|
356 | table = zeros((int(WAKTU/waktu), 2),Float) |
---|
357 | for r in range(int(WAKTU/waktu)): |
---|
358 | for c in range(2): |
---|
359 | if c==0: |
---|
360 | table[r][c] = time_instants[r] |
---|
361 | else: |
---|
362 | table[r][c] = shorelines_analytical[r] |
---|
363 | |
---|
364 | outname = "%s" %("shore_analytical.csv") |
---|
365 | outfile = open(outname, 'w') |
---|
366 | writer = csv.writer(outfile) |
---|
367 | for row in table: |
---|
368 | writer.writerow(row) |
---|
369 | outfile.close() |
---|
370 | |
---|