PoPS/doxygen/polar__condensed__schur__source__stepper_8hpp_source.html

#pragma once


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

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

#include <pops/coupling/schur/core/schur_source_kernels.hpp>  // shared geometry-free kernels + validate_krylov_params (#263)

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

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

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

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

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

#include <pops/numerics/elliptic/polar/polar_tensor_operator.hpp>  // PolarTensorKrylovSolver, apply_polar_tensor (#210)

#include <pops/numerics/linalg/lorentz_eliminator.hpp>  // B^{-1} closed form (#118)

#include <pops/parallel/comm.hpp>


#include <stdexcept>


namespace pops {


namespace detail {


struct PolarSchurOperatorCoeffKernel {

  ConstArray4 s;

  ConstArray4 bz;

  Array4 arr, att;

  Array4 art, atr;

  Real c;

  Real th_dt;

  int c_rho;


  POPS_HD void operator()(int i, int j) const {

    const Real rho = s(i, j, c_rho);

    const LorentzEliminator le(th_dt, Real(1), bz(i, j, 0));  // w = th_dt * B_z

    const Real cr = c * rho;

    arr(i, j, 0) = Real(1) + cr * le.binv_11();

    att(i, j, 0) = Real(1) + cr * le.binv_22();

    art(i, j, 0) = cr * le.binv_12();

    atr(i, j, 0) = cr * le.binv_21();

  }


};


struct PolarSchurExplicitFluxKernel {

  ConstArray4 s;

  ConstArray4 bz;

  Array4 fr, ft;

  Real th_dt;

  int c_mx, c_my;


  POPS_HD void operator()(int i, int j) const {

    const LorentzEliminator le(th_dt, Real(1), bz(i, j, 0));

    Real Fr, Ft;

    le.apply_Binv(s(i, j, c_mx), s(i, j, c_my), Fr, Ft);  // B^{-1}(mr, mtheta) = rho B^{-1} v

    fr(i, j, 0) = Fr;

    ft(i, j, 0) = Ft;

  }


};


struct PolarSchurRhsAssembleKernel {

  ConstArray4 lap;

  ConstArray4 fr, ft;

  Array4 rhs;

  Real g;

  Real half_idr, half_idth;

  Real r_min, dr;


  POPS_HD void operator()(int i, int j) const {

    const Real ri = r_min + (i + Real(0.5)) * dr;

    const Real rip = r_min + (i + Real(1.5)) * dr;  // r_cell(i+1)

    const Real rim = r_min + (i - Real(0.5)) * dr;  // r_cell(i-1)

    const Real inv_r = Real(1) / ri;

    const Real div_r =

        (rip * fr(i + 1, j, 0) - rim * fr(i - 1, j, 0)) * half_idr;      // d_r(r F_r) centered

    const Real div_t = (ft(i, j + 1, 0) - ft(i, j - 1, 0)) * half_idth;  // d_theta(F_th) centered

    const Real divF = inv_r * (div_r + div_t);

    rhs(i, j, 0) = lap(i, j, 0) + g * divF;

  }


};


struct PolarSchurReconstructKernel {

  ConstArray4 phi;

  ConstArray4 vr, vt;

  ConstArray4 bz;

  Array4 st;

  Array4 nvr, nvt;

  Real th_dt;

  Real half_idr, half_idth;

  Real r_min, dr;

  int c_rho, c_mx, c_my;


  POPS_HD void operator()(int i, int j) const {

    const Real ri = r_min + (i + Real(0.5)) * dr;

    const Real gr = (phi(i + 1, j, 0) - phi(i - 1, j, 0)) * half_idr;          // d_r phi

    const Real gt = (phi(i, j + 1, 0) - phi(i, j - 1, 0)) * (half_idth / ri);  // (1/r) d_theta phi

    const Real rhsr = vr(i, j, 0) - th_dt * gr;  // (v^n - theta dt grad_polar phi)_r

    const Real rhst = vt(i, j, 0) - th_dt * gt;

    const LorentzEliminator le(th_dt, Real(1), bz(i, j, 0));  // w = th_dt * B_z

    Real nr, nt;

    le.apply_Binv(rhsr, rhst, nr, nt);  // v^{n+theta} = B^{-1}(v^n - theta dt grad_polar phi)

    nvr(i, j, 0) = nr;

    nvt(i, j, 0) = nt;

    const Real rho = st(i, j, c_rho);  // rho^n (frozen in the source)

    st(i, j, c_mx) = rho * nr;         // mom^{n+theta} = rho^n v^{n+theta}

    st(i, j, c_my) = rho * nt;

  }


};


// The geometry-free extrapolate / energy / extract-velocity / copy-Bz kernels are shared with the

// Cartesian stepper via <pops/coupling/schur/schur_source_kernels.hpp> (#263): detail::SchurExtrapolateScalarKernel,

// SchurExtrapolateVelocityKernel, SchurEnergyKernel, ExtractVelocityKernel, CopyBzKernel. Their member

// fields are named vx/vy but hold the polar velocity (vr/vtheta): the math is frame-independent. Only the

// metric-bearing kernels above (operator-coeff, explicit-flux, RHS-assemble, reconstruct) stay local.


struct PolarSchurCopyComp0Kernel {

