Coverage for pySDC/implementations/sweeper_classes/Runge_Kutta.py: 93%

1import numpy as np 

2import logging 

3from qmat.qcoeff.butcher import RK_SCHEMES 

4 

5from pySDC.core.sweeper import Sweeper, _Pars 

6from pySDC.core.errors import ParameterError, ProblemError 

7from pySDC.core.level import Level 

8 

9 

10class ButcherTableau(object): 

11 def __init__(self, weights, nodes, matrix): 

12 """ 

13 Initialization routine to get a quadrature matrix out of a Butcher tableau 

14 

15 Args: 

16 weights (numpy.ndarray): Butcher tableau weights 

17 nodes (numpy.ndarray): Butcher tableau nodes 

18 matrix (numpy.ndarray): Butcher tableau entries 

19 """ 

20 self.check_method(weights, nodes, matrix) 

21 

22 self.tleft = 0.0 

23 self.tright = 1.0 

24 self.num_nodes = matrix.shape[0] 

25 self.weights = weights 

26 

27 self.nodes = np.append([0], nodes) 

28 self.Qmat = np.zeros([self.num_nodes + 1, self.num_nodes + 1]) 

29 self.Qmat[1:, 1:] = matrix 

30 

31 self.left_is_node = True 

32 self.right_is_node = self.nodes[-1] == self.tright 

33 

34 # compute distances between the nodes 

35 if self.num_nodes > 1: 

36 self.delta_m = self.nodes[1:] - self.nodes[:-1] 

37 else: 

38 self.delta_m = np.zeros(1) 

39 self.delta_m[0] = self.nodes[0] - self.tleft 

40 

41 # check if the RK scheme is implicit 

42 self.implicit = any(matrix[i, i] != 0 for i in range(self.num_nodes)) 

43 

44 def check_method(self, weights, nodes, matrix): 

45 """ 

46 Check that the method is entered in the correct format 

47 """ 

48 if type(matrix) != np.ndarray: 

49 raise ParameterError('Runge-Kutta matrix needs to be supplied as a numpy array!') 

50 elif len(np.unique(matrix.shape)) != 1 or len(matrix.shape) != 2: 

51 raise ParameterError('Runge-Kutta matrix needs to be a square 2D numpy array!') 

52 

53 if type(nodes) != np.ndarray: 

54 raise ParameterError('Nodes need to be supplied as a numpy array!') 

55 elif len(nodes.shape) != 1: 

56 raise ParameterError(f'Incompatible dimension of nodes! Need 1, got {len(nodes.shape)}') 

57 elif len(nodes) != matrix.shape[0]: 

58 raise ParameterError(f'Incompatible number of nodes! Need {matrix.shape[0]}, got {len(nodes)}') 

59 

60 self.check_weights(weights, nodes, matrix) 

61 

62 def check_weights(self, weights, nodes, matrix): 

63 """ 

64 Check that the weights of the method are entered in the correct format 

65 """ 

66 if type(weights) != np.ndarray: 

67 raise ParameterError('Weights need to be supplied as a numpy array!') 

68 elif len(weights.shape) != 1: 

69 raise ParameterError(f'Incompatible dimension of weights! Need 1, got {len(weights.shape)}') 

70 elif len(weights) != matrix.shape[0]: 

71 raise ParameterError(f'Incompatible number of weights! Need {matrix.shape[0]}, got {len(weights)}') 

72 

73 @property 

74 def globally_stiffly_accurate(self): 

75 return np.allclose(self.Qmat[-1, 1:], self.weights) 

76 

77 

78class ButcherTableauEmbedded(ButcherTableau): 

79 

80 def check_weights(self, weights, nodes, matrix): 

81 """ 

82 Check that the weights of the method are entered in the correct format 

83 """ 

84 if type(weights) != np.ndarray: 

85 raise ParameterError('Weights need to be supplied as a numpy array!') 

86 elif len(weights.shape) != 2: 

87 raise ParameterError(f'Incompatible dimension of weights! Need 2, got {len(weights.shape)}') 

88 elif len(weights[0]) != matrix.shape[0]: 

89 raise ParameterError(f'Incompatible number of weights! Need {matrix.shape[0]}, got {len(weights[0])}') 

