PoPS/doxygen/poisson__operator_8hpp_source.html

#pragma once


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

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

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

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

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

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


namespace pops {


// Harmonic mean of two center permittivities -> face permittivity.

// Guard against division by 0 if both centers are zero (inactive cell).


POPS_HD inline Real eps_harmonic(Real ec, Real ev) {

  const Real s = ec + ev;

  return s > Real(0) ? Real(2) * ec * ev / s : Real(0);

}


namespace detail {

// Weights of the FOUR faces (xm, xp, ym, yp) at (i,j) for the FACE permittivity path (he), with or

// without cut-cell. Each face = HARMONIC mean of the two adjacent centers (eps_x for x faces,

// eps_y for y faces); in cut-cell (hc) the face Shortley-Weller weight is multiplied by its

// permittivity, otherwise 1/h^2 (idx2 / idy2) is applied. Free POPS_HD functor (device-clean) shared

// by the three he kernels (apply / residual / smoother); body STRICTLY identical to the three

// original copies -> bit-identical output. Output by reference (wxm/wxp/wym/wyp).


POPS_HD inline void face_weights(const ConstArray4& ep, const ConstArray4& ey, int i, int j,

                                Real idx2, Real idy2, bool hc, const ConstArray4& cf, Real& wxm,

                                Real& wxp, Real& wym, Real& wyp) {

  const Real ec = ep(i, j);   // eps_x at center (x faces)

  const Real ecy = ey(i, j);  // eps_y at center (y faces); == ec when isotropic

  const Real exm = eps_harmonic(ec, ep(i - 1, j));

  const Real exp = eps_harmonic(ec, ep(i + 1, j));

  const Real eym = eps_harmonic(ecy, ey(i, j - 1));

  const Real eyp = eps_harmonic(ecy, ey(i, j + 1));

  if (hc) {  // cut-cell: eps_face multiplies each Shortley-Weller weight

    wxm = cf(i, j, 0) * exm;

    wxp = cf(i, j, 1) * exp;

    wym = cf(i, j, 2) * eym;

    wyp = cf(i, j, 3) * eyp;

  } else {  // 5-point stencil with variable face coefficient

    wxm = exm * idx2;

    wxp = exp * idx2;

    wym = eym * idy2;

    wyp = eyp * idy2;

  }

}


// Divergence of the CROSS FLUXES of the full tensor at (i,j): d_x(Axy d_y phi) + d_y(Ayx d_x phi),

// discretized by finite volumes (9-point stencil) as described in the header. Returns the

// contribution to be ADDED to div(A grad phi). Free POPS_HD functor (device-clean) shared by

// apply_laplacian and poisson_residual; face coefficient = ARITHMETIC mean. hxy/hyx

// guard each half-term: an ABSENT coefficient (field not provided) contributes 0 without deref.


POPS_HD inline Real cross_div(const ConstArray4& p, bool hxy, const ConstArray4& axy, bool hyx,

                             const ConstArray4& ayx, int i, int j, Real idx, Real idy) {

  Real out = Real(0);

  if (hxy) {  // x faces: cross flux = Axy_face * (d_y phi)_face, tangential averaged over 4 corners.

    const Real axy_xp = Real(0.5) * (axy(i, j) + axy(i + 1, j));

    const Real axy_xm = Real(0.5) * (axy(i, j) + axy(i - 1, j));

    const Real dyf_xp =

        (p(i, j + 1) + p(i + 1, j + 1) - p(i, j - 1) - p(i + 1, j - 1)) * (Real(0.25) * idy);

    const Real dyf_xm =

        (p(i - 1, j + 1) + p(i, j + 1) - p(i - 1, j - 1) - p(i, j - 1)) * (Real(0.25) * idy);

    out += (axy_xp * dyf_xp - axy_xm * dyf_xm) * idx;

  }

  if (hyx) {  // y faces: cross flux = Ayx_face * (d_x phi)_face.

    const Real ayx_yp = Real(0.5) * (ayx(i, j) + ayx(i, j + 1));

    const Real ayx_ym = Real(0.5) * (ayx(i, j) + ayx(i, j - 1));

    const Real dxf_yp =

        (p(i + 1, j) + p(i + 1, j + 1) - p(i - 1, j) - p(i - 1, j + 1)) * (Real(0.25) * idx);

    const Real dxf_ym =

        (p(i + 1, j - 1) + p(i + 1, j) - p(i - 1, j - 1) - p(i - 1, j)) * (Real(0.25) * idx);

    out += (ayx_yp * dxf_yp - ayx_ym * dxf_ym) * idy;

  }

  return out;

}


// NAMED FUNCTORS (and not POPS_HD lambdas) for the Poisson operator and Gauss-Seidel smoother kernels.

// Same reasons as the rest of the elliptic path (#93, recipe #64): these kernels are

// first-instantiated from the MG V-cycle pulled from an external TU (harness / native loader); an extended

// lambda there breaks the device kernel emission under nvcc (null kernel-stub -> Cuda segfault in

// Release -O without -g). Body STRICTLY identical to the former lambdas (same he/hc/hk branches,

// same stencil) -> bit-identical residual and potential on CPU and device.


// L = div(A grad phi) - kappa phi (apply_laplacian). cf/ep/ey/ka unused if the flag is false.

// hxy/hyx => FULL tensor: we ADD the cross fluxes d_x(Axy d_y phi) + d_y(Ayx d_x phi) (idx/idy

// = 1/dx, 1/dy; axy/ayx = off-diagonal coefficients at center). hxy=hyx=false => bit-identical.


struct ApplyLaplacianKernel {

