PoPS/doxygen/coupler_8hpp_source.html

#pragma once


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

#include <pops/coupling/base/aux_fill.hpp>  // detail::derive_aux_bc + detail::fill_bz_box (shared)

#include <pops/coupling/base/coupling_policy.hpp>

#include <pops/coupling/base/elliptic_rhs.hpp>

#include <pops/numerics/elliptic/interface/elliptic_problem.hpp>

#include <pops/numerics/elliptic/interface/elliptic_solver.hpp>

#include <pops/numerics/elliptic/mg/geometric_mg.hpp>

#include <pops/numerics/time/integrators/time_integrator.hpp>

#include <pops/numerics/time/integrators/time_steppers.hpp>  // SSPRK2Step / SSPRK3Step (shared scheme)

#include <pops/mesh/layout/box_array.hpp>

#include <pops/mesh/layout/distribution_mapping.hpp>

#include <pops/mesh/storage/fab2d.hpp>

#include <pops/mesh/execution/for_each.hpp>

#include <pops/mesh/geometry/geometry.hpp>

#include <pops/mesh/storage/mf_arith.hpp>

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

#include <pops/mesh/boundary/physical_bc.hpp>

#include <pops/numerics/fv/reconstruction.hpp>

#include <pops/numerics/fv/spatial_discretisation.hpp>

#include <pops/numerics/spatial_operator.hpp>

#include <pops/parallel/comm.hpp>


#include <functional>

#include <type_traits>

#include <utility>


namespace pops {


namespace detail {

// Namespace-scope helpers: an extended __host__ __device__ lambda CANNOT be

// defined in a private/protected method (nvcc restriction), hence the

// extraction out of the Coupler class.


// Single-model compatibility: f = model.elliptic_rhs(U) on valid cells,

// delegated to a named assembler so this responsibility is not buried in Coupler.

template <class Model>


inline void coupler_eval_rhs(const MultiFab& state, MultiFab& rhs, const Model& model) {

  SingleModelEllipticRhs<Model>{model}(state, rhs);

}


// aux = (phi, d phi/dx, d phi/dy) by centered differences. Delegates to the

// named FieldPostProcess convention with GradSign::Plus and store_phi=true: the

// coupler stores +grad phi (the physical sign E = -grad phi is carried by the

// transport drift velocity). Multiplicative form *cx / *cy kept identical

// -> bit-identical.


inline void coupler_grad_phi(const MultiFab& phi, MultiFab& aux, Real cx, Real cy) {

  field_postprocess(phi, aux, cx, cy, FieldPostProcess{FieldPostProcess::GradSign::Plus, true});

}


}  // namespace detail


template <class Model, class Elliptic = GeometricMG>


class Coupler {

  static_assert(EllipticSolver<Elliptic>, "the Coupler elliptic backend must model EllipticSolver");


 public:

  // active: optional "inside the conductor" predicate (embedded wall for

  // the Poisson solver). Empty => no internal wall.

  // bz: out-of-plane magnetic field B_z(x, y) PROVIDED by the user (constant or

  // field). Only has effect if the model declares the B_z aux component (aux_comps>3);

  // then fills aux component 3 once and for all (B_z static, external to the

  // elliptic solve: derive_aux does not touch it). Empty => no B_z. The aux channel is

  // allocated to the MODEL width: a base model (3) stays bit-identical.


  Coupler(const Model& model, const Geometry& geom, const BoxArray& ba, const BCRec& bcU,

          const BCRec& bcPhi, std::function<bool(Real, Real)> active = {},

          std::function<Real(Real, Real)> bz = {})

      : model_(model),

        geom_(geom),

        ba_(ba),

        dm_(ba.size(), n_ranks()),

        bcU_(bcU),

        bcPhi_(bcPhi),

        aux_bc_(detail::derive_aux_bc(bcPhi)),

        mg_(geom, ba, bcPhi, std::move(active)),

        aux_(ba, dm_, aux_comps<Model>(), 1),

        bz_(std::move(bz)) {

    fill_bz();  // fills the B_z component (no-op if base model or empty bz)

  }


  // Coupled SSPRK2. Three orthogonal axes, all template parameters:

  //   - Limiter: reconstruction (NoSlope / Minmod / VanLeer ...)

  //   - Policy: time coupling (PerStage = phi at each stage; OncePerStep

  //                    = a single solve per step, aux frozen)

  //   - NumericalFlux: Riemann flux (Rusanov by default, HLL, HLLC ...)

  // U must carry at least Limiter::n_ghost ghosts. The historical signature

  // advance<Limiter, Policy> stays valid (NumericalFlux default = Rusanov).

  template <class Limiter = NoSlope, class Policy = PerStageCoupling,

