Coverage for pySDC/implementations/problem_classes/acoustic_helpers/standard_integrators.py: 96%

259 statements  

« prev     ^ index     » next       coverage.py v7.6.7, created at 2024-11-16 14:51 +0000

1import math 

2from decimal import Decimal, getcontext 

3 

4import numpy as np 

5import scipy.sparse as sp 

6import scipy.sparse.linalg as LA 

7 

8 

9# 

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

11# 

12class rk_imex: 

13 def __init__(self, M_fast, M_slow, order): 

14 assert np.shape(M_fast)[0] == np.shape(M_fast)[1], "A_fast must be square" 

15 assert np.shape(M_slow)[0] == np.shape(M_slow)[1], "A_slow must be square" 

16 assert np.shape(M_fast)[0] == np.shape(M_slow)[0], "A_fast and A_slow must be of the same size" 

17 

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

19 self.order = order 

20 

21 if self.order == 2: 

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

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

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

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

26 self.nstages = 2 

27 

28 elif self.order == 3: 

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

30 alpha = 0.24169426078821 

31 beta = 0.06042356519705 

32 eta = 0.12915286960590 

33 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]]) 

34 self.A = np.array( 

35 [ 

36 [alpha, 0, 0, 0], 

37 [-alpha, alpha, 0, 0], 

38 [0, 1.0 - alpha, alpha, 0], 

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

40 ] 

41 ) 

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

43 self.b = self.b_hat 

44 self.nstages = 4 

45 

46 elif self.order == 4: 

47 self.A_hat = np.array( 

48 [ 

49 [0, 0, 0, 0, 0, 0], 

50 [1.0 / 2, 0, 0, 0, 0, 0], 

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

52 [ 

53 -116923316275.0 / 2393684061468.0, 

54 -2731218467317.0 / 15368042101831.0, 

55 9408046702089.0 / 11113171139209.0, 

56 0, 

57 0, 

58 0, 

59 ], 

60 [ 

61 -451086348788.0 / 2902428689909.0, 

62 -2682348792572.0 / 7519795681897.0, 

63 12662868775082.0 / 11960479115383.0, 

64 3355817975965.0 / 11060851509271.0, 

65 0, 

66 0, 

67 ], 

68 [ 

69 647845179188.0 / 3216320057751.0, 

70 73281519250.0 / 8382639484533.0, 

71 552539513391.0 / 3454668386233.0, 

72 3354512671639.0 / 8306763924573.0, 

73 4040.0 / 17871.0, 

74 0, 

75 ], 

76 ] 

77 ) 

78 self.A = np.array( 

79 [ 

80 [0, 0, 0, 0, 0, 0], 

81 [1.0 / 4, 1.0 / 4, 0, 0, 0, 0], 

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

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

84 [ 

85 15267082809.0 / 155376265600.0, 

86 -71443401.0 / 120774400.0, 

87 730878875.0 / 902184768.0, 

88 2285395.0 / 8070912.0, 

89 1.0 / 4, 

90 0, 

91 ], 

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

93 ] 

94 ) 

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

96 self.b_hat = np.array( 

97 [ 

98 4586570599.0 / 29645900160.0, 

99 0, 

100 178811875.0 / 945068544.0, 

101 814220225.0 / 1159782912.0, 

102 -3700637.0 / 11593932.0, 

103 61727.0 / 225920.0, 

104 ] 

105 ) 

106 self.nstages = 6 

107 

108 elif self.order == 5: 

109 # from Kennedy and Carpenter 

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

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

112 getcontext().prec = 56 

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

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

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

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

117 self.A_hat[3, 1] = 0.0 

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

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

120 self.A_hat[4, 1] = 0.0 

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

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

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

124 self.A_hat[5, 1] = 0.0 

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

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

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

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

129 self.A_hat[6, 1] = 0.0 

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

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

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

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

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

135 self.A_hat[7, 1] = 0.0 

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

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

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

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

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

141 

142 self.b_hat = np.zeros(8) 

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

144 self.b_hat[1] = 0.0 

145 self.b_hat[2] = 0.0 

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

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

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

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

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

151 

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

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

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

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

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

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

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

159 self.A[3, 1] = 0.0 

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

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

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

163 self.A[4, 1] = 0.0 

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

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

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

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

168 self.A[5, 1] = 0.0 

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

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

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

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

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

174 self.A[6, 1] = 0.0 

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

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

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

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

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

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

181 self.A[7, 1] = 0.0 

182 self.A[7, 2] = 0.0 

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

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

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

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

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

188 

189 self.b = np.zeros(8) 

190 

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

192 self.b[1] = 0.0 

193 self.b[2] = 0.0 

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

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

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

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

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

199 

200 self.nstages = 8 

201 

202 self.M_fast = sp.csc_matrix(M_fast) 

203 self.M_slow = sp.csc_matrix(M_slow) 

204 self.ndof = np.shape(M_fast)[0] 

205 

206 self.stages = np.zeros((self.nstages, self.ndof), dtype='complex') 

207 

208 def timestep(self, u0, dt): 

209 # Solve for stages 

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

211 # Construct RHS 

212 rhs = np.copy(u0) 

213 for j in range(0, i): 

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

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

216 ) 

217 

218 # Solve for stage i 

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

220 # Avoid call to spsolve with identity matrix 

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

222 else: 

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

224 

225 # Update 

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

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

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

229 ) 

230 

231 return u0 

232 

233 def f_slow(self, u): 

