PoPS/doxygen/system__block__store_8hpp_source.html

#pragma once


#include <pops/core/state/state.hpp>  // kAuxBaseComps (default aux channel of the Schur stage: B_z)

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

#include <pops/core/state/variables.hpp>   // VariableSet (role descriptor carried by each block)

#include <pops/mesh/index/box2d.hpp>       // Box2D

#include <pops/mesh/execution/for_each.hpp>  // device_fence (marshaling synchronizes the device before reading the host)

#include <pops/mesh/storage/multifab.hpp>  // MultiFab, Array4, ConstArray4


#include <functional>

#include <map>

#include <memory>

#include <stdexcept>

#include <string>

#include <vector>


namespace pops {


class CondensedSchurSourceStepper;

class PolarCondensedSchurSourceStepper;


class SystemBlockStore {

 public:

  using CellConvert = std::function<void(const double* in, double* out)>;


  struct BlockState {

    std::string name;

    MultiFab U;

    int ncomp;

    int substeps;    // static substeps (add_block)

    bool evolve;     // false = frozen species (fixed background, not advanced)

    int stride = 1;  // cadence: advance once every stride macro-steps

    double gamma;    // for the rest energy (4 var)

    std::function<void(MultiFab&, Real, int)> advance;   // (U, dt, n): n substeps of dt/n

    std::function<void(MultiFab&, MultiFab&)> rhs_into;  // R <- -div F + S (Poisson frozen)

    std::function<Real(const MultiFab&)> max_speed;      // max |wave speed| of the block

    std::function<void(const MultiFab&, MultiFab&)> add_poisson_rhs;  // += elliptic_rhs(U)

    // Descriptor of the conservative / primitive variables (names + physical ROLES) of the block.

    // The roles (provided by M::conservative_vars()) let inter-species couplings target a component

    // by its MEANING (momentum, energy) instead of a hard-coded index u[1]/u[3].

    VariableSet cons_vars, prim_vars;

    // POINTWISE cons <-> prim conversions OF THE BLOCK MODEL (one cell, ncomp doubles in/out).

    // Set at add time (install_block / push_dynamic) from the real model; empty -> identity (the

    // model exposes no conversion, e.g. pure scalar or .so generated before this work).

    // Consumed by set_primitive_state / get_primitive_state (init/diagnostic in primitive).

    CellConvert prim_to_cons, cons_to_prim;

    // Schur-CONDENSED SOURCE STAGE (OPT-IN, pops.Split(source=CondensedSchur), cf. set_source_stage).

    // nullptr (default) = no condensed source stage: the block advances EXACTLY as before

    // (bit-identical). Non null = after the hyperbolic transport, the step runs the standalone source stage

    // (CondensedSchurSourceStepper, #126) in place of the explicit / IMEX source. shared_ptr:

    // keeps BlockState MOVABLE (the stepper carries a GeometricMG, neither copyable nor trivially movable).

    std::shared_ptr<CondensedSchurSourceStepper> schur;

    // POLAR counterpart of the condensed source stage (ring (r, theta), Voie A step 2c). Exclusive with

    // schur above: set_source_stage builds ONE or the OTHER depending on the System geometry (polar

    // -> schur_polar, Cartesian -> schur). run_source_stage dispatches on whichever is non null. The

    // Cartesian path stays BIT-IDENTICAL (schur_polar == nullptr when the System is Cartesian).

    std::shared_ptr<PolarCondensedSchurSourceStepper> schur_polar;

    double schur_theta = 0.5;  // theta-scheme of the source stage (0.5 = Crank-Nicolson)

    // Component of the aux channel read as magnetic field Omega = B_z by the Schur stage (audit

    // wave 2: field TRANSPORTED in the ABI instead of the frozen literal kAuxBaseComps). Default =

    // kAuxBaseComps (canonical B_z channel), bit-identical; set_source_stage may redirect it.

    int schur_bz_comp = kAuxBaseComps;

    // GENERIC SOURCE STAGE (optional): a callable (U, dt) that advances IN PLACE the source stage of the

    // block. nullptr (default) = no generic source stage (BIT-IDENTICAL path). run_source_stage runs it

    // ONLY if no condensed Schur stage (schur / schur_polar) is wired, so it changes

    // NOTHING in the production Schur path. Used for generic splitting (pops.Strang on an arbitrary

