poisson.hpp

Go to the documentation of this file.

// core — cell-centered finite-volume Poisson on a BlockOctree, with a
// geometric multigrid built from the octree's own levels.
//
// This is the "grid + multigrid" use of the AMR octree (serial, host; the device
// and distributed paths build on the same operator later). The operator is a
// conservative two-point finite-volume Laplacian: for a leaf i,
//
//     (L u)_i = (1/V_i) * sum_faces  A_f / d_f * (u_j - u_i),
//
// with V_i the cell volume, A_f the shared face area, d_f the centre-to-centre
// normal distance. At a 2:1 interface the coarse side sums the flux over all
// 2^(Dim-1) fine neighbours (each with the fine face area), so the discretisation
// is conservative (the global integral of L u is zero on a periodic domain). The
// two-point gradient ignores the tangential centre offset, so the interface flux
// is first-order there — a documented limitation; a tangential-gradient
// correction (full 2nd-order AMR flux) is a follow-up.
//
// Multigrid: the hierarchy is the octree coarsened uniformly one level at a time
// (coarsenIf over all sibling groups). Restriction averages children -> parent;
// prolongation is piecewise-constant (correction scheme). Smoother is
// lexicographic Gauss-Seidel over the Z-order leaf slots. Periodic BCs only (the
// natural first target); the singular null space is fixed by mean removal.
//
#ifndef PECLET_AMR_POISSON_HPP
#define PECLET_AMR_POISSON_HPP
 
#include <array>
#include <cmath>
#include <functional>
#include <vector>
 
#include "peclet/amr/block_octree.hpp"
#include "peclet/amr/common.hpp"
#include "peclet/amr/face_csr.hpp"  // shared host+device FV (weight-CSR) row kernels
#include "peclet/amr/leaf_field.hpp"
#include "peclet/core/common/host_parallel.hpp"
#include "peclet/core/common/types.hpp"
 
namespace peclet::amr {
 
template <int Dim, unsigned Bits = (Dim == 2 ? 32u : (Dim == 3 ? 21u : 16u))>
class AmrPoisson {
 public:
  using Octree = BlockOctree<Dim, Bits>;
  using M = typename Octree::M;
  using Code = typename Octree::Code;
  using Coord = typename Octree::Coord;
 
  AmrPoisson() = default;
  AmrPoisson(const Octree& t, Real h0) { init(t, h0); }
  AmrPoisson(const Octree& t, const Vec<Dim>& h0) { init(t, h0); }
 
  void init(const Octree& t, Real h0) { init(t, detail::filledVec<Dim>(h0)); }
 
  void init(const Octree& t, const Vec<Dim>& h0) {
    t_ = &t;
    h0_ = h0;
    alpha_.clear();
    hasOpen_ = false;
    extResolve_ = nullptr;  // distributed seam re-installed by the owner after every init
    ghostLo_.clear();
    ghostLv_.clear();
    frameShift_ = {};
    for (int d = 0; d < Dim; ++d)
      fineExt_[d] = static_cast<Coord>(t.brick()[d] * (Index(1) << t.lmax()));
  }
 
  void setOrigin(const Vec<Dim>& o) { origin_ = o; }
 
  // ---- distributed seam (docs/amr_distributed_flow.md) ---------------------------------------
  using ExtResolver = std::function<Index(const std::array<long, Dim>&)>;
  void setResolver(ExtResolver r) { extResolve_ = std::move(r); }
  void setGhosts(std::vector<std::array<long, Dim>> lo, std::vector<unsigned> lv) {
    ghostLo_ = std::move(lo);
    ghostLv_ = std::move(lv);
  }
  Index numGhosts() const { return static_cast<Index>(ghostLv_.size()); }
  void setFrameShift(const std::array<long, Dim>& s) { frameShift_ = s; }
  unsigned levelOf(Index slot) const {
    return slot < numLeaves() ? t_->level(slot)
                              : ghostLv_[static_cast<std::size_t>(slot - numLeaves())];
  }
  std::array<long, Dim> loOf(Index slot) const {
    std::array<long, Dim> lo{};
    if (slot < numLeaves()) {
      auto b = t_->bounds(slot);
      for (int d = 0; d < Dim; ++d)
        lo[d] = static_cast<long>(b[0][d]);
    } else {
      lo = ghostLo_[static_cast<std::size_t>(slot - numLeaves())];
    }
    return lo;
  }
 
