Coverage for pySDC/implementations/problem_classes/boussinesq_helpers/standard_integrators.py: 0%

1import math 

2from decimal import Decimal, getcontext 

3 

4import numpy as np 

5import scipy.sparse as sp 

6from scipy.sparse.linalg import gmres 

7 

8from pySDC.implementations.problem_classes.Boussinesq_2D_FD_imex import boussinesq_2d_imex 

9from pySDC.implementations.problem_classes.boussinesq_helpers.helper_classes import logging, Callback 

10 

11 

12# 

13# Runge-Kutta IMEX methods of order 1 to 3 

14# 

15class rk_imex: 

16 def __init__(self, problem, order): 

17 assert order in [1, 2, 3, 4, 5], "Order must be between 1 and 5" 

18 self.order = order 

19 

20 if self.order == 1: 

21 self.A = np.array([[0, 0], [0, 1]]) 

22 self.A_hat = np.array([[0, 0], [1, 0]]) 

23 self.b = np.array([0, 1]) 

24 self.b_hat = np.array([1, 0]) 

25 self.nstages = 2 

26 

27 elif self.order == 2: 

28 self.A = np.array([[0, 0], [0, 0.5]]) 

29 self.A_hat = np.array([[0, 0], [0.5, 0]]) 

30 self.b = np.array([0, 1]) 

31 self.b_hat = np.array([0, 1]) 

32 self.nstages = 2 

33 

34 elif self.order == 3: 

35 # parameter from Pareschi and Russo, J. Sci. Comp. 2005 

36 alpha = 0.24169426078821 

37 beta = 0.06042356519705 

38 eta = 0.12915286960590 

39 self.A_hat = np.array([[0, 0, 0, 0], [0, 0, 0, 0], [0, 1.0, 0, 0], [0, 1.0 / 4.0, 1.0 / 4.0, 0]]) 

40 self.A = np.array( 

41 [ 

42 [alpha, 0, 0, 0], 

43 [-alpha, alpha, 0, 0], 

44 [0, 1.0 - alpha, alpha, 0], 

45 [beta, eta, 0.5 - beta - eta - alpha, alpha], 

46 ] 

47 ) 

48 self.b_hat = np.array([0, 1.0 / 6.0, 1.0 / 6.0, 2.0 / 3.0]) 

49 self.b = self.b_hat 

50 self.nstages = 4 

51 

52 elif self.order == 4: 

53 self.A_hat = np.array( 

54 [ 

55 [0, 0, 0, 0, 0, 0], 

56 [1.0 / 2, 0, 0, 0, 0, 0], 

57 [13861.0 / 62500.0, 6889.0 / 62500.0, 0, 0, 0, 0], 

58 [ 

59 -116923316275.0 / 2393684061468.0, 

60 -2731218467317.0 / 15368042101831.0, 

61 9408046702089.0 / 11113171139209.0, 

62 0, 

63 0, 

64 0, 

65 ], 

66 [ 

67 -451086348788.0 / 2902428689909.0, 

68 -2682348792572.0 / 7519795681897.0, 

69 12662868775082.0 / 11960479115383.0, 

70 3355817975965.0 / 11060851509271.0, 

71 0, 

72 0, 

73 ], 

74 [ 

75 647845179188.0 / 3216320057751.0, 

76 73281519250.0 / 8382639484533.0, 

77 552539513391.0 / 3454668386233.0, 

78 3354512671639.0 / 8306763924573.0, 

79 4040.0 / 17871.0, 

80 0, 

81 ], 

82 ] 

83 ) 

84 self.A = np.array( 

85 [ 

86 [0, 0, 0, 0, 0, 0], 

87 [1.0 / 4, 1.0 / 4, 0, 0, 0, 0], 

88 [8611.0 / 62500.0, -1743.0 / 31250.0, 1.0 / 4, 0, 0, 0], 

89 [5012029.0 / 34652500.0, -654441.0 / 2922500.0, 174375.0 / 388108.0, 1.0 / 4, 0, 0], 

90 [ 

91 15267082809.0 / 155376265600.0, 

92 -71443401.0 / 120774400.0, 

93 730878875.0 / 902184768.0, 

94 2285395.0 / 8070912.0, 

95 1.0 / 4, 

96 0, 

97 ], 

98 [82889.0 / 524892.0, 0, 15625.0 / 83664.0, 69875.0 / 102672.0, -2260.0 / 8211, 1.0 / 4], 

99 ] 

100 ) 