90 

91 @property 

92 def globally_stiffly_accurate(self): 

93 return np.allclose(self.Qmat[-1, 1:], self.weights[0]) 

94 

95 

96class RungeKutta(Sweeper): 

97 nodes = None 

98 weights = None 

99 matrix = None 

100 ButcherTableauClass = ButcherTableau 

101 

102 """ 

103 Runge-Kutta scheme that fits the interface of a sweeper. 

104 Actually, the sweeper idea fits the Runge-Kutta idea when using only lower triangular rules, where solutions 

105 at the nodes are successively computed from earlier nodes. However, we only perform a single iteration of this. 

106 

107 We have two choices to realise a Runge-Kutta sweeper: We can choose Q = Q_Delta = <Butcher tableau>, but in this 

108 implementation, that would lead to a lot of wasted FLOPS from integrating with Q and then with Q_Delta and 

109 subtracting the two. For that reason, we built this new sweeper, which does not have a preconditioner. 

110 

111 This class only supports lower triangular Butcher tableaux such that the system can be solved with forward 

112 substitution. In this way, we don't get the maximum order that we could for the number of stages, but computing the 

113 stages is much cheaper. In particular, if the Butcher tableaux is strictly lower triangular, we get an explicit 

114 method, which does not require us to solve a system of equations to compute the stages. 

115 

116 Please be aware that all fundamental parameters of the Sweeper are ignored. These include 

117 

118 - num_nodes 

119 - collocation_class 

120 - initial_guess 

121 - QI 

122 

123 All of these variables are either determined by the RK rule, or are not part of an RK scheme. 

124 

125 The entries of the Butcher tableau are stored as class attributes. 

126 """ 

127 

128 def __init__(self, params, level): 

129 """ 

130 Initialization routine for the custom sweeper 

131 

132 Args: 

133 params: parameters for the sweeper 

134 level (pySDC.Level.level): the level that uses this sweeper 

135 """ 

136 # set up logger 

137 self.logger = logging.getLogger('sweeper') 

138 

139 # check if some parameters are set which only apply to actual sweepers 

140 for key in ['initial_guess', 'collocation_class', 'num_nodes']: 

141 if key in params: 

142 self.logger.warning(f'"{key}" will be ignored by Runge-Kutta sweeper') 

143 

144 # set parameters to their actual values 

145 self.coll = self.get_Butcher_tableau() 

146 params['initial_guess'] = 'zero' 

147 params['collocation_class'] = type(self.ButcherTableauClass) 

148 params['num_nodes'] = self.coll.num_nodes 

149 

150 # disable residual computation by default 

151 params['skip_residual_computation'] = params.get( 

152 'skip_residual_computation', ('IT_CHECK', 'IT_FINE', 'IT_COARSE', 'IT_UP', 'IT_DOWN') 

153 ) 

154 

155 # check if we can skip some usually unnecessary right hand side evaluations 

156 params['eval_rhs_at_right_boundary'] = params.get('eval_rhs_at_right_boundary', False) 

157 

158 self.params = _Pars(params) 

159 

160 # set level using the setter in order to adapt residual tolerance if needed 

161 self.__level = None 

162 self.level = level 

163 

164 self.parallelizable = False 

165 self.QI = self.coll.Qmat 

166 

167 @classmethod 

168 def get_Q_matrix(cls): 

169 return cls.get_Butcher_tableau().Qmat 

170 

171 @classmethod 

172 def get_Butcher_tableau(cls): 

173 return cls.ButcherTableauClass(cls.weights, cls.nodes, cls.matrix) 

174 

175 @classmethod 

176 def get_update_order(cls): 

177 """ 

178 Get the order of the lower order method for doing adaptivity. Only applies to embedded methods. 

179 """ 

180 raise NotImplementedError( 

181 f"There is not an update order for RK scheme \"{cls.__name__}\" implemented. Maybe it is not an embedded scheme?" 

182 ) 

183 

184 @classmethod 

185 def is_embedded(cls): 

186 return cls.ButcherTableauClass == ButcherTableauEmbedded 

187 

188 def get_full_f(self, f): 

