PoPS/doxygen/amr__coupler__mp_8hpp_source.html

#pragma once


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

#include <pops/amr/hierarchy/refinement_ratio.hpp>

#include <pops/coupling/amr/amr_diagnostics.hpp>     // amr_mass, amr_max_drift_speed

#include <pops/coupling/amr/amr_level_storage.hpp>   // AmrLevelStack

#include <pops/coupling/amr/amr_regrid_coupler.hpp>  // amr_regrid_finest (Berger-Rigoutsos)

#include <pops/coupling/single/coupler.hpp>  // detail::coupler_eval_rhs (f = model.elliptic_rhs(U))

#include <pops/numerics/elliptic/mg/composite_fac_poisson.hpp>  // COMPOSITE FAC 2-level Poisson solver (opt-in)

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

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

#include <pops/numerics/time/amr/reflux/amr_reflux_mf.hpp>  // AmrLevelMP, amr_step_multilevel_multipatch, mf_*_mb

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

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

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

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

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

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

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

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

#include <pops/mesh/layout/refinement.hpp>  // coarsen_index

#include <pops/parallel/comm.hpp>    // all_reduce_sum / all_reduce_max (distributed mass/drift)


#include <algorithm>   // std::max

#include <cmath>       // std::hypot

#include <cstddef>     // std::size_t

#include <functional>  // std::function (conducting-wall predicate passed to the MG)

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

#include <stdexcept>  // std::runtime_error (density size guard)

#include <utility>    // std::pair, std::move

#include <vector>


namespace pops {


namespace detail {

// Coupling primitive (not policy): piecewise-constant injection of aux from

// multi-box parent -> multi-box child (valid + ghosts). DISTRIBUTED, same scheme as

// mf_fill_fine_ghosts_mb. Two parent cases:

//  - REPLICATED (level 0, replicated_parent=true): parent fully local on each rank,

//    direct read via mf_find_box. This is the replicated coarse (per-rank dmap): parallel_copy

//    would violate the replicated-metadata assumption. Bit-identical path to the historical one.

//  - DISTRIBUTED (intermediate, replicated_parent=false): the parent may be on another rank,

//    we bring its valid regions onto a LOCAL child-coarsen grid via parallel_copy, then

//    inject. A child cell outside the parent coverage (GLOBAL box_array) is left

//    intact, like the replicated path (which skipped it via mf_find_box < 0).

// In serial both paths are identical (parent local everywhere, parallel_copy = memory copy).


inline void coupler_inject_aux_mb(const MultiFab& parent, MultiFab& child,

                                  bool replicated_parent = true) {

  const int nc = child.ncomp();

  if (replicated_parent) {

    device_fence();

    for (int lc = 0; lc < child.local_size(); ++lc) {

      Array4 c = child.fab(lc).array();

      const Box2D g = child.fab(lc).grown_box();

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

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

          const int ci = coarsen_index(i, kAmrRefRatio), cj = coarsen_index(j, kAmrRefRatio);

          const int pb = mf_find_box(parent, ci, cj);

          if (pb < 0)

            continue;

          const ConstArray4 pp = parent.fab(pb).const_array();

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

            c(i, j, k) = pp(ci, cj, k);

        }

    }

    return;

  }

  const BoxArray& pba = parent.box_array();  // GLOBAL: rank-independent coverage

  auto covered = [&](int ci, int cj) {

    for (int b = 0; b < pba.size(); ++b)

      if (pba[b].contains(ci, cj))

        return true;

    return false;

  };

  const BoxArray ccoarse = coarsen_grown(child.box_array(), child.n_grow(), kAmrRefRatio);

  MultiFab Pc(ccoarse, child.dmap(), parent.ncomp(), 0);

  parallel_copy(Pc, parent);  // parent regions (from any rank) -> local grid

  device_fence();

  for (int lc = 0; lc < child.local_size(); ++lc) {

    Array4 c = child.fab(lc).array();

    const ConstArray4 pp = Pc.fab(lc).const_array();

    const Box2D g = child.fab(lc).grown_box();

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

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

        const int ci = coarsen_index(i, kAmrRefRatio), cj = coarsen_index(j, kAmrRefRatio);

        if (!covered(ci, cj))

          continue;  // outside coverage -> keep the child value

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

          c(i, j, k) = pp(ci, cj, k);

      }

  }

}


// Writes an initial density (component 0, n*n row-major in GLOBAL indices) on the coarse

// level, MULTI-BOX and DISTRIBUTION-AWARE: each rank touches only its LOCAL fabs and reads

// rho at the cell GLOBAL index (i,j). For Euler (ncomp 4) it also sets zero momentum

// + thermal energy r/(gamma-1); ncomp 1 touches only density. Replicated mono-box:

// a single fab covering the domain, global == local indices -> bit-identical to the historical

// direct write. Distributed multi-box: each local box reads its window of rho.


inline void coupler_write_coarse(MultiFab& U, const std::vector<double>& rho, int n, int ncomp,

                                 double gamma) {

  if (static_cast<int>(rho.size()) != n * n)

    throw std::runtime_error("AMR coupler: initial density of size != n*n");

  const Real gm1 = Real(gamma) - Real(1);

  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) {

        const Real r = rho[static_cast<std::size_t>(j) * n + i];

        u(i, j, 0) = r;

        if (ncomp >= 3) {

          u(i, j, 1) = 0;

          u(i, j, 2) = 0;

        }

        if (ncomp == 4)

          u(i, j, 3) = r / gm1;

      }

  }

}


