PoPS/doxygen/amr__runtime_8hpp_source.html

#pragma once


#include <pops/amr/regridding/regrid.hpp>   // tag_cells, grow_tags (per-block tags + phi for the union regrid)

#include <pops/amr/tagging/tag_box.hpp>  // TagBox, tag_union (cell-by-cell OR of the tags of all blocks)

#include <pops/core/state/state.hpp>   // kAuxBaseComps

#include <pops/core/state/variables.hpp>  // VariableSet, VariableRole, role_from_name (role -> component of coupled sources)

#include <pops/coupling/amr/amr_coupler_mp.hpp>  // detail::coupler_inject_aux_mb (aux injection coarse->fine)

#include <pops/coupling/amr/amr_regrid_coupler.hpp>  // regrid_compute_fine_layout + regrid_field_on_layout (split bricks)

#include <pops/coupling/system/amr_system_coupler.hpp>  // detail::same_layout_or_throw (shared-layout guard)

#include <pops/coupling/base/aux_fill.hpp>            // detail::derive_aux_bc (BC of the aux channel)

#include <pops/coupling/source/coupled_source_program.hpp>  // CoupledSourceKernel + CsProgram (flat ABI, P5 bytecode)

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

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

#include <pops/numerics/time/amr/reflux/amr_reflux_mf.hpp>  // AmrLevelMP, mf_average_down_mb

#include <pops/numerics/time/integrators/implicit_stepper.hpp>  // NewtonReport (OPT-IN IMEX diagnostics, aggregated per block)

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

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

#include <pops/mesh/layout/patch_box.hpp>  // PatchBox: index-space signature of a fine patch (patch_boxes())

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

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

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

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

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

#include <pops/parallel/comm.hpp>  // n_ranks() / comm_active(): MPI message+reduction counts (Spec 5 criterion 43)

#include <pops/runtime/program/profiler.hpp>  // Profiler / ProfileScope: AMR phase timings (Spec 5 criterion 43, ADC-479)


#include <algorithm>  // std::max (substeps/stride-aware CFL step)

#include <chrono>  // AmrPhaseScope wall-clock timing (Spec 5 criterion 43)

#include <cmath>      // std::isfinite (reject a degenerate dt)

#include <cstddef>

#include <functional>

#include <limits>  // std::numeric_limits (initial dt = +inf, min over the blocks)

#include <map>  // named_aux_: model-named aux fields (comp -> coarse field), re-applied each solve

#include <memory>

#include <stdexcept>

#include <string>

#include <utility>

#include <vector>


namespace pops {


struct AmrRuntimeBlock {

  std::string name;

  int ncomp = 1;

  double gamma = 1.4;

  int substeps = 1;

  int stride = 1;

  int aux_ncomp = kAuxBaseComps;


  VariableSet cons_vars;


  std::shared_ptr<std::vector<AmrLevelMP>> levels;


  std::function<void(std::vector<AmrLevelMP>&, const Box2D&, Real, Periodicity, bool)> advance;


  bool imex = false;


  std::function<void(std::vector<AmrLevelMP>&, const Box2D&, Real, Periodicity, bool)> imex_advance;


  std::function<void(std::vector<AmrLevelMP>&)> project_per_level;


  bool newton_diagnostics = false;

  std::shared_ptr<NewtonReport> newton_report;


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


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


  std::function<Real(const MultiFab&, const MultiFab&)> max_speed;


  std::function<Real(const MultiFab&, const MultiFab&)> source_frequency;

  std::function<Real(const MultiFab&, const MultiFab&)> stability_dt;


  std::function<Real()> mass;


  std::function<std::vector<double>()> density;


  std::function<std::vector<double>(const MultiFab&)> potential;

};


class AmrRuntime {

 public:


  AmrRuntime(const Geometry& geom, const BoxArray& ba_coarse, const BCRec& bcPhi,

             std::vector<AmrRuntimeBlock> blocks, Periodicity base_per = Periodicity{true, true},

             bool replicated_coarse = true, std::function<bool(Real, Real)> active = {})

      : geom_(geom),

        dom_(geom.domain),

        base_per_(base_per),

        bcPhi_(bcPhi),

        aux_bc_(detail::derive_aux_bc(bcPhi)),

        replicated_coarse_(replicated_coarse),

        mg_(geom, ba_coarse, bcPhi, active, replicated_coarse),

        ba_coarse_(ba_coarse),

        wall_active_(std::move(active)),  // copy already consumed by mg_ (earlier in decl order)

        blocks_(std::move(blocks)) {

    if (blocks_.empty())

      throw std::runtime_error("AmrRuntime : at least one block required");

    for (const auto& b : blocks_)

      if (!b.levels || b.levels->empty())

        throw std::runtime_error(

            "AmrRuntime : each block must carry at least one level "

            "(coarse) on the shared layout");

    nlev_ = static_cast<int>(blocks_[0].levels->size());


    // EXACT layout consistency between blocks (the aux is shared per level): same number of levels,

    // and per level same BoxArray (boxes AND order), same DistributionMapping, same dx/dy. SAME guard

    // as AmrSystemCoupler (detail::same_layout_or_throw): all blocks live on ALL patches of the

    // UNIQUE shared hierarchy. A single block matches itself trivially (the loop over the other blocks

    // is empty).

    {

      std::vector<std::vector<AmrLevelMP>> ref;

      ref.reserve(blocks_.size());

      for (const auto& b : blocks_)

        ref.push_back(*b.levels);

      detail::same_layout_or_throw(ref);

    }


    // Width of the SHARED aux channel: max of the blocks' aux_comps (>= kAuxBaseComps). Counterpart of

    // AmrSystemCoupler::system_aux_comps: a block reading an extra field (B_z, T_e) has the room at

    // each level, a base block ignores the extra components. PR1 does not POPULATE multi-block B_z (no

    // bz_ here), but we size the channel to the widest anyway so that load_aux<aux_comps<Model>> never

    // reads out of bounds. Without an extra-field block -> kAuxBaseComps (3) -> allocation strictly

    // identical to the base case.

    aux_ncomp_ = kAuxBaseComps;

    for (const auto& b : blocks_)

      if (b.aux_ncomp > aux_ncomp_)

        aux_ncomp_ = b.aux_ncomp;


    // SHARED aux: one MultiFab (phi, grad phi) per level, on the common grid. Sized once -> stable

    // addresses for the blocks' aux pointers. The shared layout is that of block 0

    // (same_layout_or_throw guard: identical for all).

    aux_.resize(nlev_);

    const auto& L0 = *blocks_[0].levels;

    for (int k = 0; k < nlev_; ++k)

      aux_[k] = MultiFab(L0[k].U.box_array(), L0[k].U.dmap(), aux_ncomp_, 1);

    for (auto& b : blocks_)

      for (int k = 0; k < nlev_; ++k)

        (*b.levels)[k].aux = &aux_[k];


    // Tag predicates of the union regrid: one empty slot per block (set_block_tag_predicate fills

    // them). Empty by default -> no tag -> frozen hierarchy (regrid is not called anyway as long as

    // set_regrid has not activated regrid_every_ > 0).

    block_tag_.resize(blocks_.size());

  }


  int nlev() const { return nlev_; }

  std::size_t n_blocks() const { return blocks_.size(); }


  const VariableSet& block_cons_vars(std::size_t b) const {

    if (b >= blocks_.size())

      throw std::runtime_error("AmrRuntime::block_cons_vars : block index out of bounds");

    return blocks_[b].cons_vars;

  }


  std::size_t n_coupled_sources() const { return coupled_sources_.size(); }

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

  // System Poisson right-hand side after the last solve_fields: f = Sum_b elliptic_rhs_b(U_b) on the

  // shared coarse. Exposed to check the CO-LOCATED SUM (PR1 test); same grid as the coarse (the

  // blocks' contributions are accumulated there at the same cells).

  MultiFab& poisson_rhs() { return mg_.rhs(); }

  const MultiFab& aux(int k) const { return aux_[k]; }

  std::vector<AmrLevelMP>& levels(std::size_t b) { return *blocks_[b].levels; }

  Real mass(std::size_t b) const { return blocks_[b].mass(); }

  std::vector<double> density(std::size_t b) const { return blocks_[b].density(); }

  int solve_count() const { return solve_count_; }

  int regrid_count() const { return regrid_count_; }


  using TagPredicate = std::function<bool(const ConstArray4&, int, int)>;


  void set_regrid(int every, int grow = 2, int margin = 2) {

    if (every < 0)

      throw std::runtime_error("AmrRuntime::set_regrid : regrid_every >= 0");

    regrid_every_ = every;

    regrid_grow_ = grow;

    regrid_margin_ = margin;

  }


  void set_block_tag_predicate(std::size_t b, TagPredicate crit) {

    if (b >= blocks_.size())

      throw std::runtime_error("AmrRuntime::set_block_tag_predicate : block index out of bounds");

    block_tag_[b] = std::move(crit);

  }


  void set_phi_tag_predicate(TagPredicate crit) { phi_tag_ = std::move(crit); }


  void set_named_aux(int comp, std::vector<Real> field) {

    named_aux_[comp] = std::move(field);

    if (!aux_.empty())

      apply_named_aux();  // reflect immediately if the hierarchy already exists

  }


  void set_named_aux_bc(int comp, AuxHaloPolicy policy) { named_aux_bc_[comp] = policy; }


  void register_named_field(const std::string& field, int phi_comp, int gx_comp, int gy_comp) {

    NamedField nf;

    nf.phi_comp = phi_comp;

    nf.gx_comp = gx_comp;

    nf.gy_comp = gy_comp;

    named_fields_[field] = std::move(nf);  // solver built lazily by ensure_named_elliptic

  }


  void set_block_named_elliptic_rhs(std::size_t b, const std::string& field,

                                    std::function<void(const MultiFab&, MultiFab&)> rhs) {

    if (b >= blocks_.size())

      throw std::runtime_error(

          "AmrRuntime::set_block_named_elliptic_rhs : block index out of bounds");

    blocks_[b].named_elliptic_rhs[field] = std::move(rhs);

  }


  std::size_t n_named_fields() const { return named_fields_.size(); }


  bool has_named_field(const std::string& field) const {

    return named_fields_.find(field) != named_fields_.end();

  }


  std::vector<double> named_field_values(const std::string& field) {

    auto it = named_fields_.find(field);

    if (it == named_fields_.end())

      throw std::runtime_error("AmrRuntime::named_field_values : unknown named elliptic field '" +

                               field + "' (register it via m.elliptic_field + the compiled block)");

    solve_fields();  // up-to-date phi (counterpart of potential()); named solve runs inside

    return coarse_aux_component(it->second.phi_comp);

  }


