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 """ 

76 Whether the last row of the Butcher matrix equals the weights, such that the last stage is the solution of the 

77 step. 

78 

79 Returns: 

80 bool: True if the method is stiffly accurate 

81 """ 

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

83 

84 

85class ButcherTableauEmbedded(ButcherTableau): 

86 

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

88 """ 

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

90 """ 

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

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

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

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

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

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

97 

98 @property 

99 def globally_stiffly_accurate(self): 

100 """ 

101 Whether the last row of the Butcher matrix equals the weights of the primary (higher order) method, such that 

102 the last stage is the solution of the step. 

103 

104 Returns: 

105 bool: True if the method is stiffly accurate 

106 """ 

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

108 

109 

110class RungeKutta(Sweeper): 

111 """ 

112 Base class for Runge-Kutta methods with lower triangular Butcher tableaux, wrapped in the sweeper interface. 

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

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

115 

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

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

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

119 

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

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

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

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

124 

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

126 

127 - num_nodes 

128 - collocation_class 

129 - initial_guess 

130 - QI 

131 

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

133 

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

135 """ 

136 

137 nodes = None 

138 weights = None 

139 matrix = None 

140 ButcherTableauClass = ButcherTableau 

141 

142 def __init__(self, params, level): 

143 """ 

144 Initialization routine for the custom sweeper 

145 

146 Args: 

147 params: parameters for the sweeper 

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

149 """ 

150 # set up logger 

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

152 

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

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

155 if key in params: 

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

157 

158 # set parameters to their actual values 

159 self.coll = self.get_Butcher_tableau() 

160 params['initial_guess'] = 'zero' 

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

162 

163 # disable residual computation by default 

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

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

166 ) 

167 

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

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

170 

171 self.params = _Pars(params) 

172 

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

174 self.__level = None 

175 self.level = level 

176 

177 self.parallelizable = False 

178 self.QI = self.coll.Qmat 

179 

180 @classmethod 

181 def get_Q_matrix(cls): 

182 """ 

183 Get the quadrature matrix of the scheme, i.e. the Butcher matrix padded with a leading row and column of zeros 

184 for the initial conditions. 

185 

186 Returns: 

187 numpy.ndarray: the quadrature matrix 

188 """ 

189 return cls.get_Butcher_tableau().Qmat 

190 

191 @classmethod 

192 def get_Butcher_tableau(cls): 

193 """ 

194 Build the Butcher tableau from the class attributes ``weights``, ``nodes`` and ``matrix``. 

195 

196 Returns: 

197 ButcherTableau: an instance of ``ButcherTableauClass`` 

198 """ 

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

200 

201 @classmethod 

202 def get_update_order(cls): 

203 """ 

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

205 """ 

206 raise NotImplementedError( 

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

208 ) 

209 

210 @classmethod 

211 def is_embedded(cls): 

212 """ 

213 Whether the scheme is an embedded pair, i.e. uses ``ButcherTableauEmbedded`` and has a secondary solution for 

214 error estimation. 

215 

216 Returns: 

217 bool: True if the scheme is embedded 

218 """ 

219 return cls.ButcherTableauClass == ButcherTableauEmbedded 

220 

221 def get_full_f(self, f): 

222 """ 

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

224 

225 Args: 

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

227 

228 Returns: 

229 mesh: Full right hand side as a mesh 

230 """ 

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

232 return f 

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

234 return f.impl + f.expl 

235 elif f is None: 

236 prob = self.level.prob 

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

238 else: 

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

240 

241 def integrate(self): 

242 """ 

243 Integrates the right-hand side 

244 

245 Returns: 

246 list of dtype_u: containing the integral as values 

247 """ 

248 

249 # get current level and problem 

250 lvl = self.level 

251 prob = lvl.prob 

252 

253 me = [] 

254 

255 # integrate RHS over all collocation nodes 

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

257 # new instance of dtype_u, initialize values with 0 

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

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

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

261 

262 return me 

263 

264 def update_nodes(self): 

265 """ 

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

267 

268 Returns: 

269 None 

270 """ 

271 

272 # get current level and problem 

273 lvl = self.level 

274 prob = lvl.prob 

275 

276 # only if the level has been touched before 

277 assert lvl.status.unlocked 

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

279 

280 # get number of collocation nodes for easier access 

281 M = self.coll.num_nodes 

282 

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

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

285 raise ProblemError( 

286 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.' 

287 ) 

288 

289 for m in range(0, M): 

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

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

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

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

294 

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

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

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

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

299 ) 

300 else: 

301 lvl.u[m + 1] = rhs 