189 """ 

190 Get the full right hand side as a `mesh` from the right hand side 

191 

192 Args: 

193 f (dtype_f): Right hand side at a single node 

194 

195 Returns: 

196 mesh: Full right hand side as a mesh 

197 """ 

198 if type(f).__name__ in ['mesh', 'cupy_mesh', 'firedrake_mesh']: 

199 return f 

200 elif type(f).__name__.lower() in ['imex_mesh', 'imex_cupy_mesh', 'imex_firedrake_mesh']: 

201 return f.impl + f.expl 

202 elif f is None: 

203 prob = self.level.prob 

204 return self.get_full_f(prob.dtype_f(prob.init, val=0)) 

205 else: 

206 raise NotImplementedError(f'Type \"{type(f)}\" not implemented in Runge-Kutta sweeper') 

207 

208 def integrate(self): 

209 """ 

210 Integrates the right-hand side 

211 

212 Returns: 

213 list of dtype_u: containing the integral as values 

214 """ 

215 

216 # get current level and problem 

217 lvl = self.level 

218 prob = lvl.prob 

219 

220 me = [] 

221 

222 # integrate RHS over all collocation nodes 

223 for m in range(1, self.coll.num_nodes + 1): 

224 # new instance of dtype_u, initialize values with 0 

225 me.append(prob.dtype_u(prob.init, val=0.0)) 

226 for j in range(1, self.coll.num_nodes + 1): 

227 me[-1] += lvl.dt * self.coll.Qmat[m, j] * self.get_full_f(lvl.f[j]) 

228 

229 return me 

230 

231 def update_nodes(self): 

232 """ 

233 Update the u- and f-values at the collocation nodes 

234 

235 Returns: 

236 None 

237 """ 

238 

239 # get current level and problem 

240 lvl = self.level 

241 prob = lvl.prob 

242 

243 # only if the level has been touched before 

244 assert lvl.status.unlocked 

245 assert lvl.status.sweep <= 1, "RK schemes are direct solvers. Please perform only 1 iteration!" 

246 

247 # get number of collocation nodes for easier access 

248 M = self.coll.num_nodes 

249 

250 # `solve_system` of a split problem only inverts the implicit part, so an implicit stage would drop the rest 

251 if self.coll.implicit and prob.dtype_f.__name__.lower().startswith('imex'): 

252 raise ProblemError( 

253 f'{type(self).__name__} has implicit stages, but {type(prob).__name__} splits its right hand side. Use an IMEX Runge-Kutta scheme, such as ARK548L2SA, instead.' 

254 ) 

255 

256 for m in range(0, M): 

257 # build rhs, consisting of the known values from above and new values from previous nodes (at k+1) 

258 rhs = prob.dtype_u(lvl.u[0]) 

259 for j in range(1, m + 1): 

260 rhs += lvl.dt * self.QI[m + 1, j] * self.get_full_f(lvl.f[j]) 

261 

262 # implicit solve with prefactor stemming from the diagonal of Qd, use previous stage as initial guess 

263 if self.QI[m + 1, m + 1] != 0: 

264 lvl.u[m + 1] = prob.solve_system( 

265 rhs, lvl.dt * self.QI[m + 1, m + 1], lvl.u[m], lvl.time + lvl.dt * self.coll.nodes[m + 1] 

266 ) 

267 else: 

268 lvl.u[m + 1] = rhs 

269 

270 # update function values (we don't usually need to evaluate the RHS at the solution of the step) 

271 if m < M - 1 or not self.coll.globally_stiffly_accurate or self.is_embedded(): 

272 lvl.f[m + 1] = prob.eval_f(lvl.u[m + 1], lvl.time + lvl.dt * self.coll.nodes[m + 1]) 

273 else: 

274 lvl.f[m + 1] = prob.f_init 

275 

276 # indicate presence of new values at this level 

277 lvl.status.updated = True 

278 

279 return None 

280 

281 def compute_end_point(self): 

282 """ 

283 In this Runge-Kutta implementation, the solution to the step is always stored in the last node 

284 """ 

285 lvl = self.level 

286 

287 if lvl.f[1] is None: 

288 lvl.uend = lvl.prob.dtype_u(lvl.u[0]) 