  void add_coupled_source(const std::vector<std::string>& in_blocks,

                          const std::vector<std::string>& in_roles,

                          const std::vector<double>& consts,

                          const std::vector<std::string>& out_blocks,

                          const std::vector<std::string>& out_roles,

                          const std::vector<int>& prog_ops, const std::vector<int>& prog_args,

                          const std::vector<int>& prog_lens) {

    const int n_in = static_cast<int>(in_blocks.size());

    const int n_const = static_cast<int>(consts.size());

    const int n_terms = static_cast<int>(out_blocks.size());

    // --- form validation (before any step, EXPLICIT errors); mirror of System::add_coupled_source.

    if (n_terms == 0)

      throw std::runtime_error(

          "AmrRuntime::add_coupled_source : no source term (out_blocks empty)");

    if (static_cast<int>(in_roles.size()) != n_in)

      throw std::runtime_error(

          "AmrRuntime::add_coupled_source : in_blocks / in_roles of different sizes");

    if (static_cast<int>(out_roles.size()) != n_terms ||

        static_cast<int>(prog_lens.size()) != n_terms)

      throw std::runtime_error(

          "AmrRuntime::add_coupled_source : out_blocks / out_roles / prog_lens of "

          "different sizes");

    if (prog_ops.size() != prog_args.size())

      throw std::runtime_error(

          "AmrRuntime::add_coupled_source : prog_ops / prog_args of different sizes");

    if (n_in + n_const > kCsMaxReg)

      throw std::runtime_error(

          "AmrRuntime::add_coupled_source : too many registers (inputs + constants > " +

          std::to_string(kCsMaxReg) + ")");

    if (n_terms > kCsMaxTerms)

      throw std::runtime_error("AmrRuntime::add_coupled_source : too many source terms (> " +

                               std::to_string(kCsMaxTerms) + ")");

    // Resolves (block, role) -> (block index, component) by the block CONSERVATIVE descriptor, like

    // System (#181). An unknown block throws immediately; an unknown (non-canonical) role too.

    auto resolve = [&](const std::string& block, const std::string& role) -> std::pair<int, int> {

      const int b = block_index(block);

      if (b < 0)

        throw std::runtime_error("AmrRuntime::add_coupled_source : no block named '" + block + "'");

      // STRICT (no silent fallback; mirror of System::add_coupled_source #181): a DSL coupled source

      // targets a (block, role) EXPLICITLY requested by the user. The role is addressed BY NAME: a

      // canonical role name OR a user-defined role label (index_of(string), ADC-292). If the block does

      // NOT expose this role, a fallback to component 0 would apply the source to the wrong field

      // SILENTLY (the false-positive identified at the Lot E review). We throw, listing what the block

      // exposes.

      const VariableSet& vs = blocks_[static_cast<std::size_t>(b)].cons_vars;

      const int comp = vs.index_of(role);

      if (comp < 0)

        throw std::runtime_error(

            "AmrRuntime::add_coupled_source : block '" + block + "' does not expose role '" + role +

            "' (roles: " + (vs.roles.empty() ? std::string("<none>") : roles_csv(vs)) +

            ", no silent fallback to component 0)");

      return {b, comp};

    };

    // Inputs: (block, component) read per cell. Captured by INDEX -> we rebuild the Array4 at EACH

    // application (the fabs live in the level stack, repointed per level in the splitting).

    std::vector<CsRef> ins(static_cast<std::size_t>(n_in));

    for (int c = 0; c < n_in; ++c) {

      auto [b, comp] =

          resolve(in_blocks[static_cast<std::size_t>(c)], in_roles[static_cast<std::size_t>(c)]);

      ins[static_cast<std::size_t>(c)] = {b, comp, CsProgram{}};

    }

    std::vector<CsRef> outs(static_cast<std::size_t>(n_terms));

    int off = 0;

    for (int t = 0; t < n_terms; ++t) {

      auto [b, comp] =

          resolve(out_blocks[static_cast<std::size_t>(t)], out_roles[static_cast<std::size_t>(t)]);

      const int len = prog_lens[static_cast<std::size_t>(t)];

      if (len < 0 || len > kCsMaxProg)

        throw std::runtime_error("AmrRuntime::add_coupled_source : program of term " +

                                 std::to_string(t) + " too long (> " + std::to_string(kCsMaxProg) +

                                 ")");

      if (off + len > static_cast<int>(prog_ops.size()))

        throw std::runtime_error(

            "AmrRuntime::add_coupled_source : prog_lens inconsistent with prog_ops");

      CsProgram pg;

      pg.len = len;

      for (int k = 0; k < len; ++k) {

        const int opc = prog_ops[static_cast<std::size_t>(off + k)];

        const int a = prog_args[static_cast<std::size_t>(off + k)];

        if (opc < 0 || opc > static_cast<int>(CsOp::Sqrt))

          throw std::runtime_error("AmrRuntime::add_coupled_source : invalid opcode");

        if (opc == static_cast<int>(CsOp::PushReg) && (a < 0 || a >= n_in + n_const))

          throw std::runtime_error(

              "AmrRuntime::add_coupled_source : register out of bounds in the program");

        pg.op[k] = opc;

        pg.arg[k] = a;

      }

      outs[static_cast<std::size_t>(t)] = {b, comp, pg};

      off += len;

    }

    std::vector<Real> kconsts(consts.begin(), consts.end());

    coupled_sources_.push_back(CoupledSourceSpec{std::move(ins), std::move(outs),

                                                 std::move(kconsts), n_in, n_const, n_terms});

  }


  void coupled_source_step(Real dt) {

    if (coupled_sources_.empty())

      return;        // opt-in: no source -> bit-identical path

    solve_fields();  // aux per level up to date (a term may read phi/grad via a future input)

    for (const auto& cs : coupled_sources_) {

      // PER-LEVEL application: at each level k, the blocks share EXACTLY the same layout

      // (same_layout_or_throw guard), so same local_size() and same local indexing -> we iterate in

      // parallel over the local fabs. local_size()==0 on a rank without a box -> empty loop (MPI-safe).

      for (int k = 0; k < nlev_; ++k) {

        const int sref = cs.n_in > 0 ? cs.ins[0].block : cs.outs[0].block;

        MultiFab& Uref = (*blocks_[static_cast<std::size_t>(sref)].levels)[k].U;

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

          CoupledSourceKernel kern;

          kern.dt = dt;

          kern.n_in = cs.n_in;

          kern.n_const = cs.n_const;

          kern.n_terms = cs.n_terms;

          for (int c = 0; c < cs.n_in; ++c) {

            kern.in[c] =

                (*blocks_[static_cast<std::size_t>(cs.ins[static_cast<std::size_t>(c)].block)]

                      .levels)[k]

                    .U.fab(li)

                    .array();

            kern.in_comp[c] = cs.ins[static_cast<std::size_t>(c)].comp;

          }

          for (int c = 0; c < cs.n_const; ++c)

            kern.consts[c] = cs.kconsts[static_cast<std::size_t>(c)];

          for (int t = 0; t < cs.n_terms; ++t) {

            kern.out[t] =

                (*blocks_[static_cast<std::size_t>(cs.outs[static_cast<std::size_t>(t)].block)]

                      .levels)[k]

                    .U.fab(li)

                    .array();

            kern.out_comp[t] = cs.outs[static_cast<std::size_t>(t)].comp;

            kern.prog[t] = cs.outs[static_cast<std::size_t>(t)].prog;

          }

          for_each_cell(Uref.box(li),

                        kern);  // NAMED functor (device-clean), additive forward-Euler

        }

      }

      // Restore the consistency of the covered coarse cells (cf. COVERAGE INVARIANT above).

      for (auto& b : blocks_)

        for (int k = nlev_ - 1; k >= 1; --k)

          mf_average_down_mb((*b.levels)[k].U, (*b.levels)[k - 1].U);

    }

  }


  void solve_fields() {

    ++solve_count_;

    // 1. average_down per block (fine -> coarse) over the whole hierarchy. AMR PROFILING (Spec 5

    // criterion 43): time the restriction cascade into the "average_down" scope + bump its per-solve

    // count. The scope is per-solve_fields (NOT per-cell), so a profiled run pays one clock pair here;

    // an unprofiled run constructs nothing (profiler_ null or disabled). See profile_amr_scope below.

    {

      auto _ad = profile_amr_scope("average_down");

      if (profiler_ != nullptr)

        profiler_->count("average_down");

      for (auto& b : blocks_) {

        auto& L = *b.levels;

        for (int k = nlev_ - 1; k >= 1; --k)

          mf_average_down_mb(L[k].U, L[k - 1].U);

      }

    }


    // 2. SUMMED and CO-LOCATED system RHS: f = Sum_b elliptic_rhs_b(U_b) on the coarse. We reset to

    // zero then each block ACCUMULATES (+=) its contribution on the SAME cells of the shared coarse

    // (mg_.rhs() shares the coarse layout).

    mg_.rhs().set_val(Real(0));

    for (auto& b : blocks_)

      b.add_elliptic_rhs((*b.levels)[0].U, mg_.rhs());

    mg_.solve();


    // 3. coarse aux = (phi, grad phi) via the SAME clean path as AmrSystemCoupler: fill the ghosts of

    // phi according to bcPhi_, field_postprocess (phi + grad), fill the ghosts of aux according to

    // aux_bc_ (derived from bcPhi_). Handles the non-periodic case (Foextrap).

    fill_ghosts_profiled(mg_.phi(), dom_, bcPhi_);

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

    field_postprocess(mg_.phi(), aux_[0], cx, cy,

                      FieldPostProcess{FieldPostProcess::GradSign::Plus, true});

    // 3b. model-NAMED aux (ADC-291): re-apply the static named fields onto the coarse valid cells

    // BEFORE fill_ghosts (so their ghosts are filled) and the injection (so they reach every level).

    // No-op when no named field was set; field_postprocess wrote only comps 0..2, so this never clobbers

    // phi/grad. This is what makes named aux survive a regrid (regrid re-solves -> re-applies).

    apply_named_aux();

    fill_ghosts_profiled(aux_[0], dom_, aux_bc_);

    apply_named_aux_bc();  // ADC-369: per-field halo override on the coarse physical ghosts (after the

                           // shared fill, before injection); no-op when no policy declared.

    // 4. coarse->fine injection of the aux (parent replicated only at level 1 if coarse replicated).

    for (int k = 1; k < nlev_; ++k)

      detail::coupler_inject_aux_mb(aux_[k - 1], aux_[k],

                                    /*replicated_parent=*/(k == 1) && replicated_coarse_);


    // 5. NAMED multi-elliptic fields (ADC-428): a SECOND elliptic solve per user-named field, written to

    // the field's OWN aux components and injected to the fine levels. No-op when none is registered ->

    // the default-Poisson trajectory above is strictly bit-identical.

    solve_named_fields();

  }