101 self.b = np.array([82889.0 / 524892.0, 0, 15625.0 / 83664.0, 69875.0 / 102672.0, -2260.0 / 8211, 1.0 / 4]) 

102 self.b_hat = np.array( 

103 [ 

104 4586570599.0 / 29645900160.0, 

105 0, 

106 178811875.0 / 945068544.0, 

107 814220225.0 / 1159782912.0, 

108 -3700637.0 / 11593932.0, 

109 61727.0 / 225920.0, 

110 ] 

111 ) 

112 self.nstages = 6 

113 

114 elif self.order == 5: 

115 # from Kennedy and Carpenter 

116 # copied from http://www.mcs.anl.gov/petsc/petsc-3.2/src/ts/impls/arkimex/arkimex.c 

117 self.A_hat = np.zeros((8, 8)) 

118 getcontext().prec = 56 

119 self.A_hat[1, 0] = Decimal(41.0) / Decimal(100.0) 

120 self.A_hat[2, 0] = Decimal(367902744464.0) / Decimal(2072280473677.0) 

121 self.A_hat[2, 1] = Decimal(677623207551.0) / Decimal(8224143866563.0) 

122 self.A_hat[3, 0] = Decimal(1268023523408.0) / Decimal(10340822734521.0) 

123 self.A_hat[3, 1] = 0.0 

124 self.A_hat[3, 2] = Decimal(1029933939417.0) / Decimal(13636558850479.0) 

125 self.A_hat[4, 0] = Decimal(14463281900351.0) / Decimal(6315353703477.0) 

126 self.A_hat[4, 1] = 0.0 

127 self.A_hat[4, 2] = Decimal(66114435211212.0) / Decimal(5879490589093.0) 

128 self.A_hat[4, 3] = Decimal(-54053170152839.0) / Decimal(4284798021562.0) 

129 self.A_hat[5, 0] = Decimal(14090043504691.0) / Decimal(34967701212078.0) 

130 self.A_hat[5, 1] = 0.0 

131 self.A_hat[5, 2] = Decimal(15191511035443.0) / Decimal(11219624916014.0) 

132 self.A_hat[5, 3] = Decimal(-18461159152457.0) / Decimal(12425892160975.0) 

133 self.A_hat[5, 4] = Decimal(-281667163811.0) / Decimal(9011619295870.0) 

134 self.A_hat[6, 0] = Decimal(19230459214898.0) / Decimal(13134317526959.0) 

135 self.A_hat[6, 1] = 0.0 

136 self.A_hat[6, 2] = Decimal(21275331358303.0) / Decimal(2942455364971.0) 

137 self.A_hat[6, 3] = Decimal(-38145345988419.0) / Decimal(4862620318723.0) 

138 self.A_hat[6, 4] = Decimal(-1.0) / Decimal(8.0) 

139 self.A_hat[6, 5] = Decimal(-1.0) / Decimal(8.0) 

140 self.A_hat[7, 0] = Decimal(-19977161125411.0) / Decimal(11928030595625.0) 

141 self.A_hat[7, 1] = 0.0 

142 self.A_hat[7, 2] = Decimal(-40795976796054.0) / Decimal(6384907823539.0) 

143 self.A_hat[7, 3] = Decimal(177454434618887.0) / Decimal(12078138498510.0) 

144 self.A_hat[7, 4] = Decimal(782672205425.0) / Decimal(8267701900261.0) 

145 self.A_hat[7, 5] = Decimal(-69563011059811.0) / Decimal(9646580694205.0) 

146 self.A_hat[7, 6] = Decimal(7356628210526.0) / Decimal(4942186776405.0) 

147 

148 self.b_hat = np.zeros(8) 

