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): 

129 """ 

130 Initialization routine for the custom sweeper 

131 

132 Args: 

133 params: parameters for the sweeper 

134 """ 

135 # set up logger 

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

137 

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

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

140 if key in params: 

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

142 

143 # set parameters to their actual values 

144 self.coll = self.get_Butcher_tableau() 

145 params['initial_guess'] = 'zero' 

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

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

148 

149 # disable residual computation by default 

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

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

152 ) 

153 

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

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

156 

157 self.params = _Pars(params) 

158 

159 # This will be set as soon as the sweeper is instantiated at the level 

160 self.__level = None 

161 

162 self.parallelizable = False 

163 self.QI = self.coll.Qmat 

164 

165 @classmethod 

166 def get_Q_matrix(cls): 

167 return cls.get_Butcher_tableau().Qmat 

168 

169 @classmethod 

170 def get_Butcher_tableau(cls): 

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

172 

173 @classmethod 

174 def get_update_order(cls): 

175 """ 

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

177 """ 

178 raise NotImplementedError( 

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

180 ) 

181 

182 def get_full_f(self, f): 

183 """ 

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

185 

186 Args: 

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

188 

189 Returns: 

190 mesh: Full right hand side as a mesh 

191 """ 

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

193 return f 

194 elif type(f).__name__ in ['imex_mesh', 'imex_cupy_mesh']: 

195 return f.impl + f.expl 

196 elif f is None: 

197 prob = self.level.prob 

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

199 else: 

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

201 

202 def integrate(self): 

203 """ 

204 Integrates the right-hand side 

205 

206 Returns: 

207 list of dtype_u: containing the integral as values 

208 """ 

209 

210 # get current level and problem 

211 lvl = self.level 

212 prob = lvl.prob 

213 

214 me = [] 

215 

216 # integrate RHS over all collocation nodes 

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

218 # new instance of dtype_u, initialize values with 0 

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

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

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

222 

223 return me 

224 

225 def update_nodes(self): 

226 """ 

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

228 

229 Returns: 

230 None 

231 """ 

232 

233 # get current level and problem 

234 lvl = self.level 

235 prob = lvl.prob 

236 

237 # only if the level has been touched before 

238 assert lvl.status.unlocked 

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

240 

241 # get number of collocation nodes for easier access 

242 M = self.coll.num_nodes 

243 

244 for m in range(0, M): 

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

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

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

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

249 

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

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

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

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

254 ) 

255 else: 

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

257 

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

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

260 

261 # indicate presence of new values at this level 

262 lvl.status.updated = True 

263 

264 return None 

265 

266 def compute_end_point(self): 

267 """ 

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

269 """ 

270 lvl = self.level 

271 

272 if lvl.f[1] is None: 

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

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

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

276 elif self.coll.globally_stiffly_accurate: 

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

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

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

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

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

282 else: 

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

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

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

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

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

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

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

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

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

292 

293 @property 

294 def level(self): 

295 """ 

296 Returns the current level 

297 

298 Returns: 

299 pySDC.Level.level: Current level 

300 """ 

301 return self.__level 

302 

303 @level.setter 

304 def level(self, lvl): 

305 """ 

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

307 

308 Args: 

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

310 """ 

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

312 if lvl.params.restol > 0: 

313 lvl.params.restol = -1 

314 self.logger.warning( 

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

316 ) 

317 

318 self.__level = lvl 

319 

320 def predict(self): 

321 """ 

322 Predictor to fill values at nodes before first sweep 

323 """ 

324 

325 # get current level and problem 

326 lvl = self.level 

327 prob = lvl.prob 

328 

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

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

331 

332 # indicate that this level is now ready for sweeps 

333 lvl.status.unlocked = True 

334 lvl.status.updated = True 

335 

336 

337class RungeKuttaIMEX(RungeKutta): 

338 """ 

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

340 """ 

341 

342 matrix_explicit = None 

343 weights_explicit = None 

344 ButcherTableauClass_explicit = ButcherTableau 

345 

346 def __init__(self, params): 

347 """ 

348 Initialization routine 

349 

350 Args: 

351 params: parameters for the sweeper 

352 """ 

353 super().__init__(params) 

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

355 self.coll_explicit = self.get_Butcher_tableau_explicit() 

356 self.QE = self.coll_explicit.Qmat 

357 

358 def predict(self): 

359 """ 

360 Predictor to fill values at nodes before first sweep 

361 """ 

362 

363 # get current level and problem 

364 lvl = self.level 

365 prob = lvl.prob 

366 

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

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

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

370 

371 # indicate that this level is now ready for sweeps 

372 lvl.status.unlocked = True 

373 lvl.status.updated = True 

374 

375 @classmethod 

376 def get_Butcher_tableau_explicit(cls): 

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