  void solve_named_fields() {

    if (named_fields_.empty())

      return;

    const Real dx = geom_.dx(), dy = geom_.dy();

    for (auto& [field, nf] : named_fields_) {

      if (nf.phi_comp < 0 || nf.phi_comp >= aux_ncomp_)

        throw std::runtime_error("AmrRuntime : named elliptic field '" + field +

                                 "' aux component out of the channel width (add the block that "

                                 "declares its aux fields)");

      ensure_named_elliptic(nf);

      // SUMMED + CO-LOCATED RHS on the coarse: f = Sum_b named_elliptic_rhs_b[field](U_b), exactly like

      // the default Poisson (mg_.rhs()), but reading the per-field block closures. A field with no

      // contributing block solves a zero RHS -> reject loud (mirror of the System named path).

      MultiFab& rhs = nf.mg->rhs();

      rhs.set_val(Real(0));

      bool any = false;

      for (auto& b : blocks_) {

        auto it = b.named_elliptic_rhs.find(field);

        if (it == b.named_elliptic_rhs.end() || !it->second)

          continue;

        it->second((*b.levels)[0].U, rhs);

        any = true;

      }

      if (!any)

        throw std::runtime_error(

            "AmrRuntime : named elliptic field '" + field +

            "' has no contributing block (declare m.elliptic_field on the block model)");

      nf.mg->solve();

      device_fence();  // CRITICAL: the V-cycle must finish before phi is read (same invariant as mg_)

      // Write phi (+ centered grad) into the field's OWN aux components on the coarse valid cells. The

      // default field_postprocess hardcodes comps 0..2, so we write the named comps with a dedicated

      // loop (mirror of SystemFieldSolver::solve_named_field_from_state). Per-local-fab (MPI-safe).

      MultiFab& phi_mf = nf.mg->phi();

      const int cphi = nf.phi_comp, cgx = nf.gx_comp, cgy = nf.gy_comp;

      const bool grad = (cgx >= 0 && cgx < aux_ncomp_ && cgy >= 0 && cgy < aux_ncomp_);

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

        const ConstArray4 p = phi_mf.fab(li).const_array();

        Array4 a = aux_[0].fab(li).array();

        const Box2D v = aux_[0].box(li);

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

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

            a(i, j, cphi) = p(i, j);

            if (grad) {

              a(i, j, cgx) = (p(i + 1, j) - p(i - 1, j)) / (2 * dx);

              a(i, j, cgy) = (p(i, j + 1) - p(i, j - 1)) / (2 * dy);

            }

          }

      }

    }

    // Ghost-fill the named components (shared aux fill: same routing as the default) + per-field halo

    // override (ADC-369), then inject coarse->fine so the named field reaches every level. We re-fill the

    // WHOLE aux: the default comps 0..2 were just written by the Poisson block, so their valid cells are

    // unchanged -- only ghosts are refreshed (idempotent). Cheap (one extra fill per solve_fields when a

    // named field exists; none otherwise).

    fill_ghosts_profiled(aux_[0], dom_, aux_bc_);

    apply_named_aux_bc();

    for (int k = 1; k < nlev_; ++k)

      detail::coupler_inject_aux_mb(aux_[k - 1], aux_[k],

                                    /*replicated_parent=*/(k == 1) && replicated_coarse_);

  }


  void regrid() {

    if (nlev_ < 2)

      return;  // 2 levels required (D5): nothing to re-grid in single-level

    const int fk = nlev_ - 1, pk = fk - 1;  // fine + its parent (pk == 0 in v1 with 2 levels)


    // AMR PROFILING (Spec 5 criterion 43): time the WHOLE regrid attempt (tag + cluster + prolong +

    // re-solve) into the "regrid" scope. RAII -> the scope covers EVERY early-return path below (empty

    // tags / nothing to refine), so the timing reflects the real regrid cost. The per-run "regrid"

    // COUNT is bumped at the tail (++regrid_count_) only when a regrid actually completed -- a no-op

    // attempt times itself but does not inflate the count. Null/disabled profiler -> no scope object.

    auto _rg = profile_amr_scope("regrid");


    // (R0) PRECONDITION: fields up to date (aux per level, for the |grad phi| criterion). The per-block

    // mass snapshot is NOT needed by the engine (conservation is checked test-side V1).

    solve_fields();


    // (R1)+(R2) PER-BLOCK TAGS (on the block U at the parent level) + PHI TAGS (on the shared aux).

    const int PNX = dom_.nx() << pk, PNY = dom_.ny() << pk;

    const Box2D pdom = Box2D::from_extents(PNX, PNY);

    std::vector<TagBox> parts;

    parts.reserve(blocks_.size() + 1);

    for (std::size_t b = 0; b < blocks_.size(); ++b) {

      const TagPredicate& crit = block_tag_[b];

      if (!crit)

        continue;  // block without a criterion: tags nothing on its side (re-gridded as background)

      parts.push_back(tag_cells((*blocks_[b].levels)[pk].U, pdom, crit));

    }

    if (phi_tag_)

      parts.push_back(tag_cells(aux_[pk], pdom, phi_tag_));

    if (parts.empty())

      return;  // no active criterion -> no tagged cell -> grid unchanged


    // (R3) UNION (OR) of the tags + dilation (nesting + anticipation of the structures moving).

    TagBox grown = grow_tags(tag_union(parts), regrid_grow_, pdom);


    // (R4)+(R5) cross-rank collective reduction (if coarse distributed) + UNIQUE clustering -> SHARED

    // fine layout. all_reduce_or_inplace is called INSIDE regrid_compute_fine_layout for distributed

    // pk==0: all ranks start from the SAME tag grid -> IDENTICAL fb/dmap per rank (otherwise MPI

    // desync).

    auto [fb, dmap] =

        regrid_compute_fine_layout(std::move(grown), pdom, pk, regrid_margin_, replicated_coarse_);

#ifdef POPS_HAS_MPI

    // MPI COLLECTIVE COUNT (Spec 5 criterion 43): regrid_compute_fine_layout issues ONE

    // all_reduce_or_inplace over the tag grid when the coarse is distributed (multi-rank) -- every rank

    // must cluster from the SAME gathered tags. Count it as one reduction (np>1 only; serial / single

    // rank issues no collective). np==1 is bit-identical with no count.

    if (profiler_ != nullptr && n_ranks() > 1)

      profiler_->count("mpi_reductions");

#endif

    if (fb.size() == 0)

      return;  // nothing to refine: we keep the current grid (no-op)


    // (R6) COHERENT PROLONG / RESTRICT of ALL blocks on the SAME fb/dmap (including the blocks held by

    // their stride: their frozen state is present everywhere and contributes to the Poisson, D3). The

    // ghost width is INHERITED per block (a MUSCL order-2 block carries 2 ghosts; a Minmod block and a

    // VanLeer one may differ), so the scheme does not read out of bounds at the next step (V2 / risk

    // X4).

    for (auto& b : blocks_) {

      auto& L = *b.levels;

      const int ngf = L[fk].U.n_grow();

      L[fk].U = regrid_field_on_layout(fb, dmap, L[pk].U, L[fk].U, pk, ngf, replicated_coarse_);

    }


    // (R7) REBUILD OF THE SHARED AUX (one only, width aux_ncomp_) on the new layout + RE-WIRING of the

    // aux pointer of EACH block. The address &aux_[fk] stays stable (in-place reallocation of the

    // MultiFab in the existing std::vector) -> the pointers of the other levels do not move.

    aux_[fk] = MultiFab(fb, dmap, aux_ncomp_, 1);

    for (auto& b : blocks_)

      (*b.levels)[fk].aux = &aux_[fk];


    // (V3) SHARED-LAYOUT INVARIANT: all blocks MUST live on EXACTLY the same fb/dmap (boxes, order,

    // rank per box) after the regrid. Collective guard (cross-block); catches any inconsistent

    // reconstruction before it corrupts the shared aux / the summed Poisson.

    {

      std::vector<std::vector<AmrLevelMP>> ref;

      ref.reserve(blocks_.size());

      for (const auto& b : blocks_)

        ref.push_back(*b.levels);

      detail::same_layout_or_throw(ref);

    }


    // (R8) RESTORATION OF THE COVERAGE INVARIANT: re-solve so that phi / grad phi are consistent with

    // the new grid AND to trigger the fine -> coarse cascade (mf_average_down_mb, in solve_fields) that

    // restores the covered coarse cells (otherwise a mass diagnostic, sum of the coarse only, would

    // count a phantom coarse value under the new patch, X5).

    solve_fields();

    ++regrid_count_;

    // AMR PROFILING (Spec 5 criterion 43): a regrid COMPLETED -> bump the per-run "regrid" counter

    // (parity with regrid_count_). The "regrid" TIMING scope (_rg above) already covered the whole

    // attempt; this counts only the regrids that actually rebuilt the hierarchy.

    if (profiler_ != nullptr)

      profiler_->count("regrid");

  }


  void step(Real dt) {

    solve_count_ = 0;

    // UNION-TAGS REGRID (capstone Phase 2, C.6; D2: BEFORE the macro-step's step, consistent with the

    // single-block amr_dsl_block.hpp:108). regrid_every_ cadence in MACRO-STEPS, OUTSIDE the substep

    // loops and the stride windows (macro-step granularity ONLY, D3). regrid_every_ == 0 -> FROZEN

    // hierarchy, regrid never called -> BIT-IDENTICAL trajectory to the historical one. The guard

    // macro_step_ > 0 (like the single-block) avoids a regrid at the very first step (the initial grid

    // is already the build one). The regrid sits BEFORE solve_fields below: it does its own

    // solve_fields (R0/R8), then the step's solve_fields recomputes phi on the re-gridded grid.

    if (regrid_every_ > 0 && macro_step_ > 0 && macro_step_ % regrid_every_ == 0)

      regrid();

    // System Poisson solved ONCE on the current state (OncePerStep cadence). A HELD block (stride > 1,

    // outside end-of-window) contributed with its FROZEN state since its last advance: loose coupling

    // assumed by the multirate, exactly like System::step / AmrSystemCoupler in OncePerStep. phi stays

    // frozen during the blocks' advance (no per-substep re-solve here). When reached from step_cfl this

    // re-solves an unchanged state (a second solve), kept on purpose; see the ADC-318 note in step_cfl.

    solve_fields();

    for (auto& b : blocks_) {

      // HOLD-THEN-CATCH-UP cadence (cf. AmrRuntimeBlock::stride, #140): the block is HELD as long as

      // (macro_step_+1) % stride != 0, then CATCHES UP at end-of-window by an effective step stride*dt.

      // The end-of-window catch-up keeps the block temporally consistent with the fast ones at the

      // coupling point (never in the future). stride=1: always true -> every step, bit-identical.

      if ((macro_step_ + 1) % b.stride != 0)

        continue;

      // NEWTON DIAGNOSTICS (OPT-IN): RESET of the report at the HEAD of the block advance (parity with

      // System::AdvanceImex::operator() which resets nreport before its substep loop). The report then

      // AGGREGATES over all the levels AND substeps of THIS advance (imex_advance accumulates per level

      // via backward_euler_source; step() calls imex_advance substeps times without re-resetting).

      // Placed AFTER the stride skip: a HELD block keeps the report of its LAST advance ("last advance"

      // semantics of System). No-op for a block without diagnostics (newton_report null).

      if (b.newton_diagnostics && b.newton_report)

        b.newton_report->reset();

      const Real bdt = dt * static_cast<Real>(b.stride);  // catch-up: effective step stride*dt

      // substeps equal substeps of bdt/substeps. The chosen closure does ONE advance per call;

      // substeps=1 -> a single advance of bdt (bit-identical to the single-substep case). Per-block

      // treatment SELECTION: IMEX (source-free transport + implicit stiff source, mirrors the IMEX

      // branch of AmrSystemCoupler::step) if b.imex, otherwise EXPLICIT (transport + forward-Euler

      // source). The test is PER BLOCK and stable: a single IMEX block changes nothing for the

      // neighboring explicit blocks.

      // NOTE substeps>1: the loop below calls step_block substeps times for BOTH treatments, so the

      // IMEX splitting is SUB-CYCLED (K Lie steps over bdt/K). The compile-time, for its part, applies

      // its IMEX only once over bdt (it ignores substeps on its IMEX branch): divergence INTENTIONAL

      // and sound for substeps>1 (cf. IMEX SEMANTICS UNDER substeps in the file header).

      const Real h = bdt / static_cast<Real>(b.substeps);

      auto& step_block = b.imex ? b.imex_advance : b.advance;

      for (int s = 0; s < b.substeps; ++s)

        step_block(*b.levels, dom_, h, base_per_, replicated_coarse_);

      // PROJECTION PONCTUELLE post-pas (ADC-177) : par niveau, APRES substeps + reflux/cascade.

      // Cell-local + idempotente -> conservation preservee (flux-registres deja regles). No-op si vide.

      if (b.project_per_level)

        b.project_per_level(*b.levels);

    }

    // Inter-species coupled sources AFTER the transport (same order as AmrSystemCoupler: transport then

    // coupled_source_step), by forward-Euler splitting. No-op if no source registered -> bit-identical

    // trajectory to the historical one (the feature is opt-in).

    coupled_source_step(dt);

    ++macro_step_;

  }


  Real step_cfl(Real cfl, Real h) {

    // NOTE (ADC-318): this pre-solve plus step(dt)'s own head solve below is a DOUBLE Poisson solve on

    // the SAME unchanged state (regrid_every=0 freezes the grid in between). It looks redundant but is

    // NOT, and is INTENTIONALLY kept. GeometricMG::solve() is warm-started and iterates to a RELATIVE

    // tolerance (rel_tol 1e-8; abs_tol 0 by default, so its off-step early-exit never fires here), so the

    // second solve does not recompute identical phi: starting from the first solve's iterate it

    // over-converges it by ~rel_tol. Skipping the second solve would therefore NOT be bit-identical; it

    // drifts the trajectory by ~3e-10 over 20 steps (below the solver tolerance and far below the O(dt^2)

    // scheme error, but nonzero). The de-dup was declined to preserve the exact historical bit-stream

    // (SystemStepper::step_cfl avoids the double solve by INLINING its advance, not by skipping a solve).

    solve_fields();  // aux up to date: each block's max_speed reads it on the current coarse

    Real dt = std::numeric_limits<Real>::infinity();

    last_dt_reason_ = "degenerate";

    for (auto& b : blocks_) {

      const Real w = std::max(b.max_speed((*b.levels)[0].U, aux_[0]), kCflSpeedFloor);

      Real dt_b = cfl * h * static_cast<Real>(b.substeps) / (static_cast<Real>(b.stride) * w);

      const char* why = "transport";

      // OPTIONAL block BOUNDS (AMR StabilityPolicy, audit 2026-06): same substeps/stride formulas as

      // SystemStepper::step_cfl, evaluated on the COARSE. Empty closures (model without the trait) ->

      // not queried, transport bound only (bit-identical).

      if (b.source_frequency) {

        const Real mu = b.source_frequency((*b.levels)[0].U, aux_[0]);

        if (mu > Real(0)) {

          const Real dt_src =

              cfl * static_cast<Real>(b.substeps) / (static_cast<Real>(b.stride) * mu);

          if (dt_src < dt_b) {

            dt_b = dt_src;

            why = "source_frequency";

          }

        }

      }

      if (b.stability_dt) {

        const Real db = b.stability_dt((*b.levels)[0].U, aux_[0]);

        if (db > Real(0)) {

          const Real dt_adm = db * static_cast<Real>(b.substeps) / static_cast<Real>(b.stride);

          if (dt_adm < dt_b) {

            dt_b = dt_adm;

            why = "stability_dt";

          }

        }

      }

      if (dt_b < dt) {

        dt = dt_b;

        last_dt_reason_ = std::string(why) + ":" + b.name;

      }

    }

    // Declared frequencies of the coupled sources (CoupledSource.frequency): bound on the MACRO-step

    // (the couplings apply once per macro-step), dt <= cfl / mu, without substeps/stride.

    for (const auto& cs : coupled_freqs_) {

      const Real dt_cs = cfl / cs.mu;

      if (dt_cs < dt) {

        dt = dt_cs;

        last_dt_reason_ = "coupled_source:" + cs.label;

      }

    }

    // PER-CELL frequencies (CoupledSource.frequency with an Expr): mu(U) reduced (MAX) on the COARSE

    // level of the input blocks (where the AMR CFL lives), GLOBAL all_reduce_max (ALL ranks, neutral

    // without a local box), bound dt <= cfl / max(mu). Same reason "coupled_source:<label>" as the

    // constant. No per-cell source -> empty loop (bit-identical). The Array4 are rebuilt at EACH step

    // (the hierarchy fabs are repointed by the regrid), like coupled_source_step.

    for (const auto& ce : coupled_freq_exprs_) {

      Real m = 0;

      if (ce.n_in > 0) {

        auto& Uref =

            (*blocks_[static_cast<std::size_t>(ce.ins[0].block)].levels)[0].U;  // coarse (lev 0)

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

          CoupledFreqKernel kern;

          kern.n_in = ce.n_in;

          kern.n_const = ce.n_const;

          for (int c = 0; c < ce.n_in; ++c) {

            kern.in[c] =

                (*blocks_[static_cast<std::size_t>(ce.ins[static_cast<std::size_t>(c)].block)]

                      .levels)[0]

                    .U.fab(li)

                    .array();

            kern.in_comp[c] = ce.ins[static_cast<std::size_t>(c)].comp;

          }

          for (int c = 0; c < ce.n_const; ++c)

            kern.consts[c] = ce.kconsts[static_cast<std::size_t>(c)];

          kern.prog = ce.prog;

          m = std::max(m, reduce_max_cell(Uref.box(li), kern));

        }

      } else {

        // Program WITHOUT an input field (constant in bytecode): evaluated once on the constants.

        Real reg[kCsMaxReg];

        for (int c = 0; c < ce.n_const; ++c)

          reg[c] = ce.kconsts[static_cast<std::size_t>(c)];

        const Real mu0 = ce.prog.eval(reg);

        if (mu0 > Real(0))

          m = mu0;

      }

      const double mu = all_reduce_max(static_cast<double>(m));  // ALL ranks (collective symmetry)

#ifdef POPS_HAS_MPI

      // MPI COLLECTIVE COUNT (Spec 5 criterion 43): one all_reduce_max per per-cell coupled-frequency

      // bound, multi-rank only (serial all_reduce_max is an identity, no collective).

      if (profiler_ != nullptr && n_ranks() > 1)

        profiler_->count("mpi_reductions");

#endif

      if (mu > 0.0) {

        const Real dt_cs = cfl / static_cast<Real>(mu);

        if (dt_cs < dt) {

          dt = dt_cs;

          last_dt_reason_ = "coupled_source:" + ce.label;

        }

      }

    }

    // GLOBAL bounds (AmrRuntime::add_dt_bound, parity with System::add_dt_bound): evaluated PER RANK

    // then reduced all_reduce_min (dt identical on all ranks; <= 0/non-finite = inert).

    for (const auto& g : dt_bounds_) {

      if (!g.fn)

        continue;

      double v = g.fn();

      if (!(v > 0.0) || !std::isfinite(v))

        v = std::numeric_limits<double>::infinity();

      v = all_reduce_min(v);

#ifdef POPS_HAS_MPI

      // MPI COLLECTIVE COUNT (Spec 5 criterion 43): one all_reduce_min per registered global dt bound,

      // multi-rank only (the global min keeps dt identical on all ranks). Serial -> identity, no count.

      if (profiler_ != nullptr && n_ranks() > 1)

        profiler_->count("mpi_reductions");

#endif

      if (static_cast<Real>(v) < dt) {

        dt = static_cast<Real>(v);

        last_dt_reason_ = "global:" + g.label;

      }

    }

    if (!std::isfinite(dt)) {

      dt = cfl * h / kCflSpeedFloor;  // guard (no block: impossible here)

      last_dt_reason_ = "degenerate";

    }

    step(dt);

    return dt;

  }


  int macro_step() const { return macro_step_; }

  void set_macro_step(int s) { macro_step_ = s; }


  void set_profiler(runtime::program::Profiler* prof) { profiler_ = prof; }


  void add_dt_bound(const std::string& label, std::function<double()> fn) {

    dt_bounds_.push_back(GlobalDtBound{label, std::move(fn)});

  }


  void add_coupled_frequency(const std::string& label, Real mu) {

    if (mu > Real(0))

      coupled_freqs_.push_back(CoupledFreqDecl{label, mu});

  }


  void add_coupled_frequency_expr(const std::string& label,

                                  const std::vector<std::string>& in_blocks,

                                  const std::vector<std::string>& in_roles,

                                  const std::vector<double>& consts,

                                  const std::vector<int>& freq_prog_ops,

                                  const std::vector<int>& freq_prog_args) {

    if (freq_prog_ops.empty() && freq_prog_args.empty())

      return;  // no per-cell frequency

    const int n_in = static_cast<int>(in_blocks.size());

    const int n_const = static_cast<int>(consts.size());

    if (static_cast<int>(in_roles.size()) != n_in)

      throw std::runtime_error(

          "AmrRuntime::add_coupled_frequency_expr : in_blocks / in_roles of different sizes");

    if (n_in + n_const > kCsMaxReg)

      throw std::runtime_error(

          "AmrRuntime::add_coupled_frequency_expr : too many registers (inputs + constants > " +

          std::to_string(kCsMaxReg) + ")");

    if (freq_prog_ops.size() != freq_prog_args.size())

      throw std::runtime_error(

          "AmrRuntime::add_coupled_frequency_expr : freq_prog_ops / freq_prog_args of different "

          "sizes");

    if (static_cast<int>(freq_prog_ops.size()) > kCsMaxProg)

      throw std::runtime_error(

          "AmrRuntime::add_coupled_frequency_expr : frequency program too long (> " +

          std::to_string(kCsMaxProg) + ")");

    // Resolves (block, role) -> (block index, component), STRICT (mirror of add_coupled_source).

    std::vector<CsRef> ins(static_cast<std::size_t>(n_in));

    for (int c = 0; c < n_in; ++c) {

      const std::string& block = in_blocks[static_cast<std::size_t>(c)];

      const std::string& role = in_roles[static_cast<std::size_t>(c)];

      const int b = block_index(block);

      if (b < 0)

        throw std::runtime_error("AmrRuntime::add_coupled_frequency_expr : no block named '" +

                                 block + "'");

      // Role addressed BY NAME: a canonical role name OR a user-defined role label (ADC-292), STRICT.

      const VariableSet& vs = blocks_[static_cast<std::size_t>(b)].cons_vars;

      const int comp = vs.index_of(role);

      if (comp < 0)

        throw std::runtime_error("AmrRuntime::add_coupled_frequency_expr : block '" + block +

                                 "' does not expose role '" + role + "' (roles: " +

                                 (vs.roles.empty() ? std::string("<none>") : roles_csv(vs)) +

                                 ", no silent fallback to component 0)");

      ins[static_cast<std::size_t>(c)] = {b, comp, CsProgram{}};

    }

    CsProgram pg;

    pg.len = static_cast<int>(freq_prog_ops.size());

    for (int k = 0; k < pg.len; ++k) {

      const int opc = freq_prog_ops[static_cast<std::size_t>(k)];

      const int a = freq_prog_args[static_cast<std::size_t>(k)];

      if (opc < 0 || opc > static_cast<int>(CsOp::Sqrt))

        throw std::runtime_error(

            "AmrRuntime::add_coupled_frequency_expr : invalid opcode in the frequency");

      if (opc == static_cast<int>(CsOp::PushReg) && (a < 0 || a >= n_in + n_const))

        throw std::runtime_error(

            "AmrRuntime::add_coupled_frequency_expr : register out of bounds in the frequency");

      pg.op[k] = opc;

      pg.arg[k] = a;

    }

    std::vector<Real> kconsts(consts.begin(), consts.end());

    coupled_freq_exprs_.push_back(

        CoupledFreqExprDecl{label, std::move(ins), pg, n_in, n_const, std::move(kconsts)});

  }


  const std::string& last_dt_bound() const { return last_dt_reason_; }


  const NewtonReport& newton_report(const std::string& name) const {

    const int b = block_index(name);

    if (b < 0)

      throw std::runtime_error("AmrRuntime::newton_report : no block named '" + name + "'");

    const AmrRuntimeBlock& blk = blocks_[static_cast<std::size_t>(b)];

    if (!blk.newton_diagnostics || !blk.newton_report)

      throw std::runtime_error(

          "AmrRuntime::newton_report : Newton diagnostics not enabled for block '" + name +

          "' ; add the block with newton_diagnostics=True (pops.IMEX(newton_diagnostics=True))");

    return *blk.newton_report;

  }


  std::vector<double> potential() {

    solve_fields();

    return blocks_[0].potential(aux_[0]);

  }


  Real max_speed() {

    solve_fields();

    Real w = Real(1e-12);

    for (auto& b : blocks_) {

      const Real wb = b.max_speed((*b.levels)[0].U, aux_[0]);

      if (wb > w)

        w = wb;

    }

    return w;

  }


  int n_patches() const {

    const auto& L = *blocks_[0].levels;

    return L.size() >= 2 ? static_cast<int>(L[1].U.box_array().size()) : 0;

  }


  // Index-space signatures of the fine patches (level + inclusive lo/hi corners), for ALL fine levels.

  // Read-only of the GLOBAL BoxArray (all boxes/all ranks) already stored -> rank-independent, zero

  // communication, NO hot-path cost (query between steps). Mirror of n_patches(): the same box_array()

  // that gives the COUNT gives the BOXES. Block 0 representative (SHARED layout, same_layout_or_throw

  // guard). Loop k = 1..nlev-1: a single fine level today (ratio 2), correct if a future adds levels

  // (the level field disambiguates the spacing dx = L / (n << level) Python-side).


  std::vector<PatchBox> patch_boxes() const {

    const auto& L = *blocks_[0].levels;

    std::vector<PatchBox> out;

    for (int k = 1; k < static_cast<int>(L.size()); ++k) {

      const auto& bxs = L[k].U.box_array().boxes();

      for (const Box2D& b : bxs)

        out.push_back(PatchBox{k, b.lo[0], b.lo[1], b.hi[0], b.hi[1]});

    }

    return out;

  }


  // COARSE-level (base) box counts (ADC-319, MPI ownership diagnostic). Block 0 is the SHARED layout

  // (same_layout_or_throw), so its level-0 MultiFab carries the base BoxArray + DistributionMapping

  // common to all blocks. local_size() = base boxes OWNED by this rank; box_array().size() = total base

  // boxes (all ranks). Mirror of n_patches(): a query between steps, no communication, no hot-path cost.

  int coarse_local_boxes() const { return (*blocks_[0].levels)[0].U.local_size(); }

  int coarse_total_boxes() const { return (*blocks_[0].levels)[0].U.box_array().size(); }


  // ----------------------------------------------------------------------------------------------

  // MULTI-BLOCK AMR CHECKPOINT / RESTART (ADC-509). PER-BLOCK PER-LEVEL state accessors + the

  // level-0 phi (multigrid warm-start), counterpart of AmrCouplerMP::level_state on the SHARED

  // hierarchy. The shared layout is FROZEN at build (make_shared_amr_layout: a deterministic central

  // fine patch, regrid_every==0): replaying the SAME composition reproduces the SAME hierarchy, so a

  // restart only needs to restore each block's valid cells + phi (no set_hierarchy on the runtime).

  // The _global variants all_reduce_sum the per-rank LOCAL fabs into the complete field (np>1

  // gather), MIRROR of System::state_global / gather_global; mono-rank they are the identity. @p b:

  // block index, @p k: level (0 = coarse, >= 1 = fine).

  // ----------------------------------------------------------------------------------------------


  // Conserved components of block @p b (Model::n_vars, carried by the AmrRuntimeBlock).


  int block_n_vars(std::size_t b) const {

    if (b >= blocks_.size())

      throw std::runtime_error("AmrRuntime::block_n_vars : block index out of bounds");

    return blocks_[b].ncomp;

  }


  // FULL conservative state (all components) of block @p b at level @p k, flat component-major

  // c*nf*nf + j*nf + i (nf = nx << k); zeros outside the patches at the fine level. LOCAL fabs only

  // (no gather): the facade calls this mono-rank. Mirror of AmrCouplerMP::level_state.


  std::vector<double> block_level_state(std::size_t b, int k) const {

    if (b >= blocks_.size())

      throw std::runtime_error("AmrRuntime::block_level_state : block index out of bounds");

    const std::vector<AmrLevelMP>& L = *blocks_[b].levels;

    if (k < 0 || k >= static_cast<int>(L.size()))

      throw std::runtime_error("AmrRuntime::block_level_state : level out of bounds");

    const MultiFab& U = L[k].U;

    const int nc = U.ncomp();

    const std::size_t nf = static_cast<std::size_t>(dom_.nx()) << k;

    std::vector<double> out(static_cast<std::size_t>(nc) * nf * nf, 0.0);

    device_fence();

    fill_level_state(U, nc, nf, out);

    return out;

  }


  // Same as block_level_state but all_reduce_sum the per-rank contributions -> every rank holds the

  // complete field (np>1 gather, AMR reflux pattern). COLLECTIVE: all ranks MUST call it.


  std::vector<double> block_level_state_global(std::size_t b, int k) const {

    std::vector<double> out = block_level_state(b, k);

    all_reduce_sum_inplace(out.data(), static_cast<int>(out.size()));

    return out;

  }


  // Restores block @p b at level @p k from @p s (same layout as block_level_state). Writes ONLY the

  // VALID cells of the local fabs (the ghosts are redone at the next solve_fields/step, like after a

  // regrid). NO re-prolongation: restored AS-IS. Mirror of AmrCouplerMP::set_level_state.


  void set_block_level_state(std::size_t b, int k, const std::vector<double>& s) {

    if (b >= blocks_.size())

      throw std::runtime_error("AmrRuntime::set_block_level_state : block index out of bounds");

    std::vector<AmrLevelMP>& L = *blocks_[b].levels;

    if (k < 0 || k >= static_cast<int>(L.size()))

      throw std::runtime_error("AmrRuntime::set_block_level_state : level out of bounds");

    MultiFab& U = L[k].U;

    const int nc = U.ncomp();

    const std::size_t nf = static_cast<std::size_t>(dom_.nx()) << k;

    if (s.size() != static_cast<std::size_t>(nc) * nf * nf)

      throw std::runtime_error("AmrRuntime::set_block_level_state : state size != ncomp*nf*nf");

    device_fence();

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

      Array4 u = U.fab(li).array();

      const Box2D v = U.box(li);

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

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

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

            u(i, j, c) = s[static_cast<std::size_t>(c) * nf * nf +

                           static_cast<std::size_t>(j) * nf + static_cast<std::size_t>(i)];

    }

  }


  // Potential phi of level @p k, flat nf*nf row-major, zeros outside patches. Level 0: the multigrid

  // WARM-START mg_.phi() (the state reused by the next solve -> bit-identical restart). Level >= 1:

  // shared aux comp 0 (recomputed at solve_fields). Mirror of AmrCouplerMP::level_potential; the phi

  // is SHARED by all blocks (single aux), so it carries no block index. NON-const like

  // AmrRuntime::potential() (GeometricMG::phi() returns a mutable warm-start reference).


  std::vector<double> level_potential(int k) {

    if (k < 0 || k >= nlev_)

      throw std::runtime_error("AmrRuntime::level_potential : level out of bounds");

    const std::size_t nf = static_cast<std::size_t>(dom_.nx()) << k;

    std::vector<double> out(nf * nf, 0.0);

    device_fence();

    const MultiFab& P = (k == 0) ? mg_.phi() : aux_[k];

    fill_level_phi(P, nf, out);

    return out;

  }


  // Same as level_potential but all_reduce_sum (np>1 gather). COLLECTIVE: all ranks MUST call it.


  std::vector<double> level_potential_global(int k) {

    std::vector<double> out = level_potential(k);

    all_reduce_sum_inplace(out.data(), static_cast<int>(out.size()));

    return out;

  }


  // Restores phi of level @p k. Level 0: warm-start mg_.phi() -> bit-identical restart (1st

  // post-restart solve starts from the same guess). Level >= 1: shared aux comp 0 (idempotent,

  // recomputed at solve_fields). Mirror of AmrCouplerMP::set_level_potential.


  void set_level_potential(int k, const std::vector<double>& p) {

    if (k < 0 || k >= nlev_)

      throw std::runtime_error("AmrRuntime::set_level_potential : level out of bounds");

    const std::size_t nf = static_cast<std::size_t>(dom_.nx()) << k;

    if (p.size() != nf * nf)

      throw std::runtime_error("AmrRuntime::set_level_potential : phi size != nf*nf");

    device_fence();

    MultiFab& P = (k == 0) ? mg_.phi() : aux_[k];

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

      Array4 q = P.fab(li).array();

      const Box2D v = P.box(li);

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

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

          q(i, j, 0) = p[static_cast<std::size_t>(j) * nf + static_cast<std::size_t>(i)];

    }

  }


 private:

  // Fills @p out (zero-initialized, size nc*nf*nf) from the LOCAL valid cells of @p U at GLOBAL

  // component-major indices c*nf*nf + j*nf + i. Shared by block_level_state and its _global gather

  // variant (the loop is verbatim with AmrCouplerMP::level_state -> bit-identical layout).

  static void fill_level_state(const MultiFab& U, int nc, std::size_t nf,

                               std::vector<double>& out) {

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

      const ConstArray4 u = U.fab(li).const_array();

      const Box2D v = U.box(li);

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

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

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

            out[static_cast<std::size_t>(c) * nf * nf + static_cast<std::size_t>(j) * nf +

                static_cast<std::size_t>(i)] = u(i, j, c);

    }

  }


  // Fills @p out (zero-initialized, size nf*nf) from the LOCAL valid cells of @p P (comp 0) at GLOBAL

  // row-major indices j*nf + i. Shared by level_potential and its _global gather variant.

  static void fill_level_phi(const MultiFab& P, std::size_t nf, std::vector<double>& out) {

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

      const ConstArray4 p = P.fab(li).const_array();

      const Box2D v = P.box(li);

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

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

          out[static_cast<std::size_t>(j) * nf + static_cast<std::size_t>(i)] = p(i, j, 0);

    }

  }


  // Re-applies the model-NAMED aux fields (ADC-291) onto the COARSE shared aux valid cells. Mirror of

  // SystemFieldSolver::apply_named_aux_one (cartesian System): per LOCAL fab (MPI-safe), valid cells

  // only, global flat index j*nx+i. The coarse layout is frozen across regrid (only fine levels are

  // rebuilt), so the stored coarse field stays valid; solve_fields runs the coarse->fine injection

  // right after, carrying the named comps to every level. No-op without a named field.

  void apply_named_aux() {

    if (named_aux_.empty() || aux_.empty())

      return;

    const int row = dom_.nx();

    for (const auto& [comp, field] : named_aux_) {

      if (field.empty() || comp >= aux_ncomp_)

        continue;

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

        Array4 a = aux_[0].fab(li).array();

        const Box2D v = aux_[0].box(li);

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

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

            a(i, j, comp) = field[static_cast<std::size_t>(j) * row + i];

      }

    }

  }


  // NAMED multi-elliptic field (ADC-428): a field name's aux output components + a DEDICATED coarse

  // GeometricMG (built lazily, REUSING the native solver; the operator is never reimplemented).

  // shared_ptr<GeometricMG>: GeometricMG owns Fabs (non-copyable/non-movable), and named_fields_ is a

  // std::map (stable nodes), so a heap GeometricMG gives a stable address without making NamedField

  // movable. Defined here (before ensure_named_elliptic, whose parameter type must be visible).

  struct NamedField {

    int phi_comp = -1;

    int gx_comp = -1;

    int gy_comp = -1;

    std::shared_ptr<GeometricMG>

        mg;  // dedicated coarse solver, built lazily by ensure_named_elliptic

  };


  // Builds a NAMED elliptic field's dedicated coarse GeometricMG (ADC-428), lazily, IDENTICAL to the

  // default mg_ (same coarse geometry / BoxArray / Poisson BC / wall predicate / replication). REUSES the

  // native solver -- no operator is reimplemented. The variable / anisotropic permittivity of the default

  // Poisson is NOT carried onto a named field (a named field is a plain Laplacian, like the System named

  // path). No-op if already built.

  void ensure_named_elliptic(NamedField& nf) {

    if (nf.mg)

      return;

    nf.mg =

        std::make_shared<GeometricMG>(geom_, ba_coarse_, bcPhi_, wall_active_, replicated_coarse_);

  }


  // Reads aux component @p comp of the COARSE level as a GLOBAL n*n row-major field (diagnostic /

  // read-back). Same marshaling as detail::coupler_read_coarse_phi (the default potential read-back) but

  // for an arbitrary component: local n*n buffer, all_reduce_sum_inplace when the coarse is DISTRIBUTED

  // (disjoint boxes -> exact recompose; serial / replicated is identity). Used by named_field_values.

  std::vector<double> coarse_aux_component(int comp) const {

    device_fence();

    const int nx = dom_.nx(), ny = dom_.ny();

    std::vector<double> out(static_cast<std::size_t>(nx) * ny, 0.0);

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

      const ConstArray4 a = aux_[0].fab(li).const_array();

      const Box2D v = aux_[0].box(li);

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

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

          out[static_cast<std::size_t>(j) * nx + i] = static_cast<double>(a(i, j, comp));

    }

    if (!replicated_coarse_)

      all_reduce_sum_inplace(out.data(), static_cast<int>(out.size()));

    return out;

  }


  // Per-field aux HALO override (ADC-369) on the COARSE aux, AFTER the shared fill_ghosts. Overrides

  // only each declared component's physical-face ghosts (aux_halo_override keeps periodic faces

  // periodic). No-op without a policy. Mirror of SystemFieldSolver::apply_named_aux_bc.

  void apply_named_aux_bc() {

    if (named_aux_bc_.empty() || aux_.empty())

      return;

    for (const auto& [comp, policy] : named_aux_bc_) {

      if (comp >= aux_ncomp_)

        continue;

      fill_physical_bc(aux_[0], dom_, aux_halo_override(aux_bc_, policy), comp);

    }

  }


  // Index of the block named @p name in the registry (-1 if absent). Counterpart of

  // AmrSystem::Impl::block_index (the facade names the blocks; the coupled sources target them by name,

  // resolved once at registration).

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

    for (std::size_t i = 0; i < blocks_.size(); ++i)

      if (blocks_[i].name == name)

        return static_cast<int>(i);

    return -1;

  }


  // Resolved reference of a coupled-source field: (block index, component) + the term bytecode program

  // (empty for an input). Inputs carry only block/comp; outputs carry in addition the postfix program

  // evaluated per cell. We capture the block INDEX (not a fab pointer): the Array4 are rebuilt at each

  // application, per level.

  struct CsRef {

    int block;

    int comp;

    CsProgram prog;  // outputs: term program; inputs: unused (CsProgram{})

  };

  // A registered coupled source: its inputs, its output terms and its constants, ready to be marshaled

  // into a CoupledSourceKernel per level / per fab at application.

  struct CoupledSourceSpec {

    std::vector<CsRef> ins;

    std::vector<CsRef> outs;

    std::vector<Real> kconsts;

    int n_in = 0;

    int n_const = 0;

    int n_terms = 0;

  };


  Geometry geom_;

  Box2D dom_;

  Periodicity base_per_;

  BCRec bcPhi_, aux_bc_;

  bool replicated_coarse_;

  GeometricMG mg_;

  // Coarse BoxArray + conductive-wall predicate stashed at the ctor so a NAMED elliptic field (ADC-428)

  // can build its own coarse GeometricMG identical to mg_ (ensure_named_elliptic). mg_ consumes them at

  // its construction but does not expose them, so we keep a copy here (cheap; the coarse layout is small).

  BoxArray ba_coarse_;

  std::function<bool(Real, Real)> wall_active_;

  std::vector<AmrRuntimeBlock> blocks_;

  // GLOBAL step bounds (add_dt_bound, parity with System) + ACTIVE bound of the last step_cfl.

  struct GlobalDtBound {

    std::string label;

    std::function<double()> fn;

  };

  std::vector<GlobalDtBound> dt_bounds_;

  // Declared frequencies of the coupled sources (bound dt <= cfl/mu on the macro-step, wave 3).

  struct CoupledFreqDecl {

    std::string label;

    Real mu;

  };

  std::vector<CoupledFreqDecl> coupled_freqs_;

  // PER-CELL frequencies of the coupled sources (CoupledSource.frequency with an Expr): bytecode

  // program mu(U) evaluated on the coarse at each step_cfl (MAX + all_reduce_max -> dt <= cfl/max(mu)).

  // ins = (block, comp) of the inputs (prog unused); kconsts = constants (same as the source).

  struct CoupledFreqExprDecl {

    std::string label;

    std::vector<CsRef> ins;

    CsProgram prog;

    int n_in = 0;

    int n_const = 0;

    std::vector<Real> kconsts;

  };

  std::vector<CoupledFreqExprDecl> coupled_freq_exprs_;

  std::string last_dt_reason_;

  std::vector<MultiFab> aux_;  // [level], shared by all blocks

  // Model-NAMED aux fields (ADC-291): component (>= kAuxNamedBase) -> coarse base-level field

  // (n*n row-major). STATIC user fields re-applied by solve_fields each macro-step (so they persist

  // across regrid). Empty by default -> bit-identical. cf. set_named_aux / apply_named_aux.

  std::map<int, std::vector<Real>> named_aux_;

  // Per-field aux HALO policy (ADC-369): component -> uniform boundary policy, applied to the coarse aux

  // after the shared fill (apply_named_aux_bc). Empty by default -> bit-identical.

  std::map<int, AuxHaloPolicy> named_aux_bc_;

  // NAMED multi-elliptic fields (ADC-428): field name -> its aux output components + a DEDICATED coarse

  // GeometricMG. The NamedField struct itself is defined higher up (before ensure_named_elliptic, which

  // takes it by reference: a parameter type must be visible at the function declaration, unlike a member

  // body). Empty default -> bit-identical (the solve_named_fields loop early-returns).

  std::map<std::string, NamedField> named_fields_;

  std::vector<CoupledSourceSpec>

      coupled_sources_;  // registered coupled sources (applied after transport)

  // UNION-TAGS REGRID (capstone Phase 2, C.6). regrid_every_ == 0 -> FROZEN hierarchy (default,

  // bit-identical). block_tag_: PER-BLOCK tag predicate (D1; same size as blocks_, empty = this block

  // tags nothing on its side). phi_tag_: phi tag predicate on |grad phi| (D4; empty = phi does not

  // contribute to the union).

  std::vector<TagPredicate> block_tag_;

  TagPredicate phi_tag_;

  int regrid_every_ = 0;

  int regrid_grow_ = 2;

  int regrid_margin_ = 2;

  int aux_ncomp_ = kAuxBaseComps;

  int nlev_ = 0;

  int macro_step_ = 0;

  mutable int solve_count_ = 0;

  int regrid_count_ = 0;

  // AMR / MPI PROFILING (Spec 5 criterion 43, ADC-479): non-owning pointer to the AmrSystem-owned

  // Profiler (lifetime guaranteed by the facade). Null by default -> the engine never profiles

  // (zero overhead). Set via set_profiler after build (parity with System::profiler_).

  runtime::program::Profiler* profiler_ = nullptr;


  // RAII phase-timing scope for an AMR phase (regrid / average_down / fill_boundary). Mirrors

  // runtime::program::ProfileScope but over a NULLABLE profiler pointer: it reads the clock and

  // records only when profiler_ is non-null AND enabled. A null/disabled run constructs a cheap inert

  // scope (one pointer copy + one clock read) -- the granularity is per-phase, not per-cell, so this

  // is off the hot path. Returned BY VALUE from profile_amr_scope (movable: only POD members).

  class AmrPhaseScope {

   public:

    AmrPhaseScope(runtime::program::Profiler* prof, const char* name)

        : prof_(prof != nullptr && prof->enabled() ? prof : nullptr),

          name_(name),

          t0_(std::chrono::steady_clock::now()) {}

    AmrPhaseScope(AmrPhaseScope&& o) noexcept : prof_(o.prof_), name_(o.name_), t0_(o.t0_) {

      o.prof_ = nullptr;  // the moved-from scope must not record

    }

    AmrPhaseScope(const AmrPhaseScope&) = delete;

    AmrPhaseScope& operator=(const AmrPhaseScope&) = delete;

    AmrPhaseScope& operator=(AmrPhaseScope&&) = delete;

    ~AmrPhaseScope() {

      if (prof_ == nullptr)

        return;

      const auto t1 = std::chrono::steady_clock::now();

      try {

        prof_->record(name_, std::chrono::duration<double>(t1 - t0_).count());

      } catch (...) {  // NOLINT(bugprone-empty-catch) -- a profiler never throws out of a scope

      }

    }


   private:

    runtime::program::Profiler* prof_;

    const char* name_;

    std::chrono::steady_clock::time_point t0_;

  };


  // Build an AMR phase scope (no-op when profiling is off). Used at the head of regrid /

  // average_down / fill_boundary.

  AmrPhaseScope profile_amr_scope(const char* name) { return AmrPhaseScope(profiler_, name); }


  // fill_ghosts wrapped in the "fill_boundary" timing scope + per-call count (Spec 5 criterion 43).

  // Under MPI np>1 the ghost fill is a cross-rank halo exchange -> also count one "mpi_messages" (a

  // point-to-point round, distinct from the all_reduce collectives counted as "mpi_reductions").

  // Serial / single rank: no message count (the fill is a local copy). Off-profiling: inert.

  void fill_ghosts_profiled(MultiFab& mf, const Box2D& dom, const BCRec& bc) {

    auto _fb = profile_amr_scope("fill_boundary");

    if (profiler_ != nullptr) {

      profiler_->count("fill_boundary");

#ifdef POPS_HAS_MPI

      if (n_ranks() > 1)

        profiler_->count("mpi_messages");

#endif

    }

    fill_ghosts(mf, dom, bc);

  }

};


}  // namespace pops