289 if self.is_embedded(): 

290 self.u_secondary = lvl.prob.dtype_u(lvl.u[0]) 

291 elif self.coll.globally_stiffly_accurate: 

292 lvl.uend = lvl.prob.dtype_u(lvl.u[-1]) 

293 if self.is_embedded(): 

294 self.u_secondary = lvl.prob.dtype_u(lvl.u[0]) 

295 for w2, k in zip(self.coll.weights[1], lvl.f[1:], strict=True): 

296 self.u_secondary += lvl.dt * w2 * k 

297 else: 

298 lvl.uend = lvl.prob.dtype_u(lvl.u[0]) 

299 if type(self.coll) == ButcherTableau: 

300 for w, k in zip(self.coll.weights, lvl.f[1:], strict=True): 

301 lvl.uend += lvl.dt * w * k 

302 elif self.is_embedded(): 

303 self.u_secondary = lvl.prob.dtype_u(lvl.u[0]) 

304 for w1, w2, k in zip(self.coll.weights[0], self.coll.weights[1], lvl.f[1:], strict=True): 

305 lvl.uend += lvl.dt * w1 * k 

306 self.u_secondary += lvl.dt * w2 * k 

307 

308 @property 

309 def level(self): 

310 """ 

311 Returns the current level 

312 

313 Returns: 

314 pySDC.Level.level: Current level 

315 """ 

316 return self.__level 

317 

318 @level.setter 

319 def level(self, lvl): 

320 """ 

321 Sets a reference to the current level (done in the initialization of the level) 

322 

323 Args: 

324 lvl (pySDC.Level.level): Current level 

325 """ 

326 assert isinstance(lvl, Level), f"You tried to set the sweeper's level with an instance of {type(lvl)}!" 

327 if lvl.params.restol > 0: 

328 lvl.params.restol = -1 

329 self.logger.warning( 

330 'Overwriting residual tolerance with -1 because RK methods are direct and hence may not compute a residual at all!' 

331 ) 

332 

333 self.__level = lvl 

334 

335 def predict(self): 

336 """ 

337 Predictor to fill values at nodes before first sweep 

338 """ 

339 

340 # get current level and problem 

341 lvl = self.level 

342 prob = lvl.prob 

343 

344 for m in range(1, self.coll.num_nodes + 1): 

345 lvl.u[m] = prob.dtype_u(init=prob.init, val=0.0) 

346 

347 # indicate that this level is now ready for sweeps 

348 lvl.status.unlocked = True 

349 lvl.status.updated = True 

350 

351 

352class RungeKuttaIMEX(RungeKutta): 

353 """ 

354 Implicit-explicit split Runge Kutta base class. Only supports methods that share the nodes and weights. 

355 """ 

356 

357 matrix_explicit = None 

358 weights_explicit = None 

359 ButcherTableauClass_explicit = ButcherTableau 

360 

361 def __init__(self, params, level): 

362 """ 

363 Initialization routine 

364 

365 Args: 

366 params: parameters for the sweeper 

367 level (pySDC.Level.level): the level that uses this sweeper 

368 """ 

369 super().__init__(params, level) 

370 type(self).weights_explicit = self.weights if self.weights_explicit is None else self.weights_explicit 

371 self.coll_explicit = self.get_Butcher_tableau_explicit() 

372 self.QE = self.coll_explicit.Qmat 

373 

374 def predict(self): 

375 """ 

376 Predictor to fill values at nodes before first sweep 

377 """ 

378 

379 # get current level and problem 

380 lvl = self.level 

381 prob = lvl.prob 

382 

383 for m in range(1, self.coll.num_nodes + 1): 

384 lvl.u[m] = prob.dtype_u(init=prob.init, val=0.0) 

385 lvl.f[m] = prob.dtype_f(init=prob.init, val=0.0) 

386 

387 # indicate that this level is now ready for sweeps 

388 lvl.status.unlocked = True 

389 lvl.status.updated = True 

390 

391 @classmethod 

392 def get_Butcher_tableau_explicit(cls): 

393 return cls.ButcherTableauClass_explicit(cls.weights_explicit, cls.nodes, cls.matrix_explicit) 

394 

395 def integrate(self): 

