[7837] | 1 | import os |
---|
| 2 | from scipy.special import jn |
---|
| 3 | from scipy import sin, cos, sqrt, linspace, pi, zeros |
---|
| 4 | from rootsearch import * |
---|
[7933] | 5 | from bisect_function import * |
---|
[7922] | 6 | from Numeric import Float |
---|
| 7 | from numpy import zeros,dot |
---|
[7837] | 8 | from gaussPivot import * |
---|
| 9 | from config import g |
---|
| 10 | |
---|
| 11 | def newtonRaphson2(f,q,tol=1.0e-15): ##1.0e-9 may be too large. |
---|
| 12 | for i in range(30): |
---|
| 13 | h = 1.0e-15 ##1.0e-4 may be too large. |
---|
| 14 | n = len(q) |
---|
| 15 | jac = zeros((n,n),Float) |
---|
| 16 | if 1.0+q[0]-x<0.0: |
---|
| 17 | temp1 = 1.0+q[0]-x |
---|
| 18 | q[0] = q[0]-temp1 |
---|
| 19 | q[1] = v |
---|
| 20 | return q |
---|
| 21 | f0 = f(q) |
---|
| 22 | for i in range(n): |
---|
| 23 | temp = q[i] |
---|
| 24 | q[i] = temp + h |
---|
| 25 | f1 = f(q) |
---|
| 26 | q[i] = temp |
---|
| 27 | jac[:,i] = (f1 - f0)/h |
---|
| 28 | if sqrt(dot(f0,f0)/len(q)) < tol: return q |
---|
| 29 | dq = gaussPivot(jac,-f0) |
---|
| 30 | q = q + dq |
---|
| 31 | if sqrt(dot(dq,dq)) < tol*max(max(abs(q)),1.0): return q |
---|
| 32 | print 'Too many iterations' |
---|
| 33 | |
---|
| 34 | |
---|
| 35 | def j0(x): |
---|
| 36 | return jn(0.0, x) |
---|
| 37 | |
---|
| 38 | def j1(x): |
---|
| 39 | return jn(1.0, x) |
---|
| 40 | |
---|
| 41 | def j2(x): |
---|
| 42 | return jn(2.0, x) |
---|
| 43 | |
---|
| 44 | def j3(x): |
---|
| 45 | return jn(3.0, x) |
---|
| 46 | |
---|
| 47 | def jm1(x): |
---|
| 48 | return jn(-1.0, x) |
---|
| 49 | |
---|
| 50 | def jm2(x): |
---|
| 51 | return jn(-2.0, x) |
---|
| 52 | |
---|
| 53 | def bed(x): |
---|
| 54 | return x-1.0 |
---|
| 55 | |
---|
| 56 | def prescribe(x,t): |
---|
| 57 | q = zeros(2, Float) |
---|
| 58 | def fun(q): #Here q=(z, u) |
---|
| 59 | f = zeros(2,Float) |
---|
| 60 | f[0] = q[0] + 0.5*q[1]**2.0 - A*j0(4.0*pi/T*(1.0+q[0]-x)**0.5)*cos(2.0*pi/T*(t+q[1])) |
---|
| 61 | f[1] = q[1] + A*j1(4.0*pi/T*(1.0+q[0]-x)**0.5)*sin(2.0*pi/T*(t+q[1]))/(1+q[0]-x)**0.5 |
---|
| 62 | return f |
---|
| 63 | def newtonRaphson2(f,q,tol=1.0e-15): ##1.0e-9 may be too large. |
---|
| 64 | for i in range(30): |
---|
| 65 | h = 1.0e-15 ##1.0e-4 may be too large. |
---|
| 66 | n = len(q) |
---|
| 67 | jac = zeros((n,n),Float) |
---|
| 68 | if 1.0+q[0]-x<0.0: |
---|
| 69 | temp1 = 1.0+q[0]-x |
---|
| 70 | q[0] = q[0]-temp1 |
---|
| 71 | q[1] = v |
---|
| 72 | return q |
---|
| 73 | f0 = f(q) |
---|
| 74 | for i in range(n): |
---|
| 75 | temp = q[i] |
---|
| 76 | q[i] = temp + h |
---|
| 77 | f1 = f(q) |
---|
| 78 | q[i] = temp |
---|
| 79 | jac[:,i] = (f1 - f0)/h |
---|
| 80 | if sqrt(dot(f0,f0)/len(q)) < tol: return q |
---|
| 81 | dq = gaussPivot(jac,-f0) |
---|
| 82 | q = q + dq |
---|
| 83 | if sqrt(dot(dq,dq)) < tol*max(max(abs(q)),1.0): return q |
---|
| 84 | print 'Too many iterations' |
---|
| 85 | q = newtonRaphson2(fun,q) |
---|
| 86 | return q[0], q[1] |
---|
| 87 | |
---|
| 88 | |
---|
| 89 | def root_g(a,b,t): |
---|
| 90 | dx = 0.01 |
---|
| 91 | def g(u): |
---|
| 92 | #It was u + 8.0*A*pi/T*sin(2.0*pi/T*(t+u)). See equation (10) in Johns. Use L'Hospital rule. |
---|
| 93 | #Note that there is misprint in equation (10) in Johns. |
---|
| 94 | return u + 2.0*A*pi/T*sin(2.0*pi/T*(t+u)) |
---|
| 95 | while 1: |
---|
| 96 | x1,x2 = rootsearch(g,a,b,dx) |
---|
| 97 | if x1 != None: |
---|
| 98 | a = x2 |
---|
| 99 | root = bisect(g,x1,x2,1) |
---|
| 100 | else: |
---|
| 101 | break |
---|
| 102 | return root |
---|
| 103 | |
---|
| 104 | |
---|
| 105 | def shore(t): |
---|
| 106 | a = -0.2#-1.0 |
---|
| 107 | b = 0.2#1.0 |
---|
| 108 | #dx = 0.01 |
---|
| 109 | u = root_g(a,b,t) |
---|
| 110 | xi = -0.5*u*u + A*cos(2.0*pi/T*(t+u)) |
---|
| 111 | position = 1.0 + xi |
---|
| 112 | return position, u |
---|
| 113 | |
---|
| 114 | |
---|
| 115 | def w_at_O(t): |
---|
| 116 | return eps*cos(2.0*pi*t/T) |
---|
| 117 | |
---|
| 118 | def u_at_O(t): |
---|
| 119 | a = -1.01#-1.0 |
---|
| 120 | b = 1.01#1.0 |
---|
| 121 | dx = 0.01 |
---|
| 122 | w = w_at_O(t) |
---|
| 123 | def fun(u): |
---|
| 124 | return u + A*j1(4.0*pi/T*(1.0+w)**0.5)*sin(2.0*pi/T*(t+u))/(1+w)**0.5 |
---|
| 125 | while 1: |
---|
| 126 | x1,x2 = rootsearch(fun,a,b,dx) |
---|
| 127 | if x1 != None: |
---|
| 128 | a = x2 |
---|
| 129 | root = bisect(fun,x1,x2,1) |
---|
| 130 | else: |
---|
| 131 | break |
---|
| 132 | return root |
---|
| 133 | |
---|
| 134 | |
---|
| 135 | ##==========================================================================## |
---|
| 136 | #DIMENSIONAL PARAMETERS |
---|
| 137 | L = 5e4 # Length of channel (m) |
---|
| 138 | h_0 = 5e2 # Height at origin when the water is still |
---|
| 139 | #N = 100#400 # Number of computational cells |
---|
| 140 | #cell_len = 1.1*L/N # Origin = 0.0 |
---|
| 141 | Tp = 15.0*60.0 # Period of oscillation |
---|
| 142 | a = 1.0 # Amplitude at origin |
---|
| 143 | #X_dimensionless = linspace(0.0, 1.1*L, N) # Discretized spatial domain |
---|
| 144 | ##=========================================================================## |
---|
| 145 | #DIMENSIONLESS PARAMETERS |
---|
| 146 | eps = a/h_0 |
---|
| 147 | T = Tp*sqrt(g*h_0)/L |
---|
| 148 | A = eps/j0(4.0*pi/T) |
---|
| 149 | Time = linspace(0.0,T,1000) |
---|
| 150 | #X = X_dimensionless/L |
---|
| 151 | #Z = bed(X) |
---|
| 152 | #N_X = len(X) |
---|
| 153 | N_T = len(Time) |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | Stage = zeros(N_T, Float) |
---|
| 157 | Veloc = zeros(N_T, Float) |
---|
| 158 | for i in range(N_T): |
---|
| 159 | t=Time[i] |
---|
| 160 | zet, vel = prescribe(0.0,t) |
---|
| 161 | Stage[i] = zet |
---|
| 162 | Veloc[i] = vel |
---|
| 163 | |
---|
| 164 | |
---|
| 165 | Stage_johns = zeros(N_T, Float) |
---|
| 166 | Veloc_johns = zeros(N_T, Float) |
---|
| 167 | for i in range(N_T): |
---|
| 168 | t=Time[i] |
---|
| 169 | Stage_johns[i] = w_at_O(t) |
---|
| 170 | Veloc_johns[i] = u_at_O(t) |
---|
| 171 | |
---|
| 172 | |
---|
| 173 | num=len(Stage) |
---|
| 174 | error_w=(1.0/num)*sum(abs(Stage-Stage_johns))*h_0 |
---|
| 175 | error_u=(1.0/num)*sum(abs(Veloc-Veloc_johns))*sqrt(g*h_0) |
---|
| 176 | print "error_w=", error_w |
---|
| 177 | print "error_u=", error_u |
---|
| 178 | |
---|
| 179 | """ |
---|
| 180 | |
---|
| 181 | from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot |
---|
| 182 | |
---|
| 183 | hold(False) |
---|
| 184 | clf() |
---|
| 185 | plot1 = subplot(111) |
---|
| 186 | plot(Time/T,Stage*h_0,'b-', Time/T,Stage_johns*h_0,'k--') |
---|
| 187 | xlabel('t/T') |
---|
| 188 | ylabel('Stage') |
---|
| 189 | plot1.set_xlim([0.000,0.030]) |
---|
| 190 | plot1.set_ylim([0.980,1.005]) #([-9.0e-3,9.0e-3]) |
---|
| 191 | legend(('C-G', 'Johns'), |
---|
| 192 | 'lower left', shadow=False) |
---|
| 193 | |
---|
| 194 | plot2 = subplot(212) |
---|
| 195 | plot(Time/T,Veloc*sqrt(g*h_0),'b-', Time/T,Veloc_johns*sqrt(g*h_0),'k--') |
---|
| 196 | xlabel('t/T') |
---|
| 197 | ylabel('Velocity') |
---|
| 198 | #plot1.set_xlim([0.0,1.1]) |
---|
| 199 | #plot2.set_ylim([-0.05,0.05]) #([-1.0e-12,1.0e-12]) |
---|
| 200 | legend(('C-G', 'Johns'), |
---|
| 201 | 'upper right', shadow=False) |
---|
| 202 | |
---|
| 203 | |
---|
| 204 | filename = "discrepancy-closer" |
---|
| 205 | #filename += str(i) |
---|
| 206 | filename += ".eps" |
---|
| 207 | savefig(filename) |
---|
| 208 | #show() |
---|
| 209 | |
---|
| 210 | #plot(Time,Vel_at_O) |
---|
| 211 | #show() |
---|
| 212 | """ |
---|