PoPS/doxygen/system__field__solver_8hpp_source.html

#pragma once


#include <pops/core/state/state.hpp>        // kAuxBaseComps (base component of the aux channel)

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

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

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

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

#include <pops/mesh/boundary/physical_bc.hpp>  // BCRec, fill_ghosts, fill_boundary

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

#include <pops/numerics/elliptic/poisson/poisson_fft_solver.hpp>

#include <pops/numerics/elliptic/polar/polar_poisson_solver.hpp>  // PolarPoissonSolver (direct polar Poisson)

#include <pops/parallel/comm.hpp>                           // n_ranks() (FFT MPI guard)

#include <pops/runtime/builders/block/block_builder_polar.hpp>  // derive_aux_polar (polar aux in local basis)

#include <pops/runtime/context/wall_predicate.hpp>       // detail::wall_predicate


#include <cstdio>   // POPS_TRACE_SOLVE_FIELDS: device diagnostic trace (env-gated, inert by default)

#include <cstdlib>  // getenv

#include <functional>

#include <map>  // named_aux_: NAMED aux fields (comp -> field), re-applied after channel realloc

#include <optional>

#include <stdexcept>

#include <string>

#include <type_traits>

#include <variant>

#include <vector>


namespace pops {


namespace field_solver {


inline bool pops_trace_sf() {

  static const bool on = std::getenv("POPS_TRACE_SOLVE_FIELDS") != nullptr;

  return on;

}


inline void pops_sf_mark(const char* w) {

  if (pops_trace_sf()) {

    std::fprintf(stderr, "[sf] %s\n", w);

    std::fflush(stderr);

  }

}


template <class Impl>


class SystemFieldSolver {

 public:

  explicit SystemFieldSolver(Impl* owner) : owner_(owner) {}


  static constexpr int kTeComp = kAuxBaseComps + 1;  // = 4


  bool provides_aux(const std::string& name) const {

    if (name == "B_z") {

      return !bz_field_.empty();

    }

    if (name == "T_e") {

      return te_src_ >= 0;

    }

    return true;

  }


  // --- OWNED state (elliptic solve + coefficient fields + application buffers) --------

  // Poisson configuration (elliptic solver built lazily).

  std::string p_rhs = "charge_density";

  std::string p_solver = "geometric_mg";

  std::string p_bc = "auto";

  std::string p_wall = "none";

  double p_wall_radius = 0.0;

  Real p_eps_ = 1;  // CONSTANT permittivity: div(eps grad phi) = f <=> lap phi = f/eps

  // ABSOLUTE floor of the GeometricMG V-cycle stopping criterion (same units as the residual). Default 0:

  // purely relative criterion (historical bit-identical behavior). Set > 0 (problem scale),

  // it makes the off-step solve_fields exit WITHOUT cycling on an already converged state (diagnostics, oracles,

  // restart). Inert for the FFT solver (direct, without iterative tolerance).

  Real p_abs_tol_ = 0;

  bool has_eps_field_ = false;  // VARIABLE permittivity eps(x) provided (carried by the operator)

  std::vector<double> p_eps_field_;  // field eps(x), n*n row-major (if has_eps_field_)

  bool has_eps_xy_field_ =

      false;  // ANISOTROPIC permittivity eps_x(x), eps_y(x) (operator div(diag(eps_x,eps_y) grad phi))

  std::vector<double>

      p_eps_x_field_;  // field eps_x(x), n*n row-major (faces normal to x; if has_eps_xy_field_)

  std::vector<double>

      p_eps_y_field_;  // field eps_y(x), n*n row-major (faces normal to y; if has_eps_xy_field_)

  bool has_kappa_field_ = false;  // REACTION term kappa(x) provided: div(eps grad phi) - kappa phi

  std::vector<double> p_kappa_field_;  // field kappa(x), n*n row-major (if has_kappa_field_)

  // GAUSS POLICY (restart-free Gauss-evolution option). Default "restart":

  // solve_fields re-solves -Delta phi = f (Gauss) on EVERY call -- historical behavior,

  // BIT-IDENTICAL. "evolve": after the FIRST solve (phi^0), solve_fields NO LONGER re-solves the

  // Poisson; it only DERIVES the aux (phi, grad phi) from the CURRENT phi -- the one that the condensed

  // source stage (Schur) evolves IN-PLACE in ell_phi() (cf. run_source_stage). Thus it gives a

  // restart-free evolution of -Delta phi (the Gauss constraint is only imposed at t=0).

  // INERT without a Schur stage (phi would stay frozen after t=0). The lock gauss_solved_once_ guarantees that

  // the first solve (the init of phi^0) always solves, whatever the policy.

  bool gauss_evolve_ = false;

  bool gauss_solved_once_ = false;

  std::optional<std::variant<GeometricMG, PoissonFFTSolver, RemappedFFTSolver>> ell_;

  // Direct POLAR Poisson solver (FFT-in-theta + tridiag-in-r), built lazily when

  // polar_ (cf. ensure_elliptic_polar). SEPARATE from ell_ (geom() returns a PolarGeometry, not a

  // Geometry): the cartesian path is never touched. INERT (nullopt) in cartesian.

  std::optional<PolarPoissonSolver> pell_;

  // phi buffer of the POLAR condensed SOURCE STAGE (Path A step 2c). The direct PolarPoissonSolver

  // (pell_->phi()) is WITHOUT ghost (valid box = allocation); but the PolarCondensedSchurSourceStepper

  // needs a phi WITH 1 ghost (fill_ghosts + apply_polar_tensor + centered grad + write of

  // phi^{n+1}). So we pass it this dedicated buffer (1 ghost), fed with phi^n (= aux[0] after

  // solve_fields_polar) before the source stage, which carries phi^{n+1} on output (warm start of the next

  // step). In CARTESIAN this buffer is INERT (nullopt): ell_phi() routes to ell_->phi() as

  // before, BIT-IDENTICAL.

  std::optional<MultiFab> phi_src_polar_;

  std::vector<Real> bz_field_;  // field B_z(x) n*n row-major (empty if not provided)

  int te_src_ = -1;             // index of the fluid block source of T_e (-1 = none)

  // NAMED aux fields (ADC-70 phase 1) provided by the user via System::set_aux_field: key =

  // canonical component (>= kAuxNamedBase = 5), value = field n*n (cartesian) / nr*ntheta (polar)

  // row-major. PERSISTENT like bz_field_: solve_fields touches ONLY components 0..2 (phi,

  // grad) and 4 (T_e via apply_te), so components >= 5 survive from one step to the next; but a

  // REALLOCATION of the aux channel (ensure_aux_width) starts again from a zeroed MultiFab -> we re-apply

  // them then (apply_named_aux), exactly like apply_bz / apply_te.

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

  // Per-field aux HALO policy (ADC-369): key = canonical component (>= kAuxNamedBase), value = the

  // uniform boundary policy declared via pops.AuxHalo. Applied by apply_named_aux_bc() AFTER the shared

  // aux ghost fill, overriding only that component's PHYSICAL-face ghosts (periodic faces -- Cartesian