    // source stage) and for the stepper time-order tests (non-commuting toy operators). Trailing

    // + nullptr default: the positional aggregate init of the other members stays unchanged.

    std::function<void(MultiFab&, Real)> source_step;

    // DISK TRANSPORT ADVANCES (work T5-PR3, OPT-IN). Empty (default) -> no disk routing:

    // the stepper advances via `advance` (full Cartesian path, BIT-IDENTICAL). Non empty, they MIMIC

    // `advance` (same RK / IMEX scheme, same limiter / flux) but dispatch the transport residual

    // to the disk operator, and are SELECTED only if the System is in Staircase mode (resp.

    // CutCell) AND a disk is set (set_disc_domain). Built at block add time (build_block)

    // AT THE SAME TIME as `advance`, they read the System mask / level set by pointer at step time

    // (stable address): the order add_block / set_disc_domain is irrelevant. Trailing + empty default:

    // the positional aggregate init of the other members stays unchanged.

    std::function<void(MultiFab&, Real, int)>

        advance_masked;  // residual via assemble_rhs_masked (Staircase)

    std::function<void(MultiFab&, Real, int)> advance_eb;  // residual via assemble_rhs_eb (CutCell)

    // dt_hotspot DIAGNOSTIC (ADC-182): (U, w, i, j) -> GLOBAL cell dominating the transport CFL bound

    // of the block + its speed w = max(wx, wy). ON DEMAND only (System::dt_hotspot):

    // never queried by step/step_cfl (hot path bit-identical). Trailing + empty default.

    std::function<void(const MultiFab&, Real&, int&, int&)> hotspot;

    // OPTIONAL STEP BOUNDS of the block (audit 2026-06, step_cfl work). EMPTY (default) -> the

    // stepper does not query them: STRICTLY historical step policy (transport only,

    // bit-identical). Set by set_block_dt_bounds when the model declares the traits

    // HasSourceFrequency / HasStabilityDt (cf. core/physical_model.hpp for the semantics).

    // Trailing + empty default: the positional aggregate init of the other members stays unchanged.

    std::function<Real(const MultiFab&)>

        source_frequency;  // max over cells of mu [1/s] (0 = no constraint)

    std::function<Real(const MultiFab&)>

        stability_dt;  // min over cells of the admissible step (0 = no constraint)

    // PROJECTION PONCTUELLE post-pas (ADC-177) : U <- project(U, aux) sur les cellules VALIDES,

    // appliquee par le stepper a la FIN de chaque macro-pas ENTIER (apres transport + etage source +

    // couplages ; jamais par etage RK). VIDE (defaut) -> jamais interrogee (cout nul, chemin

    // bit-identique). Trailing + defaut vide : l'init par agregat positionnel reste inchangee.

    std::function<void(MultiFab&)> project;

    // FLUX-ONLY residual R <- -div F(U) (NO default/composite source), Poisson frozen (ADC-425). The

    // SAME transport assembly as rhs_into evaluated on SourceFreeModel<Model> (zero source), so the

    // flux / ghost / geometry handling is bit-identical -- only the source is dropped (with

    // limiter='none'; the HLL wave-speed cache -- rejected on the aot/production backends compiled

    // Programs use -- is the only path where cached cell-center speeds differ from the per-face

    // reconstruction). Read by

    // System::block_neg_div_flux_into, which a compiled time Program's hyperbolic stage calls so a

    // Lie/Strang split assembles "flux but no source" (spec criterion 17). EMPTY (default) for paths

    // that do not build it (the host .so prototype loader); block_neg_div_flux_into fails loud then.

    // Trailing + empty default: the positional aggregate init of the other members stays unchanged.

    std::function<void(MultiFab&, MultiFab&)> rhs_flux_only;

    // NAMED elliptic-field RHS closures (ADC-428): field name -> (+= elliptic_field_rhs(U)). A model

    // declaring m.elliptic_field("phi2", rhs=...) carries here a SECOND Poisson right-hand side

    // (distinct from add_poisson_rhs, the default Poisson coupling), assembled the same way (host loop,

    // += per cell). EMPTY (default) -> no named elliptic field: bit-identical to the historical block.

    // The SystemFieldSolver gathers these per field (sum over blocks) into a SEPARATE elliptic solve

    // whose phi/grad are written to the field's OWN aux channel. Trailing + empty default: the

    // positional aggregate init of the other members stays unchanged.

