1 | import os |
---|
2 | from scipy.special import jn |
---|
3 | from scipy import sin, cos, sqrt, linspace, pi, zeros |
---|
4 | from rootsearch import * |
---|
5 | from bisect_function import * |
---|
6 | from Numeric import Float |
---|
7 | from numpy import zeros,dot |
---|
8 | from gaussPivot import * |
---|
9 | from config import g |
---|
10 | from analytical_prescription import * |
---|
11 | |
---|
12 | |
---|
13 | def j0(x): |
---|
14 | return jn(0.0, x) |
---|
15 | |
---|
16 | def j1(x): |
---|
17 | return jn(1.0, x) |
---|
18 | |
---|
19 | def j2(x): |
---|
20 | return jn(2.0, x) |
---|
21 | |
---|
22 | def j3(x): |
---|
23 | return jn(3.0, x) |
---|
24 | |
---|
25 | def jm1(x): |
---|
26 | return jn(-1.0, x) |
---|
27 | |
---|
28 | def jm2(x): |
---|
29 | return jn(-2.0, x) |
---|
30 | |
---|
31 | def bed(x): |
---|
32 | return x-1.0 |
---|
33 | |
---|
34 | |
---|
35 | def w_at_O(t): |
---|
36 | return eps*cos(2.0*pi*t/T) |
---|
37 | |
---|
38 | def u_at_O(t): |
---|
39 | a = -1.01#-1.0 |
---|
40 | b = 1.01#1.0 |
---|
41 | dx = 0.01 |
---|
42 | w = w_at_O(t) |
---|
43 | def fun(u): |
---|
44 | return u + A*j1(4.0*pi/T*(1.0+w)**0.5)*sin(2.0*pi/T*(t+u))/(1.0+w)**0.5 |
---|
45 | while 1: |
---|
46 | x1,x2 = rootsearch(fun,a,b,dx) |
---|
47 | if x1 != None: |
---|
48 | a = x2 |
---|
49 | root = bisect(fun,x1,x2,1) |
---|
50 | else: |
---|
51 | break |
---|
52 | return root |
---|
53 | |
---|
54 | """ |
---|
55 | ##==========================================================================## |
---|
56 | #DIMENSIONAL PARAMETERS |
---|
57 | L = 5e4 # Length of channel (m) |
---|
58 | h_0 = 5e2 # Height at origin when the water is still |
---|
59 | Tp = 15.0*60.0 # Period of oscillation |
---|
60 | a = 1.0 # Amplitude at origin |
---|
61 | ##=========================================================================## |
---|
62 | #DIMENSIONLESS PARAMETERS |
---|
63 | eps = a/h_0 |
---|
64 | T = Tp*sqrt(g*h_0)/L |
---|
65 | A = eps/j0(4.0*pi/T) |
---|
66 | """ |
---|
67 | |
---|
68 | |
---|
69 | Time = linspace(0.0,T,100) |
---|
70 | N_T = len(Time) |
---|
71 | |
---|
72 | |
---|
73 | Stage = zeros(N_T, Float) |
---|
74 | Veloc = zeros(N_T, Float) |
---|
75 | for i in range(N_T): |
---|
76 | t=Time[i] |
---|
77 | zet, vel = prescribe(0.0,t) |
---|
78 | Stage[i] = zet |
---|
79 | Veloc[i] = vel |
---|
80 | |
---|
81 | |
---|
82 | Stage_johns = zeros(N_T, Float) |
---|
83 | Veloc_johns = zeros(N_T, Float) |
---|
84 | for i in range(N_T): |
---|
85 | t=Time[i] |
---|
86 | Stage_johns[i] = w_at_O(t) |
---|
87 | Veloc_johns[i] = u_at_O(t) |
---|
88 | |
---|
89 | """ |
---|
90 | num=len(Stage) |
---|
91 | error_w=(1.0/num)*sum(abs(Stage-Stage_johns))*h_0 |
---|
92 | error_u=(1.0/num)*sum(abs(Veloc-Veloc_johns))*sqrt(g*h_0) |
---|
93 | print "error_w=", error_w |
---|
94 | print "error_u=", error_u |
---|
95 | |
---|
96 | """ |
---|
97 | |
---|
98 | from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot |
---|
99 | |
---|
100 | hold(False) |
---|
101 | clf() |
---|
102 | plot1 = subplot(211) |
---|
103 | plot(Time/T,Stage*h_0,'b-', Time/T,Stage_johns*h_0,'k--') |
---|
104 | xlabel('t/T') |
---|
105 | ylabel('Stage') |
---|
106 | #plot1.set_xlim([0.000,0.030]) |
---|
107 | #plot1.set_ylim([0.980,1.005]) #([-9.0e-3,9.0e-3]) |
---|
108 | legend(('C-G', 'Johns'), |
---|
109 | 'lower left', shadow=False) |
---|
110 | |
---|
111 | plot2 = subplot(212) |
---|
112 | plot(Time/T,Veloc*sqrt(g*h_0),'b-', Time/T,Veloc_johns*sqrt(g*h_0),'k--') |
---|
113 | xlabel('t/T') |
---|
114 | ylabel('Velocity') |
---|
115 | #plot1.set_xlim([0.0,1.1]) |
---|
116 | #plot2.set_ylim([-0.05,0.05]) #([-1.0e-12,1.0e-12]) |
---|
117 | legend(('C-G', 'Johns'), |
---|
118 | 'upper right', shadow=False) |
---|
119 | |
---|
120 | |
---|
121 | #filename = "discrepancy-closer" |
---|
122 | #filename += str(i) |
---|
123 | #filename += ".eps" |
---|
124 | #savefig(filename) |
---|
125 | #show() |
---|
126 | |
---|
127 | #plot(Time,Vel_at_O) |
---|
128 | #show() |
---|