396 """ 

397 Integrates the right-hand side 

398 

399 Returns: 

400 list of dtype_u: containing the integral as values 

401 """ 

402 

403 # get current level and problem 

404 lvl = self.level 

405 prob = lvl.prob 

406 

407 me = [] 

408 

409 # integrate RHS over all collocation nodes 

410 for m in range(1, self.coll.num_nodes + 1): 

411 # new instance of dtype_u, initialize values with 0 

412 me.append(prob.dtype_u(prob.init, val=0.0)) 

413 for j in range(1, self.coll.num_nodes + 1): 

414 me[-1] += lvl.dt * ( 

415 self.coll.Qmat[m, j] * lvl.f[j].impl + self.coll_explicit.Qmat[m, j] * lvl.f[j].expl 

416 ) 

417 

418 return me 

419 

420 def update_nodes(self): 

421 """ 

422 Update the u- and f-values at the collocation nodes 

423 

424 Returns: 

425 None 

426 """ 

427 

428 # get current level and problem 

429 lvl = self.level 

430 prob = lvl.prob 

431 

432 # only if the level has been touched before 

433 assert lvl.status.unlocked 

434 assert lvl.status.sweep <= 1, "RK schemes are direct solvers. Please perform only 1 iteration!" 

435 

436 # get number of collocation nodes for easier access 

437 M = self.coll.num_nodes 

438 

439 for m in range(0, M): 

440 # build rhs, consisting of the known values from above and new values from previous nodes (at k+1) 

441 rhs = lvl.prob.dtype_u(lvl.u[0]) 

442 for j in range(1, m + 1): 

443 rhs += lvl.dt * (self.QI[m + 1, j] * lvl.f[j].impl + self.QE[m + 1, j] * lvl.f[j].expl) 

444 

445 # implicit solve with prefactor stemming from the diagonal of Qd, use previous stage as initial guess 

446 if self.QI[m + 1, m + 1] != 0: 

447 lvl.u[m + 1] = prob.solve_system( 

448 rhs, lvl.dt * self.QI[m + 1, m + 1], lvl.u[m], lvl.time + lvl.dt * self.coll.nodes[m + 1] 

449 ) 

450 else: 

451 lvl.u[m + 1] = rhs 

452 

453 # update function values 

454 if ( 

455 m < M - 1 

456 or not (self.coll.globally_stiffly_accurate and self.coll_explicit.globally_stiffly_accurate) 

457 or self.is_embedded() 

458 ): 

459 lvl.f[m + 1] = prob.eval_f(lvl.u[m + 1], lvl.time + lvl.dt * self.coll.nodes[m + 1]) 

460 else: 

461 lvl.f[m + 1] = prob.f_init 

462 

463 # indicate presence of new values at this level 

464 lvl.status.updated = True 

465 

466 return None 

467 

468 def compute_end_point(self): 

469 """ 

470 In this Runge-Kutta implementation, the solution to the step is always stored in the last node 

471 """ 

472 lvl = self.level 

473 

474 if lvl.f[1] is None: 

475 lvl.uend = lvl.prob.dtype_u(lvl.u[0]) 

476 if self.is_embedded(): 

477 self.u_secondary = lvl.prob.dtype_u(lvl.u[0]) 

478 elif self.coll.globally_stiffly_accurate and self.coll_explicit.globally_stiffly_accurate: 

479 lvl.uend = lvl.u[-1] 

480 if self.is_embedded(): 

481 self.u_secondary = lvl.prob.dtype_u(lvl.u[0]) 

482 for w2, w2E, k in zip(self.coll.weights[1], self.coll_explicit.weights[1], lvl.f[1:], strict=True): 

483 self.u_secondary += lvl.dt * (w2 * k.impl + w2E * k.expl) 

484 else: 

485 lvl.uend = lvl.prob.dtype_u(lvl.u[0]) 

486 if type(self.coll) == ButcherTableau: 

487 for w, wE, k in zip(self.coll.weights, self.coll_explicit.weights, lvl.f[1:], strict=True): 

488 lvl.uend += lvl.dt * (w * k.impl + wE * k.expl) 

489 elif self.is_embedded(): 