    std::map<std::string, std::function<void(const MultiFab&, MultiFab&)>> named_poisson_rhs;

    // SOURCE-ONLY residual R <- S(U, aux) (the default/composite source, NO flux divergence), Poisson

    // frozen (ADC-430). The exact MIRROR of rhs_flux_only: together they split rhs_into (-div F + S).

    // SourceInto evaluates m.source per cell (the SAME source term assemble_rhs adds) with no

    // numerical-flux dispatch, so it is bit-identical to the source half of rhs_into. Read by

    // System::block_source_into, which a compiled time Program's source stage calls so a Lie/Strang

    // split assembles "the default source but no flux" (spec: rhs flux=False is source-only). EMPTY

    // (default) for paths that do not build it (the host .so prototype loader); block_source_into fails

    // loud then. Trailing + empty default: the positional aggregate init of the other members stays

    // unchanged.

    std::function<void(MultiFab&, MultiFab&)> source_only;

  };


  std::vector<BlockState> blocks;


  // --- access by NAME (0-based indexing, insertion order) ------------------------------------------


  BlockState& find(const std::string& name) {

    for (auto& s : blocks)

      if (s.name == name)

        return s;

    throw std::runtime_error("System : bloc inconnu '" + name + "'");

  }


  const BlockState& find(const std::string& name) const {

    for (auto& s : blocks)

      if (s.name == name)

        return s;

    throw std::runtime_error("System : bloc inconnu '" + name + "'");

  }


  int index(const std::string& name) const {

    for (std::size_t k = 0; k < blocks.size(); ++k)

      if (blocks[k].name == name)

        return static_cast<int>(k);

    throw std::runtime_error("System : bloc inconnu '" + name + "'");

  }


  int size() const { return static_cast<int>(blocks.size()); }


  std::vector<std::string> names() const {

    // Reads the UNIQUE block registry, populated by all add paths: a block loaded via

    // add_dynamic_block / add_compiled_block (.so) appears just like an add_block.

    std::vector<std::string> out;

    out.reserve(blocks.size());

    for (const auto& s : blocks)

      out.push_back(s.name);

    return out;

  }


  // --- state marshaling (host <-> MultiFab copy; device_fence included) ---------------------------


  std::vector<double> copy_comp0(const MultiFab& mf) const {

    device_fence();

    // MPI single-box: the box lives on the owner rank (rank 0). A rank without a box (local_size()==0)

    // has NO fab(0) -> return EMPTY rather than an OUT-OF-BOUNDS access (UB). Single-rank: local_size()

    // is always 1, behavior UNCHANGED. For the global multi-rank field: System::density_global.

    if (mf.local_size() == 0)

      return {};

    const ConstArray4 u = mf.fab(0).const_array();

    const Box2D v = mf.box(0);

    std::vector<double> out;

    out.reserve(static_cast<std::size_t>(v.nx()) * v.ny());

    for (int j = v.lo[1]; j <= v.hi[1]; ++j)

      for (int i = v.lo[0]; i <= v.hi[0]; ++i)

        out.push_back(u(i, j, 0));

    return out;

  }


  std::vector<double> copy_state(const MultiFab& mf, int ncomp) const {

    device_fence();

    // Rank without a box (MPI single-box, non-owner): return EMPTY (no fab(0)). Cf. copy_comp0;

    // the global multi-rank field goes through System::state_global (collective gather).

    if (mf.local_size() == 0)

      return {};

    const ConstArray4 u = mf.fab(0).const_array();

    const Box2D v = mf.box(0);

    std::vector<double> out;

    out.reserve(static_cast<std::size_t>(ncomp) * v.nx() * v.ny());

    for (int c = 0; c < ncomp; ++c)

      for (int j = v.lo[1]; j <= v.hi[1]; ++j)

        for (int i = v.lo[0]; i <= v.hi[0]; ++i)

          out.push_back(u(i, j, c));

    return out;

  }