149 self.b_hat[0] = Decimal(-872700587467.0) / Decimal(9133579230613.0) 

150 self.b_hat[1] = 0.0 

151 self.b_hat[2] = 0.0 

152 self.b_hat[3] = Decimal(22348218063261.0) / Decimal(9555858737531.0) 

153 self.b_hat[4] = Decimal(-1143369518992.0) / Decimal(8141816002931.0) 

154 self.b_hat[5] = Decimal(-39379526789629.0) / Decimal(19018526304540.0) 

155 self.b_hat[6] = Decimal(32727382324388.0) / Decimal(42900044865799.0) 

156 self.b_hat[7] = Decimal(41.0) / Decimal(200.0) 

157 

158 self.A = np.zeros((8, 8)) 

159 self.A[1, 0] = Decimal(41.0) / Decimal(200.0) 

160 self.A[1, 1] = Decimal(41.0) / Decimal(200.0) 

161 self.A[2, 0] = Decimal(41.0) / Decimal(400.0) 

162 self.A[2, 1] = Decimal(-567603406766.0) / Decimal(11931857230679.0) 

163 self.A[2, 2] = Decimal(41.0) / Decimal(200.0) 

164 self.A[3, 0] = Decimal(683785636431.0) / Decimal(9252920307686.0) 

165 self.A[3, 1] = 0.0 

166 self.A[3, 2] = Decimal(-110385047103.0) / Decimal(1367015193373.0) 

167 self.A[3, 3] = Decimal(41.0) / Decimal(200.0) 

168 self.A[4, 0] = Decimal(3016520224154.0) / Decimal(10081342136671.0) 

169 self.A[4, 1] = 0.0 

170 self.A[4, 2] = Decimal(30586259806659.0) / Decimal(12414158314087.0) 

171 self.A[4, 3] = Decimal(-22760509404356.0) / Decimal(11113319521817.0) 

172 self.A[4, 4] = Decimal(41.0) / Decimal(200.0) 

173 self.A[5, 0] = Decimal(218866479029.0) / Decimal(1489978393911.0) 

174 self.A[5, 1] = 0.0 

175 self.A[5, 2] = Decimal(638256894668.0) / Decimal(5436446318841.0) 

176 self.A[5, 3] = Decimal(-1179710474555.0) / Decimal(5321154724896.0) 

177 self.A[5, 4] = Decimal(-60928119172.0) / Decimal(8023461067671.0) 

178 self.A[5, 5] = Decimal(41.0) / Decimal(200.0) 

179 self.A[6, 0] = Decimal(1020004230633.0) / Decimal(5715676835656.0) 

180 self.A[6, 1] = 0.0 

181 self.A[6, 2] = Decimal(25762820946817.0) / Decimal(25263940353407.0) 

182 self.A[6, 3] = Decimal(-2161375909145.0) / Decimal(9755907335909.0) 

183 self.A[6, 4] = Decimal(-211217309593.0) / Decimal(5846859502534.0) 

184 self.A[6, 5] = Decimal(-4269925059573.0) / Decimal(7827059040749.0) 

185 self.A[6, 6] = Decimal(41.0) / Decimal(200.0) 

186 self.A[7, 0] = Decimal(-872700587467.0) / Decimal(9133579230613.0) 

187 self.A[7, 1] = 0.0 

188 self.A[7, 2] = 0.0 

189 self.A[7, 3] = Decimal(22348218063261.0) / Decimal(9555858737531.0) 

190 self.A[7, 4] = Decimal(-1143369518992.0) / Decimal(8141816002931.0) 

191 self.A[7, 5] = Decimal(-39379526789629.0) / Decimal(19018526304540.0) 

192 self.A[7, 6] = Decimal(32727382324388.0) / Decimal(42900044865799.0) 

193 self.A[7, 7] = Decimal(41.0) / Decimal(200.0) 

194 

195 self.b = np.zeros(8) 

196 

197 self.b[0] = Decimal(-975461918565.0) / Decimal(9796059967033.0) 

198 self.b[1] = 0.0 

199 self.b[2] = 0.0 

200 self.b[3] = Decimal(78070527104295.0) / Decimal(32432590147079.0) 

