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

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 

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 def get_full_f(self, f): 

185 """ 

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

187 

188 Args: 

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

190 

191 Returns: 

192 mesh: Full right hand side as a mesh 

193 """ 

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

195 return f 

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

197 return f.impl + f.expl 

198 elif f is None: 

199 prob = self.level.prob 

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

201 else: 

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

203 

204 def integrate(self): 

205 """ 

206 Integrates the right-hand side 

207 

208 Returns: 

209 list of dtype_u: containing the integral as values 

210 """ 

211 

212 # get current level and problem 

213 lvl = self.level 

214 prob = lvl.prob 

215 

216 me = [] 

217 

218 # integrate RHS over all collocation nodes 

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

220 # new instance of dtype_u, initialize values with 0 

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

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

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

224 

225 return me 

226 

227 def update_nodes(self): 

228 """ 

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

230 

231 Returns: 

232 None 

233 """ 

234 

235 # get current level and problem 

236 lvl = self.level 

237 prob = lvl.prob 

238 

239 # only if the level has been touched before 

240 assert lvl.status.unlocked 

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

242 

243 # get number of collocation nodes for easier access 

244 M = self.coll.num_nodes 

245 

246 for m in range(0, M): 

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

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

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

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

251 

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

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

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

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

256 ) 

257 else: 

258 lvl.u[m + 1] = rhs 

259 

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

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

262 

263 # indicate presence of new values at this level 

264 lvl.status.updated = True 

265 

266 return None 

267 

268 def compute_end_point(self): 

269 """ 

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

271 """ 

272 lvl = self.level 

273 

274 if lvl.f[1] is None: 

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

276 if type(self.coll) == ButcherTableauEmbedded: 

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

278 elif self.coll.globally_stiffly_accurate: 

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

280 if type(self.coll) == ButcherTableauEmbedded: 

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

282 for w2, k in zip(self.coll.weights[1], lvl.f[1:]): 

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

284 else: 

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

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

287 for w, k in zip(self.coll.weights, lvl.f[1:]): 

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

289 elif type(self.coll) == ButcherTableauEmbedded: 

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

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

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

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

294 

295 @property 

296 def level(self): 

297 """ 

298 Returns the current level 

299 

300 Returns: 

301 pySDC.Level.level: Current level 

302 """ 

303 return self.__level 

304 

305 @level.setter 

306 def level(self, lvl): 

307 """ 

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

309 

310 Args: 

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

312 """ 

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

314 if lvl.params.restol > 0: 

315 lvl.params.restol = -1 

316 self.logger.warning( 

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

318 ) 

319 

320 self.__level = lvl 

321 

322 def predict(self): 

323 """ 

324 Predictor to fill values at nodes before first sweep 

325 """ 

326 

327 # get current level and problem 

328 lvl = self.level 

329 prob = lvl.prob 

330 

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

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

333 

334 # indicate that this level is now ready for sweeps 

335 lvl.status.unlocked = True 

336 lvl.status.updated = True 

337 

338 

339class RungeKuttaIMEX(RungeKutta): 

340 """ 

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

342 """ 

343 

344 matrix_explicit = None 

345 weights_explicit = None 

346 ButcherTableauClass_explicit = ButcherTableau 

347 

348 def __init__(self, params, level): 

349 """ 

350 Initialization routine 

351 

352 Args: 

353 params: parameters for the sweeper 

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

355 """ 

356 super().__init__(params, level) 

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

358 self.coll_explicit = self.get_Butcher_tableau_explicit() 

359 self.QE = self.coll_explicit.Qmat 

360 

361 def predict(self): 

362 """ 

363 Predictor to fill values at nodes before first sweep 

364 """ 

365 

366 # get current level and problem 

367 lvl = self.level 

368 prob = lvl.prob 

369 

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

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

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

373 

374 # indicate that this level is now ready for sweeps 

375 lvl.status.unlocked = True 

376 lvl.status.updated = True 

377 

378 @classmethod 

379 def get_Butcher_tableau_explicit(cls): 

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

381 

382 def integrate(self): 

383 """ 

384 Integrates the right-hand side 

385 

386 Returns: 

387 list of dtype_u: containing the integral as values 

388 """ 

389 

390 # get current level and problem 

391 lvl = self.level 

392 prob = lvl.prob 

393 

394 me = [] 

395 

396 # integrate RHS over all collocation nodes 

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

398 # new instance of dtype_u, initialize values with 0 

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

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

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

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

403 ) 

404 

405 return me 

406 

407 def update_nodes(self): 

408 """ 

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

410 

411 Returns: 

412 None 

413 """ 

414 

415 # get current level and problem 