  Array4 d;

  ConstArray4 s;

  POPS_HD void operator()(int i, int j) const { d(i, j, 0) = s(i, j, 0); }

};


}  // namespace detail


class PolarCondensedSchurSourceStepper {

 public:


  PolarCondensedSchurSourceStepper(const VariableSet& vars, const PolarGeometry& geom,

                                   const BoxArray& ba, const BCRec& bcPhi, Real alpha,

                                   PolarPrecond precond = PolarPrecond::RadialLine)

      : PolarCondensedSchurSourceStepper(

            vars, vars.index_of(VariableRole::Density), vars.index_of(VariableRole::MomentumX),

            vars.index_of(VariableRole::MomentumY), vars.index_of(VariableRole::Energy), geom, ba,

            bcPhi, alpha, precond) {}


  PolarCondensedSchurSourceStepper(const VariableSet& vars, int c_rho, int c_mx, int c_my, int c_E,

                                   const PolarGeometry& geom, const BoxArray& ba,

                                   const BCRec& bcPhi, Real alpha,

                                   PolarPrecond precond = PolarPrecond::RadialLine)

      : vars_(vars),

        c_rho_(c_rho),

        c_mx_(c_mx),

        c_my_(c_my),

        c_E_(c_E),

        alpha_(alpha),

        geom_(geom),

        bcPhi_(bcPhi),

        precond_(precond),

        ba_(ba),

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

        // condensed tensor coefficients A = I + c rho B^{-1} (1 ghost: operator faces)

        a_rr_(ba, dm_, 1, 1),

        a_tt_(ba, dm_, 1, 1),

        a_rt_(ba, dm_, 1, 1),

        a_tr_(ba, dm_, 1, 1),

        // condensed RHS buffers

        lap_(ba, dm_, 1, 0),

        rhs_(ba, dm_, 1, 0),

        fr_(ba, dm_, 1, 1),

        ft_(ba, dm_, 1, 1),

        // a MultiFab at 1 (A=I) for the scalar Lap_polar of the RHS

        one_rr_(ba, dm_, 1, 1),

        one_tt_(ba, dm_, 1, 1),

        bz_(ba, dm_, 1, 1),

        phi_n_(ba, dm_, 1, 1),

        vr_n_(ba, dm_, 1, 0),

        vt_n_(ba, dm_, 1, 0),

        vr_t_(ba, dm_, 1, 0),

        vt_t_(ba, dm_, 1, 0) {

    // MULTI-RANK MPI (theta-only split): no more single-rank safeguard. The layout constraint

    // (each box covers the full radial range, required by the RadialLine preconditioner of the elliptic

    // solve) is checked by the PolarTensorKrylovSolver built in step() (check_radial_columns)

    // -> a clear error is raised on all ranks if the split cuts r. Single-rank / single box:

    // path unchanged (the check passes trivially, all_reduce = identity in serial).

    if (c_rho_ < 0 || c_mx_ < 0 || c_my_ < 0)

      throw std::runtime_error(

          "PolarCondensedSchurSourceStepper: the fluid block must expose the roles Density, "

          "MomentumX "

          "(radial) and MomentumY (azimuthal).");

    one_rr_.set_val(Real(1));

    one_tt_.set_val(Real(1));

  }


  bool has_energy() const { return c_E_ >= 0; }