            class NumericalFlux = RusanovFlux>


  void advance(MultiFab& U, Real dt) {

    static_assert(

        std::is_same_v<Policy, PerStageCoupling> || std::is_same_v<Policy, OncePerStepCoupling>,

        "Policy must be PerStageCoupling or OncePerStepCoupling");

    constexpr bool per = std::is_same_v<Policy, PerStageCoupling>;

    // DELEGATES the scheme to the core SSPRK2Step object (dedup, sec.8.2 A4). The residual

    // evaluator counts the stages: recompute_aux=true at stage 0, =per afterward (PerStage:

    // phi recomputed for the intermediate state; OncePerStep: aux frozen). Bit-identical.

    int stage = 0;

    SSPRK2Step{}.take_step(

        [&](MultiFab& s, MultiFab& R) {

          stage_rhs<Limiter, NumericalFlux>(s, R, (stage++ == 0) ? true : per);

        },

        U, dt);

  }


  // Coupled SSPRK3 (Shu-Osher, 3 stages). Same axes as advance.

  template <class Limiter = NoSlope, class Policy = PerStageCoupling,

            class NumericalFlux = RusanovFlux>


  void advance_ssprk3(MultiFab& U, Real dt) {

    static_assert(

        std::is_same_v<Policy, PerStageCoupling> || std::is_same_v<Policy, OncePerStepCoupling>,

        "Policy must be PerStageCoupling or OncePerStepCoupling");

    constexpr bool per = std::is_same_v<Policy, PerStageCoupling>;

    // Same as advance: delegates to SSPRK3Step, recompute_aux=true at stage 0, =per afterward.

    int stage = 0;

    SSPRK3Step{}.take_step(

        [&](MultiFab& s, MultiFab& R) {

          stage_rhs<Limiter, NumericalFlux>(s, R, (stage++ == 0) ? true : per);

        },

        U, dt);

  }


  // Unified entry point: give the coupler a SPATIAL DISCRETISATION

  // (limiter + flux) and an explicit time policy. The old SSPRK2/SSPRK3 tags

  // stay valid; the new form also enables sub-cycling:

  //   sim.step<MusclVanLeerHLLC, ExplicitTime<SSPRK3, 4>>(U, dt);

  template <class Disc = FirstOrder, class TimeInteg = SSPRK2, class Policy = PerStageCoupling>


  void step(MultiFab& U, Real dt) {

    using L = typename Disc::Limiter;

    using F = typename Disc::NumericalFlux;

    using T = typename TimePolicyTraits<TimeInteg>::Method;

    static_assert(TimePolicyTraits<TimeInteg>::treatment == TimeTreatment::Explicit,

                  "Coupler::step can only run explicit policies; "

                  "use a scheduler/system for IMEX or implicit");

    static_assert(std::is_same_v<T, SSPRK2> || std::is_same_v<T, SSPRK3>,

                  "Coupler::step supports SSPRK2 or SSPRK3");

    constexpr int n = TimePolicyTraits<TimeInteg>::substeps;

    const Real h = dt / static_cast<Real>(n);

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

      if constexpr (std::is_same_v<T, SSPRK3>)

        advance_ssprk3<L, Policy, F>(U, h);

      else

        advance<L, Policy, F>(U, h);

    }

  }


  void solve_fields(const MultiFab& U) { update_aux(U); }


  MultiFab& phi() { return mg_.phi(); }

  const MultiFab& aux() const { return aux_; }


 private:

  void update_aux(const MultiFab& state) {

    detail::coupler_eval_rhs(state, mg_.rhs(), model_);

    mg_.solve();  // EllipticSolver concept interface (backend-agnostic)

    derive_aux();

  }


  // One stage: (optional) elliptic solve, halos, hyperbolic residual into R.

  // Shared by advance (SSPRK2) and advance_ssprk3; order of operations kept

  // to stay bit-identical to the old advance.

  template <class Limiter, class NumericalFlux>

  void stage_rhs(MultiFab& s, MultiFab& R, bool recompute_aux) {

    if (recompute_aux)

      update_aux(s);

    fill_ghosts(s, geom_.domain, bcU_);

    assemble_rhs<Limiter, NumericalFlux>(model_, s, aux_, geom_, R);

  }


  void derive_aux() {

    fill_ghosts(mg_.phi(), geom_.domain, bcPhi_);

    const Real cx = Real(1) / (2 * geom_.dx());

    const Real cy = Real(1) / (2 * geom_.dy());

    detail::coupler_grad_phi(mg_.phi(), aux_, cx, cy);

    fill_ghosts(aux_, geom_.domain, aux_bc_);

  }


  // Fills the B_z aux component (index kAuxBaseComps) on valid cells from