201 self.b[4] = Decimal(-548382580838.0) / Decimal(3424219808633.0) 

202 self.b[5] = Decimal(-33438840321285.0) / Decimal(15594753105479.0) 

203 self.b[6] = Decimal(3629800801594.0) / Decimal(4656183773603.0) 

204 self.b[7] = Decimal(4035322873751.0) / Decimal(18575991585200.0) 

205 

206 self.nstages = 8 

207 

208 self.problem = problem 

209 self.ndof = np.shape(problem.M)[0] 

210 self.logger = logging() 

211 self.stages = np.zeros((self.nstages, self.ndof)) 

212 

213 def timestep(self, u0, dt): 

214 # Solve for stages 

215 for i in range(0, self.nstages): 

216 # Construct RHS 

217 rhs = np.copy(u0) 

218 for j in range(0, i): 

219 rhs += dt * self.A_hat[i, j] * (self.f_slow(self.stages[j, :])) + dt * self.A[i, j] * ( 

220 self.f_fast(self.stages[j, :]) 

221 ) 

222 

223 # Solve for stage i 

224 if self.A[i, i] == 0: 

225 # Avoid call to spsolve with identity matrix 

226 self.stages[i, :] = np.copy(rhs) 

227 else: 

228 self.stages[i, :] = self.f_fast_solve(rhs, dt * self.A[i, i], u0) 

229 

230 # Update 

231 for i in range(0, self.nstages): 

232 u0 += dt * self.b_hat[i] * (self.f_slow(self.stages[i, :])) + dt * self.b[i] * ( 

233 self.f_fast(self.stages[i, :]) 

234 ) 

235 

236 return u0 

237 

238 def f_slow(self, u): 

239 return self.problem.D_upwind.dot(u) 

240 

241 def f_fast(self, u): 

242 return self.problem.M.dot(u) 

243 

244 def f_fast_solve(self, rhs, alpha, u0): 

245 cb = Callback() 

246 sol, info = gmres( 

247 self.problem.Id - alpha * self.problem.M, 

248 rhs, 

249 x0=u0, 

250 tol=self.problem.params.gmres_tol_limit, 

251 restart=self.problem.params.gmres_restart, 

252 maxiter=self.problem.params.gmres_maxiter, 

253 atol=0, 

254 callback=cb, 

255 ) 

256 if alpha != 0.0: 

257 self.logger.add(cb.getcounter()) 

258 return sol 

259 

260 

261# 

262# Trapezoidal rule 

263# 

264class trapezoidal: 

265 def __init__(self, problem, alpha=0.5): 

266 assert isinstance(problem, boussinesq_2d_imex), "problem is wrong type of object" 

267 self.Ndof = np.shape(problem.M)[0] 

268 self.order = 2 

269 self.logger = logging() 

270 self.problem = problem 

271 self.alpha = alpha 

272 

273 def timestep(self, u0, dt): 

274 B_trap = sp.eye(self.Ndof) + self.alpha * dt * (self.problem.D_upwind + self.problem.M) 

275 b = B_trap.dot(u0) 

276 return self.f_solve(b, alpha=(1.0 - self.alpha) * dt, u0=u0) 

277 

278 # 

279 # Returns f(u) = c*u 

280 # 

281 def f(self, u): 

282 return self.problem.D_upwind.dot(u) + self.problem.M.dot(u) 

283 

284 # 

285 # Solves (Id - alpha*c)*u = b for u 

286 # 

287 def f_solve(self, b, alpha, u0): 

288 cb = Callback() 

289 sol, info = gmres( 

290 self.problem.Id - alpha * (self.problem.D_upwind + self.problem.M), 

291 b, 

292 x0=u0, 

293 tol=self.problem.params.gmres_tol_limit, 

294 restart=self.problem.params.gmres_restart, 

295 maxiter=self.problem.params.gmres_maxiter, 

296 atol=0, 

297 callback=cb, 

298 ) 

299 if alpha != 0.0: 

300 self.logger.add(cb.getcounter()) 

301 return sol 

302 

303 

304# 