// Writes the FULL INITIAL CONSERVATIVE STATE (all components) on the coarse level from a

// flat component-major field @p state (c*n*n + j*n + i), of size ncomp*n*n. Counterpart of

// coupler_write_coarse for the multi-component seed: same box traversal (replicated mono-box

// AND distributed multi-box, GLOBAL indices (i,j)), only the per-cell write differs -- here we copy

// the ncomp components positionally (no density/momentum/energy wiring; the caller already provides

// the conservative, e.g. [rho, rho*u, rho*v]). gamma omitted (no energy derived). Index computed

// in std::size_t (no int overflow at large n, unlike the int validation of

// coupler_write_coarse). Used for the drift seed (set_conservative_state).


inline void coupler_write_coarse_state(MultiFab& U, const std::vector<double>& state, int n,

                                       int ncomp) {

  const std::size_t nn = static_cast<std::size_t>(n) * static_cast<std::size_t>(n);

  if (state.size() != nn * static_cast<std::size_t>(ncomp))

    throw std::runtime_error(

        "AMR coupler: initial state of size != ncomp*n*n (full conservative "

        "state; ncomp == model n_vars)");

  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 < ncomp; ++c)

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

                             static_cast<std::size_t>(j) * static_cast<std::size_t>(n) +

                             static_cast<std::size_t>(i)];

  }

}


// Reads the coarse density (component 0) into a GLOBAL n*n row-major field, MULTI-BOX and

// DISTRIBUTION-AWARE. Each rank writes its local cells into an n*n buffer initialized to 0

// then, if distributed, all_reduce_sum_inplace recomposes the full field on ALL ranks (the

// boxes are disjoint -> the cross-rank sum reconstructs the field exactly). Replicated mono-box:

// a single fab covers everything, the buffer is already complete, all_reduce would be the identity

// -> we avoid it (bit-identical to the historical direct read fab(0)).


inline std::vector<double> coupler_read_coarse(const MultiFab& U, int n, bool replicated) {

  device_fence();

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

  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)

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

  }

  if (!replicated)

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

  return out;

}


// Reads the coarse-level potential phi (component 0 of aux(0), written by compute_aux after the

// Poisson solve) into a GLOBAL n*n row-major field, MULTI-BOX and DISTRIBUTION-AWARE. aux(0) shares

// EXACTLY the layout of the coarse U (same BoxArray + DistributionMapping, cf. amr_level_storage:

// aux_[0] is built on U.box_array()/U.dmap()), so the recomposition is identical to

// coupler_read_coarse: local n*n buffer, all_reduce_sum if distributed (disjoint boxes -> exact

// sum), avoided in replicated mono-box (field already complete). PRECONDITION: update()/compute_aux

// has run at least once (otherwise aux(0) is 0). Strict counterpart of coupler_read_coarse for phi.


inline std::vector<double> coupler_read_coarse_phi(const MultiFab& aux0, int n, bool replicated) {

  device_fence();

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

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

    const ConstArray4 a = aux0.fab(li).const_array();

    const Box2D v = aux0.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) * n + i] = a(i, j, 0);

  }

  if (!replicated)

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

  return out;

}


// Injects the coarse into the valid cells of a fine patch (piecewise constant, ratio 2),

// MULTI-BOX and DISTRIBUTION-AWARE. Makes the hierarchy consistent before the first sync_down (the

// seed patch is at 0). Replicated mono-box: coarse fully local, direct read via

// mf_find_box (always found); no collective -> bit-identical to the historical fab(0).

// Distributed multi-box: we bring the needed coarse regions onto a LOCAL child-coarsen grid

// via parallel_copy (same scheme as coupler_inject_aux_mb), then inject.


inline void coupler_inject_coarse_to_fine_mb(const MultiFab& Uc, MultiFab& Uf, bool replicated) {

  const int nc = Uf.ncomp();

  if (replicated) {

    device_fence();

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

      Array4 f = Uf.fab(li).array();

      const Box2D v = Uf.box(li);

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

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

          const int ci = coarsen_index(i, kAmrRefRatio), cj = coarsen_index(j, kAmrRefRatio);

          const int pb = mf_find_box(Uc, ci, cj);

          if (pb < 0)

            continue;

          const ConstArray4 c = Uc.fab(pb).const_array();

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

            f(i, j, k) = c(ci, cj, k);

        }

    }

    return;

  }

  const BoxArray ccoarse = coarsen(Uf.box_array(), kAmrRefRatio);  // coarse footprint (valid cells)

  MultiFab Pc(ccoarse, Uf.dmap(), Uc.ncomp(), 0);

  parallel_copy(Pc, Uc);  // coarse regions (from any rank) -> local grid

  device_fence();

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

    Array4 f = Uf.fab(li).array();

    const ConstArray4 c = Pc.fab(li).const_array();

    const Box2D v = Uf.box(li);

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

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

        const int ci = coarsen_index(i, kAmrRefRatio), cj = coarsen_index(j, kAmrRefRatio);

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

          f(i, j, k) = c(ci, cj, k);

      }

  }

}