416 lvl = self.level 

417 prob = lvl.prob 

418 

419 # only if the level has been touched before 

420 assert lvl.status.unlocked 

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

422 

423 # get number of collocation nodes for easier access 

424 M = self.coll.num_nodes 

425 

426 for m in range(0, M): 

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

428 rhs = lvl.u[0] 

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

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

431 

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

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

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

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

436 ) 

437 else: 

438 lvl.u[m + 1] = rhs[:] 

439 

440 # update function values 

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

442 

443 # indicate presence of new values at this level 

444 lvl.status.updated = True 

445 

446 return None 

447 

448 def compute_end_point(self): 

449 """ 

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

451 """ 

452 lvl = self.level 

453 

454 if lvl.f[1] is None: 

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

456 if type(self.coll) == ButcherTableauEmbedded: 

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

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

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

460 if type(self.coll) == ButcherTableauEmbedded: 

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

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

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

464 else: 

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

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

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

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

469 elif type(self.coll) == ButcherTableauEmbedded: 

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

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

472 self.coll.weights[0], 

473 self.coll.weights[1], 

474 self.coll_explicit.weights[0], 

475 self.coll_explicit.weights[1], 

476 lvl.f[1:], 

477 ): 

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

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

480 

481 

482class ForwardEuler(RungeKutta): 

483 """ 

484 Forward Euler. Still a classic. 

485 

486 Not very stable first order method. 

487 """ 

488 

489 generator = RK_SCHEMES["FE"]() 

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

491 

492 

493class BackwardEuler(RungeKutta): 

494 """ 

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

496 

497 A-stable first order method. 

498 """ 

499 

500 generator = RK_SCHEMES["BE"]() 

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

502 

503 

504class IMEXEuler(RungeKuttaIMEX): 

505 nodes = BackwardEuler.nodes 

506 weights = BackwardEuler.weights 

507 

508 matrix = BackwardEuler.matrix 

509 matrix_explicit = ForwardEuler.matrix 

510 

511 

512class CrankNicolson(RungeKutta): 

513 """ 

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

515 """ 

516 

517 generator = RK_SCHEMES["CN"]() 

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

519 

520 

521class ExplicitMidpointMethod(RungeKutta): 

522 """ 

523 Explicit Runge-Kutta method of second order. 

524 """ 

525 

526 generator = RK_SCHEMES["RK2"]() 

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

528 

529 

530class ImplicitMidpointMethod(RungeKutta): 

531 """ 

532 Implicit Runge-Kutta method of second order. 

533 """ 

534 

535 generator = RK_SCHEMES["IMP"]() 

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

537 

538 

539class RK4(RungeKutta): 

540 """ 

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

542 """ 

543 

544 generator = RK_SCHEMES["RK4"]() 

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

546 

547 

548class Heun_Euler(RungeKutta): 

549 """ 

550 Second order explicit embedded Runge-Kutta method. 

551 """ 

552 

553 ButcherTableauClass = ButcherTableauEmbedded 

554 

555 generator = RK_SCHEMES["HEUN"]() 

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

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

558 weights[0] = _weights 

559 weights[1] = matrix[-1] 

560 

561 @classmethod 

562 def get_update_order(cls): 

563 return 2 

564 

565 

566class Cash_Karp(RungeKutta): 

567 """ 

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

569 """ 

570 

571 generator = RK_SCHEMES["CashKarp"]() 

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

573 ButcherTableauClass = ButcherTableauEmbedded 

574 

575 @classmethod 

576 def get_update_order(cls): 

577 return 5 

578 

579 

580class DIRK43(RungeKutta): 

581 """ 

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

583 

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

585 """ 

586 

587 generator = RK_SCHEMES["EDIRK43"]() 

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

589 ButcherTableauClass = ButcherTableauEmbedded 

590 

591 @classmethod 

592 def get_update_order(cls): 

593 return 4 

594 

595 

596class DIRK43_2(RungeKutta): 

597 """ 

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

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

600 """ 

601 

602 generator = RK_SCHEMES["DIRK43"]() 

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

604 

605 

606class EDIRK4(RungeKutta): 

607 """ 

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

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

610 """ 

611 

612 generator = RK_SCHEMES["EDIRK4"]() 

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

614 

615 

616class ESDIRK53(RungeKutta): 

617 """ 

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

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

620 """ 

621 

622 generator = RK_SCHEMES["ESDIRK53"]() 

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

624 ButcherTableauClass = ButcherTableauEmbedded 

625 

626 @classmethod 

627 def get_update_order(cls): 

628 return 4 

629 

630 

631class ESDIRK43(RungeKutta): 

632 """ 

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

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

635 """ 

