PoPS/doxygen/composite__fac__poisson_8hpp_source.html

#pragma once


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

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

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

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

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

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

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

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

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

#include <pops/numerics/elliptic/mg/geometric_mg.hpp>  // coarse solver (geometric multigrid)

#include <pops/numerics/elliptic/poisson/poisson_operator.hpp>  // apply_laplacian (residual, reads the already-filled ghosts)

#include <pops/numerics/time/amr/levels/amr_patch_range.hpp>  // PatchRange, CoverageMask (coarse footprint of a patch)


#include <algorithm>

#include <cmath>

#include <cstdio>

#include <stdexcept>

#include <vector>


namespace pops {


namespace detail {


POPS_HD inline Real fac_bilerp_coarse(const ConstArray4& C, int i, int j, int r) {

  const Real fx = (Real(i) + Real(0.5)) / Real(r) - Real(0.5);

  const Real fy = (Real(j) + Real(0.5)) / Real(r) - Real(0.5);

  const int Ic = static_cast<int>(std::floor(fx));

  const int Jc = static_cast<int>(std::floor(fy));

  const Real tx = fx - Real(Ic), ty = fy - Real(Jc);

  const Real c00 = C(Ic, Jc, 0), c10 = C(Ic + 1, Jc, 0);

  const Real c01 = C(Ic, Jc + 1, 0), c11 = C(Ic + 1, Jc + 1, 0);

  return (Real(1) - tx) * (Real(1) - ty) * c00 + tx * (Real(1) - ty) * c10 +

         (Real(1) - tx) * ty * c01 + tx * ty * c11;

}


}  // namespace detail


class CompositeFacPoisson {

 public:


  CompositeFacPoisson(const Geometry& geom_c, const BoxArray& ba_c, const BCRec& bc,

                      const Box2D& fine_box, int ratio = 2)

      : CompositeFacPoisson(geom_c, ba_c, bc, BoxArray(std::vector<Box2D>{fine_box}), ratio) {}


  CompositeFacPoisson(const Geometry& geom_c, const BoxArray& ba_c, const BCRec& bc,

                      const BoxArray& fine_boxes, int ratio = 2)

      : geom_c_(geom_c),

        geom_f_(geom_c.refine(ratio)),

        ba_c_(ba_c),

        dm_c_(ba_c.size(), n_ranks()),

        bc_(bc),

        ratio_(ratio),

        ba_f_(fine_boxes),

        dm_f_(fine_boxes.size(), n_ranks()),

        mg_(geom_c, ba_c, bc, {}, /*replicated=*/true),

        phi_c_(ba_c, dm_c_, 1, 1),

        phi_f_(ba_f_, dm_f_, 1, 1),

        f_c_(ba_c, dm_c_, 1, 0),

        f_f_(ba_f_, dm_f_, 1, 0),

        res_c_(ba_c, dm_c_, 1, 0),

        eps_c_(ba_c, dm_c_, 1, 1),

        eps_f_(ba_f_, dm_f_, 1, 1),

        axy_c_(ba_c, dm_c_, 1, 1),

        ayx_c_(ba_c, dm_c_, 1, 1),

        axy_f_(ba_f_, dm_f_, 1, 1),

        ayx_f_(ba_f_, dm_f_, 1, 1),

        cov_(Box2D::from_extents(geom_c.domain.nx(), geom_c.domain.ny())) {

    require_separated_patches(

        fine_boxes);  // guard: NON adjacent patches (fine-fine join = Phase 4b)

    // coarse footprints (covered cells) PER PATCH: PatchRange (lo/2 .. (hi-1)/2). The global coarse

    // coverage = UNION of the footprints (any gap between disjoint patches stays NON covered).

    for (int g = 0; g < fine_boxes.size(); ++g)

      patch_coarse_.push_back(PatchRange(fine_boxes[g]).box());

    for (const Box2D& pc : patch_coarse_)

      cov_.mark(pc);

    phi_c_.set_val(Real(0));

    phi_f_.set_val(Real(0));

    eps_c_.set_val(Real(1));  // default permittivity 1 -> operator = Laplacian (scalar)

    eps_f_.set_val(Real(1));

    axy_c_.set_val(Real(0));  // default cross terms 0 -> diagonal block only

    ayx_c_.set_val(Real(0));

    axy_f_.set_val(Real(0));

    ayx_f_.set_val(Real(0));

  }