  std::pair<Index, unsigned> probeSlot(const std::array<long, Dim>& p) const {
    bool in = true;
    for (int d = 0; d < Dim; ++d)
      if (p[d] < 0 || p[d] >= static_cast<long>(fineExt_[d])) {
        in = false;
        break;
      }
    std::array<Coord, Dim> q{};
    if (!in) {
      if (extResolve_) {
        const Index s = extResolve_(p);
        return {s, s >= 0 ? levelOf(s) : 0u};
      }
      for (int d = 0; d < Dim; ++d) {
        const long e = static_cast<long>(fineExt_[d]);
        q[d] = static_cast<Coord>(((p[d] % e) + e) % e);
      }
    } else {
      for (int d = 0; d < Dim; ++d)
        q[d] = static_cast<Coord>(p[d]);
    }
    const Index j = t_->find(M::encode(q).code());
    return {j, j >= 0 ? t_->level(j) : 0u};
  }
 
  // ---- cut-cell openness (per-leaf per-face fluid fraction in [0,1]) -------
  static constexpr int kFaces = 2 * Dim;
  static int faceIndex(int axis, int dir) { return 2 * axis + (dir > 0 ? 0 : 1); }
 
  double faceOpenness(Index i, int axis, int dir) const {
    if (!hasOpen_)
      return 1.0;
    return alpha_[static_cast<std::size_t>(i) * kFaces + faceIndex(axis, dir)];
  }
 
  bool hasOpenness() const { return hasOpen_; }
  const std::vector<double>& opennessRaw() const { return alpha_; }
  void setOpennessRaw(std::vector<double> a) {
    alpha_ = std::move(a);
    hasOpen_ = true;
  }
 
  template <class OpenFn>
  void buildOpenness(OpenFn&& openFn) {
    const Index n = numLeaves();
    const Index ng = numGhosts();
    alpha_.assign(static_cast<std::size_t>(n + ng) * kFaces, 1.0);
    hasOpen_ = true;
    // Face centroids in the GLOBAL frame (localCoord + frameShift_): every rank evaluates the
    // world-coord openFn at bit-identical points, so ghost rows match the owner's exactly (the
    // symmetric-openFn contract). Single-rank: frameShift_ = 0, unchanged.
    auto fillRow = [&](Index row, const std::array<long, Dim>& lo, long s) {
      for (int axis = 0; axis < Dim; ++axis)
        for (int dir = -1; dir <= 1; dir += 2) {
          const long plane =
              (dir > 0) ? lo[axis] + frameShift_[axis] + s : lo[axis] + frameShift_[axis];
          Vec<Dim> fc{};
          for (int d = 0; d < Dim; ++d)
            fc[d] = (d == axis) ? origin_[d] + static_cast<Real>(plane) * h0_[d]
                                : origin_[d] + (static_cast<Real>(lo[d] + frameShift_[d]) +
                                                0.5 * static_cast<Real>(s)) *
                                                   h0_[d];
          double a = static_cast<double>(openFn(fc, axis));
          a = a < 0.0 ? 0.0 : (a > 1.0 ? 1.0 : a);
          alpha_[static_cast<std::size_t>(row) * kFaces + faceIndex(axis, dir)] = a;
        }
    };
    // Host-parallel (setup-parallel plan rung 2): leaf i writes only its own kFaces α slots and
    // `openFn` is a pure function of the face centroid — disjoint writes, so α is bitwise
    // identical under any thread count. Rank-local: the distributed contract is untouched.
    hostParFor(n, [&](Index i) {
      auto b = t_->bounds(i);
      std::array<long, Dim> lo{};
      for (int d = 0; d < Dim; ++d)
        lo[d] = static_cast<long>(b[0][d]);
      fillRow(i, lo, 1L << t_->level(i));
    });
    for (Index g = 0; g < ng; ++g)  // ghost α rows (distributed builds; empty single-rank)
      fillRow(n + g, ghostLo_[static_cast<std::size_t>(g)],
              1L << ghostLv_[static_cast<std::size_t>(g)]);
  }
 
  Index numLeaves() const { return t_->numLeaves(); }
 
  void setPeriodic(bool p) { periodic_ = p; }
  bool periodic() const { return periodic_; }
 
  void setImmersedWall(bool w) { immersedWall_ = w; }
  bool immersedWall() const { return immersedWall_; }
 