305# A BDF-2 implicit two-step method 

306# 

307class bdf2: 

308 def __init__(self, problem): 

309 assert isinstance(problem, boussinesq_2d_imex), "problem is wrong type of object" 

310 self.Ndof = np.shape(problem.M)[0] 

311 self.order = 2 

312 self.logger = logging() 

313 self.problem = problem 

314 

315 def firsttimestep(self, u0, dt): 

316 return self.f_solve(b=u0, alpha=dt, u0=u0) 

317 

318 def timestep(self, u0, um1, dt): 

319 b = (4.0 / 3.0) * u0 - (1.0 / 3.0) * um1 

320 return self.f_solve(b=b, alpha=(2.0 / 3.0) * dt, u0=u0) 

321 

322 # 

323 # Returns f(u) = c*u 

324 # 

325 def f(self, u): 

326 return self.problem.D_upwind.dot(u) + self.problem.M.dot(u) 

327 

328 # 

329 # Solves (Id - alpha*c)*u = b for u 

330 # 

331 def f_solve(self, b, alpha, u0): 

332 cb = Callback() 

333 sol, info = gmres( 

334 self.problem.Id - alpha * (self.problem.D_upwind + self.problem.M), 

335 b, 

336 x0=u0, 

337 tol=self.problem.params.gmres_tol_limit, 

338 restart=self.problem.params.gmres_restart, 

339 maxiter=self.problem.params.gmres_maxiter, 

340 atol=0, 

341 callback=cb, 

342 ) 

343 if alpha != 0.0: 

344 self.logger.add(cb.getcounter()) 

345 return sol 

346 

347 

348# 

349# Split-Explicit method 

350# 

351 

352 

353class SplitExplicit: 

354 def __init__(self, problem, method, pparams): 

355 assert isinstance(problem, boussinesq_2d_imex), "problem is wrong type of object" 

356 self.Ndof = np.shape(problem.M)[0] 

357 self.method = method 

358 self.logger = logging() 

359 self.problem = problem 

360 self.pparams = pparams 

361 self.NdofTher = 2 * problem.N[0] * problem.N[1] 

362 self.NdofMom = 2 * problem.N[0] * problem.N[1] 

363 

364 self.ns = None 

365 

366 # print("dx ",problem.h[0]) 

367 # print("dz ",problem.h[1]) 

368 

369 assert self.method in ["MIS4_4", "RK3"], 'Method must be MIS4_4' 

370 

371 if self.method == 'RK3': 

372 self.nstages = 3 

373 self.aRunge = np.zeros((4, 4)) 

374 self.aRunge[0, 0] = 1.0 / 3.0 

375 self.aRunge[1, 1] = 1.0 / 2.0 

376 self.aRunge[2, 2] = 1.0 

377 self.dRunge = np.zeros((4, 4)) 

378 self.gRunge = np.zeros((4, 4)) 

379 if self.method == 'MIS4_4': 

380 self.nstages = 4 

381 self.aRunge = np.zeros((4, 4)) 

382 self.aRunge[0, 0] = 0.38758444641450318 

383 self.aRunge[1, 0] = -2.5318448354142823e-002 

384 self.aRunge[1, 1] = 0.38668943087310403 

385 self.aRunge[2, 0] = 0.20899983523553325 

386 self.aRunge[2, 1] = -0.45856648476371231 

387 self.aRunge[2, 2] = 0.43423187573425748 

388 self.aRunge[3, 0] = -0.10048822195663100 

389 self.aRunge[3, 1] = -0.46186171956333327 

390 self.aRunge[3, 2] = 0.83045062122462809 

391 self.aRunge[3, 3] = 0.27014914900250392 

392 self.dRunge = np.zeros((4, 4)) 

393 self.dRunge[1, 1] = 0.52349249922385610 

394 self.dRunge[2, 1] = 1.1683374366893629 

395 self.dRunge[2, 2] = -0.75762080241712637 

396 self.dRunge[3, 1] = -3.6477233846797109e-002 

397 self.dRunge[3, 2] = 0.56936148730740477 

398 self.dRunge[3, 3] = 0.47746263002599681 

