1 | import numpy as np |
---|
2 | from scipy import sparse as sparse |
---|
3 | import math as math |
---|
4 | |
---|
5 | import getargs |
---|
6 | log_master, log_world = getargs.getLogger(__name__) |
---|
7 | INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error |
---|
8 | INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug |
---|
9 | |
---|
10 | class Struct: pass |
---|
11 | |
---|
12 | def Linf(data): |
---|
13 | return abs(data).max() |
---|
14 | def compare(msg, data1, data2): |
---|
15 | print msg, Linf(data1-data2)/Linf(data1) |
---|
16 | |
---|
17 | ############################# Tridiagonal solver ############################# |
---|
18 | |
---|
19 | class SymTriDiag: |
---|
20 | """ Solve -c(i-1)x(i-1) + b(i)x(i) - c(i)x(i+1) = d(i) |
---|
21 | Forward sweep : |
---|
22 | C(i) = -c(i) / (b(i)+c(i-1)C(i-1)), C(0)=0 |
---|
23 | D(i) = (d(i)+c(i-1)D(i-1)) / (b(i)+c(i-1)C(i-1)), D(0)=0 |
---|
24 | Back substituion : |
---|
25 | x(i) = D(i)-C(i)x(i+1), x(N+1)=0 |
---|
26 | """ |
---|
27 | def __init__(self, b,c): |
---|
28 | M, N = b.shape |
---|
29 | self.M, self.N, self.b, self.c = M,N,b,c |
---|
30 | B,C = np.zeros((M,N)), np.zeros((M,N-1)) |
---|
31 | Bi = b[:,0] |
---|
32 | B[:,0]=Bi |
---|
33 | for i in xrange(1,N): |
---|
34 | Ci = -c[:,i-1]/Bi |
---|
35 | C[:,i-1]=Ci |
---|
36 | Bi = b[:,i]+c[:,i-1]*Ci |
---|
37 | B[:,i]=Bi |
---|
38 | self.B, self.C = B,C |
---|
39 | def dot(self,x): |
---|
40 | M,N,b,c = self.M, self.N, self.b, self.c |
---|
41 | y = np.zeros((M,N)) |
---|
42 | y[:,0] = b[:,0]*x[:,0] - c[:,0]*x[:,1] |
---|
43 | for i in xrange(1,N-1): |
---|
44 | y[:,i] = -c[:,i-1]*x[:,i-1] + b[:,i]*x[:,i] - c[:,i]*x[:,i+1] |
---|
45 | y[:,N-1] = -c[:,N-2]*x[:,N-2] + b[:,N-1]*x[:,N-1] |
---|
46 | return y |
---|
47 | def solve(self,d): |
---|
48 | M,N,c,B,C = self.M, self.N, self.c, self.B, self.C |
---|
49 | D,x = np.zeros((M,N)), np.zeros((M,N)) |
---|
50 | Di = d[:,0]/B[:,0] |
---|
51 | D[:,0]=Di |
---|
52 | for i in xrange(1,N): |
---|
53 | Di = (d[:,i]+c[:,i-1]*Di)/B[:,i] |
---|
54 | D[:,i]=Di |
---|
55 | xi=D[:,N-1] |
---|
56 | x[:,N-1]=xi |
---|
57 | for i in xrange(N-2,-1,-1): |
---|
58 | xi = D[:,i]-C[:,i]*xi |
---|
59 | x[:,i]=xi |
---|
60 | return x |
---|
61 | |
---|
62 | ########################### ARK time schemes ########################## |
---|
63 | |
---|
64 | class SIRK: |
---|
65 | """ |
---|
66 | Generic SIRK scheme with fast/slow Butcher tableaux af,as |
---|
67 | Following Weller et al. (2013) |
---|
68 | sj=slow(yj), fj=fast(yj) |
---|
69 | zj = y(n) + dt*( sum(l<j)asjl*sl + sum(l<j)afjl*fl ) |
---|
70 | yj = zj + dt*afjj*fast(yj) |
---|
71 | y(n+1) = y(n)+dt*sum(wsj*sj+wsf*fj) |
---|
72 | Reformulation : remarking that |
---|
73 | y(j) = y(n) + dt*sum_(l<=j)[as(j,l)*s(l) + af(j,l)*f(l) ] |
---|
74 | z(j+1) = y(n) + dt*sum_(l<=j)[as(j+1,l)*s(l) + af(j+1,l)*f(l) ] |
---|
75 | = y(j) + dt*sum_(l<=j)[cs(j,l)*s(l) + cf(j,l)*f(l) ] |
---|
76 | where b(j,l)=a(j,l) ; b(nstage,l)=w(l) ; c(j) = b(j+1)-b(j) => c0=b1 |
---|
77 | => implementation |
---|
78 | z(0) = y(n) |
---|
79 | DO j=0, nstage-1 |
---|
80 | [y(j), fj, sj] = bwd_fast_slow(z(j),0.) |
---|
81 | z(j+1) = y(j) + dt*sum_(l<=j)[cs(j,l)*s(l) + cf(j,l)*f(l) ] |
---|
82 | END DO |
---|
83 | y(n+1) = z(nstage) |
---|
84 | |
---|
85 | where z(j) and y(j) are actually stored in a single variable zj |
---|
86 | and [y,f,s] = bwd_fast_slow(z,tau) |
---|
87 | * solves y - tau*fast(y) = z |
---|
88 | * computes f=fast(y) and s=slow(y) |
---|
89 | """ |
---|
90 | def cjl(self,b): # converts Butcher tableau b(j,l) to update coefficients c(j,l) |
---|
91 | return np.array([ b[j+1,:]-b[j,:] for j in range(self.nstage)]) |
---|
92 | def __init__(self,precision,bwd_fast_slow, dt,nstage, bsjl,bfjl): |
---|
93 | if bwd_fast_slow is None: |
---|
94 | DEBUG('SIRK with Fortran-side time stepping') |
---|
95 | precision=np.float64 |
---|
96 | else: |
---|
97 | DEBUG('SIRK with Python-side time stepping at precision %s' % precision) |
---|
98 | |
---|
99 | self.bwd_fast_slow = bwd_fast_slow |
---|
100 | self.dt, self.nstage = dt, nstage |
---|
101 | csjl = dt*self.cjl(bsjl) |
---|
102 | cfjl = dt*self.cjl(bfjl) |
---|
103 | self.csjl, self.cfjl = [x.astype(precision) for x in csjl, cfjl] |
---|
104 | self.tauj = np.array([dt*bfjl[j,j] for j in range(self.nstage)], dtype=precision) |
---|
105 | DEBUG('Types of csjl, cfjl, tauj : %s %s %s' % (self.csjl.dtype, self.cfjl.dtype, self.tauj.dtype) ) |
---|
106 | def next(self,zj): |
---|
107 | csjl, cfjl, tauj = self.csjl, self.cfjl, self.tauj |
---|
108 | fastj, slowj = [], [] |
---|
109 | for j in range(self.nstage): # here zj is z(j) |
---|
110 | zj, fj, sj = self.bwd_fast_slow(zj, tauj[j]) |
---|
111 | fastj.append(fj) |
---|
112 | slowj.append(sj) |
---|
113 | for l in range(0,j+1): # here zj is initially y(j) and finally z(j+1) |
---|
114 | zj = [z + csjl[j,l]*s + cfjl[j,l]*f for z,s,f in zip(zj, slowj[l], fastj[l]) ] |
---|
115 | # print 'SIRK : zj has type ', [ z.dtype for z in zj] |
---|
116 | return zj |
---|
117 | def advance(self, flow, N): |
---|
118 | for k in range(N): |
---|
119 | flow=self.next(flow) |
---|
120 | return flow |
---|
121 | |
---|
122 | def SIRK1(bwd_fast_slow, dt): # untested |
---|
123 | """ Forward-Backward""" |
---|
124 | wj = np.array([1]) |
---|
125 | bsjl = np.array([[0],wj]) |
---|
126 | bfjl = np.array([[1],wj]) |
---|
127 | return SIRK(bwd_fast_slow,dt,1,bsjl,bfjl) |
---|
128 | |
---|
129 | def ARK2(bwd_fast_slow, dt, precision=None, a32 = (3.+2.*np.sqrt(2.))/6.): |
---|
130 | """ ARK2 scheme by Giraldo et al. (2013), noted ARK2(2,3,2) in Weller et al. (2013) """ |
---|
131 | gamma = 1.-1./np.sqrt(2.) |
---|
132 | delta = 1./(2.*np.sqrt(2.)) |
---|
133 | wj = np.array([delta,delta,gamma]) |
---|
134 | bsjl = np.array([[0,0,0],[2.*gamma,0,0],[1.-a32,a32,0],wj]) |
---|
135 | bfjl = np.array([[0,0,0],[gamma,gamma,0],[delta,delta,gamma],wj]) |
---|
136 | return SIRK(precision,bwd_fast_slow,dt,3,bsjl,bfjl) |
---|
137 | |
---|
138 | def UJ3(bwd_fast_slow, dt): |
---|
139 | """ Strang-carryover scheme by Ullrich et al. (2013), noted UJ3(1,3,2) in Weller et al. (2013) """ |
---|
140 | zz = [0,0,0,0,0,0] |
---|
141 | wsj = [0,1/6.,1/6.,2/3.,0,0] |
---|
142 | h = [.5,0,0,0,0,0] |
---|
143 | wfj = [.5,0,0,0,0,.5] |
---|
144 | bsjl = np.array([zz,zz,[0,1,0,0,0,0], |
---|
145 | [0,.25,.25,0,0,0],wsj,wsj,wsj]) |
---|
146 | bfjl = np.array([zz,h,h,h,h,wfj,wfj]) |
---|
147 | return SIRK(bwd_fast_slow,dt,6,bsjl,bfjl) |
---|
148 | |
---|
149 | ######### Explicit RK schemes written as ARK with two identical Butcher tableaus ######## |
---|
150 | |
---|
151 | def RK4(bwd_fast_slow, dt, precision=None): |
---|
152 | """ 4-stage 4th-order RK """ |
---|
153 | bjl = np.array([[0,0,0,0], |
---|
154 | [.5,0,0,0], |
---|
155 | [0,.5,0,0], |
---|
156 | [0.,0.,1.,0], |
---|
157 | [1/6.,1/3.,1/3.,1/6.]]) |
---|
158 | return SIRK(precision, bwd_fast_slow,dt,4,bjl,bjl) |
---|
159 | |
---|
160 | def RK1(bwd_fast_slow, dt, precision=None): |
---|
161 | """ 4-stage 4th-order RK """ |
---|
162 | bjl = np.array([[0.],[1.]]) |
---|
163 | return SIRK(precision, bwd_fast_slow,dt,1,bjl,bjl) |
---|
164 | |
---|
165 | ######### Simplified explicit RK schemes (2nd-order) ######## |
---|
166 | |
---|
167 | def RK_simple(bwd_fast_slow, dt, tau): |
---|
168 | # take coefficients of simplified RK scheme |
---|
169 | # y(n) := x + tau*f(y(n-1)) |
---|
170 | # and computes the equivalent Butcher tableau |
---|
171 | tau=np.array(tau) |
---|
172 | n = tau.size # number of stages |
---|
173 | bjl = np.zeros((n+1,n)) # a[stage, coef] |
---|
174 | for i in range(n): |
---|
175 | bjl[i+1,i]=tau[i] |
---|
176 | return SIRK(bwd_fast_slow,dt,n,bjl,bjl) |
---|
177 | |
---|
178 | def RKn_simple(n, bwd_fast_slow, dt): |
---|
179 | tau = [1./(n-i) for i in range(n)] |
---|
180 | return RK_simple(bwd_fast_slow, dt, tau) |
---|
181 | |
---|
182 | def RK_Kinnmark(K, bwd_fast_slow, dt): |
---|
183 | """ RK schemes with max stability along imaginary axis |
---|
184 | Kinnmark & Gray, 1984 """ |
---|
185 | ee=np.zeros((K+1,K+1)) |
---|
186 | ee[2,0], ee[2,1], ee[2,2] = 1,1,1 |
---|
187 | ee[3,0], ee[3,1], ee[3,2], ee[3,3] = 1,2,2,2 |
---|
188 | for k in range(4,K+1): |
---|
189 | ee[k,0]=1 |
---|
190 | for n in range(1,k+1): |
---|
191 | ee[k,n]=ee[k-2,n]+2*ee[k-1,n-1] |
---|
192 | tau = [ee[K,K-n]/ee[K,K-n-1]/(K-1) for n in range(K)] |
---|
193 | return RK_simple(bwd_fast_slow, dt, tau) |
---|