Coverage for pySDC/implementations/sweeper_classes/Runge

1import numpy as np

2import logging

3from qmat.qcoeff.butcher import RK_SCHEMES

5from pySDC.core.sweeper import Sweeper, _Pars

6from pySDC.core.errors import ParameterError

7from pySDC.core.level import Level

10class ButcherTableau(object):

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

12 """

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

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)

22 self.tleft = 0.0

23 self.tright = 1.0

24 self.num_nodes = matrix.shape[0]

25 self.weights = weights

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

31 self.left_is_node = True

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

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

41 # check if the RK scheme is implicit

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

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!')

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)}')

60 self.check_weights(weights, nodes, matrix)

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)}')

73 @property

74 def globally_stiffly_accurate(self):

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

78class ButcherTableauEmbedded(ButcherTableau):

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])}')

91 @property

92 def globally_stiffly_accurate(self):

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

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 for m in range(0, M):

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

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

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

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

255

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

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

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

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

260 )

261 else:

262 lvl.u[m + 1] = rhs

263

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

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

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

267 else:

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

269

270 # indicate presence of new values at this level

271 lvl.status.updated = True

272

273 return None

274

275 def compute_end_point(self):

276 """

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

278 """

279 lvl = self.level

280

281 if lvl.f[1] is None:

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

283 if self.is_embedded():

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

285 elif self.coll.globally_stiffly_accurate:

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

287 if self.is_embedded():

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

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

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

291 else:

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

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

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

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

296 elif self.is_embedded():

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

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

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

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

301

302 @property

303 def level(self):

304 """

305 Returns the current level

306

307 Returns:

308 pySDC.Level.level: Current level

309 """

310 return self.__level

311

312 @level.setter

313 def level(self, lvl):

314 """

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

316

317 Args:

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

319 """

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

321 if lvl.params.restol > 0:

322 lvl.params.restol = -1

323 self.logger.warning(

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

325 )

326

327 self.__level = lvl

328

329 def predict(self):

330 """

331 Predictor to fill values at nodes before first sweep

332 """

333

334 # get current level and problem

335 lvl = self.level

336 prob = lvl.prob

337

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

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

340

341 # indicate that this level is now ready for sweeps

342 lvl.status.unlocked = True

343 lvl.status.updated = True

344

345

346class RungeKuttaIMEX(RungeKutta):

347 """

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

349 """

350

351 matrix_explicit = None

352 weights_explicit = None

353 ButcherTableauClass_explicit = ButcherTableau

354

355 def __init__(self, params, level):

356 """

357 Initialization routine

358

359 Args:

360 params: parameters for the sweeper

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

362 """

363 super().__init__(params, level)

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

365 self.coll_explicit = self.get_Butcher_tableau_explicit()

366 self.QE = self.coll_explicit.Qmat

367

368 def predict(self):

369 """

370 Predictor to fill values at nodes before first sweep

371 """

372

373 # get current level and problem

374 lvl = self.level

375 prob = lvl.prob

376

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

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

379 lvl.f[m] = prob.dtype_f(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 @classmethod

386 def get_Butcher_tableau_explicit(cls):

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

388

389 def integrate(self):

390 """

391 Integrates the right-hand side

392

393 Returns:

394 list of dtype_u: containing the integral as values

395 """

396

397 # get current level and problem

398 lvl = self.level

399 prob = lvl.prob

400

401 me = []

402

403 # integrate RHS over all collocation nodes

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

405 # new instance of dtype_u, initialize values with 0

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

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

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

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

410 )

411

412 return me

413

414 def update_nodes(self):

415 """

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

417

418 Returns:

419 None

420 """

421

422 # get current level and problem

423 lvl = self.level

424 prob = lvl.prob

425

426 # only if the level has been touched before

427 assert lvl.status.unlocked

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

429

430 # get number of collocation nodes for easier access

431 M = self.coll.num_nodes

432

433 for m in range(0, M):

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

435 rhs = lvl.u[0]

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

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

438

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

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

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

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

443 )

444 else:

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

446

447 # update function values

448 if (

449 m < M - 1

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

451 or self.is_embedded()

452 ):

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

454 else:

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

456

457 # indicate presence of new values at this level

458 lvl.status.updated = True

459

460 return None

461

462 def compute_end_point(self):

463 """

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