490 self.u_secondary = lvl.u[0].copy() 

491 for w1, w2, w1E, w2E, k in zip( 

492 self.coll.weights[0], 

493 self.coll.weights[1], 

494 self.coll_explicit.weights[0], 

495 self.coll_explicit.weights[1], 

496 lvl.f[1:], 

497 strict=True, 

498 ): 

499 lvl.uend += lvl.dt * (w1 * k.impl + w1E * k.expl) 

500 self.u_secondary += lvl.dt * (w2 * k.impl + w2E * k.expl) 

501 

502 

503class ForwardEuler(RungeKutta): 

504 """ 

505 Forward Euler. Still a classic. 

506 

507 Not very stable first order method. 

508 """ 

509 

510 generator = RK_SCHEMES["FE"]() 

511 nodes, weights, matrix = generator.genCoeffs() 

512 

513 

514class BackwardEuler(RungeKutta): 

515 """ 

516 Backward Euler. A favorite among true connoisseurs of the heat equation. 

517 

518 A-stable first order method. 

519 """ 

520 

521 generator = RK_SCHEMES["BE"]() 

522 nodes, weights, matrix = generator.genCoeffs() 

523 

524 

525class IMEXEuler(RungeKuttaIMEX): 

526 nodes = BackwardEuler.nodes 

527 weights = BackwardEuler.weights 

528 

529 matrix = BackwardEuler.matrix 

530 matrix_explicit = ForwardEuler.matrix 

531 

532 

533class IMEXEulerStifflyAccurate(RungeKuttaIMEX): 

534 """ 

535 This implements u = fI^-1(u0 + fE(u0)) rather than u = fI^-1(u0) + fE(u0) + u0. 

536 This implementation is slightly inefficient with two stages, but the last stage is the solution, making it stiffly 

537 accurate and suitable for some DAEs. 

538 """ 

539 

540 nodes = np.array([0, 1]) 

541 weights = np.array([0, 1]) 

542 weights_explicit = np.array([1, 0]) 

543 

544 matrix = np.array([[0, 0], [0, 1]]) 

545 matrix_explicit = np.array([[0, 0], [1, 0]]) 

546 

547 

548class CrankNicolson(RungeKutta): 

549 """ 

550 Implicit Runge-Kutta method of second order, A-stable. 

551 """ 

552 

553 generator = RK_SCHEMES["CN"]() 

554 nodes, weights, matrix = generator.genCoeffs() 

555 

556 

557class ExplicitMidpointMethod(RungeKutta): 

558 """ 

559 Explicit Runge-Kutta method of second order. 

560 """ 

561 

562 generator = RK_SCHEMES["RK2"]() 

563 nodes, weights, matrix = generator.genCoeffs() 

564 

565 

566class ImplicitMidpointMethod(RungeKutta): 

567 """ 

568 Implicit Runge-Kutta method of second order. 

569 """ 

570 

571 generator = RK_SCHEMES["IMP"]() 

572 nodes, weights, matrix = generator.genCoeffs() 

573 

574 

575class RK4(RungeKutta): 

576 """ 

577 Explicit Runge-Kutta of fourth order: Everybody's darling. 

578 """ 

579 

580 generator = RK_SCHEMES["RK4"]() 

581 nodes, weights, matrix = generator.genCoeffs() 

582 

583 

584class Heun_Euler(RungeKutta): 

585 """ 

586 Second order explicit embedded Runge-Kutta method. 

587 """ 

588 

589 ButcherTableauClass = ButcherTableauEmbedded 

590 

591 generator = RK_SCHEMES["HEUN"]() 

592 nodes, _weights, matrix = generator.genCoeffs() 

593 weights = np.zeros((2, len(_weights))) 

594 weights[0] = _weights 

595 weights[1] = matrix[-1] 

596 

597 @classmethod 

598 def get_update_order(cls): 

599 return 2 

600 

601 

602class Cash_Karp(RungeKutta): 

603 """ 

604 Fifth order explicit embedded Runge-Kutta. See [here](https://doi.org/10.1145/79505.79507). 

605 """ 

606 

607 generator = RK_SCHEMES["CashKarp"]() 

608 nodes, weights, matrix = generator.genCoeffs(embedded=True) 