  ConstArray4 p;

  Array4 L;

  Real idx2, idy2, idx, idy;

  bool hc;

  ConstArray4 cf;

  bool he;

  ConstArray4 ep, ey;

  bool hk;

  ConstArray4 ka;

  bool hxy, hyx;

  ConstArray4 axy, ayx;

  int c;  // component the matvec acts on; 0 for the scalar Poisson path (bit-identical 2-arg access)


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

    if (he) {  // face permittivity (harmonic), with or without cut-cell (coefficient path: comp 0 only)

      Real wxm, wxp, wym, wyp;

      face_weights(ep, ey, i, j, idx2, idy2, hc, cf, wxm, wxp, wym, wyp);

      L(i, j) = wxp * p(i + 1, j) + wxm * p(i - 1, j) + wyp * p(i, j + 1) + wym * p(i, j - 1) -

                (wxm + wxp + wym + wyp) * p(i, j);

    } else if (hc)  // cut-cell coefficient path (comp 0 only)

      L(i, j) = cf(i, j, 1) * p(i + 1, j) + cf(i, j, 0) * p(i - 1, j) + cf(i, j, 3) * p(i, j + 1) +

                cf(i, j, 2) * p(i, j - 1) - cf(i, j, 4) * p(i, j);

    else  // bare 5-point stencil, applied PER COMPONENT (c); c==0 => the scalar path, bit-identical.

      L(i, j, c) = (p(i + 1, j, c) - 2 * p(i, j, c) + p(i - 1, j, c)) * idx2 +

                   (p(i, j + 1, c) - 2 * p(i, j, c) + p(i, j - 1, c)) * idy2;

    // FULL block: ADDITIVE cross fluxes (after the diagonal stencil). hxy=hyx=false => +0, bit-identical.

    if (hxy || hyx)

      L(i, j) += cross_div(p, hxy, axy, hyx, ayx, i, j, idx, idy);

    // Helmholtz / screened operator: L phi = div(A grad phi) - kappa phi.

    if (hk)

      L(i, j) -= ka(i, j) * p(i, j);

  }


};


// res = f - L phi on active cells, 0 on conductor cells (poisson_residual).

// hx => FULL tensor: ADDITIVE cross fluxes (cf. ApplyLaplacianKernel). hx=false => bit-identical.


struct PoissonResidualKernel {

  ConstArray4 p, ff;

  Array4 r;

  Real idx2, idy2, idx, idy;

  bool hm;

  ConstArray4 mk;

  bool hc;

  ConstArray4 cf;

  bool he;

  ConstArray4 ep, ey;

  bool hk;

  ConstArray4 ka;

  bool hxy, hyx;

  ConstArray4 axy, ayx;


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

    if (hm && mk(i, j) == Real(0)) {

      r(i, j) = 0;

      return;

    }

    Real lap;

    if (he) {  // face permittivity (harmonic), with or without cut-cell

      Real wxm, wxp, wym, wyp;

      face_weights(ep, ey, i, j, idx2, idy2, hc, cf, wxm, wxp, wym, wyp);

      lap = wxp * p(i + 1, j) + wxm * p(i - 1, j) + wyp * p(i, j + 1) + wym * p(i, j - 1) -

            (wxm + wxp + wym + wyp) * p(i, j);

    } else if (hc)