  void step(MultiFab& state, MultiFab& phi, const MultiFab& bz_field, int c_bz, Real theta,

            Real dt) {

    const Real th_dt = theta * dt;

    const Real c = theta * theta * dt * dt * alpha_;  // c = theta^2 dt^2 alpha

    const Real g = theta * dt * alpha_;               // theta dt alpha

    const Real dr = geom_.dr();

    const Real dth = geom_.dtheta();

    const Real half_idr = Real(1) / (Real(2) * dr);

    const Real half_idth = Real(1) / (Real(2) * dth);


    // -1) freeze phi^n (final extrapolation).

    copy_comp0(phi_n_, phi);


    // 0) extract v^n = (mr, mtheta)/rho and copy B_z into the internal buffer.

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

      const ConstArray4 s = state.fab(li).const_array();

      for_each_cell(state.box(li),

                    detail::ExtractVelocityKernel{s, vr_n_.fab(li).array(), vt_n_.fab(li).array(),

                                                  c_rho_, c_mx_, c_my_});

      for_each_cell(bz_.box(li), detail::CopyBzKernel{bz_field.fab(li).const_array(),

                                                      bz_.fab(li).array(), c_bz});

    }

    const BCRec ebc = coeff_bc(bcPhi_);

    device_fence();

    fill_ghosts(bz_, geom_.domain, ebc);


    // 1a) ASSEMBLE the coefficients A = I + c rho B^{-1} at the center (4 fields).

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

      const ConstArray4 s = state.fab(li).const_array();

      const ConstArray4 b = bz_.fab(li).const_array();

      for_each_cell(a_rr_.box(li),

                    detail::PolarSchurOperatorCoeffKernel{

                        s, b, a_rr_.fab(li).array(), a_tt_.fab(li).array(), a_rt_.fab(li).array(),

                        a_tr_.fab(li).array(), c, th_dt, c_rho_});

    }

    // ghosts of the coefficients (the operator face average reads the neighbor at +-1): theta

    // periodic, radial Foextrap. PolarTensorKrylovSolver::set_coefficients also fills these ghosts, but

    // we fill them here for robustness (assemble_rhs reads fr/ft -> consistency).

    device_fence();

    fill_ghosts(a_rr_, geom_.domain, ebc);

    fill_ghosts(a_tt_, geom_.domain, ebc);

    fill_ghosts(a_rt_, geom_.domain, ebc);

    fill_ghosts(a_tr_, geom_.domain, ebc);


    // 1b) ASSEMBLE the condensed right-hand side rhs_polar = Lap_polar phi^n + g div_polar(rho B^{-1} v^n).

    //     Lap_polar phi^n: apply_polar_tensor with A = I (coefficients at 1), SAME stencil as the solve.

    device_fence();

    fill_ghosts(phi, geom_.domain, phi_bc());  // ghosts of phi^n for the boundary Laplacian

    apply_polar_tensor(phi, geom_, lap_, &one_rr_, &one_tt_, nullptr,

                       nullptr);  // lap_ = Lap_polar phi^n

    // explicit flux F = B^{-1}(mr, mtheta) at the center (1 ghost for the centered div).

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

      const ConstArray4 s = state.fab(li).const_array();

      const ConstArray4 b = bz_.fab(li).const_array();