636 

637 generator = RK_SCHEMES["ESDIRK43"]() 

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

639 ButcherTableauClass = ButcherTableauEmbedded 

640 

641 @classmethod 

642 def get_update_order(cls): 

643 return 4 

644 

645 

646class ARK548L2SAERK(RungeKutta): 

647 """ 

648 Explicit part of the ARK54 scheme. 

649 """ 

650 

651 generator = RK_SCHEMES["ARK548L2SAERK"]() 

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

653 ButcherTableauClass = ButcherTableauEmbedded 

654 

655 @classmethod 

656 def get_update_order(cls): 

657 return 5 

658 

659 

660class ARK548L2SAESDIRK(ARK548L2SAERK): 

661 """ 

662 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. 

663 """ 

664 

665 generator_IMP = RK_SCHEMES["ARK548L2SAESDIRK"]() 

666 matrix = generator_IMP.Q 

667 

668 

669class ARK54(RungeKuttaIMEX): 

670 """ 

671 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). 

672 """ 

673 

674 ButcherTableauClass = ButcherTableauEmbedded 

675 ButcherTableauClass_explicit = ButcherTableauEmbedded 

676 

677 nodes = ARK548L2SAERK.nodes 

678 weights = ARK548L2SAERK.weights 

679 

680 matrix = ARK548L2SAESDIRK.matrix 

681 matrix_explicit = ARK548L2SAERK.matrix 

682 

683 @classmethod 

684 def get_update_order(cls): 

685 return 5 

686 

687 

688class ARK548L2SAESDIRK2(RungeKutta): 

689 """ 

690 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). 

691 This method is part of the IMEX method ARK548L2SA. 

692 """ 

693 

694 generator = RK_SCHEMES["ARK548L2SAESDIRK2"]() 

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

696 ButcherTableauClass = ButcherTableauEmbedded 

697 

698 @classmethod 

699 def get_update_order(cls): 

700 return 5 

701 

702 

703class ARK548L2SAERK2(ARK548L2SAESDIRK2): 

704 """ 

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

706 This method is part of the IMEX method ARK548L2SA. 

707 """ 

708 

709 generator_EXP = RK_SCHEMES["ARK548L2SAERK2"]() 

710 matrix = generator_EXP.Q 

711 

712 

713class ARK548L2SA(RungeKuttaIMEX): 

714 """ 

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

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

717 

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

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

720 """ 

721 

722 ButcherTableauClass = ButcherTableauEmbedded 

723 ButcherTableauClass_explicit = ButcherTableauEmbedded 

724 

725 nodes = ARK548L2SAERK2.nodes 

726 weights = ARK548L2SAERK2.weights 

727 

728 matrix = ARK548L2SAESDIRK2.matrix 

729 matrix_explicit = ARK548L2SAERK2.matrix 

730 

731 @classmethod 

732 def get_update_order(cls): 

733 return 5 

734 

735 

736class ARK324L2SAERK(RungeKutta): 

737 generator = RK_SCHEMES["ARK324L2SAERK"]() 

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

739 ButcherTableauClass = ButcherTableauEmbedded 

740 

741 @classmethod 

742 def get_update_order(cls): 

743 return 3 

744 

745 

746class ARK324L2SAESDIRK(ARK324L2SAERK): 

747 generator = RK_SCHEMES["ARK324L2SAESDIRK"]() 

748 matrix = generator.Q 

749 

750 

751class ARK32(RungeKuttaIMEX): 

752 ButcherTableauClass = ButcherTableauEmbedded 

753 ButcherTableauClass_explicit = ButcherTableauEmbedded 

754 

755 nodes = ARK324L2SAESDIRK.nodes 

756 weights = ARK324L2SAESDIRK.weights 

757 

758 matrix = ARK324L2SAESDIRK.matrix 

759 matrix_explicit = ARK324L2SAERK.matrix 

760 

761 @classmethod 

762 def get_update_order(cls): 

763 return 3 

764 

765 

766class ARK2(RungeKuttaIMEX): 

767 """ 

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

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

770 """ 

771 

772 generator_IMP = RK_SCHEMES["ARK222EDIRK"]() 

773 generator_EXP = RK_SCHEMES["ARK222ERK"]() 

774 

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

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

777 

778 

779class ARK3(RungeKuttaIMEX): 

780 """ 

781 Third order four stage singly diagonally implicit globally stiffly accurate IMEX RK method with explicit first stage. 

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

783 """ 

784 

785 generator_IMP = RK_SCHEMES["ARK443ESDIRK"]() 

786 generator_EXP = RK_SCHEMES["ARK443ERK"]() 

787 

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

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