  void write_state(MultiFab& mf, int ncomp, const std::vector<double>& in) {

    // Rank without a box (MPI single-box, non-owner): NO-OP (no fab(0) to write). Lets

    // sim.set_state / sim.restart be called on ALL ranks with the GLOBAL field: only the

    // owner rank (rank 0, box = full domain) writes, the others do nothing. Single-rank:

    // local_size()==1, validation + write UNCHANGED (bit-identical).

    if (mf.local_size() == 0)

      return;

    const Box2D v = mf.box(0);

    const std::size_t need = static_cast<std::size_t>(ncomp) * v.nx() * v.ny();

    if (in.size() != need)

      throw std::runtime_error("System::set_state : taille != ncomp*n*n");

    Array4 u = mf.fab(0).array();

    std::size_t k = 0;

    for (int c = 0; c < ncomp; ++c)

      for (int j = v.lo[1]; j <= v.hi[1]; ++j)

        for (int i = v.lo[0]; i <= v.hi[0]; ++i)

          u(i, j, c) = in[k++];

  }


};


}  // namespace pops

box2d.hpp
Box2D: the integer index space of a 2D cell-centered Cartesian grid.

pops::Fab2D::const_array
ConstArray4 const_array() const
READ handle (POD device-copyable) over this Fab. Valid as long as the Fab lives.
Definition fab2d.hpp:96

pops::Fab2D::array
Array4 array()
WRITE handle (POD device-copyable) over this Fab. Valid as long as the Fab lives.
Definition fab2d.hpp:91

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

pops::SystemBlockStore
ORDERED registry of the System blocks + state marshaling helpers.
Definition system_block_store.hpp:58

pops::SystemBlockStore::write_state
void write_state(MultiFab &mf, int ncomp, const std::vector< double > &in)
Writes the ncomp components into fab(0) from a component-major buffer (same layout as copy_state).
Definition system_block_store.hpp:254

pops::SystemBlockStore::find
const BlockState & find(const std::string &name) const
Reference to block name (for reading).
Definition system_block_store.hpp:188

pops::SystemBlockStore::copy_comp0
std::vector< double > copy_comp0(const MultiFab &mf) const
Copies component 0 of fab(0) (density) in row-major (j slow, i fast).
Definition system_block_store.hpp:218

pops::SystemBlockStore::find
BlockState & find(const std::string &name)
Reference to block name (for writing).
Definition system_block_store.hpp:181

pops::SystemBlockStore::index
int index(const std::string &name) const
0-based index of block name (insertion order).
Definition system_block_store.hpp:195

pops::SystemBlockStore::blocks
std::vector< BlockState > blocks
ORDERED registry of the blocks (UNIQUE source of truth).
Definition system_block_store.hpp:177

pops::SystemBlockStore::copy_state
std::vector< double > copy_state(const MultiFab &mf, int ncomp) const
Copies the ncomp components of fab(0) in component-major layout (c slow, then j, then i).
Definition system_block_store.hpp:235

pops::SystemBlockStore::names
std::vector< std::string > names() const
Block names, in insertion order (UNIQUE source: all add paths appear here).
Definition system_block_store.hpp:205

pops::SystemBlockStore::CellConvert
std::function< void(const double *in, double *out)> CellConvert
Type-erase of the POINTWISE (one cell) cons <-> prim conversion of a block: in/out are arrays of ncom...
Definition system_block_store.hpp:63

pops::SystemBlockStore::size
int size() const
Number of registered blocks.
Definition system_block_store.hpp:203

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

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

pops
Definition amr_hierarchy.hpp:29

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

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

pops::device_fence
void device_fence()
Device barrier: waits for in-flight kernels to finish before a HOST access to unified memory.
Definition kokkos_env.hpp:43