  // periodic, polar theta -- keep their wrap). Empty -> shared aux BC for every field, bit-identical.

  std::map<int, AuxHaloPolicy> named_aux_bc_;


  // NAMED multi-elliptic fields (ADC-428): a SECOND elliptic solve (beyond the default Poisson) for a

  // user-named field m.elliptic_field("phi2", rhs=..., aux=[...]). Each named field owns:

  //   - phi_comp / gx_comp / gy_comp: the aux channel components (>= kAuxNamedBase) its solution and

  //     centered gradient are written to (the model declares them as aux_field slots; a source then

  //     reads them like any other named aux). gx_comp / gy_comp < 0 => the field declares fewer than 3

  //     aux slots and the gradient is not derived (only phi is written).

  //   - a DEDICATED elliptic solver instance (cartesian GeometricMG / FFT), built lazily the SAME way

  //     as the default ell_ (ensure_named_elliptic mirrors ensure_elliptic's poisson operator), so the

  //     native solver is REUSED, not reimplemented. The default Poisson path (ell_) is untouched.

  // The RHS = sum over blocks of s.named_poisson_rhs[name] (the per-field brick, ADC-428), assembled

  // exactly like assemble_poisson_rhs but reading the named closures.


  struct NamedField {

    int phi_comp = -1;

    int gx_comp = -1;

    int gy_comp = -1;

    std::optional<std::variant<GeometricMG, PoissonFFTSolver, RemappedFFTSolver>> ell;

  };


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


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

    if (named_aux_bc_.empty())

      return;  // hot-path fast exit (parity with the AMR counterparts)

    for (const auto& kv : named_aux_bc_) {

      if (kv.first >= owner_->aux_ncomp_)

        continue;

      fill_physical_bc(owner_->aux, owner_->dom, aux_halo_override(owner_->bc_, kv.second),

                       kv.first);

    }

  }


  void apply_bz() {

    if (bz_field_.empty() || owner_->aux_ncomp_ <= kAuxBaseComps)

      return;

    // ROW WIDTH (fast axis i) of the row-major array bz_field_: n in cartesian (square n x n,

    // BIT-IDENTICAL), nr in POLAR (ring nr x ntheta, i = r of size nr, cf. set_magnetic_field).

    // The index stays flat[j * row + i]: in cartesian row == n (unchanged); in polar row == nr.

    const int row = owner_->polar_ ? owner_->aux.box(0).nx() : owner_->cfg.n;

    // LOCAL population on the owner rank (cf. solve_fields): iteration over the local fabs of the

    // aux channel instead of fab(0) hardcoded (no-op on a rank without a local box at np>1, bit-identical to the

    // owner).

    for (int li = 0; li < owner_->aux.local_size(); ++li) {

      Array4 a = owner_->aux.fab(li).array();

      const Box2D v = owner_->aux.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, kAuxBaseComps) = bz_field_[static_cast<std::size_t>(j) * row + i];

    }

  }