  // bz_(x, y), once only (B_z static). Compile-time guard: without a B_z field in the

  // model (aux_comps == 3) the component does not exist -> no code, no out-of-bound access.

  // The B_z halos are then maintained by derive_aux (Foextrap/periodic of aux_bc_,

  // cf. grad); field_postprocess only writes phi/grad (components 0..2), B_z is preserved.

  void fill_bz() {

    if constexpr (aux_comps<Model>() > kAuxBaseComps) {

      if (!bz_)

        return;

      for (int li = 0; li < aux_.local_size(); ++li)

        detail::fill_bz_box(aux_.fab(li), aux_.box(li), geom_, bz_);  // valid box

      fill_ghosts(aux_, geom_.domain, aux_bc_);  // B_z halos before the 1st solve

    }

  }


  Model model_;

  Geometry geom_;

  BoxArray ba_;

  DistributionMapping dm_;

  BCRec bcU_, bcPhi_, aux_bc_;

  Elliptic mg_;

  MultiFab aux_;

  std::function<Real(Real, Real)> bz_;  // external B_z(x, y) (empty if not provided)

};


// The coupler elliptic backend honors the common contract: swapping

// GeometricMG for another conforming solver (FFT wrapper, PETSc) will only

// require changing the member type, not the coupling logic.

static_assert(EllipticSolver<GeometricMG>, "GeometricMG must model the EllipticSolver concept");


}  // namespace pops

aux_fill.hpp
Helpers shared by the three couplers (single-block Coupler, SystemAssembler, AmrSystemCoupler) for th...

box_array.hpp
BoxArray: the set of boxes tiling a level (disjoint, covering).

pops::BoxArray
Ordered list of boxes tiling a level.
Definition box_array.hpp:22

pops::Coupler
Single-block hyperbolic-elliptic coupler.
Definition coupler.hpp:74

pops::Coupler::Coupler
Coupler(const Model &model, const Geometry &geom, const BoxArray &ba, const BCRec &bcU, const BCRec &bcPhi, std::function< bool(Real, Real)> active={}, std::function< Real(Real, Real)> bz={})
Definition coupler.hpp:85

pops::Coupler::step
void step(MultiFab &U, Real dt)
Definition coupler.hpp:148

pops::Coupler::advance_ssprk3
void advance_ssprk3(MultiFab &U, Real dt)
Definition coupler.hpp:129

pops::Coupler::aux
const MultiFab & aux() const
Definition coupler.hpp:172

pops::Coupler::phi
MultiFab & phi()
Definition coupler.hpp:171

pops::Coupler::advance
void advance(MultiFab &U, Real dt)
Definition coupler.hpp:110