      for_each_cell(fr_.box(li),

                    detail::PolarSchurExplicitFluxKernel{s, b, fr_.fab(li).array(),

                                                         ft_.fab(li).array(), th_dt, c_mx_, c_my_});

    }

    device_fence();

    fill_ghosts(fr_, geom_.domain, ebc);  // centered div reads F(i+-1), F(j+-1)

    fill_ghosts(ft_, geom_.domain, ebc);

    // rhs_polar = Lap_polar phi^n + g div_polar F (sign L_int = -rhs_schur).

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

      for_each_cell(rhs_.box(li), detail::PolarSchurRhsAssembleKernel{

                                      lap_.fab(li).const_array(), fr_.fab(li).const_array(),

                                      ft_.fab(li).const_array(), rhs_.fab(li).array(), g, half_idr,

                                      half_idth, geom_.r_min, dr});


    // 2) SOLVE L_int(phi) = div(A grad phi) = rhs_polar via PolarTensorKrylovSolver.

    PolarTensorKrylovSolver kry(geom_, ba_, bcPhi_, precond_);

    const bool cross =

        (c != Real(0));  // non-trivial cross terms as soon as c != 0 (B_z != 0 -> binv_12 != 0)

    if (cross)

      kry.set_coefficients(&a_rr_, &a_tt_, &a_rt_, &a_tr_);

    else

      kry.set_coefficients(&a_rr_, &a_tt_);

    copy_comp0(kry.phi(), phi);  // warm start: phi^n -> kry.phi()

    copy_comp0(kry.rhs(), rhs_);

    last_result_ = kry.solve(krylov_tol_, krylov_max_iters_);

    copy_comp0(phi, kry.phi());  // phi <- phi^{n+theta}


    // 3) RECONSTRUCT v^{n+theta} = B^{-1}(v^n - theta dt grad_polar phi^{n+theta}); mom = rho v.

    device_fence();

    fill_ghosts(phi, geom_.domain, phi_bc());  // centered grad_polar reads phi(i+-1), phi(j+-1)

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

      for_each_cell(state.box(li),

                    detail::PolarSchurReconstructKernel{

                        phi.fab(li).const_array(), vr_n_.fab(li).const_array(),

                        vt_n_.fab(li).const_array(), bz_.fab(li).const_array(),

                        state.fab(li).array(), vr_t_.fab(li).array(), vt_t_.fab(li).array(), th_dt,

                        half_idr, half_idth, geom_.r_min, dr, c_rho_, c_mx_, c_my_});

    // vr_t_/vt_t_ carry v^{n+theta}.


    // 5) EXTRAPOLATE phi and v from the theta-stage to the full step: f^{n+1} = f^n + (1/theta)(f^{n+theta}-f^n).

    const Real inv_theta = Real(1) / theta;

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

      for_each_cell(phi.box(li), detail::SchurExtrapolateScalarKernel{

                                     phi_n_.fab(li).const_array(), phi.fab(li).array(), inv_theta});

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

      for_each_cell(state.box(li), detail::SchurExtrapolateVelocityKernel{

                                       vr_n_.fab(li).const_array(), vt_n_.fab(li).const_array(),

                                       vr_t_.fab(li).array(), vt_t_.fab(li).array(),

                                       state.fab(li).array(), inv_theta, c_rho_, c_mx_, c_my_});


    // 4) ENERGY (if role present): E^{n+1} = E^n + (1/2) rho (|v^{n+1}|^2 - |v^n|^2).

    if (c_E_ >= 0)

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

        for_each_cell(state.box(li), detail::SchurEnergyKernel{

                                         vr_n_.fab(li).const_array(), vt_n_.fab(li).const_array(),

                                         vr_t_.fab(li).const_array(), vt_t_.fab(li).const_array(),

                                         state.fab(li).array(), c_rho_, c_E_});


    // 6) FILL the ghosts of the state and the potential before returning.

    device_fence();

    fill_ghosts(state, geom_.domain, bcU_default());

    fill_ghosts(phi, geom_.domain, phi_bc());

  }


  const PolarKrylovResult& last_solve() const { return last_result_; }


  void set_krylov(Real tol, int max_iters) {

    detail::validate_krylov_params(tol, max_iters, "PolarCondensedSchurSourceStepper::set_krylov");

    krylov_tol_ = tol;

    krylov_max_iters_ = max_iters;

  }


  int density_comp() const { return c_rho_; }

  int momentum_x_comp() const { return c_mx_; }

  int momentum_y_comp() const { return c_my_; }

  int energy_comp() const { return c_E_; }


 private:

  BCRec bcU_default() const { return coeff_bc(bcPhi_); }


  static BCRec coeff_bc(const BCRec& bc) {

    auto fo = [](BCType t) { return t == BCType::Periodic ? t : BCType::Foextrap; };

    BCRec b;

    b.xlo = fo(bc.xlo);

    b.xhi = fo(bc.xhi);

    b.ylo = BCType::Periodic;

    b.yhi = BCType::Periodic;  // theta always periodic

    return b;

  }


  BCRec phi_bc() const {

    BCRec b = bcPhi_;

    b.ylo = BCType::Periodic;

    b.yhi = BCType::Periodic;  // theta always periodic

    b.ylo_val = Real(0);

    b.yhi_val = Real(0);

    return b;

  }


  void copy_comp0(MultiFab& dst, const MultiFab& src) {

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

      for_each_cell(dst.box(li), detail::PolarSchurCopyComp0Kernel{dst.fab(li).array(),

                                                                   src.fab(li).const_array()});

  }


  VariableSet vars_;

  int c_rho_, c_mx_, c_my_, c_E_;

  Real alpha_;

  PolarGeometry geom_;

  BCRec bcPhi_;

  PolarPrecond precond_;

  BoxArray ba_;

  DistributionMapping dm_;

  MultiFab a_rr_, a_tt_, a_rt_, a_tr_;

  MultiFab lap_, rhs_;

  MultiFab fr_, ft_;

  MultiFab one_rr_, one_tt_;

  MultiFab bz_;

  MultiFab phi_n_;

  MultiFab vr_n_, vt_n_;

  MultiFab vr_t_, vt_t_;

  PolarKrylovResult last_result_;

  Real krylov_tol_ = Real(1e-10);

  int krylov_max_iters_ = 600;

};


}  // namespace pops

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

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