465 """

466 lvl = self.level

467

468 if lvl.f[1] is None:

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

470 if self.is_embedded():

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

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

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

474 if self.is_embedded():

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

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

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

478 else:

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

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

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

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

483 elif self.is_embedded():

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

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

486 self.coll.weights[0],

487 self.coll.weights[1],

488 self.coll_explicit.weights[0],

489 self.coll_explicit.weights[1],

490 lvl.f[1:],

491 strict=True,

492 ):

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

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

495

496

497class ForwardEuler(RungeKutta):

498 """

499 Forward Euler. Still a classic.

500

501 Not very stable first order method.

502 """

503

504 generator = RK_SCHEMES["FE"]()

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

506

507

508class BackwardEuler(RungeKutta):

509 """

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

511

512 A-stable first order method.

513 """

514

515 generator = RK_SCHEMES["BE"]()

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

517

518

519class IMEXEuler(RungeKuttaIMEX):

520 nodes = BackwardEuler.nodes

521 weights = BackwardEuler.weights

522

523 matrix = BackwardEuler.matrix

524 matrix_explicit = ForwardEuler.matrix

525

526

527class IMEXEulerStifflyAccurate(RungeKuttaIMEX):

528 """

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

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

531 accurate and suitable for some DAEs.

532 """

533

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

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

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

537

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

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

540

541

542class CrankNicolson(RungeKutta):

543 """

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

545 """

546

547 generator = RK_SCHEMES["CN"]()

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

549

550

551class ExplicitMidpointMethod(RungeKutta):

552 """

553 Explicit Runge-Kutta method of second order.

554 """

555

556 generator = RK_SCHEMES["RK2"]()

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

558

559

560class ImplicitMidpointMethod(RungeKutta):

561 """

562 Implicit Runge-Kutta method of second order.

563 """

564

565 generator = RK_SCHEMES["IMP"]()

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

567

568

569class RK4(RungeKutta):

570 """

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

572 """

573

574 generator = RK_SCHEMES["RK4"]()

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

576

577

578class Heun_Euler(RungeKutta):

579 """

580 Second order explicit embedded Runge-Kutta method.

581 """

582

583 ButcherTableauClass = ButcherTableauEmbedded

584

585 generator = RK_SCHEMES["HEUN"]()

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

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

588 weights[0] = _weights

589 weights[1] = matrix[-1]

590

591 @classmethod

592 def get_update_order(cls):

593 return 2

594

595

596class Cash_Karp(RungeKutta):

597 """

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

599 """

600

601 generator = RK_SCHEMES["CashKarp"]()

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

603 ButcherTableauClass = ButcherTableauEmbedded

604

605 @classmethod

606 def get_update_order(cls):

607 return 5

608

609

610class DIRK43(RungeKutta):

611 """

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

613

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

615 """

616

617 generator = RK_SCHEMES["EDIRK43"]()

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

619 ButcherTableauClass = ButcherTableauEmbedded

620

621 @classmethod

622 def get_update_order(cls):

623 return 4

624

625

626class DIRK43_2(RungeKutta):

627 """

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

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

630 """

631

632 generator = RK_SCHEMES["DIRK43"]()

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

634

635

636class EDIRK4(RungeKutta):

637 """

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

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

640 """

641

642 generator = RK_SCHEMES["EDIRK4"]()

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

644

645

646class ESDIRK53(RungeKutta):

647 """

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

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

650 """

651

652 generator = RK_SCHEMES["ESDIRK53"]()

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

654 ButcherTableauClass = ButcherTableauEmbedded

655

656 @classmethod

657 def get_update_order(cls):

658 return 4

659

660

661class ESDIRK43(RungeKutta):

662 """

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

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