  MultiFab& rhs_coarse() {

    return f_c_;

  }


  MultiFab& rhs_fine() { return f_f_; }

  MultiFab& phi_coarse() { return phi_c_; }

  MultiFab& phi_fine() { return phi_f_; }

  MultiFab& eps_coarse() { return eps_c_; }

  MultiFab& eps_fine() { return eps_f_; }

  void use_variable_coefficient(bool v) { has_eps_ = v; }

  MultiFab& a_xy_coarse() { return axy_c_; }

  MultiFab& a_yx_coarse() { return ayx_c_; }

  MultiFab& a_xy_fine() { return axy_f_; }

  MultiFab& a_yx_fine() { return ayx_f_; }

  void use_cross_terms(bool v) { has_cross_ = v; }

  const Box2D& patch_coarse() const { return patch_coarse_[0]; }

  const Box2D& patch_coarse(int g) const { return patch_coarse_[g]; }

  int n_fine_patches() const { return ba_f_.size(); }


  void set_verbose(bool v) { verbose_ = v; }

  void set_two_way(bool v) { two_way_ = v; }


  Real solve(int max_iters = 30, int fine_sweeps = 400, Real tol = 1e-9) {

    // VARIABLE COEFFICIENT (condensed Schur operator B_z=0): sets eps on the coarse solver and

    // fills the eps ghosts PER LEVEL. eps_c ghosts = zero-gradient (coeff_bc Foextrap, like the

    // Schur builder); eps_f C-F ghosts = bilerp of eps_c (consistency of the coefficient flux across

    // the interface). Without variable coefficient -> scalar Laplacian operator (Phase 1, bit-identical).

    if (has_eps_) {

      device_fence();

      fill_ghosts(eps_c_, geom_c_.domain, coeff_bc(bc_));

      fill_cf_coarse_to_fine(eps_c_, eps_f_);

      mg_.set_epsilon(eps_c_);

    }

    if (has_cross_) {  // FULL tensor (Schur B_z != 0): cross terms on both solvers + ghosts.

      device_fence();

      fill_ghosts(axy_c_, geom_c_.domain, coeff_bc(bc_));

      fill_ghosts(ayx_c_, geom_c_.domain, coeff_bc(bc_));

      fill_cf_coarse_to_fine(axy_c_, axy_f_);

      fill_cf_coarse_to_fine(ayx_c_, ayx_f_);

      mg_.set_cross_terms(axy_c_, ayx_c_);

    }

    // 0) initial coarse solve (gives a phi_c for the 1st C-F ghost).

    copy0(mg_.rhs(), f_c_);

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

    mg_.solve(Real(1e-12), 100);

    copy0(phi_c_, mg_.phi());


    // 1) bilinear C-F ghosts + fine solve (base ONE-WAY).

    refresh_fine(fine_sweeps);


    Real rnorm = composite_coarse_residual();

    if (verbose_ && my_rank() == 0)

      std::fprintf(stderr, "[FAC] init r_c=%.4e\n", rnorm);

    if (!two_way_) {

      last_residual_ = rnorm;

      return rnorm;

    }


    // 2) FAC two-way iterations: coarse correction (C-F flux) then re-solve fine.

    for (int it = 0; it < max_iters; ++it) {

      if (rnorm < tol)

        break;

      // coarse correction: Lap e_c = r_c (homogeneous Dirichlet), phi_c += e_c (non covered).

      copy0(mg_.rhs(), res_c_);

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

      mg_.solve(Real(1e-12), 100);

      add_uncovered(phi_c_, mg_.phi());

      // re-ghost + re-solve fine on the corrected phi_c.

      refresh_fine(fine_sweeps);

      rnorm = composite_coarse_residual();

      if (verbose_ && my_rank() == 0)

        std::fprintf(stderr, "[FAC] it=%d r_c=%.4e\n", it, rnorm);

    }

    last_residual_ = rnorm;

    return rnorm;

  }