234 return self.M_slow.dot(u) 

235 

236 def f_fast(self, u): 

237 return self.M_fast.dot(u) 

238 

239 def f_fast_solve(self, rhs, alpha): 

240 L = sp.eye(self.ndof) - alpha * self.M_fast 

241 return LA.spsolve(L, rhs) 

242 

243 

244# 

245# Trapezoidal rule 

246# 

247class trapezoidal: 

248 def __init__(self, M, alpha=0.5): 

249 assert np.shape(M)[0] == np.shape(M)[1], "Matrix M must be quadratic" 

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

251 self.M = M 

252 self.alpha = alpha 

253 

254 def timestep(self, u0, dt): 

255 M_trap = sp.eye(self.Ndof) - self.alpha * dt * self.M 

256 B_trap = sp.eye(self.Ndof) + (1.0 - self.alpha) * dt * self.M 

257 b = B_trap.dot(u0) 

258 return LA.spsolve(M_trap, b) 

259 

260 

261# 

262# A BDF-2 implicit two-step method 

263# 

264class bdf2: 

265 def __init__(self, M): 

266 assert np.shape(M)[0] == np.shape(M)[1], "Matrix M must be quadratic" 

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

268 self.M = M 

269 

270 def firsttimestep(self, u0, dt): 

271 b = u0 

272 L = sp.eye(self.Ndof) - dt * self.M 

273 return LA.spsolve(L, b) 

274 

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

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

277 L = sp.eye(self.Ndof) - (2.0 / 3.0) * dt * self.M 

278 return LA.spsolve(L, b) 

279 

280 

281# 

282# A diagonally implicit Runge-Kutta method of order 2, 3 or 4 

283# 

284class dirk: 

285 def __init__(self, M, order): 

286 assert np.shape(M)[0] == np.shape(M)[1], "Matrix M must be quadratic" 

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

288 self.M = sp.csc_matrix(M) 

289 self.order = order 

290 

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

292 

293 if self.order == 2: 

294 self.nstages = 1 

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

296 self.A[0, 0] = 0.5 

297 self.tau = [0.5] 

298 self.b = [1.0] 

299 

300 if self.order == 22: 

301 self.nstages = 2 

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

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

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

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

306 

307 self.tau = np.zeros(2) 

308 self.tau[0] = 1.0 / 3.0 

309 self.tau[1] = 1.0 

310 

311 self.b = np.zeros(2) 

312 self.b[0] = 3.0 / 4.0 

313 self.b[1] = 1.0 / 4.0 

314 

315 if self.order == 3: 

316 self.nstages = 2 

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

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

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

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

321 

322 self.tau = np.zeros(2) 

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

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

325 

326 self.b = np.zeros(2) 

327 self.b[0] = 0.5 

328 self.b[1] = 0.5 

329 

330 if self.order == 4: 

331 self.nstages = 3 

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

333 

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

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

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

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

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

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

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

341 

342 self.tau = np.zeros(3) 

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

344 self.tau[1] = 1.0 / 2.0 

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

346 

347 self.b = np.zeros(3) 

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

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

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

351 

352 if self.order == 5: 

353 self.nstages = 5 

354 # From Kennedy, Carpenter "Diagonally Implicit Runge-Kutta Methods for 

355 # Ordinary Differential Equations. A Review" 

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

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

358 

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

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

361 

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

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

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

365 

366 self.A[3, 0] = -291541413000.0 / 6267936762551.0 

367 self.A[3, 1] = 226761949132.0 / 4473940808273.0 

368 self.A[3, 2] = -1282248297070.0 / 9697416712681.0 

369 self.A[3, 3] = self.A[0, 0] 

370 

371 self.A[4, 0] = -2481679516057.0 / 4626464057815.0 

372 self.A[4, 1] = -197112422687.0 / 6604378783090.0 

373 self.A[4, 2] = 3952887910906.0 / 9713059315593.0 

374 self.A[4, 3] = 4906835613583.0 / 8134926921134.0 

375 self.A[4, 4] = self.A[0, 0] 

376 

377 self.b = np.zeros(5) 

378 self.b[0] = -2522702558582.0 / 12162329469185.0 

379 self.b[1] = 1018267903655.0 / 12907234417901.0 

380 self.b[2] = 4542392826351.0 / 13702606430957.0 

381 self.b[3] = 5001116467727.0 / 12224457745473.0 

382 self.b[4] = 1509636094297.0 / 3891594770934.0 

383 

384 self.stages = np.zeros((self.nstages, self.Ndof), dtype='complex') 

385 

386 def timestep(self, u0, dt): 

387 uend = u0 

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

389 b = u0 

390 

391 # Compute right hand side for this stage's implicit step 

392 for j in range(0, i): 

393 b = b + self.A[i, j] * dt * self.f(self.stages[j, :]) 

394 

395 # Implicit solve for current stage 

396 self.stages[i, :] = self.f_solve(b, dt * self.A[i, i]) 

397 

398 # Add contribution of current stage to final value 

399 uend = uend + self.b[i] * dt * self.f(self.stages[i, :]) 

400 

401 return uend 

402 

403 # 

404 # Returns f(u) = c*u 

405 # 

406 def f(self, u): 

407 return self.M.dot(u) 

408 

409 # 

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

411 # 

412 def f_solve(self, b, alpha): 

413 L = sp.eye(self.Ndof) - alpha * self.M 

414 return LA.spsolve(L, b)