378 

379 def integrate(self): 

380 """ 

381 Integrates the right-hand side 

382 

383 Returns: 

384 list of dtype_u: containing the integral as values 

385 """ 

386 

387 # get current level and problem 

388 lvl = self.level 

389 prob = lvl.prob 

390 

391 me = [] 

392 

393 # integrate RHS over all collocation nodes 

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

395 # new instance of dtype_u, initialize values with 0 

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

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

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

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

400 ) 

401 

402 return me 

403 

404 def update_nodes(self): 

405 """ 

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

407 

408 Returns: 

409 None 

410 """ 

411 

412 # get current level and problem 

413 lvl = self.level 

414 prob = lvl.prob 

415 

416 # only if the level has been touched before 

417 assert lvl.status.unlocked 

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

419 

420 # get number of collocation nodes for easier access 

421 M = self.coll.num_nodes 

422 

423 for m in range(0, M): 

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

425 rhs = lvl.u[0] 

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

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

428 

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

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

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

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

433 ) 

434 else: 

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

436 

437 # update function values 

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

439 

440 # indicate presence of new values at this level 

441 lvl.status.updated = True 

442 

443 return None 

444 

445 def compute_end_point(self): 

446 """ 

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

448 """ 

449 lvl = self.level 

450 

451 if lvl.f[1] is None: 

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

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

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

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

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

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

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

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

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

461 else: 

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

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

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

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

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

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

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

469 self.coll.weights[0], 

470 self.coll.weights[1], 

471 self.coll_explicit.weights[0], 

472 self.coll_explicit.weights[1], 

473 lvl.f[1:], 

474 ): 

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

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

477 

478 

479class ForwardEuler(RungeKutta): 

480 """ 

481 Forward Euler. Still a classic. 

482 

483 Not very stable first order method. 

484 """ 

485 

486 generator = RK_SCHEMES["FE"]() 

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

488 

489 

490class BackwardEuler(RungeKutta): 

491 """ 

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

493 

494 A-stable first order method. 

495 """ 

496 

497 generator = RK_SCHEMES["BE"]() 

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

499 

500 

501class IMEXEuler(RungeKuttaIMEX): 

502 nodes = BackwardEuler.nodes 

503 weights = BackwardEuler.weights 

504 

505 matrix = BackwardEuler.matrix 

506 matrix_explicit = ForwardEuler.matrix 

507 

508 

509class CrankNicolson(RungeKutta): 

510 """ 

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

512 """ 

513 

514 generator = RK_SCHEMES["CN"]() 

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

516 

517 

518class ExplicitMidpointMethod(RungeKutta): 

519 """ 

520 Explicit Runge-Kutta method of second order. 

521 """ 

522 

523 generator = RK_SCHEMES["RK2"]() 

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

525 

526 

527class ImplicitMidpointMethod(RungeKutta): 

528 """ 

529 Implicit Runge-Kutta method of second order. 

530 """ 

531 

532 generator = RK_SCHEMES["IMP"]() 

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

534 

535 

536class RK4(RungeKutta): 

537 """ 

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

539 """ 

540 

541 generator = RK_SCHEMES["RK4"]() 

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

543 

544 

545class Heun_Euler(RungeKutta): 

546 """ 

547 Second order explicit embedded Runge-Kutta method. 

548 """ 

549 

550 ButcherTableauClass = ButcherTableauEmbedded 

551 

552 generator = RK_SCHEMES["HEUN"]() 

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

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

555 weights[0] = _weights 

556 weights[1] = matrix[-1] 

557 

558 @classmethod 

559 def get_update_order(cls): 

560 return 2 

561 

562 

563class Cash_Karp(RungeKutta): 

564 """ 

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

566 """ 

567 

568 generator = RK_SCHEMES["CashKarp"]() 

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

570 ButcherTableauClass = ButcherTableauEmbedded 

571 

572 @classmethod 

573 def get_update_order(cls): 

574 return 5 

575 

576 

577class DIRK43(RungeKutta): 

578 """ 

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

580 

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

582 """ 

583 

584 generator = RK_SCHEMES["EDIRK43"]() 

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

586 ButcherTableauClass = ButcherTableauEmbedded 

587 

588 @classmethod 

589 def get_update_order(cls): 

590 return 4 

591 

592 

593class DIRK43_2(RungeKutta): 

594 """ 

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

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

597 """ 

598 

599 generator = RK_SCHEMES["DIRK43"]() 

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

601 

602 

603class EDIRK4(RungeKutta): 

604 """ 

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

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

607 """ 

608 

609 generator = RK_SCHEMES["EDIRK4"]() 

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

611 

612 

613class ESDIRK53(RungeKutta): 

614 """ 

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

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

617 """ 

618 

619 generator = RK_SCHEMES["ESDIRK53"]() 

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