amr_coupler_mp.hpp
AmrCouplerMP: MULTI-PATCH E x B AMR coupler (coarse Poisson -> aux = grad phi -> fine injection -> co...

amr_reflux_mf.hpp
Umbrella for the AMR MultiFab stack: includes the numerics/time sub-headers in dependency order (flux...

amr_regrid_coupler.hpp
amr_regrid_finest: Berger-Rigoutsos regrid of the finest level (responsibility b).

amr_system_coupler.hpp
AmrSystemCoupler: multi-species SYSTEM coupler on AMR (milestone 2.3).

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

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

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

pops::AmrRuntime
AMR multi-block engine at runtime.
Definition amr_runtime.hpp:229

pops::AmrRuntime::level_potential
std::vector< double > level_potential(int k)
Definition amr_runtime.hpp:1309

pops::AmrRuntime::set_profiler
void set_profiler(runtime::program::Profiler *prof)
AMR / MPI PROFILING SEAM (Spec 5 sec.12.5, ADC-479 criterion 43).
Definition amr_runtime.hpp:1079

pops::AmrRuntime::solve_fields
void solve_fields()
sync_down (per block) + system coarse Poisson (CO-LOCATED SUMMED RHS) + coarse aux + fine injection.
Definition amr_runtime.hpp:622

pops::AmrRuntime::coarse_total_boxes
int coarse_total_boxes() const
Definition amr_runtime.hpp:1232

pops::AmrRuntime::step_cfl
Real step_cfl(Real cfl, Real h)
substeps/stride-aware CFL step (runtime counterpart of System::step_cfl, EXACT mirror of its formula)...
Definition amr_runtime.hpp:929

pops::AmrRuntime::levels
std::vector< AmrLevelMP > & levels(std::size_t b)
Definition amr_runtime.hpp:319

pops::AmrRuntime::TagPredicate
std::function< bool(const ConstArray4 &, int, int)> TagPredicate
Tag predicate of the union regrid: (ConstArray4 of the read field, i, j) -> should we refine ?...
Definition amr_runtime.hpp:330

pops::AmrRuntime::n_coupled_sources
std::size_t n_coupled_sources() const
Definition amr_runtime.hpp:312

pops::AmrRuntime::mass
Real mass(std::size_t b) const
Definition amr_runtime.hpp:320

pops::AmrRuntime::aux
const MultiFab & aux(int k) const
Definition amr_runtime.hpp:318

pops::AmrRuntime::patch_boxes
std::vector< PatchBox > patch_boxes() const
Definition amr_runtime.hpp:1216

pops::AmrRuntime::block_n_vars
int block_n_vars(std::size_t b) const
Definition amr_runtime.hpp:1246

pops::AmrRuntime::set_block_tag_predicate
void set_block_tag_predicate(std::size_t b, TagPredicate crit)
Registers the TAG PREDICATE of block b (D1: PER-BLOCK union criterion).
Definition amr_runtime.hpp:353

pops::AmrRuntime::max_speed
Real max_speed()
Max SYSTEM wave speed (max over the blocks) on the current coarse. Requires the aux up to date.
Definition amr_runtime.hpp:1194

pops::AmrRuntime::density
std::vector< double > density(std::size_t b) const
Definition amr_runtime.hpp:321

pops::AmrRuntime::AmrRuntime
AmrRuntime(const Geometry &geom, const BoxArray &ba_coarse, const BCRec &bcPhi, std::vector< AmrRuntimeBlock > blocks, Periodicity base_per=Periodicity{true, true}, bool replicated_coarse=true, std::function< bool(Real, Real)> active={})
Definition amr_runtime.hpp:238

pops::AmrRuntime::coupled_source_step
void coupled_source_step(Real dt)
Applies ALL the registered coupled sources of a step dt, by forward-Euler splitting.
Definition amr_runtime.hpp:571

pops::AmrRuntime::newton_report
const NewtonReport & newton_report(const std::string &name) const
NEWTON REPORT (OPT-IN IMEX diagnostics) of block name, AGGREGATED over the levels and substeps of its...
Definition amr_runtime.hpp:1174

pops::AmrRuntime::register_named_field
void register_named_field(const std::string &field, int phi_comp, int gx_comp, int gy_comp)
Definition amr_runtime.hpp:402

pops::AmrRuntime::potential
std::vector< double > potential()
Coarse potential (component 0 of the shared aux) as an n*n row-major field.
Definition amr_runtime.hpp:1188

pops::AmrRuntime::level_potential_global
std::vector< double > level_potential_global(int k)
Definition amr_runtime.hpp:1321

pops::AmrRuntime::set_named_aux_bc
void set_named_aux_bc(int comp, AuxHaloPolicy policy)
Registers a per-field aux HALO policy (ADC-369) for the named component comp: solve_fields applies it...
Definition amr_runtime.hpp:384

pops::AmrRuntime::add_coupled_frequency
void add_coupled_frequency(const std::string &label, Real mu)
DECLARED frequency of a coupled source (CoupledSource.frequency, wave-3 audit): step bound dt <= cfl ...
Definition amr_runtime.hpp:1090

pops::AmrRuntime::add_coupled_source
void add_coupled_source(const std::vector< std::string > &in_blocks, const std::vector< std::string > &in_roles, const std::vector< double > &consts, const std::vector< std::string > &out_blocks, const std::vector< std::string > &out_roles, const std::vector< int > &prog_ops, const std::vector< int > &prog_args, const std::vector< int > &prog_lens)
Registers an inter-species COUPLED SOURCE (DSL CoupledSource, P5 bytecode) on the runtime facade,...
Definition amr_runtime.hpp:458

pops::AmrRuntime::n_patches
int n_patches() const
Definition amr_runtime.hpp:1205

pops::AmrRuntime::poisson_rhs
MultiFab & poisson_rhs()
Definition amr_runtime.hpp:317

pops::AmrRuntime::set_block_named_elliptic_rhs
void set_block_named_elliptic_rhs(std::size_t b, const std::string &field, std::function< void(const MultiFab &, MultiFab &)> rhs)
Attaches named field's RHS contribution closure (rhs += elliptic_field_rhs(U_b)) to block b.
Definition amr_runtime.hpp:411

pops::AmrRuntime::regrid
void regrid()
UNION-TAGS REGRID (capstone Phase 2, C.6; docs/AMR_REGRID_UNION_TAGS_DESIGN.md, steps R0-R8).
Definition amr_runtime.hpp:746

pops::AmrRuntime::solve_count
int solve_count() const
Definition amr_runtime.hpp:322

pops::AmrRuntime::set_macro_step
void set_macro_step(int s)
RESTORES the macro-step counter (IO v1, reserved for restart via AmrSystem::set_clock): without it th...
Definition amr_runtime.hpp:1069

pops::AmrRuntime::set_named_aux
void set_named_aux(int comp, std::vector< Real > field)
Registers a model-NAMED aux field (ADC-291) at shared-channel component comp (= kAuxNamedBase.
Definition amr_runtime.hpp:373

pops::AmrRuntime::coarse_local_boxes
int coarse_local_boxes() const
Definition amr_runtime.hpp:1231

pops::AmrRuntime::n_blocks
std::size_t n_blocks() const
Definition amr_runtime.hpp:302

pops::AmrRuntime::block_level_state
std::vector< double > block_level_state(std::size_t b, int k) const
Definition amr_runtime.hpp:1255

pops::AmrRuntime::named_field_values
std::vector< double > named_field_values(const std::string &field)
Solved potential of named field as a COARSE n*n row-major field (diagnostic / read-back).
Definition amr_runtime.hpp:428

pops::AmrRuntime::solve_named_fields
void solve_named_fields()
Solves every registered NAMED elliptic field (ADC-428) on the coarse, writes phi (+ centered grad) in...
Definition amr_runtime.hpp:678

pops::AmrRuntime::step
void step(Real dt)
Advances the system by one macro-step dt.
Definition amr_runtime.hpp:862

pops::AmrRuntime::nlev
int nlev() const
Definition amr_runtime.hpp:301

pops::AmrRuntime::has_named_field
bool has_named_field(const std::string &field) const
True if field is a registered named elliptic field.
Definition amr_runtime.hpp:421

pops::AmrRuntime::set_phi_tag_predicate
void set_phi_tag_predicate(TagPredicate crit)
Registers the PHI TAG PREDICATE (D4: SEPARATE phi criterion, on |grad phi|).
Definition amr_runtime.hpp:363

pops::AmrRuntime::macro_step
int macro_step() const
MACRO-STEP counter of the engine (regrid + hold-then-catch-up stride cadence: regrid when macro_step_...
Definition amr_runtime.hpp:1065

pops::AmrRuntime::phi
MultiFab & phi()
Definition amr_runtime.hpp:313

pops::AmrRuntime::block_cons_vars
const VariableSet & block_cons_vars(std::size_t b) const
Conservative VariableSet (names + physical roles, Model::conservative_vars()) of block b.
Definition amr_runtime.hpp:307

pops::AmrRuntime::add_dt_bound
void add_dt_bound(const std::string &label, std::function< double()> fn)
GLOBAL step bound (AMR counterpart of System::add_dt_bound): fn() evaluated once per step_cfl,...
Definition amr_runtime.hpp:1083

pops::AmrRuntime::last_dt_bound
const std::string & last_dt_bound() const
ACTIVE bound of the last step_cfl ("transport:<block>" / "source_frequency:<block>" / "stability_dt:<...
Definition amr_runtime.hpp:1168

pops::AmrRuntime::add_coupled_frequency_expr
void add_coupled_frequency_expr(const std::string &label, const std::vector< std::string > &in_blocks, const std::vector< std::string > &in_roles, const std::vector< double > &consts, const std::vector< int > &freq_prog_ops, const std::vector< int > &freq_prog_args)
PER-CELL COUPLED frequency (CoupledSource.frequency with an Expr, refinement of the CONSTANT frequenc...
Definition amr_runtime.hpp:1103

pops::AmrRuntime::set_regrid
void set_regrid(int every, int grow=2, int margin=2)
Activates the UNION-TAGS REGRID at the cadence every (in macro-steps): every every macro-steps,...
Definition amr_runtime.hpp:339

pops::AmrRuntime::n_named_fields
std::size_t n_named_fields() const
Number of registered named elliptic fields (diagnostic / test).
Definition amr_runtime.hpp:419

pops::AmrRuntime::set_block_level_state
void set_block_level_state(std::size_t b, int k, const std::vector< double > &s)
Definition amr_runtime.hpp:1281

pops::AmrRuntime::regrid_count
int regrid_count() const
Definition amr_runtime.hpp:323

pops::AmrRuntime::block_level_state_global
std::vector< double > block_level_state_global(std::size_t b, int k) const
Definition amr_runtime.hpp:1272

pops::AmrRuntime::set_level_potential
void set_level_potential(int k, const std::vector< double > &p)
Definition amr_runtime.hpp:1330

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

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::GeometricMG::solve
int solve(Real rel_tol, int max_cycles, Real abs_tol=Real(0))
Definition geometric_mg.hpp:407

pops::GeometricMG::rhs
MultiFab & rhs()
Definition geometric_mg.hpp:222

pops::GeometricMG::phi
MultiFab & phi()
Definition geometric_mg.hpp:221

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::ncomp
int ncomp() const
Number of components.
Definition multifab.hpp:60

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

pops::MultiFab::set_val
void set_val(Real v)
Fills all cells (valid + ghosts) of every local fab with v.
Definition multifab.hpp:85

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

pops::runtime::program::Profiler
Definition profiler.hpp:27

pops::runtime::program::Profiler::count
void count(const std::string &name, std::int64_t by=1)
Definition profiler.hpp:69

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

coupled_source_program.hpp
Generic COUPLED SOURCE interpreter: postfix bytecode evaluated on device (P5 phase 1).

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

fill_boundary.hpp
fill_boundary: INTRA-level halo exchange (fills ghosts from neighbors).

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

implicit_stepper.hpp
Implicit / IMEX block step as a named CONTRACT.

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

pops::detail::coupler_inject_aux_mb
void coupler_inject_aux_mb(const MultiFab &parent, MultiFab &child, bool replicated_parent=true)
Definition amr_coupler_mp.hpp:60

pops::detail::same_layout_or_throw
void same_layout_or_throw(const std::vector< std::vector< AmrLevelMP > > &block_levels)
Definition amr_system_coupler.hpp:115

pops
Definition amr_hierarchy.hpp:29

pops::kCflSpeedFloor
constexpr Real kCflSpeedFloor
Speed FLOOR for the CFL step policies (audit 2026-06, explicit constant instead of the scattered lite...
Definition types.hpp:38

pops::for_each_cell
void for_each_cell(const Box2D &b, F f)
Applies f to EACH cell (i, j) of box b (bounds inclusive), via Kokkos::parallel_for (Serial / OpenMP ...
Definition for_each.hpp:138

pops::kCsMaxTerms
constexpr int kCsMaxTerms
Definition coupled_source_program.hpp:43

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

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

pops::mf_average_down_mb
void mf_average_down_mb(const MultiFab &Uf, MultiFab &Uc)
Definition amr_subcycling.hpp:305

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

pops::aux_halo_override
BCRec aux_halo_override(const BCRec &shared, const AuxHaloPolicy &p)
Builds the effective override BCRec for a per-field aux halo: starts from the SHARED aux BC shared (s...
Definition physical_bc.hpp:204

pops::regrid_compute_fine_layout
std::pair< BoxArray, DistributionMapping > regrid_compute_fine_layout(TagBox grown, const Box2D &pdom, int pk, int margin, bool coarse_replicated=true)
Compute the fine layout (BoxArray + DistributionMapping) of a Berger-Rigoutsos regrid from the grown ...
Definition amr_regrid_coupler.hpp:53

pops::all_reduce_sum_inplace
void all_reduce_sum_inplace(double *, int)
Definition comm.hpp:155

pops::tag_union
TagBox tag_union(const std::vector< TagBox > &parts)
Union (cell-by-cell logical OR) of several TagBox sharing EXACTLY the same box.
Definition tag_box.hpp:74

pops::kCsMaxReg
constexpr int kCsMaxReg
Definition coupled_source_program.hpp:40

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

pops::reduce_max_cell
Real reduce_max_cell(const Box2D &b, F f)
MAX reduction with a REDUCING FUNCTOR: f receives (i, j, Real& acc) and updates acc,...
Definition for_each.hpp:240

pops::fill_physical_bc
void fill_physical_bc(MultiFab &mf, const Box2D &domain, const BCRec &bc)
Fills the physical-face ghosts of ALL components per bc (historical entry point, bit-identical).
Definition physical_bc.hpp:176

pops::all_reduce_max
double all_reduce_max(double x)
Definition comm.hpp:146

pops::roles_csv
std::string roles_csv(const VariableSet &vs)
CSV of a VariableSet's roles (role_name, separator ',').
Definition variables.hpp:172

pops::all_reduce_min
double all_reduce_min(double x)
Definition comm.hpp:149

pops::CsOp::PushReg
@ PushReg

pops::CsOp::Sqrt
@ Sqrt

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

pops::grow_tags
TagBox grow_tags(const TagBox &in, int n, const Box2D &domain)
Grows the tags by n cells (square neighborhood), staying within the domain.
Definition regrid.hpp:58

pops::regrid_field_on_layout
MultiFab regrid_field_on_layout(const BoxArray &fb, const DistributionMapping &dmap, const MultiFab &par, const MultiFab &old, int pk, int ngf, bool coarse_replicated=true)
Rebuild ONE fine MultiFab on the IMPOSED layout fb / dmap (the same one for all blocks in multi-block...
Definition amr_regrid_coupler.hpp:86

pops::kCsMaxProg
constexpr int kCsMaxProg
Definition coupled_source_program.hpp:42

pops::tag_cells
TagBox tag_cells(const MultiFab &mf, const Box2D &domain, Crit crit)
Marks the valid cells where the predicate is true, on a TagBox covering the domain.
Definition regrid.hpp:40

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

patch_box.hpp
PatchBox: index-space footprint of an AMR fine patch, exposed read-only to Python.

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

profiler.hpp

regrid.hpp
Dynamic regrid: (re)builds a fine level from the tagging of a coarse level.

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

pops::AmrRuntimeBlock
Type-erased closures of ONE AMR block, placed on the shared hierarchy.
Definition amr_runtime.hpp:98

pops::AmrRuntimeBlock::potential
std::function< std::vector< double >(const MultiFab &)> potential
Coarse potential read from the shared aux (component 0) as an n*n row-major field (diagnostic).
Definition amr_runtime.hpp:223

pops::AmrRuntimeBlock::stride
int stride
HOLD-THEN-CATCH-UP cadence of the block (multirate).
Definition amr_runtime.hpp:115

pops::AmrRuntimeBlock::name
std::string name
Definition amr_runtime.hpp:99

pops::AmrRuntimeBlock::newton_diagnostics
bool newton_diagnostics
NEWTON DIAGNOSTICS (OPT-IN, wave 3: AMR counterpart of System::newton_report).
Definition amr_runtime.hpp:185

pops::AmrRuntimeBlock::advance
std::function< void(std::vector< AmrLevelMP > &, const Box2D &, Real, Periodicity, bool)> advance
Advances the block by ONE substep of size dt: AMR transport (Berger-Oliger + conservative reflux.
Definition amr_runtime.hpp:141

pops::AmrRuntimeBlock::ncomp
int ncomp
Definition amr_runtime.hpp:100

pops::AmrRuntimeBlock::gamma
double gamma
Definition amr_runtime.hpp:101

pops::AmrRuntimeBlock::density
std::function< std::vector< double >()> density
Coarse density (component 0) of the block as a global n*n row-major field (diagnostic).
Definition amr_runtime.hpp:219

pops::AmrRuntimeBlock::levels
std::shared_ptr< std::vector< AmrLevelMP > > levels
Level stack of the block (level 0 = coarse, > 0 = fine patches), ON the shared layout.
Definition amr_runtime.hpp:133

pops::AmrRuntimeBlock::substeps
int substeps
EXPLICIT substeps of the block within ITS effective macro-step: the effective step (stride * dt) is s...
Definition amr_runtime.hpp:105

pops::AmrRuntimeBlock::project_per_level
std::function< void(std::vector< AmrLevelMP > &)> project_per_level
POINTWISE PROJECTION post-pas (ADC-177) : U <- project(U, aux) appliquee PAR NIVEAU a la FIN de l'ava...
Definition amr_runtime.hpp:174

pops::AmrRuntimeBlock::newton_report
std::shared_ptr< NewtonReport > newton_report
Definition amr_runtime.hpp:186

pops::AmrRuntimeBlock::max_speed
std::function< Real(const MultiFab &, const MultiFab &)> max_speed
Speed driving the block CFL on the coarse.
Definition amr_runtime.hpp:204

pops::AmrRuntimeBlock::add_elliptic_rhs
std::function< void(const MultiFab &, MultiFab &)> add_elliptic_rhs
Contribution of the block to the Poisson right-hand side: rhs += elliptic_rhs_b(U_b) on the coarse.
Definition amr_runtime.hpp:191

pops::AmrRuntimeBlock::source_frequency
std::function< Real(const MultiFab &, const MultiFab &)> source_frequency
OPTIONAL STEP BOUNDS of the block (AMR StabilityPolicy, audit 2026-06): evaluated on the COARSE (leve...
Definition amr_runtime.hpp:212

pops::AmrRuntimeBlock::named_elliptic_rhs
std::map< std::string, std::function< void(const MultiFab &, MultiFab &)> > named_elliptic_rhs
Per-NAMED-field elliptic right-hand-side contributions of the block (ADC-428): field name -> closure ...
Definition amr_runtime.hpp:199

pops::AmrRuntimeBlock::mass
std::function< Real()> mass
Mass of component 0 of the block coarse (sum u*dV; cross-rank reduced if distributed).
Definition amr_runtime.hpp:216

pops::AmrRuntimeBlock::aux_ncomp
int aux_ncomp
Width of the aux channel READ by the block model (aux_comps<Model>(); >= kAuxBaseComps).
Definition amr_runtime.hpp:119

pops::AmrRuntimeBlock::imex_advance
std::function< void(std::vector< AmrLevelMP > &, const Box2D &, Real, Periodicity, bool)> imex_advance
IMEX advance of the block by ONE substep of size dt: (1) EXPLICIT TRANSPORT on the SOURCE-FREE model ...
Definition amr_runtime.hpp:168

pops::AmrRuntimeBlock::imex
bool imex
TEMPORAL TREATMENT of the block: false (default) = EXPLICIT (forward-Euler source,...
Definition amr_runtime.hpp:149

pops::AmrRuntimeBlock::stability_dt
std::function< Real(const MultiFab &, const MultiFab &)> stability_dt
Definition amr_runtime.hpp:213

pops::AmrRuntimeBlock::cons_vars
VariableSet cons_vars
Descriptor of the model CONSERVATIVE variables (names + physical ROLES, Model::conservative_vars()).
Definition amr_runtime.hpp:127

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::AuxHaloPolicy
Per-field aux halo policy (ADC-369): a UNIFORM boundary policy for ONE model-named aux component,...
Definition physical_bc.hpp:195

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

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

pops::Box2D::from_extents
static Box2D from_extents(int nx, int ny)
Box [0, nx-1] x [0, ny-1] covering nx*ny cells from the index origin.
Definition box2d.hpp:42

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::CoupledFreqKernel
Definition coupled_source_program.hpp:175

pops::CoupledFreqKernel::prog
CsProgram prog
Definition coupled_source_program.hpp:183

pops::CoupledFreqKernel::consts
Real consts[kCsMaxReg]
Definition coupled_source_program.hpp:180

pops::CoupledFreqKernel::in_comp
int in_comp[kCsMaxReg]
Definition coupled_source_program.hpp:177

pops::CoupledFreqKernel::n_const
int n_const
Definition coupled_source_program.hpp:181

pops::CoupledFreqKernel::n_in
int n_in
Definition coupled_source_program.hpp:178

pops::CoupledFreqKernel::in
Array4 in[kCsMaxReg]
Definition coupled_source_program.hpp:176

pops::CoupledSourceKernel
Device functor applying ONE coupled source over a box: captures the PODs by VALUE (input/output Array...
Definition coupled_source_program.hpp:133

pops::CoupledSourceKernel::out
Array4 out[kCsMaxTerms]
Definition coupled_source_program.hpp:141

pops::CoupledSourceKernel::out_comp
int out_comp[kCsMaxTerms]
Definition coupled_source_program.hpp:142

pops::CoupledSourceKernel::n_terms
int n_terms
Definition coupled_source_program.hpp:144

pops::CoupledSourceKernel::dt
Real dt
Definition coupled_source_program.hpp:146

pops::CoupledSourceKernel::consts
Real consts[kCsMaxReg]
Definition coupled_source_program.hpp:138

pops::CoupledSourceKernel::in
Array4 in[kCsMaxReg]
Definition coupled_source_program.hpp:134

pops::CoupledSourceKernel::n_const
int n_const
Definition coupled_source_program.hpp:139

pops::CoupledSourceKernel::n_in
int n_in
Definition coupled_source_program.hpp:136

pops::CoupledSourceKernel::in_comp
int in_comp[kCsMaxReg]
Definition coupled_source_program.hpp:135

pops::CoupledSourceKernel::prog
CsProgram prog[kCsMaxTerms]
Definition coupled_source_program.hpp:143

pops::CsProgram
Fixed-capacity postfix program (POD device-copyable): len opcodes, arg read only by PushReg (register...
Definition coupled_source_program.hpp:47

pops::CsProgram::len
int len
Definition coupled_source_program.hpp:48

pops::CsProgram::arg
int arg[kCsMaxProg]
Definition coupled_source_program.hpp:50

pops::CsProgram::op
int op[kCsMaxProg]
Definition coupled_source_program.hpp:49

pops::FieldPostProcess
Definition elliptic_problem.hpp:73

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::NewtonReport
AGGREGATED report (whole block, all substeps of one advance) of the implicit-source Newton.
Definition implicit_stepper.hpp:165

pops::PatchBox
INDEX-SPACE footprint of an AMR fine patch, exposed to Python by AmrSystem::patch_boxes().
Definition patch_box.hpp:30

pops::Periodicity
Per-direction periodicity: halo wrapping in x and/or y during the exchange (false = open edge,...
Definition fill_boundary.hpp:37

pops::TagBox
Dense grid of 0/1 markers over a box, input to Berger-Rigoutsos clustering.
Definition tag_box.hpp:35

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

pops::VariableSet::roles
std::vector< VariableRole > roles
parallel to names; empty = roles not provided
Definition variables.hpp:62

pops::VariableSet::index_of
int index_of(VariableRole role) const
Index of the component carrying role (first occurrence), -1 if absent.
Definition variables.hpp:67

tag_box.hpp
TagBox: dense grid of markers (0/1) over a region, input to Berger-Rigoutsos clustering.

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