// Builds the coarse level (BoxArray + DistributionMapping) of the AmrSystem path according to the

// ownership policy, in a SINGLE point for both build paths (native + compiled):

//  - replicated (distribute=false, DEFAULT): mono-box covering the domain, dmap = my_rank() everywhere

//    (the box lives on each rank). In serial my_rank()=0 -> identical to round-robin, bit for bit.

//    This is the layout GeometricMG(replicated=true) and the historical one expect.

//  - distributed (distribute=true): multi-box BoxArray::from_domain(dom, max_grid) spread round-robin

//    DistributionMapping(ba.size(), n_ranks()). Each rank carries only its tiles -> the coarse

//    Poisson and coarse transport distribute (strong-scaling). max_grid<=0 => n/2 (2x2).


inline std::pair<BoxArray, DistributionMapping> coupler_make_coarse_layout(int n, bool distribute,

                                                                           int max_grid) {

  const Box2D dom = Box2D::from_extents(n, n);

  if (!distribute) {

    BoxArray ba(std::vector<Box2D>{dom});

    return {ba, DistributionMapping(std::vector<int>{my_rank()})};

  }

  const int mg = (max_grid > 0) ? max_grid : (n / 2 > 0 ? n / 2 : n);

  BoxArray ba = BoxArray::from_domain(dom, mg);

  return {ba, DistributionMapping(ba.size(), n_ranks())};

}


}  // namespace detail


template <class Model, class Elliptic = GeometricMG>


class AmrCouplerMP {

  static_assert(EllipticSolver<Elliptic>, "Elliptic must model EllipticSolver");


 public:

  // active: optional "active cell" predicate (interior of the conductor), for the circular

  // conducting wall of the column instability (passed as-is to the multigrid). Empty

  // by default -> no wall (historical behavior unchanged). Only the coarse carries the

  // wall: the fine patches refine the ring edge, strictly inside the wall.

  // replicated_coarse: level-0 (coarse) OWNERSHIP POLICY. BOTH modes are

  // stable and their equivalence is proven bit for bit (test_mpi_decoarse, maxdiff=0):

  //   true  (performant DEFAULT): coarse mono-box REPLICATED on all ranks. Best coarse

  //          MG solve (no multigrid degeneration), zero communication for the

  //          coarse Poisson, robust reference -> the right default for small/medium cases.

  //   false (EXPLICIT scalable mode): coarse multi-box DISTRIBUTED round-robin. Lifts the

  //          O(NX*NY*nranks) memory lock of level 0, required at very large scale. But the

  //          geometric MG degenerates for a finely-split coarse (>2x2 boxes do not tile the

  //          coarsest grid): reserve for cases where the level-0 memory is the lock.

  // Criterion: set false ONLY when memory scalability requires it; otherwise keep true.

  // Removing the replicated path is DEFERRED as long as the distributed one is not strictly

  // superior. mg_ receives the same flag (otherwise, under replicated MPI, the coarse would fall on

  // the single rank 0 and compute_aux would read a phi absent elsewhere). In serial, both coincide.


  AmrCouplerMP(const Model& model, const Geometry& geom, const BoxArray& ba_coarse, const BCRec& bc,

               std::vector<AmrLevelMP> levels, std::function<bool(Real, Real)> active = {},

               bool replicated_coarse = true)

      : model_(model),

        geom_(geom),

        mg_(geom, ba_coarse, bc, std::move(active), replicated_coarse),

        stack_(geom.domain, std::move(levels), aux_comps<Model>()),

        replicated_coarse_(replicated_coarse) {}


  std::vector<AmrLevelMP>& levels() { return stack_.levels(); }

  MultiFab& coarse() { return stack_.coarse(); }

  const MultiFab& coarse() const { return stack_.coarse(); }

  // coarse-level aux: (phi, dphi/dx, dphi/dy), component 0 = phi (cf. compute_aux). Same

  // layout as coarse(). Read by the AmrSystem potential hook (coupler_read_coarse_phi).

  MultiFab& aux0() { return stack_.aux(0); }

  const MultiFab& aux0() const { return stack_.aux(0); }


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

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