pops::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::PolarCondensedSchurSourceStepper
SOURCE STAGE condensed by Schur in POLAR geometry, STANDALONE (transport frozen), GENERIC over any po...
Definition polar_condensed_schur_source_stepper.hpp:230

pops::PolarCondensedSchurSourceStepper::momentum_x_comp
int momentum_x_comp() const
Definition polar_condensed_schur_source_stepper.hpp:434

pops::PolarCondensedSchurSourceStepper::step
void step(MultiFab &state, MultiFab &phi, const MultiFab &bz_field, int c_bz, Real theta, Real dt)
POLAR condensed SOURCE STAGE, IN-PLACE on state and phi.
Definition polar_condensed_schur_source_stepper.hpp:304

pops::PolarCondensedSchurSourceStepper::set_krylov
void set_krylov(Real tol, int max_iters)
Tolerance / iteration budget of the polar Krylov solve.
Definition polar_condensed_schur_source_stepper.hpp:427

pops::PolarCondensedSchurSourceStepper::PolarCondensedSchurSourceStepper
PolarCondensedSchurSourceStepper(const VariableSet &vars, int c_rho, int c_mx, int c_my, int c_E, const PolarGeometry &geom, const BoxArray &ba, const BCRec &bcPhi, Real alpha, PolarPrecond precond=PolarPrecond::RadialLine)
EXPLICIT-COMPONENT variant (audit wave 3, parity with the Cartesian stepper): the caller DESIGNATES t...
Definition polar_condensed_schur_source_stepper.hpp:249

pops::PolarCondensedSchurSourceStepper::density_comp
int density_comp() const
Definition polar_condensed_schur_source_stepper.hpp:433

pops::PolarCondensedSchurSourceStepper::last_solve
const PolarKrylovResult & last_solve() const
Diagnostic of the last solve (BiCGStab iterations, relative residual, convergence).
Definition polar_condensed_schur_source_stepper.hpp:423

pops::PolarCondensedSchurSourceStepper::PolarCondensedSchurSourceStepper
PolarCondensedSchurSourceStepper(const VariableSet &vars, const PolarGeometry &geom, const BoxArray &ba, const BCRec &bcPhi, Real alpha, PolarPrecond precond=PolarPrecond::RadialLine)
vars: descriptor of the fluid block; MUST expose Density / MomentumX / MomentumY (Energy optional).
Definition polar_condensed_schur_source_stepper.hpp:238

pops::PolarCondensedSchurSourceStepper::has_energy
bool has_energy() const
true if the model carries an Energy role (energy update active).
Definition polar_condensed_schur_source_stepper.hpp:298

pops::PolarCondensedSchurSourceStepper::energy_comp
int energy_comp() const
Definition polar_condensed_schur_source_stepper.hpp:436

pops::PolarCondensedSchurSourceStepper::momentum_y_comp
int momentum_y_comp() const
Definition polar_condensed_schur_source_stepper.hpp:435