      lap = cf(i, j, 1) * p(i + 1, j) + cf(i, j, 0) * p(i - 1, j) + cf(i, j, 3) * p(i, j + 1) +

            cf(i, j, 2) * p(i, j - 1) - cf(i, j, 4) * p(i, j);

    else

      lap = (p(i + 1, j) - 2 * p(i, j) + p(i - 1, j)) * idx2 +

            (p(i, j + 1) - 2 * p(i, j) + p(i, j - 1)) * idy2;

    // FULL block: ADDITIVE cross fluxes (after the diagonal stencil). hxy=hyx=false => +0, bit-identical.

    if (hxy || hyx)

      lap += cross_div(p, hxy, axy, hyx, ayx, i, j, idx, idy);

    // res = f - L phi, L phi = div(A grad phi) - kappa phi = lap - kappa phi.

    r(i, j) = ff(i, j) - lap + (hk ? ka(i, j) * p(i, j) : Real(0));

  }


};


}  // namespace detail


// a_xy/a_yx: off-diagonal coefficients (FULL tensor). nullptr => cross term absent

// (bit-identical diagonal/Poisson operator). Ghosts (1 layer) assumed filled by the caller.


inline void apply_laplacian(const MultiFab& phi, const Geometry& geom, MultiFab& lap,

                            const MultiFab* coef = nullptr, const MultiFab* eps = nullptr,

                            const MultiFab* kappa = nullptr, const MultiFab* eps_y = nullptr,

                            const MultiFab* a_xy = nullptr, const MultiFab* a_yx = nullptr) {

  const Real idx2 = Real(1) / (geom.dx() * geom.dx());

  const Real idy2 = Real(1) / (geom.dy() * geom.dy());

  const Real idx = Real(1) / geom.dx();

  const Real idy = Real(1) / geom.dy();

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

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

    Array4 L = lap.fab(li).array();

    const Box2D v = lap.box(li);

    const bool hc = coef != nullptr;

    const ConstArray4 cf = hc ? coef->fab(li).const_array() : ConstArray4{};

    const bool he = eps != nullptr;

    const ConstArray4 ep = he ? eps->fab(li).const_array() : ConstArray4{};

    // eps_y==nullptr => isotropic: y faces read the same field as the x faces (eps_x).

    const ConstArray4 ey = (he && eps_y) ? eps_y->fab(li).const_array() : ep;

    const bool hk = kappa != nullptr;  // reaction term -kappa phi

    const ConstArray4 ka = hk ? kappa->fab(li).const_array() : ConstArray4{};

    const bool hxy = a_xy != nullptr;  // Axy cross half-term (x faces)

    const bool hyx = a_yx != nullptr;  // Ayx cross half-term (y faces)

    const ConstArray4 axy = hxy ? a_xy->fab(li).const_array() : ConstArray4{};

    const ConstArray4 ayx = hyx ? a_yx->fab(li).const_array() : ConstArray4{};

    // Bare 5-point matvec acts on EVERY component (a vector / state matrix-free operator: the

    // condensed-Schur block unknown). The coefficient / cross / Helmholtz branches are single-component

    // (the scalar Poisson operator), so they run for component 0 only; nc==1 reproduces the old

    // single-pass, bit-identical scalar path.

    const int nc = (he || hc || hk || hxy || hyx) ? 1 : lap.ncomp();

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

      for_each_cell(v, detail::ApplyLaplacianKernel{p, L, idx2, idy2, idx, idy, hc, cf, he, ep, ey,

                                                    hk, ka, hxy, hyx, axy, ayx, c});

  }

}


namespace detail {

// Centered finite-volume divergence of a cell-centered vector flux: the x-flux is read from

// component @c cx of fx and the y-flux from component @c cy of fy:

//   div(i,j) = (fx(i+1,j,cx) - fx(i-1,j,cx)) / (2 dx) + (fy(i,j+1,cy) - fy(i,j-1,cy)) / (2 dy),

// second order, the exact inverse of the centered gradient (field_postprocess). Named functor

// (and not an POPS_HD lambda) for the same reason as the other elliptic kernels (#93): it is

// first-instantiated from an external TU (the compiled time-program .so) and an extended lambda

// breaks the device kernel emission under nvcc. half_idx/half_idy = 1/(2 dx), 1/(2 dy).


struct DivergenceKernel {

