include/pops/numerics/time/integrators/time_steppers.hpp Source File

adc_cpp 0.3.0 Model-free C++23 core for coupled hyperbolic-elliptic systems on adaptive (AMR) meshes, with MPI and GPU (Kokkos) backends

• include
• pops
• numerics
• time
• integrators

time_steppers.hpp

Go to the documentation of this file.

1#pragma once

2

3#include <pops/core/foundation/types.hpp>

4#include <pops/mesh/storage/mf_arith.hpp> // saxpy, lincomb

5#include <pops/mesh/storage/multifab.hpp>

6

7#include <utility>

8

29

30namespace pops {

31

32// Contract: an integrator knows how to advance U by dt via a residual evaluator.

33template <class I>

34concept TimeStepper = requires(const I integ, MultiFab& U, Real dt) {

35 integ.take_step([](MultiFab&, MultiFab&) {}, U, dt);

36};

37

38// Forward Euler (order 1): U <- U + dt R(U).

39struct ForwardEuler {

40 // Reusable residual buffer, sized from a U layout. Hoist it out of a substep loop (see

41 // run_explicit_substeps) to avoid re-allocating per substep; rhs() overwrites R every call,

42 // so reuse is bit-identical.

43 struct Scratch {

44 MultiFab R;

45 explicit Scratch(const MultiFab& U) : R(U.box_array(), U.dmap(), U.ncomp(), 0) {}

46 };

47 template <class RhsEval>

48 void take_step(RhsEval&& rhs, MultiFab& U, Real dt, Scratch& s) const {

49 rhs(U, s.R);

50 saxpy(U, dt, s.R);

51 }

52 template <class RhsEval>

53 void take_step(RhsEval&& rhs, MultiFab& U, Real dt) const {

54 Scratch s(U);

55 take_step(std::forward<RhsEval>(rhs), U, dt, s);

56 }

57};

58

59// SSP-RK2 (Shu-Osher, 2 stages, order 2). Same operations as the old

60// SystemCoupler::advance_explicit_ssprk2 / ssprk.hpp::advance_ssprk2 (bit-identical).

61struct SSPRK2Step {

62 // Reusable stage buffers (residual R with 0 ghosts; stage U1 with U's ghosts, since rhs reads

63 // its ghosts). Hoist out of a substep loop via run_explicit_substeps. Reuse is bit-identical:

64 // R is overwritten by rhs and U1's valid cells are overwritten by the lincomb copy each substep,

65 // while U1's ghosts are re-derived by rhs's internal fill_ghosts before any ghost read.

66 struct Scratch {

67 MultiFab R, U1;

68 explicit Scratch(const MultiFab& U)

69 : R(U.box_array(), U.dmap(), U.ncomp(), 0),

70 U1(U.box_array(), U.dmap(), U.ncomp(), U.n_grow()) {}

71 };

72 template <class RhsEval>

73 void take_step(RhsEval&& rhs, MultiFab& U, Real dt, Scratch& s) const {

74 rhs(U, s.R);

75 lincomb(s.U1, Real(1), U, Real(0), U); // U1 = U (valid cells; ghosts re-derived by rhs)

76 saxpy(s.U1, dt, s.R);

77 rhs(s.U1, s.R);

78 saxpy(s.U1, dt, s.R);

79 lincomb(U, Real(0.5), U, Real(0.5), s.U1);

80 }

81 template <class RhsEval>

82 void take_step(RhsEval&& rhs, MultiFab& U, Real dt) const {

83 Scratch s(U);

84 take_step(std::forward<RhsEval>(rhs), U, dt, s);

85 }

86};

87

88// SSP-RK3 (Shu-Osher, 3 stages, order 3). Same operations as the old

89// SystemCoupler::advance_explicit_ssprk3 (bit-identical).

90struct SSPRK3Step {

91 // Reusable stage buffers (residual R with 0 ghosts; stages U1/U2/U3 with U's ghosts, since the

92 // first two are passed back to rhs which reads their ghosts). Hoist out of a substep loop via

93 // run_explicit_substeps. Reuse is bit-identical: every buffer is fully overwritten each substep

94 // (R by rhs, U1/U2/U3 valid cells by the lincomb copy), and stage ghosts are re-derived by rhs.

95 struct Scratch {

96 MultiFab R, U1, U2, U3;

97 explicit Scratch(const MultiFab& U)

98 : R(U.box_array(), U.dmap(), U.ncomp(), 0),

99 U1(U.box_array(), U.dmap(), U.ncomp(), U.n_grow()),

100 U2(U.box_array(), U.dmap(), U.ncomp(), U.n_grow()),

101 U3(U.box_array(), U.dmap(), U.ncomp(), U.n_grow()) {}

102 };

103 template <class RhsEval>

104 void take_step(RhsEval&& rhs, MultiFab& U, Real dt, Scratch& s) const {

105 rhs(U, s.R);

106 lincomb(s.U1, Real(1), U, Real(0), U); // U1 = U

107 saxpy(s.U1, dt, s.R);

108

109 rhs(s.U1, s.R);

110 lincomb(s.U2, Real(1), s.U1, Real(0), s.U1); // U2 = U1

111 saxpy(s.U2, dt, s.R);

112 lincomb(s.U2, Real(3) / 4, U, Real(1) / 4, s.U2);

113

114 rhs(s.U2, s.R);

115 lincomb(s.U3, Real(1), s.U2, Real(0), s.U2); // U3 = U2

116 saxpy(s.U3, dt, s.R);

117 lincomb(U, Real(1) / 3, U, Real(2) / 3, s.U3);

118 }

119 template <class RhsEval>

120 void take_step(RhsEval&& rhs, MultiFab& U, Real dt) const {

121 Scratch s(U);

122 take_step(std::forward<RhsEval>(rhs), U, dt, s);

123 }

124};

125

126// Runs @p n explicit RK substeps of @c Stepper on @p U with step @p h. When the stepper exposes a

127// reusable Scratch (the built-in ForwardEuler / SSPRK2Step / SSPRK3Step), the scratch is allocated

128// ONCE here and reused across substeps via the scratch-taking take_step overload -- removing the

129// per-substep alloc/zero/free churn (ADC-261). A custom TimeStepper without a Scratch falls back to

130// the one-shot take_step. The result is bit-identical either way (same saxpy/lincomb sequence on

131// freshly-overwritten buffers); the substep loop never changes U's layout, so one Scratch suffices.

132template <class Stepper, class RhsEval>

133inline void run_explicit_substeps(RhsEval&& rhs, MultiFab& U, Real h, int n) {

134 // Probe the actual scratch-taking overload (not just a nested type named Scratch): a custom

135 // TimeStepper that happens to expose an unrelated Scratch type but no four-arg take_step still

136 // takes the one-shot fallback instead of hitting a hard error.

137 if constexpr (requires(Stepper st, RhsEval r, MultiFab& u, Real dt,

138 typename Stepper::Scratch& sc) { st.take_step(r, u, dt, sc); }) {

139 typename Stepper::Scratch scratch(U);

140 for (int s = 0; s < n; ++s)

141 Stepper{}.take_step(rhs, U, h, scratch);

142 } else {

143 for (int s = 0; s < n; ++s)

144 Stepper{}.take_step(rhs, U, h);

145 }

146}

147

148} // namespace pops