  Real last_residual() const { return last_residual_; }


 private:

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

    device_fence();

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

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

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

      const Box2D b = dst.box(li);

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

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

          d(i, j, 0) = s(i, j, 0);

    }

  }


  void add_uncovered(MultiFab& phi, const MultiFab& e) {

    device_fence();

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

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

      const ConstArray4 ec = e.fab(li).const_array();

      const Box2D b = phi.box(li);

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

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

          if (!cov_.covered(i, j))

            p(i, j, 0) += ec(i, j, 0);

    }

  }


  static void require_separated_patches(const BoxArray& fine_boxes) {

    const int N = fine_boxes.size();

    for (int g = 0; g < N; ++g) {

      const Box2D ag = PatchRange(fine_boxes[g]).box();

      for (int h = g + 1; h < N; ++h) {

        const Box2D bh = PatchRange(fine_boxes[h]).box();

        if (!ag.grow(1).intersect(bh).empty())

          throw std::runtime_error(

              "CompositeFacPoisson: adjacent or overlapping fine patches (coarse footprints "

              "separated by less than one cell); the multi-patch fine-fine join is Phase 4b -- "

              "require disjoint patches separated by at least one coarse cell.");

      }

    }

  }


  void fill_cf_ghosts() {

    const ConstArray4 C = phi_c_.fab(0).const_array();  // replicated mono-box coarse

    const int ng = phi_f_.n_grow();

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

      Array4 F = phi_f_.fab(li).array();

      const Box2D vb = phi_f_.box(li);

      for (int j = vb.lo[1] - ng; j <= vb.hi[1] + ng; ++j)

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

          const bool inside = (i >= vb.lo[0] && i <= vb.hi[0] && j >= vb.lo[1] && j <= vb.hi[1]);

          if (inside)

            continue;  // ghosts only

          F(i, j, 0) = detail::fac_bilerp_coarse(C, i, j, ratio_);

        }

    }

  }


  void fill_cf_coarse_to_fine(const MultiFab& coarse, MultiFab& fine) {

    const ConstArray4 C = coarse.fab(0).const_array();  // replicated mono-box coarse

    const int ng = fine.n_grow();

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

      Array4 F = fine.fab(li).array();

      const Box2D vb = fine.box(li);

      for (int j = vb.lo[1] - ng; j <= vb.hi[1] + ng; ++j)

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

          const bool inside = (i >= vb.lo[0] && i <= vb.hi[0] && j >= vb.lo[1] && j <= vb.hi[1]);

          if (inside)

            continue;

          F(i, j, 0) = detail::fac_bilerp_coarse(C, i, j, ratio_);

        }

    }

  }


  static BCRec coeff_bc(const BCRec& b) {

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

    BCRec c;

    c.xlo = fo(b.xlo);

    c.xhi = fo(b.xhi);

    c.ylo = fo(b.ylo);

    c.yhi = fo(b.yhi);

    return c;

  }


  Real sor_omega(const Box2D& b) const {

    const int N = std::max(b.nx(), b.ny());

    return Real(2) / (Real(1) + std::sin(Real(kPi_) / Real(N)));

  }


  void refresh_fine(int sweeps) {

    device_fence();

    fill_ghosts(phi_c_, geom_c_.domain,

                bc_);  // phi_c physical ghosts (the bilerp reads up to the border)

    fill_cf_ghosts();

    fine_sor(sweeps);

    average_down(phi_f_, phi_c_,

                 ratio_);  // consistency: coarse covered = fine average (multi-box OK)

  }


  void fine_sor(int sweeps) {

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

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

    const bool he = has_eps_;

    const bool hc = has_cross_;

    const Real idx = Real(1) / geom_f_.dx(), idy = Real(1) / geom_f_.dy();  // cross_div: 1/dx, 1/dy

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

      const Box2D vb = phi_f_.box(li);

      const Real omega = sor_omega(vb);

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

      const ConstArray4 Pc =

          phi_f_.fab(li).const_array();  // const view (same memory) for cross stencil

      const ConstArray4 F = f_f_.fab(li).const_array();

      const ConstArray4 E = eps_f_.fab(li).const_array();

      const ConstArray4 AXY = axy_f_.fab(li).const_array();

      const ConstArray4 AYX = ayx_f_.fab(li).const_array();

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

        for (int color = 0; color < 2; ++color)

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

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

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

                continue;

              // FACE permittivities (harmonic mean of the 2 centers); eps==1 -> faces == 1.

              const Real exm = he ? eps_harmonic(E(i, j, 0), E(i - 1, j, 0)) : Real(1);

              const Real exp = he ? eps_harmonic(E(i, j, 0), E(i + 1, j, 0)) : Real(1);

              const Real eym = he ? eps_harmonic(E(i, j, 0), E(i, j - 1, 0)) : Real(1);

              const Real eyp = he ? eps_harmonic(E(i, j, 0), E(i, j + 1, 0)) : Real(1);

              const Real diag = (exm + exp) * idx2 + (eym + eyp) * idy2;

              const Real nb = (exm * P(i - 1, j, 0) + exp * P(i + 1, j, 0)) * idx2 +

                              (eym * P(i, j - 1, 0) + eyp * P(i, j + 1, 0)) * idy2;

              // EXPLICIT cross terms (9 points, read from the current P): div(A grad phi) =

              // diag_block + cross. We solve diag_block(P) + cross(P) = f -> P = (nb + cross - f)/diag.

              const Real cross =

                  hc ? detail::cross_div(Pc, true, AXY, true, AYX, i, j, idx, idy) : Real(0);

              const Real pgs = (nb + cross - F(i, j, 0)) / diag;

              P(i, j, 0) = (Real(1) - omega) * P(i, j, 0) +

                           omega * pgs;  // over-relax (under-relax if strong)

            }

    }

  }


  Real composite_coarse_residual() {

    device_fence();

    fill_ghosts(phi_c_, geom_c_.domain, bc_);

    // r_c = f_c - div(A grad phi_c) (apply_laplacian reads the already-filled ghosts; eps + cross if active).

    // The cross terms are read also on the COVERED cells (= fine average after average_down) -> the

    // 9-point stencil stays consistent at the interface; only the NORMAL flux is explicitly joined C-F

    // (the cross flux, tangential and small for the Schur step, is carried by the volume stencil).

    MultiFab lap(ba_c_, dm_c_, 1, 0);

    apply_laplacian(phi_c_, geom_c_, lap, /*coef=*/nullptr, has_eps_ ? &eps_c_ : nullptr,

                    /*kappa=*/nullptr, /*eps_y=*/nullptr, has_cross_ ? &axy_c_ : nullptr,

                    has_cross_ ? &ayx_c_ : nullptr);

    device_fence();

    Array4 R = res_c_.fab(0).array();

    const ConstArray4 LAP = lap.fab(0).const_array();

    const ConstArray4 FC = f_c_.fab(0).const_array();

    const Box2D b = res_c_.box(0);

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

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

        R(i, j, 0) = cov_.covered(i, j) ? Real(0) : (FC(i, j, 0) - LAP(i, j, 0));


    // C-F FLUX CORRECTION, PER FINE PATCH. On each coarse cell BORDERING a patch (non covered,

    // covered neighbor), we REPLACE the contribution of the C-F face in div(eps grad phi_c) by the

    // FINE contribution (conservative sum of the r fine faces, harmonic face eps): r_c += (coarse

    // - fine). Since the patches are separated by at least one coarse cell, each border is a TRUE

    // coarse-fine join; the test !cov_.covered(I, J) defensively skips a bordering cell that would be

    // covered by ANOTHER patch (impossible under the guard, but robust: a covered bordering

    // cell is already interior to another patch, its residual stays 0). A cell SEPARATING two

    // patches (right border of one, left border of the other) gets TWO corrections, one per face: correct.

    const ConstArray4 PC = phi_c_.fab(0).const_array();

    const ConstArray4 EC = eps_c_.fab(0).const_array();

    const bool he = has_eps_;

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

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

    const int r = ratio_;

    for (int g = 0; g < phi_f_.local_size(); ++g) {

      const ConstArray4 PF = phi_f_.fab(g).const_array();

      const ConstArray4 EF = eps_f_.fab(g).const_array();

      const int Ic0 = patch_coarse_[g].lo[0], Ic1 = patch_coarse_[g].hi[0];

      const int Jc0 = patch_coarse_[g].lo[1], Jc1 = patch_coarse_[g].hi[1];

      // Faces NORMAL TO X: bordering columns I = Ic0-1 (covered +x face) and I = Ic1+1 (-x face).

      for (int J = Jc0; J <= Jc1; ++J) {

        if (!cov_.covered(Ic0 - 1, J)) {  // left: cell (Ic0-1, J), fine face at i = r*Ic0.

          const int I = Ic0 - 1;

          const Real efc = he ? eps_harmonic(EC(I, J, 0), EC(I + 1, J, 0)) : Real(1);

          const Real coarse_c = efc * (PC(I + 1, J, 0) - PC(I, J, 0)) * idx2;

          Real fine_sum = Real(0);

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

            const int jf = r * J + t;

            const Real eff =

                he ? eps_harmonic(EF(r * Ic0 - 1, jf, 0), EF(r * Ic0, jf, 0)) : Real(1);

            fine_sum += eff * (PF(r * Ic0, jf, 0) - PF(r * Ic0 - 1, jf, 0));  // interior - ghost

          }

          R(I, J, 0) += coarse_c - fine_sum * idx2;

        }

        if (!cov_.covered(Ic1 + 1, J)) {  // right: cell (Ic1+1, J), fine faces at i = r*Ic1+r.

          const int I = Ic1 + 1;

          const Real efc = he ? eps_harmonic(EC(I, J, 0), EC(I - 1, J, 0)) : Real(1);

          const Real coarse_c = efc * (PC(I - 1, J, 0) - PC(I, J, 0)) * idx2;

          Real fine_sum = Real(0);

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

            const int jf = r * J + t;

            const Real eff =

                he ? eps_harmonic(EF(r * Ic1 + r - 1, jf, 0), EF(r * Ic1 + r, jf, 0)) : Real(1);

            fine_sum += eff * (PF(r * Ic1 + r - 1, jf, 0) - PF(r * Ic1 + r, jf, 0));

          }

          R(I, J, 0) += coarse_c - fine_sum * idx2;

        }

      }

      // Faces NORMAL TO Y: bordering rows J = Jc0-1 (+y face) and J = Jc1+1 (-y face).

      for (int I = Ic0; I <= Ic1; ++I) {

        if (!cov_.covered(I, Jc0 - 1)) {

          const int J = Jc0 - 1;

          const Real efc = he ? eps_harmonic(EC(I, J, 0), EC(I, J + 1, 0)) : Real(1);

          const Real coarse_c = efc * (PC(I, J + 1, 0) - PC(I, J, 0)) * idy2;

          Real fine_sum = Real(0);

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

            const int iff = r * I + t;

            const Real eff =

                he ? eps_harmonic(EF(iff, r * Jc0 - 1, 0), EF(iff, r * Jc0, 0)) : Real(1);

            fine_sum += eff * (PF(iff, r * Jc0, 0) - PF(iff, r * Jc0 - 1, 0));

          }

          R(I, J, 0) += coarse_c - fine_sum * idy2;

        }

        if (!cov_.covered(I, Jc1 + 1)) {

          const int J = Jc1 + 1;

          const Real efc = he ? eps_harmonic(EC(I, J, 0), EC(I, J - 1, 0)) : Real(1);

          const Real coarse_c = efc * (PC(I, J - 1, 0) - PC(I, J, 0)) * idy2;

          Real fine_sum = Real(0);

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

            const int iff = r * I + t;

            const Real eff =

                he ? eps_harmonic(EF(iff, r * Jc1 + r - 1, 0), EF(iff, r * Jc1 + r, 0)) : Real(1);

            fine_sum += eff * (PF(iff, r * Jc1 + r - 1, 0) - PF(iff, r * Jc1 + r, 0));

          }

          R(I, J, 0) += coarse_c - fine_sum * idy2;

        }

      }

    }


    // inf norm of the residual over the NON covered cells.

    Real nrm = Real(0);

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

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

        if (!cov_.covered(i, j))

          nrm = std::fmax(nrm, std::fabs(R(i, j, 0)));

    return nrm;

  }


  Geometry geom_c_, geom_f_;

  BoxArray ba_c_;

  DistributionMapping dm_c_;

  BCRec bc_;

  int ratio_;

  BoxArray ba_f_;

  DistributionMapping dm_f_;

  GeometricMG mg_;

  MultiFab phi_c_, phi_f_, f_c_, f_f_, res_c_;

  MultiFab eps_c_, eps_f_;

  MultiFab axy_c_, ayx_c_, axy_f_, ayx_f_;

  std::vector<Box2D> patch_coarse_;

  CoverageMask cov_;

  Real last_residual_ = 0;

  bool has_eps_ = false;

  bool has_cross_ = false;

  bool verbose_ = false;

  bool two_way_ = true;

  static constexpr Real kPi_ = Real(3.14159265358979323846);

};


}  // namespace pops