  ConstArray4 fx, fy;

  Array4 div;

  Real half_idx, half_idy;

  int cx, cy;


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

    div(i, j) = (fx(i + 1, j, cx) - fx(i - 1, j, cx)) * half_idx +

                (fy(i, j + 1, cy) - fy(i, j - 1, cy)) * half_idy;

  }


};


}  // namespace detail


// Centered FV divergence of a cell-centered vector flux into div_out (component 0):

//   div = d fx/dx + d fy/dy, centered (the matching inverse of the centered gradient). The x-flux is

// read from component @p cx of @p fx and the y-flux from component @p cy of @p fy, so the 2-component

// output of the centered gradient (field_postprocess: d/dx in component 0, d/dy in component 1) can be

// fed back as a SINGLE field passed for both arguments (apply_divergence(g, g, geom, out, 0, 1)) to

// recover the 5-point Laplacian, or two distinct single-component fluxes with cx = cy = 0. Ghosts (1

// layer) assumed filled by the caller. Mirrors apply_laplacian's structure (one local-fab loop, named

// functor). The Schur condensation builder uses the same stencil inline (coupling/schur/core/

// schur_condensation.hpp SchurRhsAssembleKernel, which fuses it with -Lap phi^n -- not refactored to

// call this so the native source path stays bit-identical).


inline void apply_divergence(const MultiFab& fx, const MultiFab& fy, const Geometry& geom,

                             MultiFab& div_out, int cx = 0, int cy = 0) {

  const Real half_idx = Real(1) / (Real(2) * geom.dx());

  const Real half_idy = Real(1) / (Real(2) * geom.dy());

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

    const ConstArray4 fxv = fx.fab(li).const_array();

    const ConstArray4 fyv = fy.fab(li).const_array();

    Array4 d = div_out.fab(li).array();

    for_each_cell(div_out.box(li),

                  detail::DivergenceKernel{fxv, fyv, d, half_idx, half_idy, cx, cy});

  }

}


// res = f - div(A grad phi) on active cells, 0 on conductor cells.

// a_xy/a_yx: off-diagonal coefficients (cf. apply_laplacian). nullptr => bit-identical.


inline void poisson_residual(MultiFab& phi, const MultiFab& f, const Geometry& geom,

                             const BCRec& bc, MultiFab& res, const MultiFab* mask = nullptr,

                             const MultiFab* coef = nullptr, const MultiFab* eps = nullptr,

                             const MultiFab* kappa = nullptr, const MultiFab* eps_y = nullptr,

                             const MultiFab* a_xy = nullptr, const MultiFab* a_yx = nullptr) {

  device_fence();  // GPU: phi may have been written by a kernel (smoother); we

                   // wait before the host read in fill_ghosts.

  fill_ghosts(phi, geom.domain, bc);

  const Real idx2 = Real(1) / (geom.dx() * geom.dx());

  const Real idy2 = Real(1) / (geom.dy() * geom.dy());

  const Real idx = Real(1) / geom.dx();

  const Real idy = Real(1) / geom.dy();

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

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

    const ConstArray4 ff = f.fab(li).const_array();

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

    const Box2D v = res.box(li);

    const bool hm = mask != nullptr;

    const ConstArray4 mk = hm ? mask->fab(li).const_array() : ConstArray4{};

    const bool hc = coef != nullptr;

    const ConstArray4 cf = hc ? coef->fab(li).const_array() : ConstArray4{};

    const bool he = eps != nullptr;

    const ConstArray4 ep = he ? eps->fab(li).const_array() : ConstArray4{};

    // eps_y==nullptr => isotropic: y faces read the same field as the x faces (eps_x).

    const ConstArray4 ey = (he && eps_y) ? eps_y->fab(li).const_array() : ep;

    const bool hk = kappa != nullptr;  // reaction term -kappa phi

    const ConstArray4 ka = hk ? kappa->fab(li).const_array() : ConstArray4{};

    const bool hxy = a_xy != nullptr;  // Axy cross half-term (x faces)

    const bool hyx = a_yx != nullptr;  // Ayx cross half-term (y faces)

    const ConstArray4 axy = hxy ? a_xy->fab(li).const_array() : ConstArray4{};

    const ConstArray4 ayx = hyx ? a_yx->fab(li).const_array() : ConstArray4{};

    for_each_cell(v,

                  detail::PoissonResidualKernel{p,  ff, r,  idx2, idy2, idx, idy, hm,  mk,  hc,

                                                cf, he, ep, ey,   hk,   ka,  hxy, hyx, axy, ayx});

  }

}