    apply_named_aux();  // stack_ exists at ctor: reflect onto the coarse aux right away

  }


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

  const Box2D& domain() const { return stack_.domain(); }

  int nlev() const { return stack_.nlev(); }


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

  // SINGLE-RANK AMR CHECKPOINT / RESTART (ADC-65). The mono-block coupler carries the FULL

  // CONSERVATIVE STATE per level (all components) + the phi (multigrid warm-start), and can IMPOSE

  // a SAVED fine hierarchy (instead of Berger-Rigoutsos clustering on tags). SINGLE-RANK: the

  // accessors loop over local_size() + device_fence(), WITHOUT MPI gather (the facade rejects np>1

  // AND multi-block upstream; multi-rank/multi-block restart is a documented follow-up).

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


  // Reads the FULL conservative state (all components) of level @p k into a flat

  // component-major field c*nf*nf + j*nf + i, nf = n << k (n = coarse side). The cells OUTSIDE

  // patches (uncovered fine level) stay at 0: a fine level is only defined within its patches

  // (at restart we rewrite ONLY the patch cells, cf. set_level_state).


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

    std::vector<AmrLevelMP>& L = stack_.L();

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

      throw std::runtime_error("AmrCouplerMP::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>(stack_.domain().nx()) << k;

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

    device_fence();

    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);

    }

    return out;

  }


  // Restores the full conservative state of level @p k from @p s (same layout as level_state).

  // Writes ONLY the VALID cells of the local fabs (the patches): the ghosts are redone at the

  // next update()/advance (exactly like after a regrid), and a fine cell outside a patch

  // does not exist. NO RE-PROLONGATION: the state is restored AS-IS (no coarse->fine injection).


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

    std::vector<AmrLevelMP>& L = stack_.L();

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

      throw std::runtime_error("AmrCouplerMP::set_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>(stack_.domain().nx()) << k;

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

      throw std::runtime_error("AmrCouplerMP::set_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)];

    }

  }


  // Reads the potential phi of level @p k, flat nf*nf row-major field (nf = n << k), zeros outside patches.

  // Level 0: the multigrid WARM-START -- mg_.phi() (VALID cells), the state actually

  // reused by the NEXT solve (GeometricMG::solve keeps phi between calls). Level >= 1:

  // aux(k) component 0 (informational; recomputed at update). It is mg_.phi() level 0 that makes the

  // restart BIT-IDENTICAL (the 1st post-restart solve starts from the same guess as the continuous run).


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

    if (k < 0 || k >= stack_.nlev())

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

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

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

    device_fence();

    const MultiFab& P = (k == 0) ? mg_.phi() : stack_.aux(k);

    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);

    }

    return out;

  }


  // Restores the potential of level @p k. Level 0: warm-start mg_.phi() (valid cells) -> the

  // multigrid restart is BIT-IDENTICAL (the 1st post-restart solve starts from the same guess). Level

  // >= 1: aux(k) comp 0 (recomputed at update; idempotent restore, no effect on the dynamics).


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

    if (k < 0 || k >= stack_.nlev())

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

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

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

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

    device_fence();

    MultiFab& P = (k == 0) ? mg_.phi() : stack_.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)];

    }

  }


  // GLOBAL (np>1 GATHER) variants of level_state / level_potential (ADC-509). The base accessors

  // above fill the GLOBAL buffer from the LOCAL fabs only (zeros where this rank owns no box); these

  // add an all_reduce_sum_inplace so EACH rank holds the complete field (AMR reflux pattern, comm.hpp;

  // MIRROR of System::state_global / gather_global). COLLECTIVE: all ranks MUST call them. Mono-rank

  // the reduce is the identity -> bit-identical to level_state / level_potential. The base distributed

  // coarse + round-robin fine patches are disjoint, so the sum double-counts nothing.


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

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

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

    return out;

  }


  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;

  }


  // Imposes the fine-level hierarchy (restart): rebuilds level 1 on the SAVED @p

  // fine_boxes BoxArray (instead of Berger-Rigoutsos clustering on tags), via the SAME mechanism as

  // regrid (regrid_field_on_layout: parent interp + fine carry-over), then reattaches the level-1 aux.

  // The rebuilt valid content is OVERWRITTEN afterwards by set_level_state (restore as-is): here we

  // rely only on the IMPOSED LAYOUT. SINGLE-RANK, 2-level mono-block hierarchy (so we impose

  // ONLY level 1). Clear rejection if the hierarchy has no fine level or if no box was saved.


  void set_hierarchy(const std::vector<Box2D>& fine_boxes) {

    std::vector<AmrLevelMP>& L = stack_.L();

    if (L.size() < 2)

      throw std::runtime_error(

          "AmrCouplerMP::set_hierarchy: mono-level hierarchy (no fine patch "

          "to impose)");

    if (fine_boxes.empty())

      throw std::runtime_error(

          "AmrCouplerMP::set_hierarchy: no saved fine box (restart of a "

          "fine-patch hierarchy required)");

    const int ngf = L[1].U.n_grow();  // inherit the ghost width of the current fine (scheme parity)

    BoxArray fb(fine_boxes);

    DistributionMapping dmap(static_cast<int>(fb.size()),

                             n_ranks());  // single-rank -> all on rank 0

    L[1].U = regrid_field_on_layout(fb, dmap, L[0].U, L[1].U, /*pk=*/0, ngf, replicated_coarse_);

    stack_.reattach_aux(1);  // realloc aux[1] on the new layout + rewire L[1].aux

  }


  void sync_down() {  // average fine -> coarse over the whole hierarchy (multi-box)

    auto& L = stack_.L();

    for (int k = stack_.nlev() - 1; k >= 1; --k)

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

  }


  void set_composite_poisson(bool v) { composite_poisson_ = v; }

  bool composite_poisson() const { return composite_poisson_; }


  void compute_aux() {  // coarse Poisson + grad phi + injection to the fine levels

    auto& L = stack_.L();

    const Box2D& dom = stack_.domain();

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

    // COMPOSITE path (opt-in): the fine patch TRULY refines the elliptic. Supported scope = 2 levels,

    // ONE mono-box fine patch, replicated coarse (Phase 2). Otherwise Option A below (bit-identical).

    if (composite_poisson_ && replicated_coarse_ && stack_.nlev() == 2 &&

        L[1].U.box_array().size() == 1) {

      compute_aux_composite();

      return;

    }

    // right-hand side via the model (no copied formula): f = elliptic_rhs(U)

    detail::coupler_eval_rhs(L[0].U, mg_.rhs(), model_);

    mg_.solve();  // leaves phi with its ghosts filled (last gs_rb_sweep -> fill_ghosts)

    device_fence();

    // aux = (phi, grad phi) per LOCAL coarse fab: covers the replicated mono-box (1 fab) as well as

    // the distributed multi-box (de-replication). The box-edge derivatives read the ghosts of

    // phi filled by the solve (distributed inter-box exchange via fill_boundary). mg_.phi() and

    // aux(0) share the same layout (same BoxArray + DistributionMapping) -> fab(li) <-> box(li).

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

      const ConstArray4 p = mg_.phi().fab(li).const_array();

      Array4 a = stack_.aux(0).fab(li).array();

      const Box2D b = stack_.aux(0).box(li);

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

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

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

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

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

        }

    }

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

    // fill_boundary (ghosts) and the injection (lines below), so they reach every level and survive a

    // regrid. The loop above wrote only comps 0..2 (phi/grad), so this never clobbers them. No-op

    // without a named field (bit-identical).

    apply_named_aux();

    fill_boundary(stack_.aux(0), dom, Periodicity{true, true});

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

                           // shared fill); no-op on a periodic domain / without a policy.

    // parent aux(k-1) replicated only if level 0 is: otherwise it is DISTRIBUTED (multi-box)

    // and the injection goes through parallel_copy. Beyond level 1, the parent is always distributed.

    for (int k = 1; k < stack_.nlev(); ++k)

      detail::coupler_inject_aux_mb(stack_.aux(k - 1), stack_.aux(k),

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

  }


  void update() {

    sync_down();

    compute_aux();

  }


  // Selectable spatial discretization (default FirstOrder = NoSlope + Rusanov,

  // strictly identical to the old step()). recon_prim selects the primitive

  // reconstruction (same parameter as assemble_rhs / System); false (default) -> conservative.

  // imex: treats the stiff source IMPLICITLY (backward_euler) rather than forward Euler;

  // false (default) -> historical explicit treatment, bit-identical. The source being

  // cell-local (outside reflux registers), the implicit split preserves conservation.

  template <class Disc = FirstOrder>


  void step(Real dt, bool recon_prim = false, bool imex = false, const NewtonOptions& nopts = {},

            AmrTimeMethod tmethod = AmrTimeMethod::kEuler, Real pos_floor = Real(0)) {

    update();

    advance_amr<typename Disc::Limiter, typename Disc::NumericalFlux>(

        model_, stack_.L(), stack_.domain(), dt, Periodicity{true, true}, replicated_coarse_,

        recon_prim, imex, nopts, tmethod, pos_floor);

  }


  template <class Disc = FirstOrder>


  void advance_transport(Real dt, bool recon_prim = false, Real pos_floor = Real(0)) {

    advance_amr<typename Disc::Limiter, typename Disc::NumericalFlux>(

        model_, stack_.L(), stack_.domain(), dt, Periodicity{true, true}, replicated_coarse_,

        recon_prim, /*imex=*/false, NewtonOptions{}, AmrTimeMethod::kEuler, pos_floor);

  }


  // Regrid of the FINE level by Berger-Rigoutsos (delegated to amr_regrid_finest):

  // rebuilds the patches (carry over fine data, otherwise parent interp) + the aux.

  // margin = nesting. The coupler only orders the call.

  template <class Crit>


  void regrid(Crit crit, int grow = 2, int margin = 2) {

    amr_regrid_finest(stack_.L(), stack_.aux(), stack_.domain(), crit, grow, margin,

                      aux_comps<Model>(), replicated_coarse_);

  }


  // coarse mass via the shared diagnostic amr_mass_mb (replicated mono-box as well as

  // distributed multi-box). Replicated coarse: the local sum IS already the total mass

  // (each rank holds everything) -> no all_reduce. Distributed: local part -> all_reduce_sum.


  Real mass() const {

    const Real M = amr_mass_mb(stack_.coarse(), geom_.dx(), geom_.dy());

    return replicated_coarse_ ? M : all_reduce_sum(M);

  }


  // max drift speed via amr_max_drift_speed_mb + floor. all_reduce_max correct

  // in BOTH cases: under replication the local max is already global (idempotent);

  // distributed, we take the max of the parts.


  Real max_drift_speed() const {

    const Real v = amr_max_drift_speed_mb(stack_.aux(0), model_.B0);

    return all_reduce_max(std::max(v, Real(1e-12)));

  }


  Real max_wave_speed() {

    Real w = Real(1e-12);

    MultiFab& U = stack_.coarse();

    MultiFab& A = stack_.aux(0);

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

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

      const ConstArray4 a = A.fab(li).const_array();

      const Box2D b = U.box(li);

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

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

          const auto us = load_state<Model>(u, i, j);

          const Aux ax = load_aux<aux_comps<Model>()>(a, i, j);

          w = std::max(

              w, std::max(model_.max_wave_speed(us, ax, 0), model_.max_wave_speed(us, ax, 1)));

        }

    }

    return all_reduce_max(w);

  }


 private:

  void compute_aux_composite() {

    // ADC-291 NOTE: unlike compute_aux (Option A), this opt-in composite-FAC path does NOT re-apply

    // named aux onto the fine level (it derives each level's aux from the FAC phi, with no coarse->fine

    // aux injection). The coarse named comp survives (the grad writes touch only comps 0..2), but a

    // fine-level model reading extra_field(k) would read 0. set_composite_poisson is C++-only and not

    // facade-reachable, so named aux cannot hit this path today; carrying named aux to the composite

    // fine level is a documented follow-up (cf. pops.capabilities()['aux']['followups']).

    auto& L = stack_.L();

    const Box2D& dom = stack_.domain();

    const Box2D fine_box = L[1].U.box_array()[0];

    if (!fac_built_ || !same_box(fac_fine_box_, fine_box)) {

      fac_ = std::make_shared<CompositeFacPoisson>(geom_, mg_.box_array(), mg_.bc(), fine_box, 2);

      fac_fine_box_ = fine_box;

      fac_built_ = true;

    }

    // f = elliptic_rhs(U) PER LEVEL: the fine has its OWN refined right-hand side (not an injection).

    detail::coupler_eval_rhs(L[0].U, fac_->rhs_coarse(), model_);

    detail::coupler_eval_rhs(L[1].U, fac_->rhs_fine(), model_);

    fac_->solve();

    device_fence();

    // level-0 aux (coarse): phi + grad from phi_coarse (same centered stencils as the Option A path).

    fill_ghosts(fac_->phi_coarse(), dom, mg_.bc());

    detail::coupler_grad_phi(fac_->phi_coarse(), stack_.aux(0), Real(1) / (Real(2) * geom_.dx()),

                             Real(1) / (Real(2) * geom_.dy()));

    fill_boundary(stack_.aux(0), dom, Periodicity{true, true});

    // level-1 aux (fine): phi + grad from phi_fine -> FINE grad (fine centered diff, reads the C-F

    // bilinear ghosts) = the fidelity gain vs the constant coarse grad injected by Option A.

    detail::coupler_grad_phi(fac_->phi_fine(), stack_.aux(1), Real(1) / (Real(2) * L[1].dx),

                             Real(1) / (Real(2) * L[1].dy));

  }


  static bool same_box(const Box2D& a, const Box2D& b) {

    return a.lo[0] == b.lo[0] && a.lo[1] == b.lo[1] && a.hi[0] == b.hi[0] && a.hi[1] == b.hi[1];

  }


  Model model_;

  Geometry geom_;

  Elliptic mg_;

  AmrLevelStack<AmrLevelMP> stack_;

  bool

      replicated_coarse_;  // level 0 replicated (true) or distributed multi-box (false, de-replication)

  // COMPOSITE FAC Poisson path (opt-in, set_composite_poisson). fac_ built lazily on the

  // current fine patch (rebuilt if the patch changes after regrid). Default OFF -> Option A bit-identical.

  bool composite_poisson_ = false;

  bool fac_built_ = false;

  std::shared_ptr<CompositeFacPoisson> fac_;

  Box2D fac_fine_box_{};

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

  // (n*n row-major). STATIC user fields re-applied by compute_aux each update (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_;


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

  // SystemFieldSolver::apply_named_aux_one and AmrRuntime::apply_named_aux: per local fab (MPI-safe),

  // valid cells only, global flat index j*nx+i. compute_aux runs the coarse->fine injection right

  // after, carrying the named comps to the fine levels. No-op without a named field.

  void apply_named_aux() {

    if (named_aux_.empty())

      return;

    const int row = stack_.domain().nx();

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

      if (field.empty() || comp >= stack_.aux(0).ncomp())

        continue;

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

        Array4 a = stack_.aux(0).fab(li).array();

        const Box2D v = stack_.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];

      }

    }

  }


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

  // declared component's physical-face ghosts; aux_halo_override(mg_.bc(), policy) keeps periodic faces

  // periodic (so on a periodic domain this is a no-op). Mirror of SystemFieldSolver::apply_named_aux_bc.

  void apply_named_aux_bc() {

    if (named_aux_bc_.empty())

      return;

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

      if (comp >= stack_.aux(0).ncomp())

        continue;

      fill_physical_bc(stack_.aux(0), stack_.domain(), aux_halo_override(mg_.bc(), policy), comp);

    }

  }

};


}  // namespace pops