  void apply_te() {

    if (te_src_ < 0 || owner_->aux_ncomp_ <= kTeComp)

      return;

    const auto& s = owner_->sp[static_cast<std::size_t>(te_src_)];

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

    // LOCAL population on the owner rank (cf. solve_fields): we iterate over the local fabs of the

    // aux channel instead of fab(0) hardcoded (no-op on a rank without a local box at np>1, bit-identical to the

    // owner). s.U and aux share the same DistributionMapping -> same local indexing.

    for (int li = 0; li < owner_->aux.local_size(); ++li) {

      const ConstArray4 us = s.U.fab(li).const_array();

      Array4 a = owner_->aux.fab(li).array();

      const Box2D v = owner_->aux.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 rho = us(i, j, 0), mx = us(i, j, 1), my = us(i, j, 2), E = us(i, j, 3);

          const Real p = gm1 * (E - Real(0.5) * (mx * mx + my * my) / rho);

          a(i, j, kTeComp) = p / rho;  // T = p / rho

        }

    }

  }


  void apply_named_aux_one(int comp, const std::vector<Real>& field) {

    if (field.empty() || owner_->aux_ncomp_ <= comp)

      return;

    // ROW WIDTH (fast axis i): n in cartesian (square n x n), nr in polar (ring nr x

    // ntheta). Index flat[j * row + i], identical to apply_bz / set_density.

    const int row = owner_->polar_ ? owner_->aux.box(0).nx() : owner_->cfg.n;

    for (int li = 0; li < owner_->aux.local_size(); ++li) {

      Array4 a = owner_->aux.fab(li).array();

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

    }

  }


  void apply_named_aux() {

    for (const auto& kv : named_aux_)

      apply_named_aux_one(kv.first, kv.second);

  }


  // --- elliptic solver (system Poisson) -----------------------------


  BCRec poisson_bc() {

    std::string mode = p_bc;

    if (mode == "auto")

      mode = (p_wall == "circle" || !owner_->cfg.periodic) ? "dirichlet" : "periodic";

    BCRec b;

    if (mode == "periodic")

      return b;

    if (mode == "dirichlet") {

      b.xlo = b.xhi = b.ylo = b.yhi = BCType::Dirichlet;

      return b;

    }

    if (mode == "neumann") {

      b.xlo = b.xhi = b.ylo = b.yhi = BCType::Foextrap;

      return b;

    }

    throw std::runtime_error("System::set_poisson: unknown bc '" + mode + "'");

  }


  std::function<bool(Real, Real)> wall_active() {

    return detail::wall_predicate(p_wall, p_wall_radius, owner_->cfg.L, "System::set_poisson");

  }


  void ensure_elliptic() {

    if (ell_)

      return;

    // The system right-hand side is ALWAYS f = Sum_s elliptic_rhs_s(u_s), assembled by

    // solve_fields from the elliptic brick of EACH block (charge q n, background alpha (n-n0),

    // gravity coupling 4piG (rho-rho0)). The token is thus NOT a computation mode but a LABEL

    // of this composite right-hand side. "composite" names this behavior honestly; "charge_density"

    // stays the historical alias (default, bit-identical) since the usual case is a charge block.

    if (p_rhs != "charge_density" && p_rhs != "composite")

      throw std::runtime_error("System::set_poisson: unknown rhs '" + p_rhs +

                               "' (charge_density|composite; the right-hand side = sum of the "

                               "per-block elliptic bricks)");

    const BCRec pbc = poisson_bc();

    std::function<bool(Real, Real)> active = wall_active();

    if (p_solver == "fft" || p_solver == "fft_spectral") {

      // The FFT path is CONSTANT-coefficient and PERIODIC only: reject walls, variable / anisotropic

      // permittivity and reaction kappa FIRST (each guard fires identically on all ranks -> no

      // deadlock), so only the pure periodic constant-coefficient case reaches the layout switch below.

      if (active)

        throw std::runtime_error("System: solver '" + p_solver +

                                 "' incompatible with a wall -> 'geometric_mg'");

      if (has_eps_field_)

        throw std::runtime_error("System: solver '" + p_solver +

                                 "' has a CONSTANT coefficient, incompatible with a "

                                 "variable eps(x) field -> use solver='geometric_mg'");

      if (has_eps_xy_field_)

        throw std::runtime_error(

            "System: solver '" + p_solver +

            "' has a CONSTANT coefficient, incompatible with an "

            "ANISOTROPIC permittivity eps_x(x), eps_y(x) -> use solver='geometric_mg'");

      if (has_kappa_field_)

        throw std::runtime_error("System: solver '" + p_solver +

                                 "' (pure Poisson) incompatible with a "

                                 "reaction term kappa(x) -> use solver='geometric_mg'");

      // 'fft_spectral': same plumbing, CONTINUOUS symbol -(kx^2+ky^2) (fidelity to the spectral

      // references, e.g. poisson_fft.m of RIEMOM2D); 'fft' keeps the discrete stencil (bit-identical).

      const bool spectral = (p_solver == "fft_spectral");

      if (n_ranks() > 1) {

        // ADC-287: distributed periodic FFT via a box<->slab remap. System distributes ONE box

        // round-robin (DistributionMapping(1, n_ranks())), so a direct PoissonFFTSolver would

        // dereference a nonexistent fab(0) on a rank without a box (SIGSEGV). RemappedFFTSolver presents

        // the SAME single-box layout outward (rhs()/phi() on owner_->ba/dm, aligned with the aux) and

        // hides a scatter/gather around PoissonFFT inside solve() -> the field-solve path is untouched.

        // COLLECTIVE: every rank constructs the same type; its Ny % n_ranks() guard throws on all ranks

        // (no deadlock). NOTE: pending the ADC-273 design vote (structural change to the ell_ variant).

        ell_.emplace(std::in_place_type<RemappedFFTSolver>, owner_->geom, owner_->ba, pbc, active,

                     spectral);

      } else {

        // Single-rank: the proven direct FFT on the System box. PoissonFFTSolver keeps a hard guard in

        // its constructor (rejects n_ranks()>1 / ba.size()!=1).

        ell_.emplace(std::in_place_type<PoissonFFTSolver>, owner_->geom, owner_->ba, pbc, active,

                     spectral);

      }

    } else if (p_solver == "geometric_mg") {

      ell_.emplace(std::in_place_type<GeometricMG>, owner_->geom, owner_->ba, pbc,

                   std::move(active));

      std::get<GeometricMG>(*ell_).set_abs_tol(

          p_abs_tol_);  // absolute floor of the V-cycle (0 = relative only)

      if (has_eps_field_)

        apply_epsilon_field();  // operator div(eps grad phi) with variable eps(x)

      if (has_eps_xy_field_)

        apply_epsilon_anisotropic_field();  // div(diag(eps_x, eps_y) grad phi)

      if (has_kappa_field_)

        apply_reaction_field();  // term - kappa phi (screened Poisson / Helmholtz)

      // Guard: with kappa and a CONSTANT permittivity eps != 1 (without an eps(x) field), the rhs

      // would be scaled by 1/eps (shortcut lap phi = f/eps) -- inconsistent with the term -kappa phi.

      // We then require eps = 1 or an eps(x) field (carried by the operator, without scaling).

      if (has_kappa_field_ && !has_eps_field_ && !has_eps_xy_field_ && p_eps_ != Real(1))

        throw std::runtime_error(

            "System: reaction term kappa(x) + CONSTANT permittivity eps != 1 "

            "unsupported; use eps = 1 or an eps(x) field (set_epsilon_field)");

    } else {

      throw std::runtime_error("System::set_poisson: unknown solver '" + p_solver +

                               "' (geometric_mg|fft|fft_spectral)");

    }

  }


  void apply_epsilon_field() {

    GeometricMG& mg = std::get<GeometricMG>(*ell_);

    MultiFab eps_fine(owner_->ba, owner_->dm, 1, 0);

    const int n = owner_->cfg.n;

    // Filling of the source field LOCAL to the owner rank (iteration over the local fabs, never

    // fab(0) hardcoded): no-op on a rank without a local box at np>1, identical to before on the

    // owner. mg.set_epsilon is then COLLECTIVE (local copy + MPI-safe restriction).

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

      Array4 e = eps_fine.fab(li).array();

      const Box2D v = eps_fine.box(li);

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

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

          e(i, j, 0) = static_cast<Real>(p_eps_field_[static_cast<std::size_t>(j) * n + i]);

    }

    mg.set_epsilon(

        eps_fine);  // copy on the fine level + restriction (average_down) to the coarse ones

  }


  void apply_epsilon_anisotropic_field() {

    GeometricMG& mg = std::get<GeometricMG>(*ell_);

    MultiFab eps_x_fine(owner_->ba, owner_->dm, 1, 0), eps_y_fine(owner_->ba, owner_->dm, 1, 0);

    const int n = owner_->cfg.n;

    // LOCAL filling on the owner rank (cf. apply_epsilon_field): iteration over the local fabs

    // (no-op on an empty rank at np>1). eps_x_fine and eps_y_fine share ba/dm -> same

    // local indexing. mg.set_epsilon_anisotropic is then COLLECTIVE (copy + restriction).

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

      Array4 ex = eps_x_fine.fab(li).array();

      Array4 ey = eps_y_fine.fab(li).array();

      const Box2D v = eps_x_fine.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 std::size_t k = static_cast<std::size_t>(j) * n + i;

          ex(i, j, 0) = static_cast<Real>(p_eps_x_field_[k]);

          ey(i, j, 0) = static_cast<Real>(p_eps_y_field_[k]);

        }

    }

    mg.set_epsilon_anisotropic(eps_x_fine,

                               eps_y_fine);  // faces x <- eps_x, faces y <- eps_y (+ restriction)

  }


  void apply_reaction_field() {

    GeometricMG& mg = std::get<GeometricMG>(*ell_);

    MultiFab kappa_fine(owner_->ba, owner_->dm, 1, 0);

    const int n = owner_->cfg.n;

    // LOCAL filling on the owner rank (cf. apply_epsilon_field): iteration over the local fabs

    // (no-op on an empty rank at np>1). mg.set_reaction is then COLLECTIVE.

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

      Array4 k = kappa_fine.fab(li).array();

      const Box2D v = kappa_fine.box(li);

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

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

          k(i, j, 0) = static_cast<Real>(p_kappa_field_[static_cast<std::size_t>(j) * n + i]);

    }

    mg.set_reaction(kappa_fine);

  }


  MultiFab& ell_rhs() {

    return std::visit([](auto& e) -> MultiFab& { return e.rhs(); }, *ell_);

  }


  MultiFab& ell_phi() {

    if (owner_->polar_) {

      // Allocate lazily (1 ghost) on the System layout, then copy phi^n from aux[0].

      if (!phi_src_polar_)

        phi_src_polar_.emplace(owner_->ba, owner_->dm, 1, 1);

      for (int li = 0; li < phi_src_polar_->local_size(); ++li) {

        const ConstArray4 a = owner_->aux.fab(li).const_array();

        Array4 p = phi_src_polar_->fab(li).array();

        const Box2D v = phi_src_polar_->box(li);

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

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

            p(i, j, 0) = a(i, j, 0);  // aux[0] = phi^n

      }

      return *phi_src_polar_;

    }

    return std::visit([](auto& e) -> MultiFab& { return e.phi(); }, *ell_);

  }


  void ell_solve() {

    pops_sf_mark("ell_solve: before std::visit");

    std::visit(

        [](auto& e) {

          using T = std::decay_t<decltype(e)>;

          if (pops_trace_sf())

            pops_sf_mark(std::is_same_v<T, GeometricMG> ? "ell_solve: GeometricMG::solve() start"

                                                       : "ell_solve: FFT solver::solve() start");

          e.solve();

          pops_sf_mark("ell_solve: solve() return");

        },

        *ell_);

    pops_sf_mark("ell_solve: after std::visit");

  }


  // --- ELLIPTIC-SOLVER PROFILING COUNTERS (Spec 5 sec.13.11.1, ADC-479 criteria 42/43) ----------

  // Read-back accessors for the per-solve stats of the ACTIVE cartesian elliptic solver, queried at

  // the System solve_fields seam AFTER ell_solve()+device_fence() to populate the native profiler

  // counters (mg_cycles / krylov_iters / mg_levels / elliptic_bottom). LOW-INVASIVE by design: the

  // deep numerics never sees the profiler -- GeometricMG caches the stats per solve (chrono only) and

  // we visit the variant here. Solvers that lack a notion (FFT: direct, no cycles/levels/iters)

  // return 0 HONESTLY -- not every solver has these. nullopt ell_ (never solved) -> 0.

  //

  // krylov_iters is 0 on the cartesian field-solve path: the default Poisson uses GeometricMG (or a

  // direct FFT), never a Krylov solver (the Krylov path lives in the condensed Schur SOURCE stage,

  // which is not the ell_ elliptic solve). The counter is emitted for completeness / future Krylov

  // elliptic backends; it stays an honest 0 here.


  int last_mg_cycles() const {

    if (!ell_)

      return 0;

    return std::visit(

        [](const auto& e) -> int {

          using T = std::decay_t<decltype(e)>;

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

            return e.last_cycles();

          else

            return 0;  // direct FFT solver: no V-cycles

        },

        *ell_);

  }


  int last_krylov_iters() const {

    return 0;  // cartesian ell_ = GeometricMG / FFT, never a Krylov elliptic solver (see note above)

  }


  int last_num_levels() const {

    if (!ell_)

      return 0;

    return std::visit(

        [](const auto& e) -> int {

          using T = std::decay_t<decltype(e)>;

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

            return e.num_levels();

          else

            return 0;  // direct FFT solver: no multigrid hierarchy

        },

        *ell_);

  }


  double last_bottom_seconds() const {

    if (!ell_)

      return 0.0;

    return std::visit(

        [](const auto& e) -> double {

          using T = std::decay_t<decltype(e)>;

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

            return e.last_bottom_seconds();

          else

            return 0.0;  // direct FFT solver: no coarsest-grid bottom solve

        },

        *ell_);

  }


  // --- direct POLAR Poisson (PolarPoissonSolver) -------------------------


  void ensure_elliptic_polar() {

    if (pell_)

      return;

    if (p_rhs != "charge_density" && p_rhs != "composite")

      throw std::runtime_error("System::set_poisson (polar): unknown rhs '" + p_rhs +

                               "' (charge_density|composite)");

    if (p_solver != "geometric_mg" && p_solver != "polar")

      throw std::runtime_error(

          "System::set_poisson (polar): solver '" + p_solver +

          "' unsupported on a ring; the polar Poisson is direct (FFT-in-theta + tridiag-in-r). "

          "Leave the default ('geometric_mg') or request 'polar'.");

    if (has_eps_field_ || has_eps_xy_field_ || has_kappa_field_)

      throw std::runtime_error(

          "System::set_poisson (polar): variable / anisotropic permittivity / reaction unsupported "

          "by the direct polar Poisson (Phase 2b; operator (1/r) d_r(r d_r) + (1/r^2) d_theta^2)");

    // MULTI-BOX GUARD (ADC-67): the DIRECT polar Poisson (FFT-in-theta + tridiag-in-r) requires

    // complete theta ROWS and radial COLUMNS on ONE box (PolarPoissonSolver already rejects

    // ba.size()!=1). We raise a clear UPSTREAM message HERE -- before the construction of the solver, so from the

    // 1st solve_fields / potential / set_potential of a System with theta_boxes > 1 -- rather than letting

    // the low-level rejection of PolarPoissonSolver bubble up. The theta SPLITTING (theta_boxes > 1) only serves

    // the TRANSPORT; for a multi-box electrostatic field, go through the tensor Schur stage.

    if (owner_->ba.size() != 1)

      throw std::runtime_error(

          "System: DIRECT polar Poisson incompatible with the theta splitting (theta_boxes=" +

          std::to_string(owner_->ba.size()) +

          " boxes); it requires a single-box grid (theta_boxes=1). For a multi-box theta "

          "splitting, "

          "use the polar tensor Schur stage (pops.Split + pops.CondensedSchur), multi-box, or "

          "go back to theta_boxes=1.");

    // Radial BC: Dirichlet/Neumann from poisson_bc() (xlo/xhi). theta always periodic.

    const BCRec pbc = poisson_bc();

    pell_.emplace(owner_->pgeom_, owner_->ba, pbc);

  }


  void assemble_poisson_rhs(MultiFab& rhs, int target_block = -1,

                            const MultiFab* U_stage = nullptr) {

    rhs.set_val(Real(0));

    for (std::size_t b = 0; b < owner_->sp.size(); ++b) {

      auto& s = owner_->sp[b];

      const bool override_here = (U_stage != nullptr && static_cast<int>(b) == target_block);

      s.add_poisson_rhs(override_here ? *U_stage : s.U, rhs);

    }

  }


  void assemble_poisson_rhs_from_blocks(MultiFab& rhs,

                                        const std::vector<const MultiFab*>& U_stages) {

    if (U_stages.size() != owner_->sp.size())

      throw std::invalid_argument(

          "assemble_poisson_rhs_from_blocks: U_stages size " + std::to_string(U_stages.size()) +

          " != block count " + std::to_string(owner_->sp.size()) +

          " (index U_stages by block index; nullptr = use the block's live state)");

    rhs.set_val(Real(0));

    for (std::size_t b = 0; b < owner_->sp.size(); ++b) {

      auto& s = owner_->sp[b];

      s.add_poisson_rhs(U_stages[b] != nullptr ? *U_stages[b] : s.U, rhs);

    }

  }


  void solve_fields_polar(int target_block = -1, const MultiFab* U_stage = nullptr) {

    ensure_elliptic_polar();

    MultiFab& rhs = pell_->rhs();

    assemble_poisson_rhs(rhs, target_block, U_stage);

    // CONSTANT permittivity eps != 1: lap phi = f/eps (1/eps scaling of the rhs), like the

    // cartesian. (variable/aniso eps(x) is refused by ensure_elliptic_polar.)

    if (p_eps_ != Real(1)) {

      const Real inv = Real(1) / p_eps_;

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

        Array4 r = rhs.fab(li).array();

        const Box2D v = rhs.box(li);

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

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

            r(i, j, 0) *= inv;

      }

    }

    pell_->solve();

    device_fence();

    // Derivation (phi, grad_r, grad_theta) in the local basis (e_r, e_theta) via the SAME helper as the C++

    // test (derive_aux_polar of block_builder_polar.hpp). phi is WITHOUT ghost (direct single-box solver):

    // the helper thus never reads an out-of-domain index (radial OFFSET at the walls, theta WRAPPED in

    // periodic) -- that was the bug: the centered difference read phi(lo-1)/phi(hi+1)/phi(.,jlo-1) out of

    // allocation -> spurious gradient -> divergent velocity -> nan.

    derive_aux_polar(pell_->phi(), owner_->aux, owner_->pgeom_);

    apply_te();  // inert in polar ExB (no fluid block source of T_e), kept by symmetry

    // Aux ghosts: theta PERIODIC (joint 0/2pi), r PHYSICAL (extrapolation at the boundary). fill_ghosts

    // already routes through bc_ (xlo/xhi Foextrap, ylo/yhi Periodic) -> correct periodic azimuthal halo.

    fill_ghosts(owner_->aux, owner_->dom, owner_->bc_);

    apply_named_aux_bc();  // ADC-369: per-field halo override on the RADIAL faces (theta stays periodic)

  }


  void solve_fields(int target_block = -1, const MultiFab* U_stage = nullptr) {

    if (owner_->polar_)

      // ring: polar Poisson + aux in local basis (e_r, e_theta)

      return solve_fields_polar(target_block, U_stage);

    pops_sf_mark("solve_fields: start");

    ensure_elliptic();

    pops_sf_mark("solve_fields: after ensure_elliptic");

    // GAUSS POLICY: "restart" (default, gauss_evolve_==false) re-solves -Delta phi = f on

    // EVERY call (bit-identical to history). "evolve": after the FIRST solve (phi^0), we SKIP

    // the RHS assembly + the elliptic solve -> ell_phi() keeps the current phi (the one that the

    // Schur source stage evolves in-place), giving a restart-free -Delta phi evolution (Gauss imposed only at t=0).

    // The derivation of the aux (phi, grad phi) below ALWAYS runs, on the current phi.

    if (!(gauss_evolve_ && gauss_solved_once_)) {

      MultiFab& rhs = ell_rhs();

      // f = Sum_s elliptic_rhs_s(U_s). By default the CURRENT state of each block; with a

      // target_block / U_stage override the target block reads U_stage (per-stage field solve, ADC-409).

      assemble_poisson_rhs(rhs, target_block, U_stage);

      pops_sf_mark("solve_fields: after add_poisson_rhs");

      // CONSTANT permittivity: div(eps grad phi) = f <=> lap phi = f/eps, so we scale the rhs by

      // 1/eps. With a VARIABLE or ANISOTROPIC eps(x) field we DO NOT do it: the GeometricMG

      // operator carries eps directly (apply_epsilon_field / apply_epsilon_anisotropic_field), the

      // rhs stays f as is.

      if (!has_eps_field_ && !has_eps_xy_field_ && p_eps_ != Real(1)) {

        const Real inv = Real(1) / p_eps_;

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

          Array4 r = rhs.fab(li).array();

          const Box2D v = rhs.box(li);

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

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

              r(i, j, 0) *= inv;

        }

      }

      pops_sf_mark("solve_fields: before ell_solve");

      ell_solve();

      gauss_solved_once_ = true;

      pops_sf_mark("solve_fields: after ell_solve, before device_fence");

    }

    device_fence();

    pops_sf_mark("solve_fields: after device_fence (aux derivation)");

    const Real dx = owner_->geom.dx(), dy = owner_->geom.dy();

    // Per-cell derivation (phi, grad phi) -> aux channel: LOCAL to the owner rank. System

    // distributes ONE box (round-robin DistributionMapping(1, n_ranks())), so at np>1 a single rank

    // owns it; the others have local_size()==0 and HAVE NO fab to derive. We iterate over the LOCAL

    // fabs (never fab(0) hardcoded): no-op on an empty rank, identical to before on the owner

    // (loop executed once, bit-identical to np=1). ell_phi() and aux share the same

    // DistributionMapping -> same local indexing.

    MultiFab& phi_mf = ell_phi();

    for (int li = 0; li < owner_->aux.local_size(); ++li) {

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

      Array4 a = owner_->aux.fab(li).array();

      const Box2D v = owner_->aux.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, 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);

        }

    }

    pops_sf_mark("solve_fields: after aux derivation (phi, grad phi)");

    apply_te();  // T_e = p/rho of the fluid block source, recomputed on each solve (B_z, comp 3, preserved)

    pops_sf_mark("solve_fields: after apply_te");

    if (owner_->periodic_)

      fill_boundary(owner_->aux, owner_->dom, owner_->per_);

    else

      fill_ghosts(owner_->aux, owner_->dom,

                  owner_->bc_);  // extrapolation at the boundary (wall / free outflow)

    apply_named_aux_bc();  // ADC-369: per-field halo override (after the shared fill; no-op if none)

    pops_sf_mark("solve_fields: end (fill ghosts aux)");

  }


  void solve_fields_from_state(int block_idx, const MultiFab& U_stage) {

    if (block_idx < 0 || block_idx >= static_cast<int>(owner_->sp.size()))

      throw std::out_of_range("solve_fields_from_state: block index " + std::to_string(block_idx) +

                              " out of range (" + std::to_string(owner_->sp.size()) + " blocks)");

    solve_fields(block_idx, &U_stage);

  }


  void solve_fields_polar_from_blocks(const std::vector<const MultiFab*>& U_stages) {

    ensure_elliptic_polar();

    MultiFab& rhs = pell_->rhs();

    assemble_poisson_rhs_from_blocks(rhs, U_stages);

    if (p_eps_ != Real(1)) {

      const Real inv = Real(1) / p_eps_;

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

        Array4 r = rhs.fab(li).array();

        const Box2D v = rhs.box(li);

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

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

            r(i, j, 0) *= inv;

      }

    }

    pell_->solve();

    device_fence();

    derive_aux_polar(pell_->phi(), owner_->aux, owner_->pgeom_);

    apply_te();

    fill_ghosts(owner_->aux, owner_->dom, owner_->bc_);

    apply_named_aux_bc();

  }


  void solve_fields_from_blocks(const std::vector<const MultiFab*>& U_stages) {

    if (owner_->polar_)

      return solve_fields_polar_from_blocks(U_stages);

    ensure_elliptic();

    // Coupled multi-block solve always re-solves the Poisson (it is requested explicitly with the

    // stage states); it is NOT subject to the gauss_evolve_ skip (that policy is for the single

    // default head solve). The aux derivation below runs on the freshly solved phi.

    {

      MultiFab& rhs = ell_rhs();

      assemble_poisson_rhs_from_blocks(rhs, U_stages);

      if (!has_eps_field_ && !has_eps_xy_field_ && p_eps_ != Real(1)) {

        const Real inv = Real(1) / p_eps_;

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

          Array4 r = rhs.fab(li).array();

          const Box2D v = rhs.box(li);

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

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

              r(i, j, 0) *= inv;

        }

      }

      ell_solve();

      gauss_solved_once_ = true;

    }

    device_fence();

    const Real dx = owner_->geom.dx(), dy = owner_->geom.dy();

    MultiFab& phi_mf = ell_phi();

    for (int li = 0; li < owner_->aux.local_size(); ++li) {

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

      Array4 a = owner_->aux.fab(li).array();

      const Box2D v = owner_->aux.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, 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);

        }

    }

    apply_te();

    if (owner_->periodic_)

      fill_boundary(owner_->aux, owner_->dom, owner_->per_);

    else

      fill_ghosts(owner_->aux, owner_->dom, owner_->bc_);

    apply_named_aux_bc();

  }


  // --- NAMED multi-elliptic field (ADC-428) ----------------------------------


  void ensure_named_elliptic(NamedField& nf) {

    if (nf.ell)

      return;

    const BCRec pbc = poisson_bc();

    std::function<bool(Real, Real)> active = wall_active();

    if (p_solver == "fft" || p_solver == "fft_spectral") {

      if (active)

        throw std::runtime_error("System: named elliptic field solver '" + p_solver +

                                 "' incompatible with a wall -> 'geometric_mg'");

      const bool spectral = (p_solver == "fft_spectral");

      if (n_ranks() > 1)

        nf.ell.emplace(std::in_place_type<RemappedFFTSolver>, owner_->geom, owner_->ba, pbc, active,

                       spectral);

      else

        nf.ell.emplace(std::in_place_type<PoissonFFTSolver>, owner_->geom, owner_->ba, pbc, active,

                       spectral);

    } else if (p_solver == "geometric_mg") {

      nf.ell.emplace(std::in_place_type<GeometricMG>, owner_->geom, owner_->ba, pbc,

                     std::move(active));

      std::get<GeometricMG>(*nf.ell).set_abs_tol(p_abs_tol_);

    } else {

      throw std::runtime_error("System: named elliptic field solver '" + p_solver +

                               "' unsupported (geometric_mg|fft|fft_spectral)");

    }

  }


  void assemble_named_poisson_rhs(const std::string& field, MultiFab& rhs, int target_block,

                                  const MultiFab* U_stage) {

    rhs.set_val(Real(0));

    bool any = false;

    for (std::size_t b = 0; b < owner_->sp.size(); ++b) {

      auto& s = owner_->sp[b];

      auto it = s.named_poisson_rhs.find(field);

      if (it == s.named_poisson_rhs.end() || !it->second)

        continue;

      const bool override_here = (U_stage != nullptr && static_cast<int>(b) == target_block);

      it->second(override_here ? *U_stage : s.U, rhs);

      any = true;

    }

    if (!any)

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

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

                               "model)");

  }


  void solve_named_field_from_state(const std::string& field, int block_idx,

                                    const MultiFab& U_stage) {

    if (block_idx < 0 || block_idx >= static_cast<int>(owner_->sp.size()))

      throw std::out_of_range("solve_fields_from_state (named): block index " +

                              std::to_string(block_idx) + " out of range (" +

                              std::to_string(owner_->sp.size()) + " blocks)");

    auto it = named_fields_.find(field);

    if (it == named_fields_.end())

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

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

    if (owner_->polar_)

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

                               "' on a polar (ring) grid is not supported yet (cartesian only)");

    NamedField& nf = it->second;

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

      throw std::runtime_error(

          "System: named elliptic field '" + field +

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

          "its aux fields before solving)");

    ensure_named_elliptic(nf);

    MultiFab& rhs = std::visit([](auto& e) -> MultiFab& { return e.rhs(); }, *nf.ell);

    assemble_named_poisson_rhs(field, rhs, block_idx, &U_stage);

    // CONSTANT permittivity eps != 1: lap phi = f/eps (1/eps scaling), like the default Poisson. The

    // variable / anisotropic eps(x) and reaction kappa of the DEFAULT Poisson are NOT applied to a

    // named field (plain Laplacian operator; see ensure_named_elliptic).

    if (p_eps_ != Real(1)) {

      const Real inv = Real(1) / p_eps_;

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

        Array4 r = rhs.fab(li).array();

        const Box2D v = rhs.box(li);

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

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

            r(i, j, 0) *= inv;

      }

    }

    std::visit([](auto& e) { e.solve(); }, *nf.ell);

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

    MultiFab& phi_mf = std::visit([](auto& e) -> MultiFab& { return e.phi(); }, *nf.ell);

    const Real dx = owner_->geom.dx(), dy = owner_->geom.dy();

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

    const bool grad =

        (cgx >= 0 && cgx < owner_->aux_ncomp_ && cgy >= 0 && cgy < owner_->aux_ncomp_);

    for (int li = 0; li < owner_->aux.local_size(); ++li) {

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

      Array4 a = owner_->aux.fab(li).array();

      const Box2D v = owner_->aux.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 field's aux components: the shared aux fill (same routing as solve_fields)

    // then the per-field halo override (ADC-369). This re-fills ALL components -- the shared phi/grad

    // were last written by the default solve_fields, so their valid cells are unchanged and only the

    // halos are recomputed (idempotent for those components).

    if (owner_->periodic_)

      fill_boundary(owner_->aux, owner_->dom, owner_->per_);

    else

      fill_ghosts(owner_->aux, owner_->dom, owner_->bc_);

    apply_named_aux_bc();

  }


 private:

  Impl* owner_;

};


}  // namespace field_solver


}  // namespace pops