namespace detail {

// Red-black Gauss-Seidel smoother on one color (gs_color). p is WRITTEN in place. Body identical to

// the former lambda -> bit-identical. See the comment of the other kernels (#93) for the motivation

// of the named functor.


struct GsColorKernel {

  Array4 p;

  ConstArray4 ff;

  Real idx2, idy2, diag0;

  int color;

  bool hm;

  ConstArray4 mk;

  bool hc;

  ConstArray4 cf;

  bool he;

  ConstArray4 ep, ey;

  bool hk;

  ConstArray4 ka;


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

    if (((i + j) & 1) != color)

      return;

    if (hm && mk(i, j) == Real(0))

      return;  // conductor: pins phi=0

    Real off, diag;

    if (he) {  // face permittivity (harmonic), with or without cut-cell

      Real wxm, wxp, wym, wyp;

      face_weights(ep, ey, i, j, idx2, idy2, hc, cf, wxm, wxp, wym, wyp);

      off = wxp * p(i + 1, j) + wxm * p(i - 1, j) + wyp * p(i, j + 1) + wym * p(i, j - 1);

      diag = wxm + wxp + wym + wyp;

    } else if (

        hc) {  // cut-cell stencil (Shortley-Weller); conductor neighbor = phi=0 on the circle

      off = cf(i, j, 1) * p(i + 1, j) + cf(i, j, 0) * p(i - 1, j) + cf(i, j, 3) * p(i, j + 1) +

            cf(i, j, 2) * p(i, j - 1);

      diag = cf(i, j, 4);

    } else {

      off = (p(i + 1, j) + p(i - 1, j)) * idx2 + (p(i, j + 1) + p(i, j - 1)) * idy2;

      diag = diag0;

    }

    // Reaction term: the operator becomes div(eps grad phi) - kappa phi, so the

    // diagonal gains +kappa (kappa >= 0 => more diagonally dominant, MG converges better).

    p(i, j) = (off - ff(i, j)) / (diag + (hk ? ka(i, j) : Real(0)));

  }


};


inline void gs_color(MultiFab& phi, const MultiFab& f, const Geometry& geom, int color,

                     const MultiFab* mask, const MultiFab* coef, const MultiFab* eps,

                     const MultiFab* kappa = nullptr, const MultiFab* eps_y = nullptr) {

  const Real idx2 = Real(1) / (geom.dx() * geom.dx());

  const Real idy2 = Real(1) / (geom.dy() * geom.dy());

  const Real diag0 = 2 * idx2 + 2 * idy2;

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

    Array4 p = phi.fab(li).array();

    const ConstArray4 ff = f.fab(li).const_array();

    const Box2D v = phi.box(li);

    const bool hm = mask != nullptr;

    const ConstArray4 mk = hm ? mask->fab(li).const_array() : ConstArray4{};

    const bool hc = coef != nullptr;

    const ConstArray4 cf = hc ? coef->fab(li).const_array() : ConstArray4{};

    const bool he = eps != nullptr;

    const ConstArray4 ep = he ? eps->fab(li).const_array() : ConstArray4{};

    // eps_y==nullptr => isotropic: y faces read the same field as the x faces (eps_x).

    const ConstArray4 ey = (he && eps_y) ? eps_y->fab(li).const_array() : ep;

    const bool hk = kappa != nullptr;  // reaction term -kappa phi (Helmholtz / screened)

    const ConstArray4 ka = hk ? kappa->fab(li).const_array() : ConstArray4{};

    for_each_cell(

        v, GsColorKernel{p, ff, idx2, idy2, diag0, color, hm, mk, hc, cf, he, ep, ey, hk, ka});

  }

}


}  // namespace detail


inline void gs_rb_sweep(MultiFab& phi, const MultiFab& f, const Geometry& geom, const BCRec& bc,

                        const MultiFab* mask = nullptr, const MultiFab* coef = nullptr,

                        const MultiFab* eps = nullptr, const MultiFab* kappa = nullptr,

                        const MultiFab* eps_y = nullptr) {

  device_fence();  // wait for the previous kernel before the host read of the halos

  fill_ghosts(phi, geom.domain, bc);

  detail::gs_color(phi, f, geom, 0, mask, coef, eps, kappa, eps_y);  // red (GPU kernel)

  device_fence();  // the black sweep reads the red values via host fill_ghosts

  fill_ghosts(phi, geom.domain, bc);

  detail::gs_color(phi, f, geom, 1, mask, coef, eps, kappa, eps_y);  // black

}