  double boundaryDiag(Index i) const {
    if (periodic_ && !immersedWall_)
      return 0.0;
    auto b = t_->bounds(i);
    const auto& lo = b[0];
    const Coord si = Coord(Coord(1) << t_->level(i));
    double s = 0.0;
    for (int axis = 0; axis < Dim; ++axis) {
      // Phase 3: A_f/d_f of a BOX cell — the face area normal to `axis` over the half width along
      // `axis`. Identical to the old single `wall` when the three spacings are equal.
      const double wall = areaOf(si, axis) / (0.5 * static_cast<Real>(si) * h0_[axis]);
      for (int dir = -1; dir <= 1; dir += 2) {
        const long pc = (dir > 0) ? static_cast<long>(lo[axis]) + static_cast<long>(si)
                                  : static_cast<long>(lo[axis]) - 1;
        const bool domainBoundary =
            !periodic_ && (pc < 0 || pc >= static_cast<long>(fineExt_[axis]));
        if (domainBoundary)
          s += faceOpenness(i, axis, dir) * wall;  // domain Dirichlet wall (open part)
        else if (immersedWall_)
          s += (1.0 - faceOpenness(i, axis, dir)) * wall;  // immersed no-slip wall (solid part)
      }
    }
    return s;
  }
 
  Real cellWidth(Index i, int axis) const {
    return h0_[axis] * static_cast<Real>(Index(1) << levelOf(i));
  }
  Real cellWidth(Index i) const { return h0_[0] * static_cast<Real>(Index(1) << levelOf(i)); }
  Real cellVolume(Index i) const {
    const Real f = static_cast<Real>(Index(1) << levelOf(i));
    Real v = 1;
    for (int d = 0; d < Dim; ++d)
      v *= h0_[d] * f;
    return v;
  }
 
  template <class Fn>
  void forEachFaceNeighbor(Index i, Fn&& fn) const {
    auto b = t_->bounds(i);
    const auto& lo = b[0];
    const unsigned Li = t_->level(i);
    const Coord si = Coord(Coord(1) << Li);
    for (int axis = 0; axis < Dim; ++axis)
      for (int dir = -1; dir <= 1; dir += 2) {
        const long pc = (dir > 0) ? static_cast<long>(lo[axis]) + static_cast<long>(si)
                                  : static_cast<long>(lo[axis]) - 1;
        // Non-periodic: a domain-boundary face has no neighbour cell (it is a Dirichlet
        // wall handled by boundaryDiag) — skip it.
        if (!periodic_ && (pc < 0 || pc >= static_cast<long>(fineExt_[axis])))
          continue;
        std::array<long, Dim> p{};
        for (int d = 0; d < Dim; ++d)
          p[d] = static_cast<long>(lo[d]);
        p[axis] = pc;
        const auto [j, Lj] = probeSlot(p);
        if (j < 0)
          continue;  // unresolved during the distributed discovery fixpoint only
        if (Lj >= Li) {
          // same level or coarser: one neighbour, shared face = this cell's face.
          // Openness lives on the finer side (here, this cell i).
          fn(j, coeff(si, Coord(Coord(1) << Lj), axis), axis, faceOpenness(i, axis, dir));
        } else {
          // finer neighbour: 2^(Dim-1) sub-faces, each the fine face area.
          const Coord sj = Coord(si >> 1);
          const int nsub = 1 << (Dim - 1);
          for (int k = 0; k < nsub; ++k) {
            std::array<long, Dim> q = p;
            int bit = 0;
            for (int t = 0; t < Dim; ++t) {
              if (t == axis)
                continue;
              q[t] = static_cast<long>(lo[t]) + (((k >> bit) & 1) ? static_cast<long>(sj) : 0L);
              ++bit;
            }
            const Index jj = probeSlot(q).first;
            if (jj < 0)
              continue;
            // Fine face area (min = sj) but the true centre-to-centre distance
            // (si+sj)/2 — same value the fine side computes, so the operator is
            // symmetric / conservative across the 2:1 interface. Openness lives on
            // the finer side (the neighbour jj), its face toward i is -dir.
            fn(jj, coeff(si, sj, axis), axis, faceOpenness(jj, axis, -dir));
          }
        }
      }
  }
 