609 ButcherTableauClass = ButcherTableauEmbedded 

610 

611 @classmethod 

612 def get_update_order(cls): 

613 return 5 

614 

615 

616class DIRK43(RungeKutta): 

617 """ 

618 Embedded A-stable diagonally implicit RK pair of order 3 and 4. 

619 

620 Taken from [here](https://doi.org/10.1007/BF01934920). 

621 """ 

622 

623 generator = RK_SCHEMES["EDIRK43"]() 

624 nodes, weights, matrix = generator.genCoeffs(embedded=True) 

625 ButcherTableauClass = ButcherTableauEmbedded 

626 

627 @classmethod 

628 def get_update_order(cls): 

629 return 4 

630 

631 

632class DIRK43_2(RungeKutta): 

633 """ 

634 L-stable Diagonally Implicit RK method with four stages of order 3. 

635 Taken from [here](https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods). 

636 """ 

637 

638 generator = RK_SCHEMES["DIRK43"]() 

639 nodes, weights, matrix = generator.genCoeffs() 

640 

641 

642class EDIRK4(RungeKutta): 

643 """ 

644 Stiffly accurate, fourth-order EDIRK with four stages. Taken from 

645 [here](https://ntrs.nasa.gov/citations/20160005923), second one in eq. (216). 

646 """ 

647 

648 generator = RK_SCHEMES["EDIRK4"]() 

649 nodes, weights, matrix = generator.genCoeffs() 

650 

651 

652class ESDIRK53(RungeKutta): 

653 """ 

654 A-stable embedded RK pair of orders 5 and 3, ESDIRK5(3)6L[2]SA. 

655 Taken from [here](https://ntrs.nasa.gov/citations/20160005923) 

656 """ 

657 

658 generator = RK_SCHEMES["ESDIRK53"]() 

659 nodes, weights, matrix = generator.genCoeffs(embedded=True) 

660 ButcherTableauClass = ButcherTableauEmbedded 

661 

662 @classmethod 

663 def get_update_order(cls): 

664 return 4 

665 

666 

667class ESDIRK43(RungeKutta): 

668 """ 

669 A-stable embedded RK pair of orders 4 and 3, ESDIRK4(3)6L[2]SA. 

670 Taken from [here](https://ntrs.nasa.gov/citations/20160005923) 

671 """ 

672 

673 generator = RK_SCHEMES["ESDIRK43"]() 

674 nodes, weights, matrix = generator.genCoeffs(embedded=True) 

675 ButcherTableauClass = ButcherTableauEmbedded 

676 

677 @classmethod 

678 def get_update_order(cls): 

679 return 4 

680 

681 

682class ARK548L2SAERK(RungeKutta): 

683 """ 

684 Explicit part of the ARK54 scheme. 

685 """ 

686 

687 generator = RK_SCHEMES["ARK548L2SAERK"]() 

688 nodes, weights, matrix = generator.genCoeffs(embedded=True) 

689 ButcherTableauClass = ButcherTableauEmbedded 

690 

691 @classmethod 

692 def get_update_order(cls): 

693 return 5 

694 

695 

696class ARK548L2SAESDIRK(ARK548L2SAERK): 

697 """ 

698 Implicit part of the ARK54 scheme. Be careful with the embedded scheme. It seems that both schemes are order 5 as opposed to 5 and 4 as claimed. This may cause issues when doing adaptive time-stepping. 

699 """ 

700 

701 generator_IMP = RK_SCHEMES["ARK548L2SAESDIRK"]() 

702 matrix = generator_IMP.Q 

703 

704 

705class ARK54(RungeKuttaIMEX): 

706 """ 

707 Pair of pairs of ARK5(4)8L[2]SA-ERK and ARK5(4)8L[2]SA-ESDIRK from [here](https://doi.org/10.1016/S0168-9274(02)00138-1). 

708 """ 

709 

710 ButcherTableauClass = ButcherTableauEmbedded 

711 ButcherTableauClass_explicit = ButcherTableauEmbedded 

712 

713 nodes = ARK548L2SAERK.nodes 

714 weights = ARK548L2SAERK.weights 

715 

716 matrix = ARK548L2SAESDIRK.matrix 