inline void gs_smooth(MultiFab& phi, const MultiFab& f, const Geometry& geom, const BCRec& bc,

                      int nsweeps, const MultiFab* mask = nullptr, const MultiFab* coef = nullptr,

                      const MultiFab* eps = nullptr, const MultiFab* kappa = nullptr,

                      const MultiFab* eps_y = nullptr) {

  for (int s = 0; s < nsweeps; ++s)

    gs_rb_sweep(phi, f, geom, bc, mask, coef, eps, kappa, eps_y);

}


namespace detail {

// Pins phi=0 in the conductor cells (mask==0). Named functor (#93); body identical.


struct ZeroConductorKernel {

  Array4 p;

  ConstArray4 mk;


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

    if (mk(i, j) == Real(0))

      p(i, j) = 0;

  }


};


}  // namespace detail


// Forces phi=0 in the conductor cells (mask==0).


inline void zero_conductor(MultiFab& phi, const MultiFab& mask) {

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

    Array4 p = phi.fab(li).array();

    const ConstArray4 mk = mask.fab(li).const_array();

    const Box2D v = phi.box(li);

    for_each_cell(v, detail::ZeroConductorKernel{p, mk});

  }

}


}  // namespace pops

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

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

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

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

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

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

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

fab2d.hpp
Fab2D: single-grid data on a Box2D (in-house equivalent of AMReX's FArrayBox); Array4 / ConstArray4: ...

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

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

pops::detail::gs_color
void gs_color(MultiFab &phi, const MultiFab &f, const Geometry &geom, int color, const MultiFab *mask, const MultiFab *coef, const MultiFab *eps, const MultiFab *kappa=nullptr, const MultiFab *eps_y=nullptr)
Definition poisson_operator.hpp:351

pops::detail::face_weights
POPS_HD void face_weights(const ConstArray4 &ep, const ConstArray4 &ey, int i, int j, Real idx2, Real idy2, bool hc, const ConstArray4 &cf, Real &wxm, Real &wxp, Real &wym, Real &wyp)
Definition poisson_operator.hpp:52

pops::detail::cross_div
POPS_HD Real cross_div(const ConstArray4 &p, bool hxy, const ConstArray4 &axy, bool hyx, const ConstArray4 &ayx, int i, int j, Real idx, Real idy)
Definition poisson_operator.hpp:79

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::zero_conductor
void zero_conductor(MultiFab &phi, const MultiFab &mask)
Definition poisson_operator.hpp:410

pops::apply_laplacian
void apply_laplacian(const MultiFab &phi, const Geometry &geom, MultiFab &lap, const MultiFab *coef=nullptr, const MultiFab *eps=nullptr, const MultiFab *kappa=nullptr, const MultiFab *eps_y=nullptr, const MultiFab *a_xy=nullptr, const MultiFab *a_yx=nullptr)
Definition poisson_operator.hpp:191

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::gs_smooth
void gs_smooth(MultiFab &phi, const MultiFab &f, const Geometry &geom, const BCRec &bc, int nsweeps, const MultiFab *mask=nullptr, const MultiFab *coef=nullptr, const MultiFab *eps=nullptr, const MultiFab *kappa=nullptr, const MultiFab *eps_y=nullptr)
Definition poisson_operator.hpp:389

pops::gs_rb_sweep
void gs_rb_sweep(MultiFab &phi, const MultiFab &f, const Geometry &geom, const BCRec &bc, const MultiFab *mask=nullptr, const MultiFab *coef=nullptr, const MultiFab *eps=nullptr, const MultiFab *kappa=nullptr, const MultiFab *eps_y=nullptr)
Definition poisson_operator.hpp:377

pops::poisson_residual
void poisson_residual(MultiFab &phi, const MultiFab &f, const Geometry &geom, const BCRec &bc, MultiFab &res, const MultiFab *mask=nullptr, const MultiFab *coef=nullptr, const MultiFab *eps=nullptr, const MultiFab *kappa=nullptr, const MultiFab *eps_y=nullptr, const MultiFab *a_xy=nullptr, const MultiFab *a_yx=nullptr)
Definition poisson_operator.hpp:271

pops::eps_harmonic
POPS_HD Real eps_harmonic(Real ec, Real ev)
Definition poisson_operator.hpp:40

pops::apply_divergence
void apply_divergence(const MultiFab &fx, const MultiFab &fy, const Geometry &geom, MultiFab &div_out, int cx=0, int cy=0)
Definition poisson_operator.hpp:256

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

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::Box2D
2D integer index space, cell-centered.
Definition box2d.hpp:37

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

pops::Geometry
Cartesian geometry of a level: index domain + physical bounds [xlo, xhi] x [ylo, yhi].
Definition geometry.hpp:20

pops::Geometry::dy
POPS_HD Real dy() const
Grid spacing in y (= (yhi - ylo) / domain.ny()). POPS_HD.
Definition geometry.hpp:33

pops::Geometry::dx
POPS_HD Real dx() const
Grid spacing in x (= (xhi - xlo) / domain.nx()). POPS_HD.
Definition geometry.hpp:31

pops::Geometry::domain
Box2D domain
Definition geometry.hpp:21

pops::detail::ApplyLaplacianKernel
Definition poisson_operator.hpp:113

pops::detail::ApplyLaplacianKernel::cf
ConstArray4 cf
Definition poisson_operator.hpp:118

pops::detail::ApplyLaplacianKernel::ep
ConstArray4 ep
Definition poisson_operator.hpp:120

pops::detail::ApplyLaplacianKernel::axy
ConstArray4 axy
Definition poisson_operator.hpp:124

pops::detail::ApplyLaplacianKernel::he
bool he
Definition poisson_operator.hpp:119

pops::detail::ApplyLaplacianKernel::hyx
bool hyx
Definition poisson_operator.hpp:123

pops::detail::ApplyLaplacianKernel::idy2
Real idy2
Definition poisson_operator.hpp:116

pops::detail::ApplyLaplacianKernel::hk
bool hk
Definition poisson_operator.hpp:121

pops::detail::ApplyLaplacianKernel::idx
Real idx
Definition poisson_operator.hpp:116

pops::detail::ApplyLaplacianKernel::idx2
Real idx2
Definition poisson_operator.hpp:116

pops::detail::ApplyLaplacianKernel::L
Array4 L
Definition poisson_operator.hpp:115

pops::detail::ApplyLaplacianKernel::idy
Real idy
Definition poisson_operator.hpp:116

pops::detail::ApplyLaplacianKernel::hc
bool hc
Definition poisson_operator.hpp:117

pops::detail::ApplyLaplacianKernel::ayx
ConstArray4 ayx
Definition poisson_operator.hpp:124

pops::detail::ApplyLaplacianKernel::c
int c
Definition poisson_operator.hpp:125

pops::detail::ApplyLaplacianKernel::ey
ConstArray4 ey
Definition poisson_operator.hpp:120

pops::detail::ApplyLaplacianKernel::p
ConstArray4 p
Definition poisson_operator.hpp:114

pops::detail::ApplyLaplacianKernel::hxy
bool hxy
Definition poisson_operator.hpp:123

pops::detail::ApplyLaplacianKernel::ka
ConstArray4 ka
Definition poisson_operator.hpp:122

pops::detail::ApplyLaplacianKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition poisson_operator.hpp:126

pops::detail::DivergenceKernel
Definition poisson_operator.hpp:234

pops::detail::DivergenceKernel::cx
int cx
Definition poisson_operator.hpp:238

pops::detail::DivergenceKernel::fy
ConstArray4 fy
Definition poisson_operator.hpp:235

pops::detail::DivergenceKernel::fx
ConstArray4 fx
Definition poisson_operator.hpp:235

pops::detail::DivergenceKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition poisson_operator.hpp:239

pops::detail::DivergenceKernel::cy
int cy
Definition poisson_operator.hpp:238

pops::detail::DivergenceKernel::div
Array4 div
Definition poisson_operator.hpp:236

pops::detail::DivergenceKernel::half_idy
Real half_idy
Definition poisson_operator.hpp:237

pops::detail::DivergenceKernel::half_idx
Real half_idx
Definition poisson_operator.hpp:237

pops::detail::GsColorKernel
Definition poisson_operator.hpp:312

pops::detail::GsColorKernel::p
Array4 p
Definition poisson_operator.hpp:313

pops::detail::GsColorKernel::cf
ConstArray4 cf
Definition poisson_operator.hpp:320

pops::detail::GsColorKernel::hm
bool hm
Definition poisson_operator.hpp:317

pops::detail::GsColorKernel::idy2
Real idy2
Definition poisson_operator.hpp:315

pops::detail::GsColorKernel::ey
ConstArray4 ey
Definition poisson_operator.hpp:322

pops::detail::GsColorKernel::idx2
Real idx2
Definition poisson_operator.hpp:315

pops::detail::GsColorKernel::mk
ConstArray4 mk
Definition poisson_operator.hpp:318

pops::detail::GsColorKernel::hk
bool hk
Definition poisson_operator.hpp:323

pops::detail::GsColorKernel::hc
bool hc
Definition poisson_operator.hpp:319

pops::detail::GsColorKernel::ep
ConstArray4 ep
Definition poisson_operator.hpp:322

pops::detail::GsColorKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition poisson_operator.hpp:325

pops::detail::GsColorKernel::diag0
Real diag0
Definition poisson_operator.hpp:315

pops::detail::GsColorKernel::ka
ConstArray4 ka
Definition poisson_operator.hpp:324

pops::detail::GsColorKernel::color
int color
Definition poisson_operator.hpp:316

pops::detail::GsColorKernel::he
bool he
Definition poisson_operator.hpp:321

pops::detail::GsColorKernel::ff
ConstArray4 ff
Definition poisson_operator.hpp:314

pops::detail::PoissonResidualKernel
Definition poisson_operator.hpp:149

pops::detail::PoissonResidualKernel::idx2
Real idx2
Definition poisson_operator.hpp:152

pops::detail::PoissonResidualKernel::idy2
Real idy2
Definition poisson_operator.hpp:152

pops::detail::PoissonResidualKernel::ep
ConstArray4 ep
Definition poisson_operator.hpp:158

pops::detail::PoissonResidualKernel::r
Array4 r
Definition poisson_operator.hpp:151

pops::detail::PoissonResidualKernel::ka
ConstArray4 ka
Definition poisson_operator.hpp:160

pops::detail::PoissonResidualKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition poisson_operator.hpp:163

pops::detail::PoissonResidualKernel::ey
ConstArray4 ey
Definition poisson_operator.hpp:158

pops::detail::PoissonResidualKernel::cf
ConstArray4 cf
Definition poisson_operator.hpp:156

pops::detail::PoissonResidualKernel::idx
Real idx
Definition poisson_operator.hpp:152

pops::detail::PoissonResidualKernel::mk
ConstArray4 mk
Definition poisson_operator.hpp:154

pops::detail::PoissonResidualKernel::ff
ConstArray4 ff
Definition poisson_operator.hpp:150

pops::detail::PoissonResidualKernel::axy
ConstArray4 axy
Definition poisson_operator.hpp:162

pops::detail::PoissonResidualKernel::idy
Real idy
Definition poisson_operator.hpp:152

pops::detail::PoissonResidualKernel::hk
bool hk
Definition poisson_operator.hpp:159

pops::detail::PoissonResidualKernel::p
ConstArray4 p
Definition poisson_operator.hpp:150

pops::detail::PoissonResidualKernel::hc
bool hc
Definition poisson_operator.hpp:155

pops::detail::PoissonResidualKernel::hyx
bool hyx
Definition poisson_operator.hpp:161

pops::detail::PoissonResidualKernel::ayx
ConstArray4 ayx
Definition poisson_operator.hpp:162

pops::detail::PoissonResidualKernel::he
bool he
Definition poisson_operator.hpp:157

pops::detail::PoissonResidualKernel::hxy
bool hxy
Definition poisson_operator.hpp:161

pops::detail::PoissonResidualKernel::hm
bool hm
Definition poisson_operator.hpp:153

pops::detail::ZeroConductorKernel
Definition poisson_operator.hpp:399

pops::detail::ZeroConductorKernel::operator()
POPS_HD void operator()(int i, int j) const
Definition poisson_operator.hpp:402

pops::detail::ZeroConductorKernel::p
Array4 p
Definition poisson_operator.hpp:400

pops::detail::ZeroConductorKernel::mk
ConstArray4 mk
Definition poisson_operator.hpp:401

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

POPS_HD
#define POPS_HD
Definition types.hpp:25