block_builder_polar.hpp
POLAR counterpart of block_builder.hpp: builds a block's closures (time advance + residual + Poisson ...

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

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

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

pops::GeometricMG
Definition geometric_mg.hpp:79

pops::GeometricMG::set_epsilon
void set_epsilon(std::function< Real(Real, Real)> eps_fn)
Definition geometric_mg.hpp:249

pops::GeometricMG::set_reaction
void set_reaction(std::function< Real(Real, Real)> kappa_fn)
Definition geometric_mg.hpp:304

pops::GeometricMG::set_epsilon_anisotropic
void set_epsilon_anisotropic(std::function< Real(Real, Real)> eps_x_fn, std::function< Real(Real, Real)> eps_y_fn)
Definition geometric_mg.hpp:274

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

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

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

pops::MultiFab::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::field_solver::SystemFieldSolver
SystemFieldSolver<Impl>: see contract above.
Definition system_field_solver.hpp:83

pops::field_solver::SystemFieldSolver::solve_fields
void solve_fields(int target_block=-1, const MultiFab *U_stage=nullptr)
Solves the system Poisson then DERIVES the aux = (phi, grad phi[, B_z, T_e]).
Definition system_field_solver.hpp:689

pops::field_solver::SystemFieldSolver::te_src_
int te_src_
Definition system_field_solver.hpp:154