amr_diagnostics.hpp
Diagnostics extracted from the AMR couplers: mass and max drift speed (responsibility c).

amr_level_storage.hpp
AmrLevelStack: AMR hierarchy storage (levels + aux) extracted from the couplers.

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

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::AmrCouplerMP
Multi-patch E x B AMR coupler.
Definition amr_coupler_mp.hpp:276

pops::AmrCouplerMP::level_potential
std::vector< double > level_potential(int k)
Definition amr_coupler_mp.hpp:394

pops::AmrCouplerMP::domain
const Box2D & domain() const
Definition amr_coupler_mp.hpp:328

pops::AmrCouplerMP::coarse
MultiFab & coarse()
Definition amr_coupler_mp.hpp:307

pops::AmrCouplerMP::set_level_state
void set_level_state(int k, const std::vector< double > &s)
Definition amr_coupler_mp.hpp:368

pops::AmrCouplerMP::set_composite_poisson
void set_composite_poisson(bool v)
OPT-IN: replaces the Option A AMR Poisson (coarse solve + piecewise-constant gradient injection) with...
Definition amr_coupler_mp.hpp:482

pops::AmrCouplerMP::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.
Definition amr_coupler_mp.hpp:327

pops::AmrCouplerMP::levels
std::vector< AmrLevelMP > & levels()
Definition amr_coupler_mp.hpp:306