pops::Coupler::solve_fields
void solve_fields(const MultiFab &U)
Solve phi and derive aux = (phi, grad phi) for U WITHOUT advancing in time (useful to estimate the E ...
Definition coupler.hpp:169

pops::MultiFab
Field distributed over a level: decomposition (BoxArray) + distribution (DistributionMapping) + ncomp...
Definition multifab.hpp:33

pops::MultiFab::fab
Fab2D & fab(int li)
Local fab at index li (0 <= li < local_size()), for writing.
Definition multifab.hpp:67

pops::MultiFab::box
const Box2D & box(int li) const
VALID box of local fab li.
Definition multifab.hpp:71

pops::MultiFab::local_size
int local_size() const
Number of fabs OWNED by this rank (bound on local indices).
Definition multifab.hpp:65

comm.hpp
Parallel seam: minimal MPI abstraction (rank/size + collectives) with serial fallback.

pops::EllipticSolver
Definition elliptic_solver.hpp:30

coupling_policy.hpp
Hyperbolic-elliptic temporal coupling policies (compile-time tag types).

distribution_mapping.hpp
DistributionMapping: maps each box (by global index) to its owning MPI rank.

elliptic_problem.hpp
DESCRIPTIVE types of the elliptic stage: EllipticProblem (problem definition) and FieldPostProcess (f...

elliptic_rhs.hpp
Elliptic (Poisson) RIGHT-HAND-SIDE assemblers: single-model and N-species.

elliptic_solver.hpp
EllipticSolver concept: common contract for elliptic solvers at the MultiFab level (solve D phi = f),...

fab2d.hpp
Fab2D: single-grid data on a Box2D (in-house equivalent of AMReX's FArrayBox); Array4 / ConstArray4: ...

for_each.hpp
for_each_cell and reductions: the parallelism SEAM over the cells of a Box2D; sync_host / sync_device...

geometric_mg.hpp
GeometricMG: in-house geometric multigrid (V-cycle) for the elliptic operator, Gauss-Seidel smoother ...

geometry.hpp
Geometry: index-space (Box2D) <-> Cartesian physical-space mapping; PolarGeometry: SIBLING for a glob...

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::detail::coupler_eval_rhs
void coupler_eval_rhs(const MultiFab &state, MultiFab &rhs, const Model &model)
Assemble the single-model elliptic RHS: rhs = model.elliptic_rhs(U) on valid cells (delegated to Sing...
Definition coupler.hpp:53

pops::detail::fill_bz_box
void fill_bz_box(Fab2D &f, const Box2D &box, const Geometry &g, const Bz &bz)
Writes B_z(x, y) at component kAuxBaseComps on box box of fab f, sampling bz at the cell centers of g...
Definition aux_fill.hpp:42

pops::detail::coupler_grad_phi
void coupler_grad_phi(const MultiFab &phi, MultiFab &aux, Real cx, Real cy)
Set aux = (phi, d phi/dx, d phi/dy) by centered differences (factors cx, cy = 1/(2 dx),...
Definition coupler.hpp:64

pops
Definition amr_hierarchy.hpp:29

pops::kAuxBaseComps
constexpr int kAuxBaseComps
Definition state.hpp:147

pops::Real
double Real
Definition types.hpp:30

pops::n_ranks
int n_ranks()
Definition comm.hpp:139

pops::TimeTreatment::Explicit
@ Explicit

pops::aux_comps
POPS_HD constexpr int aux_comps()
Width of the aux channel a model CONSUMES.
Definition physical_model.hpp:67

pops::field_postprocess
void field_postprocess(const MultiFab &phi, MultiFab &out, Real cx, Real cy, FieldPostProcess spec)
Definition elliptic_problem.hpp:104

pops::fill_ghosts
void fill_ghosts(MultiFab &mf, const Box2D &domain, const BCRec &bc)
COMPLETE ghost filling: fill_boundary (interior + periodic, periodicity deduced from bc) THEN fill_ph...
Definition physical_bc.hpp:227

physical_bc.hpp
PHYSICAL boundary conditions at the domain edge (BCType, BCRec, fill_physical_bc, fill_ghosts).

reconstruction.hpp
Interface reconstruction policies: MUSCL limiters and WENO5-Z.

spatial_discretisation.hpp
"Spatial method" aggregate: reconstruction (limiter) + numerical flux in one named type.

spatial_operator.hpp
Cartesian spatial operator: assembles R(U, aux) = -div F + S over the cells of a level.

pops::BCRec
Boundary conditions for the FOUR faces of the domain (type + associated Dirichlet value).
Definition physical_bc.hpp:29

pops::FieldPostProcess
Definition elliptic_problem.hpp:73

pops::FieldPostProcess::GradSign::Plus
@ Plus

pops::Geometry
Cartesian geometry of a level: index domain + physical bounds [xlo, xhi] x [ylo, yhi].
Definition geometry.hpp:20

pops::Geometry::dy
POPS_HD Real dy() const
Grid spacing in y (= (yhi - ylo) / domain.ny()). POPS_HD.
Definition geometry.hpp:33

pops::Geometry::dx
POPS_HD Real dx() const
Grid spacing in x (= (xhi - xlo) / domain.nx()). POPS_HD.
Definition geometry.hpp:31

pops::Geometry::domain
Box2D domain
Definition geometry.hpp:21

pops::NoSlope
First-order reconstruction (piecewise constant): zero slope, 1 ghost.
Definition reconstruction.hpp:26

pops::PerStageCoupling
Tag: solves the elliptic problem at EVERY RK stage (aux follows the intermediate state,...
Definition coupling_policy.hpp:17

pops::RusanovFlux
RusanovFlux (local Lax-Friedrichs): robust flux, compatible with any minimal PhysicalModel.
Definition numerical_flux.hpp:55

pops::SSPRK2Step
Definition time_steppers.hpp:61

pops::SSPRK2Step::take_step
void take_step(RhsEval &&rhs, MultiFab &U, Real dt, Scratch &s) const
Definition time_steppers.hpp:73

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::SingleModelEllipticRhs
SINGLE-model RHS assembler: rhs(.,.,0) = model.elliptic_rhs(U) over the valid cells.
Definition elliptic_rhs.hpp:43

pops::TimePolicyTraits
Definition time_integrator.hpp:46

pops::TimePolicyTraits::Method
T Method
Definition time_integrator.hpp:47

time_integrator.hpp
Scheme tags (SSPRK2, SSPRK3, UserTimeIntegrator), TimeTreatment enum and per-block time policies: tem...

time_steppers.hpp
Time integrators as first-class OBJECTS with a take_step method: TimeStepper concept,...

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