pops::field_solver::SystemFieldSolver::has_kappa_field_
bool has_kappa_field_
Definition system_field_solver.hpp:128

pops::field_solver::SystemFieldSolver::poisson_bc
BCRec poisson_bc()
Resolves the BC mode into a BCRec: "auto" -> dirichlet if wall/non-periodic, otherwise periodic; "per...
Definition system_field_solver.hpp:292

pops::field_solver::SystemFieldSolver::has_eps_xy_field_
bool has_eps_xy_field_
Definition system_field_solver.hpp:122

pops::field_solver::SystemFieldSolver::p_eps_x_field_
std::vector< double > p_eps_x_field_
Definition system_field_solver.hpp:125

pops::field_solver::SystemFieldSolver::solve_fields_from_blocks
void solve_fields_from_blocks(const std::vector< const MultiFab * > &U_stages)
Coupled multi-block field solve (Spec 3 criterion 24, ADC-457): SAME elliptic solve + aux derivation ...
Definition system_field_solver.hpp:813

pops::field_solver::SystemFieldSolver::ell_
std::optional< std::variant< GeometricMG, PoissonFFTSolver, RemappedFFTSolver > > ell_
Definition system_field_solver.hpp:140

pops::field_solver::SystemFieldSolver::phi_src_polar_
std::optional< MultiFab > phi_src_polar_
Definition system_field_solver.hpp:152