  template <class Fn>
  void forEachFaceFull(Index i, Fn&& fn) const {
    auto b = t_->bounds(i);
    const auto& lo = b[0];
    const unsigned Li = t_->level(i);
    const Coord si = Coord(Coord(1) << Li);
    for (int axis = 0; axis < Dim; ++axis)
      for (int dir = -1; dir <= 1; dir += 2) {
        const long pc = (dir > 0) ? static_cast<long>(lo[axis]) + static_cast<long>(si)
                                  : static_cast<long>(lo[axis]) - 1;
        std::array<long, Dim> p{};
        for (int d = 0; d < Dim; ++d)
          p[d] = static_cast<long>(lo[d]);
        p[axis] = pc;
        const auto [j, Lj] = probeSlot(p);
        if (j < 0)
          continue;  // unresolved during the distributed discovery fixpoint only
        if (Lj >= Li) {
          const Coord sj = Coord(Coord(1) << Lj);
          fn(j, axis, dir, areaOf(si, axis),
             0.5 * (static_cast<Real>(si) + static_cast<Real>(sj)) * h0_[axis],
             faceOpenness(i, axis, dir));
        } else {
          const Coord sj = Coord(si >> 1);
          const int nsub = 1 << (Dim - 1);
          for (int k = 0; k < nsub; ++k) {
            std::array<long, Dim> q = p;
            int bit = 0;
            for (int t = 0; t < Dim; ++t) {
              if (t == axis)
                continue;
              q[t] = static_cast<long>(lo[t]) + (((k >> bit) & 1) ? static_cast<long>(sj) : 0L);
              ++bit;
            }
            const Index jj = probeSlot(q).first;
            if (jj < 0)
              continue;
            fn(jj, axis, dir, areaOf(sj, axis),
               0.5 * (static_cast<Real>(si) + static_cast<Real>(sj)) * h0_[axis],
               faceOpenness(jj, axis, -dir));
          }
        }
      }
  }
 
  Index periodicNeighbor(Index i, int axis, int dir) const {
    std::array<long, Dim> p = loOf(i);
    const long si = 1L << levelOf(i);
    p[axis] = (dir > 0) ? p[axis] + si : p[axis] - 1;
    return probeSlot(p).first;
  }
 
  double coarseStar(const std::vector<double>& u, Index coarse, Index fine, int axis) const {
    const double uc = u[static_cast<std::size_t>(coarse)];
    const std::array<long, Dim> bc = loOf(coarse);  // ghost-safe (block-local longs)
    const std::array<long, Dim> bf = loOf(fine);
    const double sc = static_cast<double>(Index(1) << levelOf(coarse));
    const double sf = static_cast<double>(Index(1) << levelOf(fine));
    double val = uc;
    for (int t = 0; t < Dim; ++t) {
      if (t == axis)
        continue;
      // Phase 3: the tangential offset AND the differencing width are the spacings of the
      // TANGENTIAL axis `t` (they were one `H` when the cells were cubes).
      const double H = cellWidth(coarse, t);
      const double dt =
          ((static_cast<double>(bf[t]) + 0.5 * sf) - (static_cast<double>(bc[t]) + 0.5 * sc)) *
          h0_[t];
      Index cp = periodicNeighbor(coarse, t, +1);
      Index cm = periodicNeighbor(coarse, t, -1);
      if (cp < 0 || cm < 0)
        continue;
      if (levelOf(cp) != levelOf(coarse) || levelOf(cm) != levelOf(coarse))
        continue;
      // Skip the correction near a solid: a nearly-closed tangential face means
      // the quadratic stencil would lean on a solid-side value. Drop to the raw
      // coarse value on this axis (locally lower order, but robust).
      if (faceOpenness(coarse, t, +1) < 0.5 || faceOpenness(coarse, t, -1) < 0.5)
        continue;
      const double up = u[static_cast<std::size_t>(cp)];
      const double um = u[static_cast<std::size_t>(cm)];
      const double Dt = (up - um) / (2.0 * H);
      const double Dtt = (up - 2.0 * uc + um) / (H * H);
      val += dt * Dt + 0.5 * dt * dt * Dtt;
    }
    return val;
  }
 
  void applyLaplacianQuad(const std::vector<double>& u, std::vector<double>& out) const {
    const Index n = numLeaves();
    out.assign(static_cast<std::size_t>(n), 0.0);
    for (Index i = 0; i < n; ++i) {
      const double ui = u[static_cast<std::size_t>(i)];
      const unsigned Li = t_->level(i);
      double acc = 0.0;
      forEachFaceNeighbor(i, [&](Index j, Real c, int axis, double a) {
        const unsigned Lj = levelOf(j);  // ghost-safe
        double uj = u[static_cast<std::size_t>(j)];
        double uii = ui;
        if (Lj > Li)
          uj = coarseStar(u, j, i, axis);  // j coarser: correct its value
        else if (Lj < Li)
          uii = coarseStar(u, i, j, axis);  // i coarser: correct our value for this sub-face
        acc += a * c * (uj - uii);
      });
      out[static_cast<std::size_t>(i)] = acc / cellVolume(i);
    }
  }
 