399 self.gRunge = np.zeros((4, 4)) 

400 self.gRunge[1, 1] = 0.13145089796226542 

401 self.gRunge[2, 1] = -0.36855857648747881 

402 self.gRunge[2, 2] = 0.33159232636600550 

403 self.gRunge[3, 1] = -6.5767130537473045e-002 

404 self.gRunge[3, 2] = 4.0591093109036858e-002 

405 self.gRunge[3, 3] = 6.4902111640806712e-002 

406 self.dtRunge = np.zeros(self.nstages) 

407 for i in range(0, self.nstages): 

408 self.dtRunge[i] = 0 

409 temp = 1.0 

410 for j in range(0, i + 1): 

411 self.dtRunge[i] = self.dtRunge[i] + self.aRunge[i, j] 

412 temp = temp - self.dRunge[i, j] 

413 self.dRunge[i, 0] = temp 

414 for j in range(0, i + 1): 

415 self.aRunge[i, j] = self.aRunge[i, j] / self.dtRunge[i] 

416 self.gRunge[i, j] = self.gRunge[i, j] / self.dtRunge[i] 

417 

418 self.U = np.zeros((self.Ndof, self.nstages + 1)) 

419 self.F = np.zeros((self.Ndof, self.nstages)) 

420 self.FSlow = np.zeros(self.Ndof) 

421 self.nsMin = 8 

422 self.logger.nsmall = 0 

423 

424 def NumSmallTimeSteps(self, dx, dz, dt): 

425 cs = self.pparams['c_s'] 

426 ns = dt / (0.9 / np.sqrt(1 / (dx * dx) + 1 / (dz * dz)) / cs) 

427 ns = max(np.int(np.ceil(ns)), self.nsMin) 

428 return ns 

429 

430 def timestep(self, u0, dt): 

431 self.U[:, 0] = u0 

432 

433 self.ns = self.NumSmallTimeSteps(self.problem.h[0], self.problem.h[1], dt) 

434 

435 for i in range(0, self.nstages): 

436 self.F[:, i] = self.f_slow(self.U[:, i]) 

437 self.FSlow[:] = 0.0 

438 for j in range(0, i + 1): 

439 self.FSlow += self.aRunge[i, j] * self.F[:, j] + self.gRunge[i, j] / dt * (self.U[:, j] - u0) 

440 self.U[:, i + 1] = 0 

441 for j in range(0, i + 1): 

442 self.U[:, i + 1] += self.dRunge[i, j] * self.U[:, j] 

443 nsLoc = np.int(np.ceil(self.ns * self.dtRunge[i])) 

444 self.logger.nsmall += nsLoc 

445 dtLoc = dt * self.dtRunge[i] 

446 dTau = dtLoc / nsLoc 

447 self.U[:, i + 1] = self.VerletLin(self.U[:, i + 1], self.FSlow, nsLoc, dTau) 

448 u0 = self.U[:, self.nstages] 

449 return u0 

450 

451 def VerletLin(self, u0, FSlow, ns, dTau): 

452 for _ in range(0, ns): 

453 u0[0 : self.NdofMom] += dTau * (self.f_fastMom(u0) + FSlow[0 : self.NdofMom]) 

454 u0[self.NdofMom : self.Ndof] += dTau * (self.f_fastTher(u0) + FSlow[self.NdofMom : self.Ndof]) 

455 

456 return u0 

457 

458 def RK3Lin(self, u0, FSlow, ns, dTau): 

459 u = u0 

460 for _ in range(0, ns): 

461 u = u0 + dTau / 3.0 * (self.f_fast(u) + FSlow) 

462 u = u0 + dTau / 2.0 * (self.f_fast(u) + FSlow) 

463 u = u0 + dTau * (self.f_fast(u) + FSlow) 

464 u0 = u 

465 

466 return u0 

467 

468 def f_slow(self, u): 

469 return self.problem.D_upwind.dot(u) 

470 

471 def f_fast(self, u): 

472 return self.problem.M.dot(u) 

473 

474 def f_fastMom(self, u): 