717 matrix_explicit = ARK548L2SAERK.matrix 

718 

719 @classmethod 

720 def get_update_order(cls): 

721 return 5 

722 

723 

724class ARK548L2SAESDIRK2(RungeKutta): 

725 """ 

726 Stiffly accurate singly diagonally L-stable implicit embedded Runge-Kutta pair of orders 5 and 4 with explicit first stage from [here](https://doi.org/10.1016/j.apnum.2018.10.007). 

727 This method is part of the IMEX method ARK548L2SA. 

728 """ 

729 

730 generator = RK_SCHEMES["ARK548L2SAESDIRK2"]() 

731 nodes, weights, matrix = generator.genCoeffs(embedded=True) 

732 ButcherTableauClass = ButcherTableauEmbedded 

733 

734 @classmethod 

735 def get_update_order(cls): 

736 return 5 

737 

738 

739class ARK548L2SAERK2(ARK548L2SAESDIRK2): 

740 """ 

741 Explicit embedded pair of Runge-Kutta methods of orders 5 and 4 from [here](https://doi.org/10.1016/j.apnum.2018.10.007). 

742 This method is part of the IMEX method ARK548L2SA. 

743 """ 

744 

745 generator_EXP = RK_SCHEMES["ARK548L2SAERK2"]() 

746 matrix = generator_EXP.Q 

747 

748 

749class ARK548L2SA(RungeKuttaIMEX): 

750 """ 

751 IMEX Runge-Kutta method of order 5 based on the explicit method ARK548L2SAERK2 and the implicit method 

752 ARK548L2SAESDIRK2 from [here](https://doi.org/10.1016/j.apnum.2018.10.007). 

753 

754 According to Kennedy and Carpenter (see reference), the two IMEX RK methods of order 5 are the only ones available 

755 as of now. And we are not aware of higher order ones. This one is newer then the other one and apparently better. 

756 """ 

757 

758 ButcherTableauClass = ButcherTableauEmbedded 

759 ButcherTableauClass_explicit = ButcherTableauEmbedded 

760 

761 nodes = ARK548L2SAERK2.nodes 

762 weights = ARK548L2SAERK2.weights 

763 

764 matrix = ARK548L2SAESDIRK2.matrix 

765 matrix_explicit = ARK548L2SAERK2.matrix 

766 

767 @classmethod 

768 def get_update_order(cls): 

769 return 5 

770 

771 

772class ARK324L2SAERK(RungeKutta): 

773 generator = RK_SCHEMES["ARK324L2SAERK"]() 

774 nodes, weights, matrix = generator.genCoeffs(embedded=True) 

775 ButcherTableauClass = ButcherTableauEmbedded 

776 

777 @classmethod 

778 def get_update_order(cls): 

779 return 3 

780 

781 

782class ARK324L2SAESDIRK(ARK324L2SAERK): 

783 generator = RK_SCHEMES["ARK324L2SAESDIRK"]() 

784 matrix = generator.Q 

785 

786 

787class ARK32(RungeKuttaIMEX): 

788 ButcherTableauClass = ButcherTableauEmbedded 

789 ButcherTableauClass_explicit = ButcherTableauEmbedded 

790 

791 nodes = ARK324L2SAESDIRK.nodes 

792 weights = ARK324L2SAESDIRK.weights 

793 

794 matrix = ARK324L2SAESDIRK.matrix 

795 matrix_explicit = ARK324L2SAERK.matrix 

796 

797 @classmethod 

798 def get_update_order(cls): 

799 return 3 

800 

801 

802class ARK2(RungeKuttaIMEX): 

803 """ 

804 Second order two stage singly diagonally implicit globally stiffly accurate IMEX RK method with explicit first stage. 

805 Can be used to integrate simple DAEs because explicit and implicit part are both stiffly accurate. 

806 """ 

807 

808 generator_IMP = RK_SCHEMES["ARK222EDIRK"]() 

809 generator_EXP = RK_SCHEMES["ARK222ERK"]() 

810 

811 nodes, weights, matrix = generator_IMP.genCoeffs() 

812 _, weights_explicit, matrix_explicit = generator_EXP.genCoeffs() 

813