  double residualQuad(const std::vector<double>& u, const std::vector<double>& rhs,
                      std::vector<double>& res) const {
    std::vector<double> lu;
    applyLaplacianQuad(u, lu);
    const Index n = numLeaves();
    res.assign(static_cast<std::size_t>(n), 0.0);
    double s = 0.0;
    for (Index i = 0; i < n; ++i) {
      double r = rhs[static_cast<std::size_t>(i)] - lu[static_cast<std::size_t>(i)];
      res[static_cast<std::size_t>(i)] = r;
      s += cellVolume(i) * r * r;
    }
    return std::sqrt(s);
  }
 
  struct FvAssembled {
    std::vector<double> invVol;  
    std::vector<Index> start;    
    std::vector<Index> nbr;      
    std::vector<double> coef;    
    std::vector<double> bcDiag;  
  };
  FvAssembled assembleFv() const {
    const Index n = numLeaves();
    FvAssembled A;
    A.start.assign(static_cast<std::size_t>(n) + 1, 0);
    for (Index i = 0; i < n; ++i) {
      Index cnt = 0;
      forEachFaceNeighbor(i, [&](Index, Real, int, double) { ++cnt; });
      A.start[static_cast<std::size_t>(i) + 1] = A.start[static_cast<std::size_t>(i)] + cnt;
    }
    const Index nf = A.start[static_cast<std::size_t>(n)];
    A.nbr.resize(static_cast<std::size_t>(nf));
    A.coef.resize(static_cast<std::size_t>(nf));
    A.invVol.resize(static_cast<std::size_t>(n));
    A.bcDiag.resize(static_cast<std::size_t>(n));
    for (Index i = 0; i < n; ++i) {
      A.invVol[static_cast<std::size_t>(i)] = 1.0 / cellVolume(i);
      A.bcDiag[static_cast<std::size_t>(i)] = boundaryDiag(i);
      Index k = A.start[static_cast<std::size_t>(i)];
      forEachFaceNeighbor(i, [&](Index j, Real c, int, double a) {
        A.nbr[static_cast<std::size_t>(k)] = j;
        A.coef[static_cast<std::size_t>(k)] = a * c;
        ++k;
      });
    }
    return A;
  }
  FvCsrOpT<HostArr<double>, HostArr<Index>> hostFvOp(const FvAssembled& A) const {
    FvCsrOpT<HostArr<double>, HostArr<Index>> v;
    v.n = numLeaves();
    v.invVol = HostArr<double>(A.invVol.data());
    v.coef = HostArr<double>(A.coef.data());
    v.start = HostArr<Index>(A.start.data());
    v.nbr = HostArr<Index>(A.nbr.data());
    v.bcDiag = HostArr<double>(A.bcDiag.data());
    return v;
  }
  void applyFvShared(const std::vector<double>& u, std::vector<double>& out) const {
    const FvAssembled A = assembleFv();
    const auto op = hostFvOp(A);
    const Index n = numLeaves();
    out.assign(static_cast<std::size_t>(n), 0.0);
    const HostArr<double> uacc(u.data());
    for (Index i = 0; i < n; ++i)
      out[static_cast<std::size_t>(i)] = fvApplyRow(op, i, uacc);
  }
 
  void applyLaplacian(const std::vector<double>& u, std::vector<double>& out) const {
    const Index n = numLeaves();
    out.assign(static_cast<std::size_t>(n), 0.0);
    for (Index i = 0; i < n; ++i) {
      const double ui = u[static_cast<std::size_t>(i)];
      double acc = 0.0;
      forEachFaceNeighbor(i, [&](Index j, Real c, int, double a) {
        acc += a * c * (u[static_cast<std::size_t>(j)] - ui);
      });
      out[static_cast<std::size_t>(i)] = acc / cellVolume(i);
    }
  }
 