302 

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

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

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

306 else: 

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

308 

309 # indicate presence of new values at this level 

310 lvl.status.updated = True 

311 

312 return None 

313 

314 def compute_end_point(self): 

315 """ 

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

317 """ 

318 lvl = self.level 

319 

320 if lvl.f[1] is None: 

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

322 if self.is_embedded(): 

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

324 elif self.coll.globally_stiffly_accurate: 

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

326 if self.is_embedded(): 

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

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

329 self.u_secondary += lvl.dt * w2 * self.get_full_f(k) 

330 else: 

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

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

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

334 lvl.uend += lvl.dt * w * self.get_full_f(k) 

335 elif self.is_embedded(): 

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

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

338 k = self.get_full_f(k) 

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

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

341 

342 @property 

343 def level(self): 

344 """ 

345 Returns the current level 

346 

347 Returns: 

348 pySDC.Level.level: Current level 

349 """ 

350 return self.__level 

351 

352 @level.setter 

353 def level(self, lvl): 

354 """ 

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

356 

357 Args: 

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

359 """ 

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

361 if lvl.params.restol > 0: 

362 lvl.params.restol = -1 

363 self.logger.warning( 

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

365 ) 

366 

367 self.__level = lvl 

368 

369 def predict(self): 

370 """ 

371 Predictor to fill values at nodes before first sweep 

372 """ 

373 

374 # get current level and problem 

375 lvl = self.level 

376 prob = lvl.prob 

377 

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

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

380 

381 # indicate that this level is now ready for sweeps 

382 lvl.status.unlocked = True 

383 lvl.status.updated = True 

384 

385 

386class RungeKuttaIMEX(RungeKutta): 

387 """ 

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

389 """ 

390 

391 matrix_explicit = None 

392 weights_explicit = None 

393 ButcherTableauClass_explicit = ButcherTableau 

394 

395 def __init__(self, params, level): 

396 """ 

397 Initialization routine 

398 

399 Args: 

400 params: parameters for the sweeper 

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

402 """ 

403 super().__init__(params, level) 

404 self.coll_explicit = self.get_Butcher_tableau_explicit() 

405 self.QE = self.coll_explicit.Qmat 

406 

407 def predict(self): 

408 """ 

409 Predictor to fill values at nodes before first sweep 

410 """ 

411 

412 # get current level and problem 

413 lvl = self.level 

414 prob = lvl.prob 

415 

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

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

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

419 

420 # indicate that this level is now ready for sweeps 

421 lvl.status.unlocked = True 

422 lvl.status.updated = True 

423 

424 @classmethod 

425 def get_Butcher_tableau_explicit(cls): 

426 """ 

427 Build the Butcher tableau of the explicit part from ``weights_explicit``, ``nodes`` and ``matrix_explicit``. 

428 Without ``weights_explicit``, the explicit part uses the implicit ``weights``. 

429 

430 Returns: 

431 ButcherTableau: an instance of ``ButcherTableauClass_explicit`` 

432 """ 

433 weights = cls.weights if cls.weights_explicit is None else cls.weights_explicit 

434 return cls.ButcherTableauClass_explicit(weights, cls.nodes, cls.matrix_explicit) 

435 

436 def integrate(self): 

437 """ 

438 Integrates the right-hand side 

439 

440 Returns: 

441 list of dtype_u: containing the integral as values 

442 """ 

443 

444 # get current level and problem 

445 lvl = self.level 

446 prob = lvl.prob 

447 

448 me = [] 

449 

450 # integrate RHS over all collocation nodes 

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

452 # new instance of dtype_u, initialize values with 0 

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

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

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

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

457 ) 

458 

459 return me 

460 

461 def update_nodes(self): 

462 """ 

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

464 

465 Returns: 

466 None 

467 """ 

468 

469 # get current level and problem 

470 lvl = self.level 

471 prob = lvl.prob 

472 

473 # only if the level has been touched before 

474 assert lvl.status.unlocked 

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

476 

477 # get number of collocation nodes for easier access 

478 M = self.coll.num_nodes 

479 

480 for m in range(0, M): 

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

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

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

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

485 

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

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

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

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

490 ) 

491 else: 

492 lvl.u[m + 1] = rhs 

493 

494 # update function values 

495 if ( 

496 m < M - 1 

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

498 or self.is_embedded() 

499 ): 

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

501 else: 

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

503 

504 # indicate presence of new values at this level 

505 lvl.status.updated = True 

506 

507 return None 

508 

509 def compute_end_point(self): 

510 """ 

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

512 """ 

513 lvl = self.level 

514 