pops::AmrCouplerMP::composite_poisson
bool composite_poisson() const
Definition amr_coupler_mp.hpp:483

pops::AmrCouplerMP::set_hierarchy
void set_hierarchy(const std::vector< Box2D > &fine_boxes)
Definition amr_coupler_mp.hpp:454

pops::AmrCouplerMP::level_potential_global
std::vector< double > level_potential_global(int k)
Definition amr_coupler_mp.hpp:442

pops::AmrCouplerMP::AmrCouplerMP
AmrCouplerMP(const Model &model, const Geometry &geom, const BoxArray &ba_coarse, const BCRec &bc, std::vector< AmrLevelMP > levels, std::function< bool(Real, Real)> active={}, bool replicated_coarse=true)
Definition amr_coupler_mp.hpp:297

pops::AmrCouplerMP::step
void step(Real dt, bool recon_prim=false, bool imex=false, const NewtonOptions &nopts={}, AmrTimeMethod tmethod=AmrTimeMethod::kEuler, Real pos_floor=Real(0))
Advances the hierarchy by one step dt: update() then advance_amr (Berger-Oliger subcycling + reflux +...
Definition amr_coupler_mp.hpp:556

pops::AmrCouplerMP::set_level_potential
void set_level_potential(int k, const std::vector< double > &p)
Definition amr_coupler_mp.hpp:414

pops::AmrCouplerMP::nlev
int nlev() const
Definition amr_coupler_mp.hpp:329

pops::AmrCouplerMP::advance_transport
void advance_transport(Real dt, bool recon_prim=false, Real pos_floor=Real(0))
TRANSPORT-ONLY ADVANCE (hyperbolic), WITHOUT update() or source.
Definition amr_coupler_mp.hpp:573

pops::AmrCouplerMP::aux0
const MultiFab & aux0() const
Definition amr_coupler_mp.hpp:312

pops::AmrCouplerMP::compute_aux
void compute_aux()
Definition amr_coupler_mp.hpp:485

pops::AmrCouplerMP::sync_down
void sync_down()
Definition amr_coupler_mp.hpp:472

pops::AmrCouplerMP::max_drift_speed
Real max_drift_speed() const
Definition amr_coupler_mp.hpp:599

pops::AmrCouplerMP::mass
Real mass() const
Definition amr_coupler_mp.hpp:591

pops::AmrCouplerMP::aux0
MultiFab & aux0()
Definition amr_coupler_mp.hpp:311

pops::AmrCouplerMP::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_coupler_mp.hpp:319

pops::AmrCouplerMP::level_state
std::vector< double > level_state(int k)
Definition amr_coupler_mp.hpp:343

pops::AmrCouplerMP::coarse
const MultiFab & coarse() const
Definition amr_coupler_mp.hpp:308

pops::AmrCouplerMP::max_wave_speed
Real max_wave_speed()
Max wave speed on the coarse level via model.max_wave_speed.
Definition amr_coupler_mp.hpp:612

pops::AmrCouplerMP::regrid
void regrid(Crit crit, int grow=2, int margin=2)
Definition amr_coupler_mp.hpp:583

pops::AmrCouplerMP::update
void update()
Updates the hierarchy before a step: sync_down (fine -> coarse) then compute_aux (coarse Poisson + gr...
Definition amr_coupler_mp.hpp:532

pops::AmrCouplerMP::level_state_global
std::vector< double > level_state_global(int k)
Definition amr_coupler_mp.hpp:437

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

pops::BoxArray::size
int size() const
Number of boxes in the tiling.
Definition box_array.hpp:42

pops::BoxArray::from_domain
static BoxArray from_domain(const Box2D &domain, int max_grid_size)
Tile the domain into tiles of at most max_grid_size per direction, distributed evenly.
Definition box_array.hpp:30

pops::DistributionMapping
Owning MPI rank of each box, indexed by GLOBAL box index (parallel to a BoxArray).
Definition distribution_mapping.hpp:19

pops::Fab2D::grown_box
const Box2D & grown_box() const
Grown box (valid + ng ghosts) = actual memory footprint.
Definition fab2d.hpp:76

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

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

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

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

pops::MultiFab::ncomp
int ncomp() const
Number of components.
Definition multifab.hpp:60

pops::MultiFab::dmap
const DistributionMapping & dmap() const
GLOBAL distribution (owner rank per box).
Definition multifab.hpp:58

pops::MultiFab::n_grow
int n_grow() const
Number of ghost layers.
Definition multifab.hpp:62

pops::MultiFab::box_array
const BoxArray & box_array() const
GLOBAL decomposition of the level (all boxes, all ranks).
Definition multifab.hpp:56

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

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

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

composite_fac_poisson.hpp
CompositeFacPoisson: 2-level AMR COMPOSITE elliptic solver (Fast Adaptive Composite,...

pops::EllipticSolver
Definition elliptic_solver.hpp:30

coupler.hpp
Coupler: single-block hyperbolic-elliptic coupler (Poisson -> aux -> advance loop).

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

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

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

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

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

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

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::coupler_read_coarse
std::vector< double > coupler_read_coarse(const MultiFab &U, int n, bool replicated)
Definition amr_coupler_mp.hpp:169

pops::detail::coupler_make_coarse_layout
std::pair< BoxArray, DistributionMapping > coupler_make_coarse_layout(int n, bool distribute, int max_grid)
Definition amr_coupler_mp.hpp:257

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

pops::detail::coupler_read_coarse_phi
std::vector< double > coupler_read_coarse_phi(const MultiFab &aux0, int n, bool replicated)
Definition amr_coupler_mp.hpp:191

pops::detail::coupler_write_coarse
void coupler_write_coarse(MultiFab &U, const std::vector< double > &rho, int n, int ncomp, double gamma)
Definition amr_coupler_mp.hpp:112

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

pops::detail::coupler_inject_coarse_to_fine_mb
void coupler_inject_coarse_to_fine_mb(const MultiFab &Uc, MultiFab &Uf, bool replicated)
Definition amr_coupler_mp.hpp:212

pops::detail::coupler_write_coarse_state
void coupler_write_coarse_state(MultiFab &U, const std::vector< double > &state, int n, int ncomp)
Definition amr_coupler_mp.hpp:143

pops
Definition amr_hierarchy.hpp:29

pops::mf_find_box
int mf_find_box(const MultiFab &mf, int I, int J)
Definition amr_subcycling.hpp:231

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::all_reduce_sum
double all_reduce_sum(double x)
Definition comm.hpp:143

pops::coarsen
BoxArray coarsen(const BoxArray &ba, int r)
Coarsens each box of the BoxArray by a ratio r (coarsen box by box, order preserved).
Definition refinement.hpp:39

pops::coarsen_grown
BoxArray coarsen_grown(const BoxArray &ba, int ngrow, int r)
Definition amr_subcycling.hpp:241

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

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

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::fill_boundary
void fill_boundary(MultiFab &mf, const Box2D &domain, Periodicity per={})
BLOCKING halo exchange: begin then end immediately (no overlap).
Definition fill_boundary.hpp:333

pops::amr_max_drift_speed_mb
Real amr_max_drift_speed_mb(const MultiFab &aux0, Real B0)
LOCAL max drift speed: max of |grad phi| / B0 (aux comp 1, 2 = grad phi) over the valid cells,...
Definition amr_diagnostics.hpp:48

pops::amr_mass_mb
Real amr_mass_mb(const MultiFab &coarse, Real dx, Real dy)
LOCAL mass: sum of u(.,.,0) * dx * dy over the valid cells of ALL local fabs, WITHOUT MPI reduction (...
Definition amr_diagnostics.hpp:33

pops::my_rank
int my_rank()
Definition comm.hpp:136

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::AmrTimeMethod
AmrTimeMethod
Definition amr_flux_helpers.hpp:46

pops::AmrTimeMethod::kEuler
@ kEuler

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

pops::kAmrRefRatio
constexpr int kAmrRefRatio
The native AMR refinement ratio between two consecutive levels.
Definition refinement_ratio.hpp:30

pops::amr_regrid_finest
void amr_regrid_finest(std::vector< AmrLevelMP > &L, std::vector< MultiFab > &aux, const Box2D &dom, Crit crit, int grow, int margin, int aux_ncomp=kAuxBaseComps, bool coarse_replicated=true)
Regrid the finest level (L.back()) by Berger-Rigoutsos on the criterion crit applied to the parent: r...
Definition amr_regrid_coupler.hpp:151

pops::parallel_copy
void parallel_copy(MultiFab &dst, const MultiFab &src)
Copies the valid regions that OVERLAP from src to dst (same indices, no shift).
Definition refinement.hpp:50

pops::coarsen_index
POPS_HD int coarsen_index(int a, int r)
Index of the coarse cell containing the fine cell a (FLOOR division by r, handles a < 0).
Definition refinement.hpp:34

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::fill_ghosts
void fill_ghosts(MultiFab &mf, const Box2D &domain, const BCRec &bc)
COMPLETE ghost filling: fill_boundary (interior + periodic, periodicity deduced from bc) THEN fill_ph...
Definition physical_bc.hpp:227

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

refinement.hpp
AMR inter-level transfer operators (integer ratio r) + parallel_copy.

refinement_ratio.hpp
Single source of truth for the native AMR refinement ratio.

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::Aux
POINTWISE auxiliary fields shared with the physics: single coupling channel.
Definition state.hpp:123

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::hi
int hi[2]
Definition box2d.hpp:39

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::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::NewtonOptions
Options of the local Newton of the implicit source (backward-Euler).
Definition implicit_stepper.hpp:114

pops::Periodicity
Per-direction periodicity: halo wrapping in x and/or y during the exchange (false = open edge,...
Definition fill_boundary.hpp:37

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