  double residual(const std::vector<double>& u, const std::vector<double>& rhs,
                  std::vector<double>& res) const {
    const Index n = numLeaves();
    res.assign(static_cast<std::size_t>(n), 0.0);
    double s = 0.0;
    for (Index i = 0; i < n; ++i) {
      const double ui = u[static_cast<std::size_t>(i)];
      double acc = 0.0;
      forEachFaceNeighbor(i, [&](Index j, Real c, int, double a) {
        acc += a * c * (u[static_cast<std::size_t>(j)] - ui);
      });
      double r = rhs[static_cast<std::size_t>(i)] - acc / cellVolume(i);
      res[static_cast<std::size_t>(i)] = r;
      s += cellVolume(i) * r * r;
    }
    return std::sqrt(s);
  }
 
  void gaussSeidel(std::vector<double>& u, const std::vector<double>& rhs, int sweeps) const {
    const Index n = numLeaves();
    for (int s = 0; s < sweeps; ++s)
      for (Index i = 0; i < n; ++i) {
        double sumOff = 0.0, diag = 0.0;
        forEachFaceNeighbor(i, [&](Index j, Real c, int, double a) {
          sumOff += a * c * u[static_cast<std::size_t>(j)];
          diag += a * c;
        });
        if (diag != 0.0)
          u[static_cast<std::size_t>(i)] =
              (sumOff - cellVolume(i) * rhs[static_cast<std::size_t>(i)]) / diag;
      }
  }
 
  void removeMean(std::vector<double>& u) const {
    const Index n = numLeaves();
    double sum = 0.0, vol = 0.0;
    for (Index i = 0; i < n; ++i) {
      sum += cellVolume(i) * u[static_cast<std::size_t>(i)];
      vol += cellVolume(i);
    }
    double m = sum / vol;
    for (Index i = 0; i < n; ++i)
      u[static_cast<std::size_t>(i)] -= m;
  }
 
  const Octree& octree() const { return *t_; }
  const Vec<Dim>& h0() const { return h0_; }
  Real h0Scalar() const { return h0_[0]; }
  const std::array<Coord, Dim>& fineExt() const { return fineExt_; }
  const Vec<Dim>& origin() const { return origin_; }
 
 private:
  Coord wrap(long c, int axis) const {
    long e = static_cast<long>(fineExt_[axis]);
    return static_cast<Coord>(((c % e) + e) % e);
  }
  // A_f / d_f (physical) for a cell of width-units `si` next to one of `sj`.
  Real coeff(Coord si, Coord sj, int axis) const {
    Real dist = 0.5 * (static_cast<Real>(si) + static_cast<Real>(sj)) * h0_[axis];
    return areaOf(si < sj ? si : sj, axis) / dist;
  }
  Real coeff(Coord si, Coord sj) const { return coeff(si, sj, 0); }
  // Physical area of a face NORMAL TO `axis` of a cell of width-units `s`: prod_{d != axis} s*h0_d.
  Real areaOf(Coord s, int axis) const {
    Real area = 1;
    for (int d = 0; d < Dim; ++d)
      if (d != axis)
        area *= static_cast<Real>(s) * h0_[d];
    return area;
  }
  Real areaOf(Coord s) const { return areaOf(s, 0); }
 
  const Octree* t_ = nullptr;
  Vec<Dim> h0_ = detail::filledVec<Dim>(1.0);
  std::array<Coord, Dim> fineExt_{};
  Vec<Dim> origin_{};
  ExtResolver extResolve_;                      // distributed seam: out-of-block probe resolver
  std::vector<std::array<long, Dim>> ghostLo_;  // ghost slot → block-local lo (longs)
  std::vector<unsigned> ghostLv_;               // ghost slot → covering-leaf level
  std::array<long, Dim> frameShift_{};          // block global fine origin (0 single-rank)
  std::vector<double> alpha_;  // per-(leaf+ghost) per-face openness (kFaces per row), or empty
  bool hasOpen_ = false;
  bool periodic_ = true;       // false ⇒ homogeneous Dirichlet domain walls (boundaryDiag)
  bool immersedWall_ = false;  // true ⇒ velocity operator: (1−α) interior faces are no-slip walls
};
 
template <int Dim, unsigned Bits = (Dim == 2 ? 32u : (Dim == 3 ? 21u : 16u))>
class AmrMultigrid {
 public:
  using Octree = BlockOctree<Dim, Bits>;
  using M = typename Octree::M;
  using Code = typename Octree::Code;
 
  static bool canLiftLevel(const Octree& t, Index bottomExtent) {
    const IVec<Dim>& G = t.brick();
    Index mx = 0;
    for (int d = 0; d < Dim; ++d)
      mx = std::max(mx, G[d]);
    if (mx <= bottomExtent)
      return false;
    for (int d = 0; d < Dim; ++d)
      if ((G[d] % 2) != 0 || (G[d] / 2) < 2)
        return false;
    return t.canLiftRoot();
  }
 
