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

Last change on this file since 7837 was 7837, checked in by mungkasi, 14 years ago

Again, adding some codes for 1d problems on debris avalanche and periodic waves.

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