amr_patch_range.hpp
Named types of the multi-patch coarse-fine interface: PatchRange (coarse footprint of a fine patch),...

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

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

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

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

pops::CompositeFacPoisson
2-level COMPOSITE FAC Poisson solver (scalar).
Definition composite_fac_poisson.hpp:87

pops::CompositeFacPoisson::a_yx_fine
MultiFab & a_yx_fine()
Definition composite_fac_poisson.hpp:163

pops::CompositeFacPoisson::use_variable_coefficient
void use_variable_coefficient(bool v)
Definition composite_fac_poisson.hpp:154

pops::CompositeFacPoisson::set_two_way
void set_two_way(bool v)
true: iterate the FAC two-way coupling (C-F flux correction + coarse correction).
Definition composite_fac_poisson.hpp:175

pops::CompositeFacPoisson::n_fine_patches
int n_fine_patches() const
Number of fine patches (size of the fine BoxArray).
Definition composite_fac_poisson.hpp:170

pops::CompositeFacPoisson::CompositeFacPoisson
CompositeFacPoisson(const Geometry &geom_c, const BoxArray &ba_c, const BCRec &bc, const Box2D &fine_box, int ratio=2)
MONO-PATCH CTOR (Phase 1): DELEGATES to the multi-patch ctor with a fine BoxArray of a single box,...
Definition composite_fac_poisson.hpp:95