  void build(const Octree& finest, Real h0, bool liftRoot = true, Index bottomExtent = 4) {
    build(finest, detail::filledVec<Dim>(h0), liftRoot, bottomExtent);
  }
  void build(const Octree& finest, const Vec<Dim>& h0, bool liftRoot = true,
             Index bottomExtent = 4) {
    levels_.clear();
    levels_.push_back(finest);
    for (;;) {
      Octree c = levels_.back();
      Index merged = c.coarsenIf([](Code, unsigned) { return true; });
      if (merged == 0 || c.numLeaves() == levels_.back().numLeaves())
        break;
      levels_.push_back(c);
      if (c.numLeaves() == 1)
        break;
    }
    // Below the root brick: lift and merge while §6.2 allows it. Every leaf is a root cell here
    // (the octree loop ran to exhaustion), and an even brick with an even origin tiles complete
    // octets, so the coarsenIf after a lift merges ALL of them.
    if (liftRoot) {
      while (canLiftLevel(levels_.back(), bottomExtent)) {
        Octree c = levels_.back();
        c.liftRoot();
        const Index merged = c.coarsenIf([](Code, unsigned) { return true; });
        if (merged == 0 || c.numLeaves() == levels_.back().numLeaves())
          break;  // defensive: a level that did not shrink would stall the ladder
        levels_.push_back(std::move(c));
      }
    }
    ops_.resize(levels_.size());
    // All levels share the finest h0: a coarse octree's leaves carry a higher
    // `level`, and cellWidth = h0 * 2^level already encodes the doubled width.
    for (std::size_t L = 0; L < levels_.size(); ++L)
      ops_[L].init(levels_[L], h0);
    // child(fine slot) -> parent(coarse slot) for each fine/coarse pair.
    c2p_.assign(levels_.size() ? levels_.size() - 1 : 0, {});
    for (std::size_t L = 0; L + 1 < levels_.size(); ++L) {
      const Octree& f = levels_[L];
      const Octree& c = levels_[L + 1];
      c2p_[L].resize(static_cast<std::size_t>(f.numLeaves()));
      // Host-parallel (rung 2): `find` is a pure binary search over the frozen coarse octree and
      // each fine row owns one c2p slot — disjoint writes.
      hostParFor(f.numLeaves(), [&](Index i) {
        // Covering-leaf c2p (see multigrid.hpp): == ancestor+find for merged children, correct
        // (identity) for root-level rows in mixed-depth ladders, block-alignment-independent.
        c2p_[L][static_cast<std::size_t>(i)] = c.find(f.code(i));
      });
    }
  }
 
  std::size_t numLevels() const { return levels_.size(); }
  const AmrPoisson<Dim, Bits>& op(std::size_t L = 0) const { return ops_[L]; }
 
  void setPeriodic(bool p) {
    for (auto& o : ops_)
      o.setPeriodic(p);
  }
 
  void setImmersedWall(bool w) {
    for (auto& o : ops_)
      o.setImmersedWall(w);
  }
 
  void vcycle(std::size_t L, std::vector<double>& u, const std::vector<double>& rhs, int pre = 2,
              int post = 2) {
    if (L + 1 == levels_.size()) {
      ops_[L].gaussSeidel(u, rhs, 40);  // coarsest: solve hard
      ops_[L].removeMean(u);
      return;
    }
    ops_[L].gaussSeidel(u, rhs, pre);
    std::vector<double> res;
    ops_[L].residual(u, rhs, res);
 
    // Restrict residual: volume-weighted average of children -> parent.
    const Octree& f = levels_[L];
    const Octree& c = levels_[L + 1];
    std::vector<double> crhs(static_cast<std::size_t>(c.numLeaves()), 0.0);
    std::vector<double> cvol(static_cast<std::size_t>(c.numLeaves()), 0.0);
    for (Index i = 0; i < f.numLeaves(); ++i) {
      Index p = c2p_[L][static_cast<std::size_t>(i)];
      if (p < 0)
        continue;
      Real v = ops_[L].cellVolume(i);
      crhs[static_cast<std::size_t>(p)] += v * res[static_cast<std::size_t>(i)];
      cvol[static_cast<std::size_t>(p)] += v;
    }
    for (Index p = 0; p < c.numLeaves(); ++p)
      if (cvol[static_cast<std::size_t>(p)] > 0)
        crhs[static_cast<std::size_t>(p)] /= cvol[static_cast<std::size_t>(p)];
 
    std::vector<double> ccorr(static_cast<std::size_t>(c.numLeaves()), 0.0);
    vcycle(L + 1, ccorr, crhs, pre, post);
 
    // Prolong correction (piecewise constant) and add.
    for (Index i = 0; i < f.numLeaves(); ++i) {
      Index p = c2p_[L][static_cast<std::size_t>(i)];
      if (p >= 0)
        u[static_cast<std::size_t>(i)] += ccorr[static_cast<std::size_t>(p)];
    }
    ops_[L].gaussSeidel(u, rhs, post);
    ops_[L].removeMean(u);
  }
 
