source: trunk/anuga_work/development/sudi/sw_1d/periodic_waves/johns/run-discrepancy.py @ 7922

Last change on this file since 7922 was 7922, checked in by mungkasi, 13 years ago

Modifying codes so that arrays are compatible with numpy instead of Numeric.

File size: 5.5 KB
Line 
1import os
2from scipy.special import jn
3from scipy import sin, cos, sqrt, linspace, pi, zeros
4from rootsearch import *
5from bisect import *
6from Numeric import Float
7from numpy import zeros,dot
8from gaussPivot import *
9from config import g
10
11def 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
35def j0(x):
36    return jn(0.0, x)
37
38def j1(x):
39    return jn(1.0, x)
40
41def j2(x):
42    return jn(2.0, x)
43
44def j3(x):
45    return jn(3.0, x)
46
47def jm1(x):
48    return jn(-1.0, x)
49
50def jm2(x):
51    return jn(-2.0, x)
52
53def bed(x):
54    return x-1.0
55
56def 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
89def 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
105def 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
115def w_at_O(t):
116    return eps*cos(2.0*pi*t/T)
117
118def 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
137L = 5e4                 # Length of channel (m)
138h_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
141Tp = 15.0*60.0                           # Period of oscillation
142a = 1.0                                 # Amplitude at origin
143#X_dimensionless = linspace(0.0, 1.1*L, N)             # Discretized spatial domain
144##=========================================================================##
145#DIMENSIONLESS PARAMETERS
146eps = a/h_0
147T = Tp*sqrt(g*h_0)/L
148A = eps/j0(4.0*pi/T)
149Time = linspace(0.0,T,1000)
150#X = X_dimensionless/L
151#Z = bed(X)
152#N_X = len(X)
153N_T = len(Time)
154
155
156Stage = zeros(N_T, Float)   
157Veloc = zeros(N_T, Float)
158for 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
165Stage_johns = zeros(N_T, Float)   
166Veloc_johns = zeros(N_T, Float)
167for 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
173num=len(Stage)
174error_w=(1.0/num)*sum(abs(Stage-Stage_johns))*h_0
175error_u=(1.0/num)*sum(abs(Veloc-Veloc_johns))*sqrt(g*h_0)
176print "error_w=", error_w
177print "error_u=", error_u
178
179"""
180
181from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
182
183hold(False)
184clf()
185plot1 = subplot(111)
186plot(Time/T,Stage*h_0,'b-', Time/T,Stage_johns*h_0,'k--')
187xlabel('t/T')
188ylabel('Stage')
189plot1.set_xlim([0.000,0.030])
190plot1.set_ylim([0.980,1.005]) #([-9.0e-3,9.0e-3])
191legend(('C-G', 'Johns'),
192       'lower left', shadow=False)
193
194plot2 = subplot(212)
195plot(Time/T,Veloc*sqrt(g*h_0),'b-', Time/T,Veloc_johns*sqrt(g*h_0),'k--')
196xlabel('t/T')
197ylabel('Velocity')
198#plot1.set_xlim([0.0,1.1])
199#plot2.set_ylim([-0.05,0.05]) #([-1.0e-12,1.0e-12])
200legend(('C-G', 'Johns'),
201       'upper right', shadow=False)
202
203
204filename = "discrepancy-closer"
205#filename += str(i)
206filename += ".eps"
207savefig(filename)
208#show()
209
210#plot(Time,Vel_at_O)
211#show()
212"""
Note: See TracBrowser for help on using the repository browser.