pops::PolarTensorKrylovSolver
MATRIX-FREE BiCGStab Krylov solver for the FULL-tensor POLAR elliptic operator L_int(phi) = div(A gra...
Definition polar_tensor_operator.hpp:332

pops::PolarTensorKrylovSolver::set_coefficients
void set_coefficients(MultiFab *a_rr, MultiFab *a_tt, MultiFab *a_rt=nullptr, MultiFab *a_tr=nullptr)
Sets the tensor A coefficients.
Definition polar_tensor_operator.hpp:391

pops::PolarTensorKrylovSolver::rhs
MultiFab & rhs()
Definition polar_tensor_operator.hpp:383

pops::PolarTensorKrylovSolver::phi
MultiFab & phi()
Definition polar_tensor_operator.hpp:384

pops::PolarTensorKrylovSolver::solve
void solve()
Definition polar_tensor_operator.hpp:409

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

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

lorentz_eliminator.hpp
LorentzEliminator: 2x2 operator B of the Schur scheme for implicit velocity elimination.

mf_arith.hpp
MultiFab arithmetic (saxpy, lincomb, norm_inf, dot) over VALID cells.

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

pops::detail::validate_krylov_params
void validate_krylov_params(Real tol, int max_iters, const char *who)
Validates the Krylov tolerance / iteration budget shared by the Schur source steppers (historical con...
Definition schur_source_kernels.hpp:103

pops
Definition amr_hierarchy.hpp:29

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

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

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

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::VariableRole
VariableRole
PHYSICAL role of a component.
Definition variables.hpp:27

pops::VariableRole::Energy
@ Energy

pops::VariableRole::Density
@ Density

pops::VariableRole::MomentumX
@ MomentumX

pops::VariableRole::MomentumY
@ MomentumY

pops::apply_polar_tensor
void apply_polar_tensor(const MultiFab &phi, const PolarGeometry &geom, MultiFab &lap, const MultiFab *a_rr, const MultiFab *a_tt, const MultiFab *a_rt, const MultiFab *a_tr)
Applies L_int(phi) = div(A grad phi) in polar over the whole MultiFab.
Definition polar_tensor_operator.hpp:259

pops::PolarPrecond
PolarPrecond
Choice of the SIMPLE BiCGStab PRECONDITIONER (NO MG V-cycle – stagnation on polar 1/r^2,...
Definition polar_tensor_operator.hpp:319

pops::PolarPrecond::RadialLine
@ RadialLine

pops::BCType
BCType
Boundary condition type for a face: Periodic (handled by fill_boundary), Foextrap (zero gradient,...
Definition physical_bc.hpp:25

pops::BCType::Foextrap
@ Foextrap

pops::BCType::Periodic
@ Periodic

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

polar_tensor_operator.hpp
Iterative POLAR elliptic operator with anisotropic TENSOR coefficient (cross terms).

schur_source_kernels.hpp
Geometry-INDEPENDENT device kernels shared by the Schur SOURCE-STAGE steppers: the Cartesian condense...

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::ylo
BCType ylo
Definition physical_bc.hpp:31

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

pops::LorentzEliminator
LorentzEliminator: operator B = [[1,-w],[w,1]] and its analytic inverse.
Definition lorentz_eliminator.hpp:55

pops::LorentzEliminator::binv_11
POPS_HD Real binv_11() const
Definition lorentz_eliminator.hpp:81

pops::LorentzEliminator::binv_21
POPS_HD Real binv_21() const
Definition lorentz_eliminator.hpp:83

pops::LorentzEliminator::binv_12
POPS_HD Real binv_12() const
Definition lorentz_eliminator.hpp:82

pops::LorentzEliminator::apply_Binv
POPS_HD void apply_Binv(Real vx, Real vy, Real &vxp, Real &vyp) const
apply_Binv: applies B^{-1} = (1/det)*[[1,w],[-w,1]] to (vx, vy), writes (vxp, vyp)....
Definition lorentz_eliminator.hpp:73

pops::LorentzEliminator::binv_22
POPS_HD Real binv_22() const
Definition lorentz_eliminator.hpp:84

pops::PolarGeometry
Definition geometry.hpp:59

pops::PolarGeometry::dr
POPS_HD Real dr() const
Definition geometry.hpp:70

pops::PolarGeometry::domain
Box2D domain
nx() = nr (radial cells), ny() = ntheta (azimuthal cells)
Definition geometry.hpp:60

pops::PolarGeometry::dtheta
POPS_HD Real dtheta() const
Definition geometry.hpp:71

pops::PolarKrylovResult
Result of a polar BiCGStab solve: iterations, relative residual, convergence.
Definition polar_tensor_operator.hpp:304

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

pops::detail::CopyBzKernel
Copies the B_z field (aux channel) into an internal scalar MultiFab (0 ghost is enough,...
Definition schur_source_kernels.hpp:93

pops::detail::ExtractVelocityKernel
Extracts the velocity v = (mx, my) / rho from the state (Density / MomentumX / MomentumY roles) into ...
Definition schur_source_kernels.hpp:79

pops::detail::PolarSchurCopyComp0Kernel
dst <- src (component 0).
Definition polar_condensed_schur_source_stepper.hpp:211

pops::detail::PolarSchurCopyComp0Kernel::d
Array4 d
Definition polar_condensed_schur_source_stepper.hpp:212

pops::detail::PolarSchurCopyComp0Kernel::operator()
POPS_HD void operator()(int i, int j) const
Definition polar_condensed_schur_source_stepper.hpp:214

pops::detail::PolarSchurCopyComp0Kernel::s
ConstArray4 s
Definition polar_condensed_schur_source_stepper.hpp:213

pops::detail::PolarSchurExplicitFluxKernel
EXPLICIT flux F = rho B^{-1} v^n = B^{-1}(mr, mtheta) at the center (physical components (e_r,...
Definition polar_condensed_schur_source_stepper.hpp:132

pops::detail::PolarSchurExplicitFluxKernel::c_my
int c_my
MomentumX (= radial), MomentumY (= azimuthal) components.
Definition polar_condensed_schur_source_stepper.hpp:137

pops::detail::PolarSchurExplicitFluxKernel::th_dt
Real th_dt
theta * dt (w = th_dt * B_z)
Definition polar_condensed_schur_source_stepper.hpp:136

pops::detail::PolarSchurExplicitFluxKernel::bz
ConstArray4 bz
B_z field at the center.
Definition polar_condensed_schur_source_stepper.hpp:134

pops::detail::PolarSchurExplicitFluxKernel::s
ConstArray4 s
fluid state (mr, mtheta read at components c_mx, c_my)
Definition polar_condensed_schur_source_stepper.hpp:133

pops::detail::PolarSchurExplicitFluxKernel::c_mx
int c_mx
Definition polar_condensed_schur_source_stepper.hpp:137

pops::detail::PolarSchurExplicitFluxKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition polar_condensed_schur_source_stepper.hpp:138

pops::detail::PolarSchurExplicitFluxKernel::ft
Array4 ft
output: F_r, F_theta = B^{-1}(mr, mtheta)
Definition polar_condensed_schur_source_stepper.hpp:135

pops::detail::PolarSchurExplicitFluxKernel::fr
Array4 fr
Definition polar_condensed_schur_source_stepper.hpp:135

pops::detail::PolarSchurOperatorCoeffKernel
POLAR condensed tensor coefficients A = I + c rho B^{-1} at cell centers.
Definition polar_condensed_schur_source_stepper.hpp:110

pops::detail::PolarSchurOperatorCoeffKernel::s
ConstArray4 s
fluid state (read rho)
Definition polar_condensed_schur_source_stepper.hpp:111

pops::detail::PolarSchurOperatorCoeffKernel::c_rho
int c_rho
Density component.
Definition polar_condensed_schur_source_stepper.hpp:117

pops::detail::PolarSchurOperatorCoeffKernel::att
Array4 att
output: a_rr, a_tt (diagonal of A)
Definition polar_condensed_schur_source_stepper.hpp:113

pops::detail::PolarSchurOperatorCoeffKernel::art
Array4 art
Definition polar_condensed_schur_source_stepper.hpp:114

pops::detail::PolarSchurOperatorCoeffKernel::bz
ConstArray4 bz
B_z field at the center.
Definition polar_condensed_schur_source_stepper.hpp:112

pops::detail::PolarSchurOperatorCoeffKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition polar_condensed_schur_source_stepper.hpp:118

pops::detail::PolarSchurOperatorCoeffKernel::arr
Array4 arr
Definition polar_condensed_schur_source_stepper.hpp:113

pops::detail::PolarSchurOperatorCoeffKernel::c
Real c
c = theta^2 dt^2 alpha
Definition polar_condensed_schur_source_stepper.hpp:115

pops::detail::PolarSchurOperatorCoeffKernel::atr
Array4 atr
output: cross terms a_rt, a_tr
Definition polar_condensed_schur_source_stepper.hpp:114

pops::detail::PolarSchurOperatorCoeffKernel::th_dt
Real th_dt
theta * dt (w = th_dt * B_z, binv depends only on w)
Definition polar_condensed_schur_source_stepper.hpp:116

pops::detail::PolarSchurReconstructKernel
Reconstructs v^{n+theta} = B^{-1}(v^n - theta dt grad_polar phi^{n+theta}) and writes mom = rho^n v^{...
Definition polar_condensed_schur_source_stepper.hpp:176

pops::detail::PolarSchurReconstructKernel::vt
ConstArray4 vt
v^n (components 0: velocity, NOT momentum)
Definition polar_condensed_schur_source_stepper.hpp:178

pops::detail::PolarSchurReconstructKernel::th_dt
Real th_dt
theta * dt (w = th_dt * B_z, and gradient factor)
Definition polar_condensed_schur_source_stepper.hpp:182

pops::detail::PolarSchurReconstructKernel::c_my
int c_my
Density / MomentumX (radial) / MomentumY (azimuthal) components.
Definition polar_condensed_schur_source_stepper.hpp:185

pops::detail::PolarSchurReconstructKernel::r_min
Real r_min
Definition polar_condensed_schur_source_stepper.hpp:184

pops::detail::PolarSchurReconstructKernel::half_idr
Real half_idr
Definition polar_condensed_schur_source_stepper.hpp:183

pops::detail::PolarSchurReconstructKernel::nvt
Array4 nvt
output: v^{n+theta} (component 0) for the energy / the diagnostic
Definition polar_condensed_schur_source_stepper.hpp:181

pops::detail::PolarSchurReconstructKernel::vr
ConstArray4 vr
Definition polar_condensed_schur_source_stepper.hpp:178

pops::detail::PolarSchurReconstructKernel::c_mx
int c_mx
Definition polar_condensed_schur_source_stepper.hpp:185

pops::detail::PolarSchurReconstructKernel::phi
ConstArray4 phi
phi^{n+theta} (ghosts filled: centered grad reads i+-1, j+-1)
Definition polar_condensed_schur_source_stepper.hpp:177

pops::detail::PolarSchurReconstructKernel::half_idth
Real half_idth
1/(2 dr), 1/(2 dtheta)
Definition polar_condensed_schur_source_stepper.hpp:183

pops::detail::PolarSchurReconstructKernel::bz
ConstArray4 bz
B_z field at the center.
Definition polar_condensed_schur_source_stepper.hpp:179

pops::detail::PolarSchurReconstructKernel::nvr
Array4 nvr
Definition polar_condensed_schur_source_stepper.hpp:181

pops::detail::PolarSchurReconstructKernel::st
Array4 st
fluid state (WRITE mr, mtheta; READ rho)
Definition polar_condensed_schur_source_stepper.hpp:180

pops::detail::PolarSchurReconstructKernel::dr
Real dr
for r_cell(i) (azimuthal metric 1/r)
Definition polar_condensed_schur_source_stepper.hpp:184

pops::detail::PolarSchurReconstructKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition polar_condensed_schur_source_stepper.hpp:186

pops::detail::PolarSchurReconstructKernel::c_rho
int c_rho
Definition polar_condensed_schur_source_stepper.hpp:185

pops::detail::PolarSchurRhsAssembleKernel
rhs_polar(i,j) = lap_polar(i,j) (= Lap_polar phi^n) + g * div_polar F, second-order centered POLAR di...
Definition polar_condensed_schur_source_stepper.hpp:153

pops::detail::PolarSchurRhsAssembleKernel::lap
ConstArray4 lap
Lap_polar phi^n (positive sign, A=I)
Definition polar_condensed_schur_source_stepper.hpp:154

pops::detail::PolarSchurRhsAssembleKernel::half_idth
Real half_idth
1/(2 dr), 1/(2 dtheta)
Definition polar_condensed_schur_source_stepper.hpp:158

pops::detail::PolarSchurRhsAssembleKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition polar_condensed_schur_source_stepper.hpp:160

pops::detail::PolarSchurRhsAssembleKernel::fr
ConstArray4 fr
Definition polar_condensed_schur_source_stepper.hpp:155

pops::detail::PolarSchurRhsAssembleKernel::r_min
Real r_min
Definition polar_condensed_schur_source_stepper.hpp:159

pops::detail::PolarSchurRhsAssembleKernel::dr
Real dr
for r_cell(i), r_cell(i+-1)
Definition polar_condensed_schur_source_stepper.hpp:159

pops::detail::PolarSchurRhsAssembleKernel::half_idr
Real half_idr
Definition polar_condensed_schur_source_stepper.hpp:158

pops::detail::PolarSchurRhsAssembleKernel::g
Real g
theta dt alpha
Definition polar_condensed_schur_source_stepper.hpp:157

pops::detail::PolarSchurRhsAssembleKernel::ft
ConstArray4 ft
flux F at the center (ghosts filled)
Definition polar_condensed_schur_source_stepper.hpp:155

pops::detail::PolarSchurRhsAssembleKernel::rhs
Array4 rhs
output: condensed right-hand side (L_int sign)
Definition polar_condensed_schur_source_stepper.hpp:156

pops::detail::SchurEnergyKernel
Energy update: E^{n+1} = E^n + (1/2) rho^n (|v^{n+1}|^2 - |v^n|^2).
Definition schur_source_kernels.hpp:62

pops::detail::SchurExtrapolateScalarKernel
Linear extrapolation of a SCALAR field from the theta-stage to the full step: f^{n+1} = f^n + (1/thet...
Definition schur_source_kernels.hpp:32

pops::detail::SchurExtrapolateVelocityKernel
Linear extrapolation of the VELOCITY (vx, vy) from the theta-stage to the full step,...
Definition schur_source_kernels.hpp:43

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

POPS_HD
#define POPS_HD
Definition types.hpp:25

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