pops::MultiFab

Field distributed over a level: decomposition (BoxArray) + distribution (DistributionMapping) + ncomp...

Definition multifab.hpp:33

pops::TimeStepper

Definition time_steppers.hpp:34

mf_arith.hpp

MultiFab arithmetic (saxpy, lincomb, norm_inf, dot) over VALID cells.

multifab.hpp

MultiFab: a field DISTRIBUTED over a level (equivalent of AMReX's MultiFab).

pops

Definition amr_hierarchy.hpp:29

pops::Real

double Real

Definition types.hpp:30

pops::saxpy

void saxpy(MultiFab &y, Real a, const MultiFab &x)

y <- y + a x over ALL components of the valid cells. Identical layouts required.

Definition mf_arith.hpp:102

pops::run_explicit_substeps

void run_explicit_substeps(RhsEval &&rhs, MultiFab &U, Real h, int n)

Definition time_steppers.hpp:133

pops::lincomb

void lincomb(MultiFab &z, Real a, const MultiFab &x, Real b, const MultiFab &y)

z <- a x + b y over ALL components of the valid cells. Identical layouts; aliasing safe.

Definition mf_arith.hpp:133

pops::ForwardEuler::Scratch

Definition time_steppers.hpp:43

pops::ForwardEuler::Scratch::R

MultiFab R

Definition time_steppers.hpp:44

pops::ForwardEuler::Scratch::Scratch

Scratch(const MultiFab &U)

Definition time_steppers.hpp:45

pops::ForwardEuler

Definition time_steppers.hpp:39

pops::ForwardEuler::take_step

void take_step(RhsEval &&rhs, MultiFab &U, Real dt, Scratch &s) const

Definition time_steppers.hpp:48

pops::ForwardEuler::take_step

void take_step(RhsEval &&rhs, MultiFab &U, Real dt) const

Definition time_steppers.hpp:53

pops::SSPRK2Step::Scratch

Definition time_steppers.hpp:66

pops::SSPRK2Step::Scratch::U1

MultiFab U1

Definition time_steppers.hpp:67

pops::SSPRK2Step::Scratch::R

MultiFab R

Definition time_steppers.hpp:67

pops::SSPRK2Step::Scratch::Scratch

Scratch(const MultiFab &U)

Definition time_steppers.hpp:68

pops::SSPRK2Step

Definition time_steppers.hpp:61

pops::SSPRK2Step::take_step

void take_step(RhsEval &&rhs, MultiFab &U, Real dt) const

Definition time_steppers.hpp:82

pops::SSPRK2Step::take_step

void take_step(RhsEval &&rhs, MultiFab &U, Real dt, Scratch &s) const

Definition time_steppers.hpp:73

pops::SSPRK3Step::Scratch

Definition time_steppers.hpp:95

pops::SSPRK3Step::Scratch::U1

MultiFab U1

Definition time_steppers.hpp:96

pops::SSPRK3Step::Scratch::Scratch

Scratch(const MultiFab &U)

Definition time_steppers.hpp:97

pops::SSPRK3Step::Scratch::R

MultiFab R

Definition time_steppers.hpp:96

pops::SSPRK3Step::Scratch::U2

MultiFab U2

Definition time_steppers.hpp:96

pops::SSPRK3Step::Scratch::U3

MultiFab U3

Definition time_steppers.hpp:96

pops::SSPRK3Step

Definition time_steppers.hpp:90

pops::SSPRK3Step::take_step

void take_step(RhsEval &&rhs, MultiFab &U, Real dt, Scratch &s) const

Definition time_steppers.hpp:104

pops::SSPRK3Step::take_step

void take_step(RhsEval &&rhs, MultiFab &U, Real dt) const

Definition time_steppers.hpp:120

types.hpp

Base scalar types and the POPS_HD macro (host+device portability).

---

Generated by  1.9.8