state.hpp
Pointwise types of the physics layer: StateVec<N> (conserved state) and Aux (auxiliary fields from th...

pops::Array4
WRITE POD handle (raw pointer + strides) over a Fab2D buffer, indexed by (i, j, c) IN GLOBAL INDICES ...
Definition fab2d.hpp:29

pops::Box2D
2D integer index space, cell-centered.
Definition box2d.hpp:37

pops::Box2D::nx
POPS_HD int nx() const
Width (direction 0). POPS_HD (called from Geometry::dx() in a device kernel).
Definition box2d.hpp:50

pops::Box2D::hi
int hi[2]
Definition box2d.hpp:39

pops::Box2D::ny
POPS_HD int ny() const
Height (direction 1). POPS_HD (called from Geometry::dy() in a device kernel).
Definition box2d.hpp:52

pops::Box2D::lo
int lo[2]
Definition box2d.hpp:38

pops::ConstArray4
READ-only handle (const counterpart of Array4): same layout and same contract (POD device-copyable,...
Definition fab2d.hpp:44

pops::SystemBlockStore::BlockState
Compiled closures frozen at block add time (composite model + spatial scheme + time).
Definition system_block_store.hpp:71

pops::SystemBlockStore::BlockState::project
std::function< void(MultiFab &)> project
Definition system_block_store.hpp:143

pops::SystemBlockStore::BlockState::source_frequency
std::function< Real(const MultiFab &)> source_frequency
Definition system_block_store.hpp:136

pops::SystemBlockStore::BlockState::stride
int stride
Definition system_block_store.hpp:77

pops::SystemBlockStore::BlockState::prim_to_cons
CellConvert prim_to_cons
Definition system_block_store.hpp:91

pops::SystemBlockStore::BlockState::cons_vars
VariableSet cons_vars
Definition system_block_store.hpp:86

pops::SystemBlockStore::BlockState::source_step
std::function< void(MultiFab &, Real)> source_step
Definition system_block_store.hpp:114

pops::SystemBlockStore::BlockState::schur_theta
double schur_theta
Definition system_block_store.hpp:103

pops::SystemBlockStore::BlockState::U
MultiFab U
Definition system_block_store.hpp:73

pops::SystemBlockStore::BlockState::advance_eb
std::function< void(MultiFab &, Real, int)> advance_eb
Definition system_block_store.hpp:125

pops::SystemBlockStore::BlockState::stability_dt
std::function< Real(const MultiFab &)> stability_dt
Definition system_block_store.hpp:138

pops::SystemBlockStore::BlockState::named_poisson_rhs
std::map< std::string, std::function< void(const MultiFab &, MultiFab &)> > named_poisson_rhs
Definition system_block_store.hpp:162

pops::SystemBlockStore::BlockState::schur_polar
std::shared_ptr< PolarCondensedSchurSourceStepper > schur_polar
Definition system_block_store.hpp:102

pops::SystemBlockStore::BlockState::rhs_flux_only
std::function< void(MultiFab &, MultiFab &)> rhs_flux_only
Definition system_block_store.hpp:154

pops::SystemBlockStore::BlockState::source_only
std::function< void(MultiFab &, MultiFab &)> source_only
Definition system_block_store.hpp:172

pops::SystemBlockStore::BlockState::advance
std::function< void(MultiFab &, Real, int)> advance
Definition system_block_store.hpp:79

pops::SystemBlockStore::BlockState::prim_vars
VariableSet prim_vars
Definition system_block_store.hpp:86

pops::SystemBlockStore::BlockState::schur_bz_comp
int schur_bz_comp
Definition system_block_store.hpp:107

pops::SystemBlockStore::BlockState::substeps
int substeps
Definition system_block_store.hpp:75

pops::SystemBlockStore::BlockState::evolve
bool evolve
Definition system_block_store.hpp:76

pops::SystemBlockStore::BlockState::name
std::string name
Definition system_block_store.hpp:72

pops::SystemBlockStore::BlockState::rhs_into
std::function< void(MultiFab &, MultiFab &)> rhs_into
Definition system_block_store.hpp:80

pops::SystemBlockStore::BlockState::advance_masked
std::function< void(MultiFab &, Real, int)> advance_masked
Definition system_block_store.hpp:124

pops::SystemBlockStore::BlockState::cons_to_prim
CellConvert cons_to_prim
Definition system_block_store.hpp:91

pops::SystemBlockStore::BlockState::add_poisson_rhs
std::function< void(const MultiFab &, MultiFab &)> add_poisson_rhs
Definition system_block_store.hpp:82

pops::SystemBlockStore::BlockState::schur
std::shared_ptr< CondensedSchurSourceStepper > schur
Definition system_block_store.hpp:97

pops::SystemBlockStore::BlockState::hotspot
std::function< void(const MultiFab &, Real &, int &, int &)> hotspot
Definition system_block_store.hpp:129

pops::SystemBlockStore::BlockState::gamma
double gamma
Definition system_block_store.hpp:78

pops::SystemBlockStore::BlockState::max_speed
std::function< Real(const MultiFab &)> max_speed
Definition system_block_store.hpp:81

pops::SystemBlockStore::BlockState::ncomp
int ncomp
Definition system_block_store.hpp:74

pops::VariableSet
A model's variable set: kind (cons/prim), names, size, canonical roles (optional, parallel to names; ...
Definition variables.hpp:58

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

variables.hpp
Descriptor of a model's variables (Vars).