pops::CompositeFacPoisson::rhs_fine
MultiFab & rhs_fine()
fine right-hand side f_f (div(eps grad phi_f) = f_f)
Definition composite_fac_poisson.hpp:146

pops::CompositeFacPoisson::set_verbose
void set_verbose(bool v)
Definition composite_fac_poisson.hpp:172

pops::CompositeFacPoisson::use_cross_terms
void use_cross_terms(bool v)
Definition composite_fac_poisson.hpp:164

pops::CompositeFacPoisson::a_xy_coarse
MultiFab & a_xy_coarse()
Cross terms a_xy / a_yx (at cell centers) PER LEVEL: FULL tensor A = diag(eps,eps) + [[0,...
Definition composite_fac_poisson.hpp:160

pops::CompositeFacPoisson::phi_fine
MultiFab & phi_fine()
Definition composite_fac_poisson.hpp:148

pops::CompositeFacPoisson::CompositeFacPoisson
CompositeFacPoisson(const Geometry &geom_c, const BoxArray &ba_c, const BCRec &bc, const BoxArray &fine_boxes, int ratio=2)
MULTI-PATCH CTOR (Phase 4a).
Definition composite_fac_poisson.hpp:102

pops::CompositeFacPoisson::a_yx_coarse
MultiFab & a_yx_coarse()
Definition composite_fac_poisson.hpp:161

pops::CompositeFacPoisson::solve
Real solve(int max_iters=30, int fine_sweeps=400, Real tol=1e-9)
Solves the composite system.
Definition composite_fac_poisson.hpp:179

pops::CompositeFacPoisson::phi_coarse
MultiFab & phi_coarse()
Definition composite_fac_poisson.hpp:147

pops::CompositeFacPoisson::rhs_coarse
MultiFab & rhs_coarse()
coarse right-hand side f_c (div(eps grad phi_c) = f_c)
Definition composite_fac_poisson.hpp:143

pops::CompositeFacPoisson::a_xy_fine
MultiFab & a_xy_fine()
Definition composite_fac_poisson.hpp:162

pops::CompositeFacPoisson::last_residual
Real last_residual() const
Definition composite_fac_poisson.hpp:234

pops::CompositeFacPoisson::patch_coarse
const Box2D & patch_coarse(int g) const
Coarse footprint of fine patch g (0 <= g < n_fine_patches()).
Definition composite_fac_poisson.hpp:168

pops::CompositeFacPoisson::patch_coarse
const Box2D & patch_coarse() const
Coarse footprint of the FIRST fine patch (mono-patch compat). Multi-patch: see patch_coarse(g).
Definition composite_fac_poisson.hpp:166

pops::CompositeFacPoisson::eps_fine
MultiFab & eps_fine()
Definition composite_fac_poisson.hpp:153

pops::CompositeFacPoisson::eps_coarse
MultiFab & eps_coarse()
VARIABLE permittivity eps (at cell centers) PER LEVEL.
Definition composite_fac_poisson.hpp:152

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

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

pops::GeometricMG::solve
int solve(Real rel_tol, int max_cycles, Real abs_tol=Real(0))
Definition geometric_mg.hpp:407

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

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

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

pops::GeometricMG::set_cross_terms
void set_cross_terms(std::function< Real(Real, Real)> a_xy_fn, std::function< Real(Real, Real)> a_yx_fn)
Definition geometric_mg.hpp:329

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::n_grow
int n_grow() const
Number of ghost layers.
Definition multifab.hpp:62

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

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

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

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

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

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

pops::detail::fac_bilerp_coarse
POPS_HD Real fac_bilerp_coarse(const ConstArray4 &C, int i, int j, int r)
BILINEAR interpolation of the coarse potential (cell-centered, C with ghosts) at the CENTER of the fi...
Definition composite_fac_poisson.hpp:70

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::average_down
void average_down(const MultiFab &fine, MultiFab &coarse, int r, MultiFab &cfine)
CONSERVATIVE average fine -> coarse (ratio r): coarse(I, J) = average of the r^2 fine cells of the bl...
Definition refinement.hpp:181

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

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

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::my_rank
int my_rank()
Definition comm.hpp:136

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

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

poisson_operator.hpp
Free functions of the elliptic operator: apply_laplacian (matvec), poisson_residual (residual),...

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

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::Box2D::from_extents
static Box2D from_extents(int nx, int ny)
Box [0, nx-1] x [0, ny-1] covering nx*ny cells from the index origin.
Definition box2d.hpp:42

pops::Box2D::hi
int hi[2]
Definition box2d.hpp:39

pops::Box2D::lo
int lo[2]
Definition box2d.hpp:38

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

pops::CoverageMask::covered
bool covered(int I, int J) const
Definition amr_patch_range.hpp:168

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::PatchRange
Definition amr_patch_range.hpp:42

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

POPS_HD
#define POPS_HD
Definition types.hpp:25