pops::field_solver::SystemFieldSolver::has_eps_field_
bool has_eps_field_
Definition system_field_solver.hpp:120

pops::field_solver::SystemFieldSolver::p_wall
std::string p_wall
Definition system_field_solver.hpp:112

pops::field_solver::SystemFieldSolver::solve_fields_from_state
void solve_fields_from_state(int block_idx, const MultiFab &U_stage)
Per-stage field solve (ADC-409): SAME solve + derive-aux as solve_fields(), but the target block bloc...
Definition system_field_solver.hpp:767

pops::field_solver::SystemFieldSolver::p_wall_radius
double p_wall_radius
Definition system_field_solver.hpp:113

pops::field_solver::SystemFieldSolver::apply_te
void apply_te()
Populates the T_e component (electron temperature) = p/rho of the fluid block source te_src_.
Definition system_field_solver.hpp:240

pops::field_solver::SystemFieldSolver::apply_bz
void apply_bz()
Populates the B_z component (index kAuxBaseComps) of the shared channel from bz_field_,...
Definition system_field_solver.hpp:217

pops::field_solver::SystemFieldSolver::ell_phi
MultiFab & ell_phi()
Potential phi read (and rewritten) by the condensed source stage.
Definition system_field_solver.hpp:468

pops::field_solver::SystemFieldSolver::p_eps_y_field_
std::vector< double > p_eps_y_field_
Definition system_field_solver.hpp:127