665 """

666

667 generator = RK_SCHEMES["ESDIRK43"]()

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

669 ButcherTableauClass = ButcherTableauEmbedded

670

671 @classmethod

672 def get_update_order(cls):

673 return 4

674

675

676class ARK548L2SAERK(RungeKutta):

677 """

678 Explicit part of the ARK54 scheme.

679 """

680

681 generator = RK_SCHEMES["ARK548L2SAERK"]()

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

683 ButcherTableauClass = ButcherTableauEmbedded

684

685 @classmethod

686 def get_update_order(cls):

687 return 5

688

689

690class ARK548L2SAESDIRK(ARK548L2SAERK):

691 """

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

693 """

694

695 generator_IMP = RK_SCHEMES["ARK548L2SAESDIRK"]()

696 matrix = generator_IMP.Q

697

698

699class ARK54(RungeKuttaIMEX):

700 """

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

702 """

703

704 ButcherTableauClass = ButcherTableauEmbedded

705 ButcherTableauClass_explicit = ButcherTableauEmbedded

706

707 nodes = ARK548L2SAERK.nodes

708 weights = ARK548L2SAERK.weights

709

710 matrix = ARK548L2SAESDIRK.matrix

711 matrix_explicit = ARK548L2SAERK.matrix

712

713 @classmethod

714 def get_update_order(cls):

715 return 5

716

717

718class ARK548L2SAESDIRK2(RungeKutta):

719 """

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

721 This method is part of the IMEX method ARK548L2SA.

722 """

723

724 generator = RK_SCHEMES["ARK548L2SAESDIRK2"]()

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

726 ButcherTableauClass = ButcherTableauEmbedded

727

728 @classmethod

729 def get_update_order(cls):

730 return 5

731

732

733class ARK548L2SAERK2(ARK548L2SAESDIRK2):

734 """

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

736 This method is part of the IMEX method ARK548L2SA.

737 """

738

739 generator_EXP = RK_SCHEMES["ARK548L2SAERK2"]()

740 matrix = generator_EXP.Q

741

742

743class ARK548L2SA(RungeKuttaIMEX):

744 """

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

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

747

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

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

750 """

751

752 ButcherTableauClass = ButcherTableauEmbedded

753 ButcherTableauClass_explicit = ButcherTableauEmbedded

754

755 nodes = ARK548L2SAERK2.nodes

756 weights = ARK548L2SAERK2.weights

757

758 matrix = ARK548L2SAESDIRK2.matrix

759 matrix_explicit = ARK548L2SAERK2.matrix

760

761 @classmethod

762 def get_update_order(cls):

763 return 5

764

765

766class ARK324L2SAERK(RungeKutta):

767 generator = RK_SCHEMES["ARK324L2SAERK"]()

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

769 ButcherTableauClass = ButcherTableauEmbedded

770

771 @classmethod

772 def get_update_order(cls):

773 return 3

774

775

776class ARK324L2SAESDIRK(ARK324L2SAERK):

777 generator = RK_SCHEMES["ARK324L2SAESDIRK"]()

778 matrix = generator.Q

779

780

781class ARK32(RungeKuttaIMEX):

782 ButcherTableauClass = ButcherTableauEmbedded

783 ButcherTableauClass_explicit = ButcherTableauEmbedded

784

785 nodes = ARK324L2SAESDIRK.nodes

786 weights = ARK324L2SAESDIRK.weights

787

788 matrix = ARK324L2SAESDIRK.matrix

789 matrix_explicit = ARK324L2SAERK.matrix

790

791 @classmethod

792 def get_update_order(cls):

793 return 3

794

795

796class ARK2(RungeKuttaIMEX):

797 """

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

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

800 """

801

802 generator_IMP = RK_SCHEMES["ARK222EDIRK"]()

803 generator_EXP = RK_SCHEMES["ARK222ERK"]()

804

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

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

807

808

809class ARK3(RungeKuttaIMEX):

810 """

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

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

813 """

814

815 generator_IMP = RK_SCHEMES["ARK443ESDIRK"]()

816 generator_EXP = RK_SCHEMES["ARK443ERK"]()

817

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

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

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

362 statements