475 return self.problem.M[0 : self.NdofMom, self.NdofMom : self.Ndof].dot(u[self.NdofMom : self.Ndof]) 

476 

477 def f_fastTher(self, u): 

478 return self.problem.M[self.NdofMom : self.Ndof, 0 : self.NdofMom].dot(u[0 : self.NdofMom]) 

479 

480 

481class dirk: 

482 def __init__(self, problem, order): 

483 assert isinstance(problem, boussinesq_2d_imex), "problem is wrong type of object" 

484 self.Ndof = np.shape(problem.M)[0] 

485 self.order = order 

486 self.logger = logging() 

487 self.problem = problem 

488 

489 assert self.order in [2, 22, 3, 4, 5], 'Order must be 2,22,3,4' 

490 

491 if self.order == 2: 

492 self.nstages = 1 

493 self.A = np.zeros((1, 1)) 

494 self.A[0, 0] = 0.5 

495 self.tau = [0.5] 

496 self.b = [1.0] 

497 

498 if self.order == 22: 

499 self.nstages = 2 

500 self.A = np.zeros((2, 2)) 

501 self.A[0, 0] = 1.0 / 3.0 

502 self.A[1, 0] = 1.0 / 2.0 

503 self.A[1, 1] = 1.0 / 2.0 

504 

505 self.tau = np.zeros(2) 

506 self.tau[0] = 1.0 / 3.0 

507 self.tau[1] = 1.0 

508 

509 self.b = np.zeros(2) 

510 self.b[0] = 3.0 / 4.0 

511 self.b[1] = 1.0 / 4.0 

512 

513 if self.order == 3: 

514 self.nstages = 2 

515 self.A = np.zeros((2, 2)) 

516 self.A[0, 0] = 0.5 + 1.0 / (2.0 * math.sqrt(3.0)) 

517 self.A[1, 0] = -1.0 / math.sqrt(3.0) 

518 self.A[1, 1] = self.A[0, 0] 

519 

520 self.tau = np.zeros(2) 

521 self.tau[0] = 0.5 + 1.0 / (2.0 * math.sqrt(3.0)) 

522 self.tau[1] = 0.5 - 1.0 / (2.0 * math.sqrt(3.0)) 

523 

524 self.b = np.zeros(2) 

525 self.b[0] = 0.5 

526 self.b[1] = 0.5 

527 

528 if self.order == 4: 

529 self.nstages = 3 

530 alpha = 2.0 * math.cos(math.pi / 18.0) / math.sqrt(3.0) 

531 

532 self.A = np.zeros((3, 3)) 

533 self.A[0, 0] = (1.0 + alpha) / 2.0 

534 self.A[1, 0] = -alpha / 2.0 

535 self.A[1, 1] = self.A[0, 0] 

536 self.A[2, 0] = 1.0 + alpha 

537 self.A[2, 1] = -(1.0 + 2.0 * alpha) 

538 self.A[2, 2] = self.A[0, 0] 

539 

540 self.tau = np.zeros(3) 

541 self.tau[0] = (1.0 + alpha) / 2.0 

542 self.tau[1] = 1.0 / 2.0 

543 self.tau[2] = (1.0 - alpha) / 2.0 

544 

545 self.b = np.zeros(3) 

546 self.b[0] = 1.0 / (6.0 * alpha * alpha) 

547 self.b[1] = 1.0 - 1.0 / (3.0 * alpha * alpha) 

548 self.b[2] = 1.0 / (6.0 * alpha * alpha) 

549 

550 if self.order == 5: 

551 self.nstages = 5 

552 # From Kennedy, Carpenter "Diagonally Implicit Runge-Kutta Methods for Ordinary Differential Equations. 

553 # A Review" 

554 self.A = np.zeros((5, 5)) 

555 self.A[0, 0] = 4024571134387.0 / 14474071345096.0 

556 

557 self.A[1, 0] = 9365021263232.0 / 12572342979331.0 

558 self.A[1, 1] = self.A[0, 0] 

559 

560 self.A[2, 0] = 2144716224527.0 / 9320917548702.0 

561 self.A[2, 1] = -397905335951.0 / 4008788611757.0