515 if lvl.f[1] is None: 

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

517 if self.is_embedded(): 

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

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

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

521 if self.is_embedded(): 

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

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

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

525 else: 

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

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

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

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

530 elif self.is_embedded(): 

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

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

533 self.coll.weights[0], 

534 self.coll.weights[1], 

535 self.coll_explicit.weights[0], 

536 self.coll_explicit.weights[1], 

537 lvl.f[1:], 

538 strict=True, 

539 ): 

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

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

542 

543 

544class ForwardEuler(RungeKutta): 

545 """ 

546 Forward Euler. Still a classic. 

547 

548 Not very stable first order method. 

549 """ 

550 

551 generator = RK_SCHEMES["FE"]() 

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

553 

554 

555class BackwardEuler(RungeKutta): 

556 """ 

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

558 

559 A-stable first order method. 

560 """ 

561 

562 generator = RK_SCHEMES["BE"]() 

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

564 

565 

566class IMEXEuler(RungeKuttaIMEX): 

567 """ 

568 First-order IMEX Euler with one backward Euler stage, at which both parts of the right hand side are evaluated. 

569 

570 This is a Lie splitting of a backward Euler step for fI and a forward Euler step for fE from its result: 

571 u1 = u0 + dt fI(u1), u = u1 + dt fE(u1, t + dt). The explicit tableau has node 1 and matrix entry 0, so fE is 

572 evaluated at the new time and the implicit stage, not at the start of the step. For the variant with fE at the 

573 start of the step, see IMEXEulerStifflyAccurate. 

574 """ 

575 

576 nodes = BackwardEuler.nodes 

577 weights = BackwardEuler.weights 

578 

579 matrix = BackwardEuler.matrix 

580 matrix_explicit = ForwardEuler.matrix 

581 

582 

583class IMEXEulerStifflyAccurate(RungeKuttaIMEX): 

584 """ 

585 IMEX Euler with the explicit part evaluated at the start of the step, as a stiffly accurate two-stage method. 

586 This implements u = fI^-1(u0 + fE(u0)) rather than IMEXEuler's u = fI^-1(u0) + fE(fI^-1(u0)). 

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

588 accurate and suitable for some DAEs. 

589 """ 

590 

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

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

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

594 

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

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

597 

598 

599class CrankNicolson(RungeKutta): 

600 """ 

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

602 """ 

603 

604 generator = RK_SCHEMES["CN"]() 

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

606 

607 

608class ExplicitMidpointMethod(RungeKutta): 

609 """ 

610 Explicit Runge-Kutta method of second order. 

611 """ 

612 

613 generator = RK_SCHEMES["RK2"]() 

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

615 

616 

617class ImplicitMidpointMethod(RungeKutta): 

618 """ 

619 Implicit Runge-Kutta method of second order. 

620 """ 

621 

622 generator = RK_SCHEMES["IMP"]() 

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

624 

625 

626class RK4(RungeKutta): 

627 """ 

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

629 """ 

630 

631 generator = RK_SCHEMES["RK4"]() 

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

633 

634 

635class Heun_Euler(RungeKutta): 

636 """ 

637 Second order explicit embedded Runge-Kutta method. 

638 """ 

639 

640 ButcherTableauClass = ButcherTableauEmbedded 

641 

642 generator = RK_SCHEMES["HEUN"]() 

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

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

645 weights[0] = _weights 

646 weights[1] = matrix[-1] 

647 

648 @classmethod 

649 def get_update_order(cls): 

650 """ 

651 Order in dt of the embedded error estimate, i.e. of the local error of forward Euler. 

652 

653 Returns: 

654 int: 2 

655 """ 

656 return 2 

657 

658 

659class Cash_Karp(RungeKutta): 

660 """ 

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

662 """ 

663 

664 generator = RK_SCHEMES["CashKarp"]() 

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

666 ButcherTableauClass = ButcherTableauEmbedded 

667 

668 @classmethod 

669 def get_update_order(cls): 

670 """ 

671 Order in dt of the embedded error estimate, i.e. of the local error of the fourth order secondary method. 

672 

673 Returns: 

674 int: 5 

675 """ 

676 return 5 

677 

678 

679class DIRK43(RungeKutta): 

680 """ 

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

682 

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

684 """ 

685 

686 generator = RK_SCHEMES["EDIRK43"]() 

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

688 ButcherTableauClass = ButcherTableauEmbedded 

689 

690 @classmethod 

691 def get_update_order(cls): 

692 """ 

693 Order in dt of the embedded error estimate, i.e. of the local error of the third order secondary method. 

694 

695 Returns: 

696 int: 4 

697 """ 