621 ButcherTableauClass = ButcherTableauEmbedded 

622 

623 @classmethod 

624 def get_update_order(cls): 

625 return 4 

626 

627 

628class ESDIRK43(RungeKutta): 

629 """ 

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

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

632 """ 

633 

634 generator = RK_SCHEMES["ESDIRK43"]() 

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

636 ButcherTableauClass = ButcherTableauEmbedded 

637 

638 @classmethod 

639 def get_update_order(cls): 

640 return 4 

641 

642 

643class ARK548L2SAERK(RungeKutta): 

644 """ 

645 Explicit part of the ARK54 scheme. 

646 """ 

647 

648 generator = RK_SCHEMES["ARK548L2SAERK"]() 

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

650 ButcherTableauClass = ButcherTableauEmbedded 

651 

652 @classmethod 

653 def get_update_order(cls): 

654 return 5 

655 

656 

657class ARK548L2SAESDIRK(ARK548L2SAERK): 

658 """ 

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

660 """ 

661 

662 generator_IMP = RK_SCHEMES["ARK548L2SAESDIRK"]() 

663 matrix = generator_IMP.Q 

664 

665 

666class ARK54(RungeKuttaIMEX): 

667 """ 

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

669 """ 

670 

671 ButcherTableauClass = ButcherTableauEmbedded 

672 ButcherTableauClass_explicit = ButcherTableauEmbedded 

673 

674 nodes = ARK548L2SAERK.nodes 

675 weights = ARK548L2SAERK.weights 

676 

677 matrix = ARK548L2SAESDIRK.matrix 

678 matrix_explicit = ARK548L2SAERK.matrix 

679 

680 @classmethod 

681 def get_update_order(cls): 

682 return 5 

683 

684 

685class ARK548L2SAESDIRK2(RungeKutta): 

686 """ 

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

688 This method is part of the IMEX method ARK548L2SA. 

689 """ 

690 

691 generator = RK_SCHEMES["ARK548L2SAESDIRK2"]() 

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

693 ButcherTableauClass = ButcherTableauEmbedded 

694 

695 @classmethod 

696 def get_update_order(cls): 

697 return 5 

698 

699 

700class ARK548L2SAERK2(ARK548L2SAESDIRK2): 

701 """ 

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

703 This method is part of the IMEX method ARK548L2SA. 

704 """ 

705 

706 generator_EXP = RK_SCHEMES["ARK548L2SAERK2"]() 

707 matrix = generator_EXP.Q 

708 

709 

710class ARK548L2SA(RungeKuttaIMEX): 

711 """ 

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

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

714 

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

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

717 """ 

718 

719 ButcherTableauClass = ButcherTableauEmbedded 

720 ButcherTableauClass_explicit = ButcherTableauEmbedded 

721 

722 nodes = ARK548L2SAERK2.nodes 

723 weights = ARK548L2SAERK2.weights 

724 

725 matrix = ARK548L2SAESDIRK2.matrix 

726 matrix_explicit = ARK548L2SAERK2.matrix 

727 

728 @classmethod 

729 def get_update_order(cls): 

730 return 5 

731 

732 

733class ARK324L2SAERK(RungeKutta): 

734 generator = RK_SCHEMES["ARK324L2SAERK"]() 

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

736 ButcherTableauClass = ButcherTableauEmbedded 

737 

738 @classmethod 

739 def get_update_order(cls): 

740 return 3 

741 

742 

743class ARK324L2SAESDIRK(ARK324L2SAERK): 

744 generator = RK_SCHEMES["ARK324L2SAESDIRK"]() 

745 matrix = generator.Q 

746 

747 

748class ARK32(RungeKuttaIMEX): 

749 ButcherTableauClass = ButcherTableauEmbedded 

750 ButcherTableauClass_explicit = ButcherTableauEmbedded 

751 

752 nodes = ARK324L2SAESDIRK.nodes 

753 weights = ARK324L2SAESDIRK.weights 

754 

755 matrix = ARK324L2SAESDIRK.matrix 

756 matrix_explicit = ARK324L2SAERK.matrix 

757 

758 @classmethod 

759 def get_update_order(cls): 

760 return 3 

761 

762 

763class ARK2(RungeKuttaIMEX): 

764 """ 

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

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

767 """ 

768 

769 generator_IMP = RK_SCHEMES["ARK222EDIRK"]() 

770 generator_EXP = RK_SCHEMES["ARK222ERK"]() 

771 

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

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

774 

775 

776class ARK3(RungeKuttaIMEX): 

777 """ 

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

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

780 """ 

781 

782 generator_IMP = RK_SCHEMES["ARK443ESDIRK"]() 

783 generator_EXP = RK_SCHEMES["ARK443ERK"]() 

784 

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

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