  double solveQuad(std::vector<double>& u, const std::vector<double>& rhs, int outer = 30,
                   int cyclesPerOuter = 1) {
    AmrPoisson<Dim, Bits>& P = ops_[0];
    const std::size_t n = u.size();
    std::vector<double> lq, ls, rhsp(n), res;
    double r = 0.0;
    for (int o = 0; o < outer; ++o) {
      P.applyLaplacianQuad(u, lq);
      P.applyLaplacian(u, ls);
      for (std::size_t i = 0; i < n; ++i)
        rhsp[i] = rhs[i] - (lq[i] - ls[i]);
      for (int c = 0; c < cyclesPerOuter; ++c)
        vcycle(0, u, rhsp);
      r = P.residualQuad(u, rhs, res);
    }
    return r;
  }
 
  template <class OpenFn>
  void setOpenness(OpenFn&& openFn) {
    if (levels_.empty())
      return;
    ops_[0].buildOpenness(openFn);
    for (std::size_t L = 0; L + 1 < levels_.size(); ++L)
      coarsenOpenness(L);
  }
 
  void setOpennessRaw(std::vector<double> alpha) {
    if (levels_.empty())
      return;
    ops_[0].setOpennessRaw(std::move(alpha));
    for (std::size_t L = 0; L + 1 < levels_.size(); ++L)
      coarsenOpenness(L);
  }
 
 private:
  static int faceIdx(int axis, int dir) { return 2 * axis + (dir > 0 ? 0 : 1); }
 
  // Area-average the level-L face openness onto level L+1 (each coarse face is the
  // mean of the 2^(Dim-1) fine sub-faces covering it; equal areas => plain mean).
  void coarsenOpenness(std::size_t L) {
    const Octree& f = levels_[L];
    const Octree& c = levels_[L + 1];
    const int F = 2 * Dim;
    std::vector<double> ca(static_cast<std::size_t>(c.numLeaves()) * F, 0.0);
    std::vector<int> cnt(static_cast<std::size_t>(c.numLeaves()) * F, 0);
    for (Index i = 0; i < f.numLeaves(); ++i) {
      Index p = c2p_[L][static_cast<std::size_t>(i)];
      if (p < 0)
        continue;
      const std::size_t base = static_cast<std::size_t>(p) * F;
      if (c.level(p) == f.level(i)) {
        // identity (cell not coarsened this round): copy every face.
        for (int axis = 0; axis < Dim; ++axis)
          for (int dir = -1; dir <= 1; dir += 2) {
            int fi = faceIdx(axis, dir);
            ca[base + fi] += ops_[L].faceOpenness(i, axis, dir);
            cnt[base + fi] += 1;
          }
      } else {
        // merged child: each axis contributes its outward face to the parent face.
        unsigned oct = M::from_code(f.code(i)).child_index(f.level(i));
        for (int axis = 0; axis < Dim; ++axis) {
          int dir = ((oct >> axis) & 1) ? +1 : -1;
          int fi = faceIdx(axis, dir);
          ca[base + fi] += ops_[L].faceOpenness(i, axis, dir);
          cnt[base + fi] += 1;
        }
      }
    }
    for (std::size_t k = 0; k < ca.size(); ++k)
      ca[k] = cnt[k] ? ca[k] / cnt[k] : 1.0;
    ops_[L + 1].setOpennessRaw(std::move(ca));
  }
 
  std::vector<Octree> levels_;
  std::vector<AmrPoisson<Dim, Bits>> ops_;
  std::vector<std::vector<Index>> c2p_;
};
 
}  // namespace peclet::amr
 
#endif  // PECLET_AMR_POISSON_HPP