698 return 4 

699 

700 

701class DIRK43_2(RungeKutta): 

702 """ 

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

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

705 """ 

706 

707 generator = RK_SCHEMES["DIRK43"]() 

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

709 

710 

711class EDIRK4(RungeKutta): 

712 """ 

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

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

715 """ 

716 

717 generator = RK_SCHEMES["EDIRK4"]() 

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

719 

720 

721class ESDIRK53(RungeKutta): 

722 """ 

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

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

725 """ 

726 

727 generator = RK_SCHEMES["ESDIRK53"]() 

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

729 ButcherTableauClass = ButcherTableauEmbedded 

730 

731 @classmethod 

732 def get_update_order(cls): 

733 """ 

734 Order in dt of the embedded error estimate, i.e. of the local error of the third order secondary method. 

735 

736 Returns: 

737 int: 4 

738 """ 

739 return 4 

740 

741 

742class ESDIRK43(RungeKutta): 

743 """ 

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

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

746 """ 

747 

748 generator = RK_SCHEMES["ESDIRK43"]() 

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

750 ButcherTableauClass = ButcherTableauEmbedded 

751 

752 @classmethod 

753 def get_update_order(cls): 

754 """ 

755 Order in dt of the embedded error estimate, i.e. of the local error of the third order secondary method. 

756 

757 Returns: 

758 int: 4 

759 """ 

760 return 4 

761 

762 

763class ARK548L2SAERK(RungeKutta): 

764 """ 

765 Explicit part of the ARK54 scheme. 

766 """ 

767 

768 generator = RK_SCHEMES["ARK548L2SAERK"]() 

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

770 ButcherTableauClass = ButcherTableauEmbedded 

771 

772 @classmethod 

773 def get_update_order(cls): 

774 """ 

775 Order in dt of the embedded error estimate, i.e. of the local error of the fourth order secondary method. 

776 

777 Returns: 

778 int: 5 

779 """ 

780 return 5 

781 

782 

783class ARK548L2SAESDIRK(ARK548L2SAERK): 

784 """ 

785 Implicit part of the ARK54 scheme. As Kennedy and Carpenter state, the embedded method has order 4 and the main 

786 method order 5. The fifth-order conditions of the embedded method are violated only by about 1e-4, however, so it 

787 can look like order 5 in convergence tests at moderate step sizes. 

788 """ 

789 

790 generator_IMP = RK_SCHEMES["ARK548L2SAESDIRK"]() 

791 matrix = generator_IMP.Q 

792 

793 

794class ARK54(RungeKuttaIMEX): 

795 """ 

796 IMEX Runge-Kutta method ARK5(4)8L[2]SA of Kennedy and Carpenter, combining ARK548L2SAERK and ARK548L2SAESDIRK. 

797 Pair of pairs of ARK5(4)8L[2]SA-ERK and ARK5(4)8L[2]SA-ESDIRK from 

798 [here](https://doi.org/10.1016/S0168-9274(02)00138-1). 

799 """ 

800 

801 ButcherTableauClass = ButcherTableauEmbedded 

802 ButcherTableauClass_explicit = ButcherTableauEmbedded 

803 

804 nodes = ARK548L2SAERK.nodes 

805 weights = ARK548L2SAERK.weights 

806 

807 matrix = ARK548L2SAESDIRK.matrix 

808 matrix_explicit = ARK548L2SAERK.matrix 

809 

810 @classmethod 

811 def get_update_order(cls): 

812 """ 

813 Order in dt of the embedded error estimate, i.e. of the local error of the fourth order secondary method. 

814 

815 Returns: 

816 int: 5 

817 """ 

818 return 5 

819 

820 

821class ARK548L2SAESDIRK2(RungeKutta): 

822 """ 

823 Implicit part of ARK548L2SA: an L-stable, stiffly accurate ESDIRK pair of orders 5 and 4. 

824 Stiffly accurate singly diagonally L-stable implicit embedded Runge-Kutta pair of orders 5 and 4 with explicit 

825 first stage from [here](https://doi.org/10.1016/j.apnum.2018.10.007). 

826 This method is part of the IMEX method ARK548L2SA. 

827 """ 

828 

829 generator = RK_SCHEMES["ARK548L2SAESDIRK2"]() 

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

831 ButcherTableauClass = ButcherTableauEmbedded 

832 

833 @classmethod 

834 def get_update_order(cls): 

835 """ 

836 Order in dt of the embedded error estimate, i.e. of the local error of the fourth order secondary method. 

837 

838 Returns: 

839 int: 5 

840 """