pops::field_solver::SystemFieldSolver::p_bc
std::string p_bc
Definition system_field_solver.hpp:111

pops::field_solver::SystemFieldSolver::ensure_elliptic_polar
void ensure_elliptic_polar()
Builds the direct POLAR Poisson (PolarPoissonSolver, single-rank, single box covering the ring) LAZIL...
Definition system_field_solver.hpp:564

pops::field_solver::SystemFieldSolver::SystemFieldSolver
SystemFieldSolver(Impl *owner)
Definition system_field_solver.hpp:86

pops::field_solver::SystemFieldSolver::ell_solve
void ell_solve()
Solves the active cartesian Poisson (GeometricMG V-cycle or direct FFT).
Definition system_field_solver.hpp:487

pops::field_solver::SystemFieldSolver::register_named_field
void register_named_field(const std::string &field, int phi_comp, int gx_comp, int gy_comp)
Register a named elliptic field (ADC-428): records the aux output components (where the field's solve...
Definition system_field_solver.hpp:192

pops::field_solver::SystemFieldSolver::kTeComp
static constexpr int kTeComp
Canonical component of T_e (after phi/grad/B_z); cf. pops::Aux and AUX_CANONICAL on the DSL side.
Definition system_field_solver.hpp:89

pops::field_solver::SystemFieldSolver::apply_named_aux_bc
void apply_named_aux_bc()
Re-applies the per-field aux HALO policies (ADC-369) onto the shared channel, AFTER the shared fill_g...
Definition system_field_solver.hpp:203

pops::field_solver::SystemFieldSolver::gauss_evolve_
bool gauss_evolve_
Definition system_field_solver.hpp:138

pops::field_solver::SystemFieldSolver::p_abs_tol_
Real p_abs_tol_
Definition system_field_solver.hpp:119

pops::field_solver::SystemFieldSolver::pell_
std::optional< PolarPoissonSolver > pell_
Definition system_field_solver.hpp:144

pops::field_solver::SystemFieldSolver::apply_named_aux_one
void apply_named_aux_one(int comp, const std::vector< Real > &field)
Populates ONE NAMED aux component (canonical index comp >= kAuxNamedBase) of the shared channel from ...
Definition system_field_solver.hpp:267

pops::field_solver::SystemFieldSolver::bz_field_
std::vector< Real > bz_field_
Definition system_field_solver.hpp:153

pops::field_solver::SystemFieldSolver::solve_fields_polar_from_blocks
void solve_fields_polar_from_blocks(const std::vector< const MultiFab * > &U_stages)
POLAR coupled multi-block solve (Spec 3 criterion 24, ADC-457): same solve + aux derivation as solve_...
Definition system_field_solver.hpp:780

pops::field_solver::SystemFieldSolver::last_mg_cycles
int last_mg_cycles() const
Definition system_field_solver.hpp:514

pops::field_solver::SystemFieldSolver::named_fields_
std::map< std::string, NamedField > named_fields_
Definition system_field_solver.hpp:185

pops::field_solver::SystemFieldSolver::p_rhs
std::string p_rhs
Definition system_field_solver.hpp:109

pops::field_solver::SystemFieldSolver::p_eps_
Real p_eps_
Definition system_field_solver.hpp:114

pops::field_solver::SystemFieldSolver::last_num_levels
int last_num_levels() const
Definition system_field_solver.hpp:530

pops::field_solver::SystemFieldSolver::apply_reaction_field
void apply_reaction_field()
Installs the reaction term kappa(x) (n*n row-major) on the GeometricMG: the operator becomes div(eps ...
Definition system_field_solver.hpp:443

pops::field_solver::SystemFieldSolver::assemble_named_poisson_rhs
void assemble_named_poisson_rhs(const std::string &field, MultiFab &rhs, int target_block, const MultiFab *U_stage)
Assembles the RIGHT-HAND SIDE of named field field into rhs: f = Sum_s named_poisson_rhs_s[field](U_s...
Definition system_field_solver.hpp:896

pops::field_solver::SystemFieldSolver::last_bottom_seconds
double last_bottom_seconds() const
Definition system_field_solver.hpp:543

pops::field_solver::SystemFieldSolver::assemble_poisson_rhs_from_blocks
void assemble_poisson_rhs_from_blocks(MultiFab &rhs, const std::vector< const MultiFab * > &U_stages)
Assembles the system Poisson RIGHT-HAND SIDE into rhs for a SIMULTANEOUS coupled multi-block solve (S...
Definition system_field_solver.hpp:627

pops::field_solver::SystemFieldSolver::ensure_named_elliptic
void ensure_named_elliptic(NamedField &nf)
Builds the DEDICATED cartesian elliptic solver of named field nf, the SAME poisson operator as the de...
Definition system_field_solver.hpp:865

pops::field_solver::SystemFieldSolver::provides_aux
bool provides_aux(const std::string &name) const
True if the named aux field is provided at install time, for Spec-2 criterion 24 install-time require...
Definition system_field_solver.hpp:97

pops::field_solver::SystemFieldSolver::p_kappa_field_
std::vector< double > p_kappa_field_
Definition system_field_solver.hpp:129

pops::field_solver::SystemFieldSolver::p_eps_field_
std::vector< double > p_eps_field_
Definition system_field_solver.hpp:121

pops::field_solver::SystemFieldSolver::named_aux_bc_
std::map< int, AuxHaloPolicy > named_aux_bc_
Definition system_field_solver.hpp:166

pops::field_solver::SystemFieldSolver::apply_epsilon_field
void apply_epsilon_field()
Installs the eps(x) field (n*n row-major) on the GeometricMG: the operator becomes div(eps grad phi) ...
Definition system_field_solver.hpp:398

pops::field_solver::SystemFieldSolver::apply_named_aux
void apply_named_aux()
Re-applies ALL the stored named aux fields (cf.
Definition system_field_solver.hpp:284

pops::field_solver::SystemFieldSolver::gauss_solved_once_
bool gauss_solved_once_
Definition system_field_solver.hpp:139

pops::field_solver::SystemFieldSolver::named_aux_
std::map< int, std::vector< Real > > named_aux_
Definition system_field_solver.hpp:161

pops::field_solver::SystemFieldSolver::apply_epsilon_anisotropic_field
void apply_epsilon_anisotropic_field()
Installs the eps_x(x), eps_y(x) fields (n*n row-major each) on the GeometricMG: the operator becomes ...
Definition system_field_solver.hpp:419

pops::field_solver::SystemFieldSolver::assemble_poisson_rhs
void assemble_poisson_rhs(MultiFab &rhs, int target_block=-1, const MultiFab *U_stage=nullptr)
Assembles the system Poisson RIGHT-HAND SIDE into rhs: f = Sum_s elliptic_rhs_s(U_s),...
Definition system_field_solver.hpp:606

pops::field_solver::SystemFieldSolver::wall_active
std::function< bool(Real, Real)> wall_active()
"Conductor interior" predicate from p_wall / p_wall_radius / cfg.L (cf.
Definition system_field_solver.hpp:311

pops::field_solver::SystemFieldSolver::ensure_elliptic
void ensure_elliptic()
Builds the cartesian elliptic solver (ell_) LAZILY according to p_solver: GeometricMG (carries eps(x)...
Definition system_field_solver.hpp:319

pops::field_solver::SystemFieldSolver::last_krylov_iters
int last_krylov_iters() const
Definition system_field_solver.hpp:527

pops::field_solver::SystemFieldSolver::p_solver
std::string p_solver
Definition system_field_solver.hpp:110

pops::field_solver::SystemFieldSolver::solve_named_field_from_state
void solve_named_field_from_state(const std::string &field, int block_idx, const MultiFab &U_stage)
Solves named field field's SECOND elliptic problem from block block_idx's stage state U_stage and wri...
Definition system_field_solver.hpp:922

pops::field_solver::SystemFieldSolver::solve_fields_polar
void solve_fields_polar(int target_block=-1, const MultiFab *U_stage=nullptr)
POLAR solve_fields: assembles f = Sum_s elliptic_rhs_s(U_s) (host loop per cell), solves the polar Po...
Definition system_field_solver.hpp:647

pops::field_solver::SystemFieldSolver::ell_rhs
MultiFab & ell_rhs()
Right-hand side f of the active cartesian elliptic solver (GeometricMG or FFT), via std::visit.
Definition system_field_solver.hpp:459

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

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

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

pops::detail::wall_predicate
std::function< bool(Real, Real)> wall_predicate(const std::string &wall, double wall_radius, double L, const std::string &err_context)
Builds the "inside the conductor" predicate (embedded wall for the Poisson solver) from the wall mode...
Definition wall_predicate.hpp:30

pops::field_solver::pops_trace_sf
bool pops_trace_sf()
True if the DIAGNOSTIC trace of the solve_fields path is active (environment variable POPS_TRACE_SOLV...
Definition system_field_solver.hpp:66

pops::field_solver::pops_sf_mark
void pops_sf_mark(const char *w)
Writes the marker w to stderr (with flush) ONLY if pops_trace_sf(); no-op otherwise.
Definition system_field_solver.hpp:71

pops
Definition amr_hierarchy.hpp:29

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

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

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

pops::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::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::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::derive_aux_polar
void derive_aux_polar(const MultiFab &phi, MultiFab &aux, const PolarGeometry &g)
Derives the POLAR aux in the local basis (e_r, e_theta) from the potential phi resolved by PolarPoiss...
Definition block_builder_polar.hpp:231

pops::BCType::Foextrap
@ Foextrap

pops::BCType::Dirichlet
@ Dirichlet

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

poisson_fft_solver.hpp
DIRECT EllipticSolver backends via spectral FFT (periodic BCs), MultiFab wrappers around PoissonFFT: ...

polar_poisson_solver.hpp
Direct POLAR Poisson solver on an annulus (PolarGeometry), Phase 2a.

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

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

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

pops::BCRec::yhi
BCType yhi
Definition physical_bc.hpp:31

pops::BCRec::xlo
BCType xlo
Definition physical_bc.hpp:30

pops::BCRec::ylo
BCType ylo
Definition physical_bc.hpp:31

pops::BCRec::xhi
BCType xhi
Definition physical_bc.hpp:30

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

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::field_solver::SystemFieldSolver::NamedField
Definition system_field_solver.hpp:179

pops::field_solver::SystemFieldSolver::NamedField::gy_comp
int gy_comp
Definition system_field_solver.hpp:182

pops::field_solver::SystemFieldSolver::NamedField::phi_comp
int phi_comp
Definition system_field_solver.hpp:180

pops::field_solver::SystemFieldSolver::NamedField::ell
std::optional< std::variant< GeometricMG, PoissonFFTSolver, RemappedFFTSolver > > ell
Definition system_field_solver.hpp:183

pops::field_solver::SystemFieldSolver::NamedField::gx_comp
int gx_comp
Definition system_field_solver.hpp:181

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

wall_predicate.hpp
ELLIPTIC conductor-wall predicate shared by the System and AmrSystem runtimes (the wall acts on the P...