block_exchange.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_VOF_BLOCK_EXCHANGE_HPP
#define PECLET_FLOW_VOF_BLOCK_EXCHANGE_HPP
 
#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <Kokkos_Core.hpp>
#include <map>
#include <stdexcept>
#include <type_traits>
#include <vector>
 
#ifdef PECLET_FLOW_MPI
#include <mpi.h>
#endif
 
#include "mac_cutcell.hpp"  // ccFor3 (rule H)
#include "vof/block_container.hpp"
#include "policy.hpp"
 
namespace peclet::flow::vof {
 
struct VofPiece {
  int rank = 0;
  VofBox g;             
  int loc[3]{0, 0, 0};  
  long cells() const { return g.cells(); }
};
 
struct VofMoveJob {
  const double* src;
  double* dst;
  I3 e, be;  
  I3 s, q;   
  I3 n;      
};
 
struct VofMoveTable {
  VofMoveJob job[kVofBlockBatch];
  int nj;
};
 
inline bool vofMoveOverlap(const VofMoveJob& a, const VofMoveJob& b) {
  if (a.dst != b.dst)
    return false;
  return a.q.x < b.q.x + b.n.x && b.q.x < a.q.x + a.n.x && a.q.y < b.q.y + b.n.y &&
         b.q.y < a.q.y + a.n.y && a.q.z < b.q.z + b.n.z && b.q.z < a.q.z + a.n.z;
}
 
inline void vofBuildPieces(const VofBox& box, I3 gs, const std::array<bool, 3>& per,
                           const std::vector<VofBox>& rankBox, std::vector<VofPiece>& out,
                           int outside = 0) {
  out.clear();
  VofRun rx[5], ry[5], rz[5];
  const int nx = vofAxisRuns(box.lo[0], box.hi[0], gs.x, per[0], rx, outside);
  const int ny = vofAxisRuns(box.lo[1], box.hi[1], gs.y, per[1], ry, outside);
  const int nz = vofAxisRuns(box.lo[2], box.hi[2], gs.z, per[2], rz, outside);
  const int gsa[3] = {gs.x, gs.y, gs.z};
  for (int iz = 0; iz < nz; ++iz)
    for (int iy = 0; iy < ny; ++iy)
      for (int ix = 0; ix < nx; ++ix) {
        VofBox s;
        s.lo[0] = rx[ix].g0;
        s.hi[0] = rx[ix].g0 + rx[ix].len;
        s.lo[1] = ry[iy].g0;
        s.hi[1] = ry[iy].g0 + ry[iy].len;
        s.lo[2] = rz[iz].g0;
        s.hi[2] = rz[iz].g0 + rz[iz].len;
        const int base[3] = {rx[ix].loc0, ry[iy].loc0, rz[iz].loc0};
        for (std::size_t r = 0; r < rankBox.size(); ++r) {
          VofBox own = rankBox[r];
          if (outside > 0)
            for (int d = 0; d < 3; ++d) {
              if (per[d])
                continue;
              if (own.lo[d] == 0)
                own.lo[d] = -outside;
              if (own.hi[d] == gsa[d])
                own.hi[d] = gsa[d] + outside;
            }
          const VofBox o = VofBox::intersect(s, own);
          if (o.empty())
            continue;
          VofPiece p;
          p.rank = static_cast<int>(r);
          p.g = o;
          for (int d = 0; d < 3; ++d)
            p.loc[d] = base[d] + (o.lo[d] - s.lo[d]);
          out.push_back(p);
        }
      }
}
 
class VofBlockExchange : public VofBlockExchangeBase {
 public:
  struct Patch {
    I3 e{0, 0, 0};  
    I3 n{0, 0, 0};  
    I3 o{0, 0, 0};  
    int g = 0;      
  };
 
  void init(I3 gs, std::array<bool, 3> per, std::vector<VofBox> rankBox, int rank) {
    gs_ = gs;
    per_ = per;
    rankBox_ = std::move(rankBox);
    rank_ = rank;
    size_ = static_cast<int>(rankBox_.size());
    if (rank_ < 0 || rank_ >= size_)
      throw std::runtime_error("peclet::flow::vof::VofBlockExchange: rank out of range");
  }
 
  void setRankBoxes(std::vector<VofBox> rankBox) {
    if (rankBox.size() != rankBox_.size())
      throw std::runtime_error(
          "peclet::flow::vof::VofBlockExchange::setRankBoxes: rank count changed");
    rankBox_ = std::move(rankBox);
  }
 
  void setPatch(Patch p, SField uf, SField vf, SField wf) {
    patch_ = p;
    lu_ = uf;
    lv_ = vf;
    lw_ = wf;
  }
 
#ifdef PECLET_FLOW_MPI
  void setComm(MPI_Comm c) {
    comm_ = c;
    hasComm_ = true;
  }
#endif
 
  long gatherBytes() const override { return gBytes_; }
  long scatterBytes() const override { return sBytes_; }
  long gatherMessages() const { return gMsgs_; }
  long scatterMessages() const { return sMsgs_; }
 
  // ================================================================== rung W1 item (b):
  // DEVICE-RESIDENT packing.
  //
  // The pack/unpack kernels run in the block's own memory space, so the only host traffic left is
  // the MPI staging copy of each MESSAGE (one `deep_copy` per peer, not one mirror of the whole
  // local patch per step) -- and for the master's OWN cells there is no host traffic at all: the
  // copy is device-to-device. That is exactly the order `core`'s grid halo grew its device-resident
  // variant (host-staged first, proven, then the device pack with a host staging buffer at the MPI
  // boundary), and the reason the staging survives is that the MPI here is not guaranteed
  // CUDA-aware. Every step is a copy of a double, so the device path is BITWISE the host path;
  // `deviceStaging = false` selects the W0 arm for the measurement.
  //
  // The four transfers of the container are one pattern with four instantiations:
  //
  //   gatherFaceVel     EXTENDED box, 3 components, patch -> block          (blockBase = 0)
  //   gatherColour      INNER box,    1 component,  patch -> block          (blockBase = ghost)
  //   scatterColourMax  INNER box,    1 component,  block -> patch, MAX     (blockBase = ghost)
  //   scatterForceSum   INNER box,    3 components, block -> patch, SUM     (blockBase = ghost)
  //
  // so `gatherImpl` / `scatterImpl` take the box selector, the component count, the block-array
  // base offset and (for the scatter) the combine rule, and nothing else differs.
  bool deviceStaging = true;
 
  void gatherImpl(std::vector<VofBlock>& blocks, int ghost, bool useExtended, int nc,
                  const SField* loc, SField (*blockView)(VofBlock&, int), int blockBase,
                  int outside = 0) {
    gBytes_ = 0;
    gMsgs_ = 0;
    bufSeq_ = 0;
    struct Recv {
      std::size_t bi;
      int from;
      long cells;
      SField dbuf;
      HostView hbuf;
    };
    std::vector<Recv> recvs;
    std::vector<SField> sendD(blocks.size());
    std::vector<HostView> sendH(blocks.size());
    std::vector<VofMoveJob> moves;  // H-4(e), host only: the master's own pieces, every block
#ifdef PECLET_FLOW_MPI
    std::vector<MPI_Request> reqs;
#endif
    for (std::size_t bi = 0; bi < blocks.size(); ++bi) {
      VofBlock& b = blocks[bi];
      const VofBox box = useExtended ? b.extended(ghost) : b.box;
      vofBuildPieces(box, gs_, per_, rankBox_, pieces_, outside);
      const bool master = (b.master == rank_);
      if (master)
        for (int c = 0; c < nc; ++c)
          Kokkos::deep_copy(SExec(), blockView(b, c), 0.0);  // async: device work follows
      std::map<int, long> counts;
      for (const auto& p : pieces_)
        counts[p.rank] += p.cells();
      if (master)
        for (const auto& kv : counts) {
          if (kv.first == rank_)
            continue;
          Recv r;
          r.bi = bi;
          r.from = kv.first;
          r.cells = kv.second;
          r.dbuf = buffer(nc * kv.second);
          r.hbuf = Kokkos::create_mirror_view(r.dbuf);
          recvs.push_back(std::move(r));
        }
      if (!counts.count(rank_))
        continue;
      if (master) {  // device-to-device, no host traffic at all
        for (const auto& p : pieces_)
          if (p.rank == rank_) {
            if (kVofHostExec<SExec>)  // H-4(e): collected, moved in one launch after the loop
              appendMoveJobs(p, nc, loc, blockView, b, blockBase, moves);
            else
              movePiece(p, nc, loc, blockView, b, blockBase, /*toBlock=*/true);
          }
      } else {
        sendD[bi] = buffer(nc * counts[rank_]);
        long off = 0;
        for (const auto& p : pieces_)
          if (p.rank == rank_)
            packFromPatch(p, nc, loc, sendD[bi], off);
        Kokkos::fence();
        sendH[bi] = Kokkos::create_mirror_view(sendD[bi]);
        Kokkos::deep_copy(sendH[bi], sendD[bi]);
      }
    }
    moveLocalPieces(moves);  // after every block's zero fill; the pieces are disjoint
    Kokkos::fence();
#ifdef PECLET_FLOW_MPI
    if (size_ > 1) {
      requireComm();
      for (auto& r : recvs) {
        reqs.push_back(MPI_REQUEST_NULL);
        MPI_Irecv(r.hbuf.data(), static_cast<int>(r.hbuf.extent(0)), MPI_DOUBLE, r.from,
                  tagOf(blocks[r.bi].id, 0), comm_, &reqs.back());
        gBytes_ += static_cast<long>(r.hbuf.extent(0)) * 8;
        ++gMsgs_;
      }
      for (std::size_t bi = 0; bi < blocks.size(); ++bi) {
        if (sendH[bi].extent(0) == 0)
          continue;
        reqs.push_back(MPI_REQUEST_NULL);
        MPI_Isend(sendH[bi].data(), static_cast<int>(sendH[bi].extent(0)), MPI_DOUBLE,
                  blocks[bi].master, tagOf(blocks[bi].id, 0), comm_, &reqs.back());
        gBytes_ += static_cast<long>(sendH[bi].extent(0)) * 8;
        ++gMsgs_;
      }
      if (!reqs.empty())
        MPI_Waitall(static_cast<int>(reqs.size()), reqs.data(), MPI_STATUSES_IGNORE);
    }
#endif
    for (auto& r : recvs) {
      VofBlock& b = blocks[r.bi];
      Kokkos::deep_copy(r.dbuf, r.hbuf);
      vofBuildPieces(useExtended ? b.extended(ghost) : b.box, gs_, per_, rankBox_, pieces_,
                     outside);
      long off = 0;
      for (const auto& p : pieces_)
        if (p.rank == r.from)
          unpackToBlock(p, nc, r.dbuf, off, blockView, b, blockBase);
    }
  }
 
  void scatterImpl(std::vector<VofBlock>& blocks, int nc, SField* loc,
                   SField (*blockView)(VofBlock&, int), int blockBase, int op,
                   bool gatherLocal = false) {
    sBytes_ = 0;
    sMsgs_ = 0;
    bufSeq_ = 0;
    if (!gatherLocal)
      zeroPatchInner(nc, loc);
    struct Recv {
      std::size_t bi;
      int from;
      SField dbuf;
      HostView hbuf;
    };
    std::vector<Recv> recvs;
    std::vector<SField> sendD;
    std::vector<HostView> sendH;
    std::vector<int> sendTo;
    std::vector<long> sendTag;
#ifdef PECLET_FLOW_MPI
    std::vector<MPI_Request> reqs;
#endif
    for (std::size_t bi = 0; bi < blocks.size(); ++bi) {
      VofBlock& b = blocks[bi];
      vofBuildPieces(b.box, gs_, per_, rankBox_, pieces_);
      const bool master = (b.master == rank_);
      std::map<int, long> counts;
      for (const auto& p : pieces_)
        counts[p.rank] += p.cells();
      if (master) {
        for (const auto& kv : counts) {
          if (kv.first == rank_) {
            if (!gatherLocal)
              for (const auto& p : pieces_)
                if (p.rank == rank_)
                  combineBlockIntoPatch(p, nc, blockView, b, blockBase, loc, op);
            continue;
          }
          SField d = buffer(nc * kv.second);
          long off = 0;
          for (const auto& p : pieces_)
            if (p.rank == kv.first)
              packFromBlock(p, nc, blockView, b, blockBase, d, off);
          Kokkos::fence();
          HostView h = Kokkos::create_mirror_view(d);
          Kokkos::deep_copy(h, d);
          sendD.push_back(d);
          sendH.push_back(h);
          sendTo.push_back(kv.first);
          sendTag.push_back(b.id);
        }
      } else if (counts.count(rank_)) {
        Recv r;
        r.bi = bi;
        r.from = b.master;
        r.dbuf = buffer(nc * counts[rank_]);
        r.hbuf = Kokkos::create_mirror_view(r.dbuf);
        recvs.push_back(std::move(r));
      }
    }
    if (gatherLocal)
      gatherLocalSum(blocks, nc, loc, blockView, blockBase);
    Kokkos::fence();
#ifdef PECLET_FLOW_MPI
    if (size_ > 1) {
      requireComm();
      for (auto& r : recvs) {
        reqs.push_back(MPI_REQUEST_NULL);
        MPI_Irecv(r.hbuf.data(), static_cast<int>(r.hbuf.extent(0)), MPI_DOUBLE, r.from,
                  tagOf(blocks[r.bi].id, 1), comm_, &reqs.back());
        sBytes_ += static_cast<long>(r.hbuf.extent(0)) * 8;
        ++sMsgs_;
      }
      for (std::size_t k = 0; k < sendH.size(); ++k) {
        reqs.push_back(MPI_REQUEST_NULL);
        MPI_Isend(sendH[k].data(), static_cast<int>(sendH[k].extent(0)), MPI_DOUBLE, sendTo[k],
                  tagOf(sendTag[k], 1), comm_, &reqs.back());
        sBytes_ += static_cast<long>(sendH[k].extent(0)) * 8;
        ++sMsgs_;
      }
      if (!reqs.empty())
        MPI_Waitall(static_cast<int>(reqs.size()), reqs.data(), MPI_STATUSES_IGNORE);
    }
#endif
    for (auto& r : recvs) {
      VofBlock& b = blocks[r.bi];
      Kokkos::deep_copy(r.dbuf, r.hbuf);
      vofBuildPieces(b.box, gs_, per_, rankBox_, pieces_);
      long off = 0;
      for (const auto& p : pieces_)
        if (p.rank == rank_)
          combineBufIntoPatch(p, nc, r.dbuf, off, loc, op);
    }
  }
 
  void gatherLocalSum(std::vector<VofBlock>& blocks, int nc, SField* loc,
                      SField (*blockView)(VofBlock&, int), int blockBase) {
    struct Job {
      const double* f[3];
      int lo[3], n[3];
      I3 be;
    };
    std::vector<Job> jobs;
    for (auto& b : blocks) {
      if (b.master != rank_)
        continue;
      Job j;
      for (int c = 0; c < 3; ++c)
        j.f[c] = blockView(b, c < nc ? c : 0).data();
      for (int d = 0; d < 3; ++d) {
        j.lo[d] = b.box.lo[d];
        j.n[d] = b.box.n(d);
      }
      j.be = b.advector().extent();
      jobs.push_back(j);
    }
    GatherTable T;
    T.nc = nc;
    T.gh = blockBase;
    for (int d = 0; d < 3; ++d)
      T.per[d] = per_[d];
    T.L[0] = gs_.x;
    T.L[1] = gs_.y;
    T.L[2] = gs_.z;
    std::size_t j0 = 0;
    do {
      T.nj = static_cast<int>(std::min<std::size_t>(kVofBlockBatch, jobs.size() - j0));
      T.first = (j0 == 0);
      for (int k = 0; k < T.nj; ++k) {
        for (int c = 0; c < 3; ++c)
          T.f[k][c] = jobs[j0 + k].f[c];
        for (int d = 0; d < 3; ++d) {
          T.lo[k][d] = jobs[j0 + k].lo[d];
          T.n[k][d] = jobs[j0 + k].n[d];
        }
        T.be[k] = jobs[j0 + k].be;
      }
      gatherLocalLaunch(T, loc);
      j0 += kVofBlockBatch;
    } while (j0 < jobs.size());
  }
  struct GatherTable {
    const double* f[kVofBlockBatch][3];
    int lo[kVofBlockBatch][3], n[kVofBlockBatch][3];
    I3 be[kVofBlockBatch];
    int nj, nc, gh, L[3];
    bool per[3], first;
  };
  void gatherLocalLaunch(const GatherTable& T, SField* loc) {
    const I3 e = patch_.e, n = patch_.n, o = patch_.o;
    const int g = patch_.g;
    SField d0 = loc[0], d1 = T.nc > 1 ? loc[1] : loc[0], d2 = T.nc > 2 ? loc[2] : loc[0];
    if constexpr (std::is_same_v<typename SField::memory_space, Kokkos::HostSpace>) {
      // §14 H-3(b), host: rule H (ccFor3: rows, x inside) and per-axis job masks. m_d[x_d] holds
      // bit k iff job k's box contains patch coordinate x_d on axis d under the periodic predicate
      // below, so a cell visits exactly the jobs in mx[x] & my[y] & mz[z], in ascending k, with
      // l[d] computed as there: the same jobs, order and additions -- bitwise -- without the 16
      // jobs x 3 axes of ((v % L) + L) % L per cell. The device keeps the launch below.
      static_assert(kVofBlockBatch <= 32, "the job masks are 32-bit");
      const int nn[3] = {n.x, n.y, n.z}, oo[3] = {o.x, o.y, o.z};
      std::vector<std::uint32_t> mk((std::size_t)n.x + n.y + n.z, 0u);
      std::uint32_t* md[3] = {mk.data(), mk.data() + n.x, mk.data() + n.x + n.y};
      for (int d = 0; d < 3; ++d)
        for (int xd = 0; xd < nn[d]; ++xd)
          for (int k = 0; k < T.nj; ++k) {
            int v = xd + oo[d] - T.lo[k][d];
            if (T.per[d])
              v = ((v % T.L[d]) + T.L[d]) % T.L[d];
            if (v >= 0 && v < T.n[k][d])
              md[d][xd] |= (1u << k);
          }
      const std::uint32_t *mx = md[0], *my = md[1], *mz = md[2];
      ccFor3(
          "vof::block::gather_local_sum", C3{0, 0, 0}, C3{n.x, n.y, n.z},
          KOKKOS_LAMBDA(int x, int y, int z) {
            const int gp[3] = {x + o.x, y + o.y, z + o.z};
            const long i = L3(x + g, y + g, z + g, e);
            double acc[3] = {0.0, 0.0, 0.0};
            if (!T.first) {
              acc[0] = d0(i);
              acc[1] = T.nc > 1 ? d1(i) : 0.0;
              acc[2] = T.nc > 2 ? d2(i) : 0.0;
            }
            const std::uint32_t m = mx[x] & my[y] & mz[z];
            for (int k = 0; k < T.nj; ++k) {
              if (!((m >> k) & 1u))
                continue;
              int l[3];
              for (int d = 0; d < 3; ++d) {
                int v = gp[d] - T.lo[k][d];
                if (T.per[d])
                  v = ((v % T.L[d]) + T.L[d]) % T.L[d];
                l[d] = v;
              }
              const long q = L3(l[0] + T.gh, l[1] + T.gh, l[2] + T.gh, T.be[k]);
              for (int c = 0; c < T.nc; ++c)
                acc[c] += T.f[k][c][q];
            }
            d0(i) = acc[0];
            if (T.nc > 1)
              d1(i) = acc[1];
            if (T.nc > 2)
              d2(i) = acc[2];
          });
      Kokkos::fence();  // the masks live on this stack frame
      return;
    }
    Kokkos::parallel_for(
        "vof::block::gather_local_sum", MDRange3<SExec>(SExec(), {0, 0, 0}, {n.x, n.y, n.z}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          const int gp[3] = {x + o.x, y + o.y, z + o.z};
          const long i = L3(x + g, y + g, z + g, e);
          double acc[3] = {0.0, 0.0, 0.0};
          if (!T.first) {
            acc[0] = d0(i);
            acc[1] = T.nc > 1 ? d1(i) : 0.0;
            acc[2] = T.nc > 2 ? d2(i) : 0.0;
          }
          for (int k = 0; k < T.nj; ++k) {
            int l[3];
            bool in = true;
            for (int d = 0; d < 3 && in; ++d) {
              int v = gp[d] - T.lo[k][d];
              if (T.per[d])
                v = ((v % T.L[d]) + T.L[d]) % T.L[d];
              in = v >= 0 && v < T.n[k][d];
              l[d] = v;
            }
            if (!in)
              continue;
            const long q = L3(l[0] + T.gh, l[1] + T.gh, l[2] + T.gh, T.be[k]);
            for (int c = 0; c < T.nc; ++c)
              acc[c] += T.f[k][c][q];
          }
          d0(i) = acc[0];
          if (T.nc > 1)
            d1(i) = acc[1];
          if (T.nc > 2)
            d2(i) = acc[2];
        });
  }
 
  // ---- the four device kernels (public for nvcc's extended-lambda rule) -----------------------
  void zeroPatchInner(int nc, SField* loc) {
    const I3 e = patch_.e, n = patch_.n;
    const int g = patch_.g;
    for (int c = 0; c < nc; ++c) {
      SField f = loc[c];
      Kokkos::parallel_for(
          "vof::block::zero_patch",
          MDRange3<SExec>(SExec(), {0, 0, 0}, {n.x, n.y, n.z}),
          KOKKOS_LAMBDA(int x, int y, int z) { f(L3(x + g, y + g, z + g, e)) = 0.0; });
    }
  }
 
  void packFromPatch(const VofPiece& p, int nc, const SField* loc, SField buf, long& off) {
    const I3 e = patch_.e, o = patch_.o;
    const int g = patch_.g;
    const int lo0 = p.g.lo[0], lo1 = p.g.lo[1], lo2 = p.g.lo[2];
    const int n0 = p.g.n(0), n1 = p.g.n(1), n2 = p.g.n(2);
    const long base = off;
    SField s0 = loc[0], s1 = nc > 1 ? loc[1] : loc[0], s2 = nc > 2 ? loc[2] : loc[0];
    Kokkos::parallel_for(
        "vof::block::pack_patch",
        MDRange3<SExec>(SExec(), {0, 0, 0}, {n0, n1, n2}),
        KOKKOS_LAMBDA(int x, int y, int z) {
          const long i = L3(lo0 + x - o.x + g, lo1 + y - o.y + g, lo2 + z - o.z + g, e);
          const long k = base + x + static_cast<long>(y) * n0 + static_cast<long>(z) * n0 * n1;
          buf(nc * k + 0) = s0(i);
          if (nc > 1)
            buf(nc * k + 1) = s1(i);
          if (nc > 2)
            buf(nc * k + 2) = s2(i);
        });
    off += static_cast<long>(n0) * n1 * n2;
  }
 
  void unpackToBlock(const VofPiece& p, int nc, SField buf, long& off,
                     SField (*blockView)(VofBlock&, int), VofBlock& b, int blockBase) {
    const I3 be = b.advector().extent();
    const int lo0 = p.g.lo[0], lo1 = p.g.lo[1], lo2 = p.g.lo[2];
    const int n0 = p.g.n(0), n1 = p.g.n(1), n2 = p.g.n(2);
    const int q0 = p.loc[0] + blockBase, q1 = p.loc[1] + blockBase, q2 = p.loc[2] + blockBase;
    const long base = off;
    (void)lo0;
    (void)lo1;
    (void)lo2;
    for (int c = 0; c < nc; ++c) {
      SField d = blockView(b, c);
      const int cc = c;
      Kokkos::parallel_for(
          "vof::block::unpack_block",
          MDRange3<SExec>(SExec(), {0, 0, 0}, {n0, n1, n2}),
          KOKKOS_LAMBDA(int x, int y, int z) {
            const long k = base + x + static_cast<long>(y) * n0 + static_cast<long>(z) * n0 * n1;
            d(L3(q0 + x, q1 + y, q2 + z, be)) = buf(nc * k + cc);
          });
    }
    off += static_cast<long>(n0) * n1 * n2;
  }
 
  void movePiece(const VofPiece& p, int nc, const SField* loc,
                 SField (*blockView)(VofBlock&, int), VofBlock& b, int blockBase, bool toBlock) {
    const I3 e = patch_.e, o = patch_.o, be = b.advector().extent();
    const int g = patch_.g;
    const int lo0 = p.g.lo[0], lo1 = p.g.lo[1], lo2 = p.g.lo[2];
    const int n0 = p.g.n(0), n1 = p.g.n(1), n2 = p.g.n(2);
    const int q0 = p.loc[0] + blockBase, q1 = p.loc[1] + blockBase, q2 = p.loc[2] + blockBase;
    for (int c = 0; c < nc; ++c) {
      SField src = loc[c], dst = blockView(b, c);
      if (!toBlock)
        throw std::runtime_error("vof::block: movePiece is patch->block only");
      Kokkos::parallel_for(
          "vof::block::move_local",
          MDRange3<SExec>(SExec(), {0, 0, 0}, {n0, n1, n2}),
          KOKKOS_LAMBDA(int x, int y, int z) {
            dst(L3(q0 + x, q1 + y, q2 + z, be)) =
                src(L3(lo0 + x - o.x + g, lo1 + y - o.y + g, lo2 + z - o.z + g, e));
          });
    }
  }
 
  void appendMoveJobs(const VofPiece& p, int nc, const SField* loc,
                      SField (*blockView)(VofBlock&, int), VofBlock& b, int blockBase,
                      std::vector<VofMoveJob>& out) {
    const I3 o = patch_.o;
    const int g = patch_.g;
    for (int c = 0; c < nc; ++c) {
      VofMoveJob J;
      J.src = loc[c].data();
      J.dst = blockView(b, c).data();
      J.e = patch_.e;
      J.be = b.advector().extent();
      J.s = I3{p.g.lo[0] - o.x + g, p.g.lo[1] - o.y + g, p.g.lo[2] - o.z + g};
      J.q = I3{p.loc[0] + blockBase, p.loc[1] + blockBase, p.loc[2] + blockBase};
      J.n = I3{p.g.n(0), p.g.n(1), p.g.n(2)};
      out.push_back(J);
    }
  }
 
  template <class Exec = SExec>
  void moveLocalPieces(const std::vector<VofMoveJob>& mv) {
    if constexpr (kVofHostExec<Exec>) {
      for (std::size_t j0 = 0; j0 < mv.size(); j0 += kVofBlockBatch) {
        VofMoveTable T;
        T.nj = static_cast<int>(std::min<std::size_t>(kVofBlockBatch, mv.size() - j0));
        for (int k = 0; k < T.nj; ++k)
          T.job[k] = mv[j0 + k];
#ifndef NDEBUG
        for (int k = 0; k < T.nj; ++k)
          for (int l = 0; l < k; ++l)
            assert(!vofMoveOverlap(T.job[k], T.job[l]) && "vof::block: movePiece pieces overlap");
#endif
        vofHostRowFor<Exec>(
            "vof::block::move_local", T.nj,
            [&](const int k) { return static_cast<long>(T.job[k].n.y) * T.job[k].n.z; },
            [=](const int k, const long r) {
              const VofMoveJob& J = T.job[k];
              const int y = static_cast<int>(r % J.n.y), z = static_cast<int>(r / J.n.y);
              const double* src = J.src + L3(J.s.x, J.s.y + y, J.s.z + z, J.e);
              double* dst = J.dst + L3(J.q.x, J.q.y + y, J.q.z + z, J.be);
              PECLET_FLOW_OMP_SIMD
              for (int x = 0; x < J.n.x; ++x)
                dst[x] = src[x];
            });
      }
    }
  }
 
  void packFromBlock(const VofPiece& p, int nc, SField (*blockView)(VofBlock&, int), VofBlock& b,
                     int blockBase, SField buf, long& off) {
    const I3 be = b.advector().extent();
    const int n0 = p.g.n(0), n1 = p.g.n(1), n2 = p.g.n(2);
    const int q0 = p.loc[0] + blockBase, q1 = p.loc[1] + blockBase, q2 = p.loc[2] + blockBase;
    const long base = off;
    for (int c = 0; c < nc; ++c) {
      SField src = blockView(b, c);
      const int cc = c;
      Kokkos::parallel_for(
          "vof::block::pack_block",
          MDRange3<SExec>(SExec(), {0, 0, 0}, {n0, n1, n2}),
          KOKKOS_LAMBDA(int x, int y, int z) {
            const long k = base + x + static_cast<long>(y) * n0 + static_cast<long>(z) * n0 * n1;
            buf(nc * k + cc) = src(L3(q0 + x, q1 + y, q2 + z, be));
          });
    }
    off += static_cast<long>(n0) * n1 * n2;
  }
 
  void combineBlockIntoPatch(const VofPiece& p, int nc, SField (*blockView)(VofBlock&, int),
                             VofBlock& b, int blockBase, SField* loc, int op) {
    const I3 e = patch_.e, o = patch_.o, be = b.advector().extent();
    const int g = patch_.g;
    const int lo0 = p.g.lo[0], lo1 = p.g.lo[1], lo2 = p.g.lo[2];
    const int n0 = p.g.n(0), n1 = p.g.n(1), n2 = p.g.n(2);
    const int q0 = p.loc[0] + blockBase, q1 = p.loc[1] + blockBase, q2 = p.loc[2] + blockBase;
    for (int c = 0; c < nc; ++c) {
      SField src = blockView(b, c), dst = loc[c];
      const int opc = op;
      Kokkos::parallel_for(
          "vof::block::combine_block",
          MDRange3<SExec>(SExec(), {0, 0, 0}, {n0, n1, n2}),
          KOKKOS_LAMBDA(int x, int y, int z) {
            const double v = src(L3(q0 + x, q1 + y, q2 + z, be));
            const long i = L3(lo0 + x - o.x + g, lo1 + y - o.y + g, lo2 + z - o.z + g, e);
            if (opc == 0) {
              if (v > dst(i))
                dst(i) = v;
            } else {
              dst(i) += v;
            }
          });
    }
  }
 
  void combineBufIntoPatch(const VofPiece& p, int nc, SField buf, long& off, SField* loc, int op) {
    const I3 e = patch_.e, o = patch_.o;
    const int g = patch_.g;
    const int lo0 = p.g.lo[0], lo1 = p.g.lo[1], lo2 = p.g.lo[2];
    const int n0 = p.g.n(0), n1 = p.g.n(1), n2 = p.g.n(2);
    const long base = off;
    for (int c = 0; c < nc; ++c) {
      SField dst = loc[c];
      const int cc = c, opc = op;
      Kokkos::parallel_for(
          "vof::block::combine_buf",
          MDRange3<SExec>(SExec(), {0, 0, 0}, {n0, n1, n2}),
          KOKKOS_LAMBDA(int x, int y, int z) {
            const long k = base + x + static_cast<long>(y) * n0 + static_cast<long>(z) * n0 * n1;
            const double v = buf(nc * k + cc);
            const long i = L3(lo0 + x - o.x + g, lo1 + y - o.y + g, lo2 + z - o.z + g, e);
            if (opc == 0) {
              if (v > dst(i))
                dst(i) = v;
            } else {
              dst(i) += v;
            }
          });
    }
    off += static_cast<long>(n0) * n1 * n2;
  }
 
  // ---- the four public transfers ----------------------------------------------------------
  static SField faceViewOf(VofBlock& b, int c) { return b.advector().faceVel(c); }
  static SField colourViewOf(VofBlock& b, int c) {
    (void)c;
    return b.advector().colour();
  }
  static SField forceViewOf(VofBlock& b, int c) { return b.csfForce(c); }
 
  void gatherFaceVel(std::vector<VofBlock>& blocks, int ghost) override {
    // The block's extended box reaches `ghost` cells beyond a non-periodic domain face; those
    // cells are read from the boundary rank's patch ghosts, so the patch must be at least as wide.
    if (ghost > patch_.g)
      throw std::runtime_error(
          "peclet::flow::vof::VofBlockExchange: block ghost wider than the face-velocity patch");
    if (!deviceStaging) {
      gatherFaceVelHost(blocks, ghost);
      return;
    }
    SField loc[3] = {lu_, lv_, lw_};
    gatherImpl(blocks, ghost, /*useExtended=*/true, 3, loc, &faceViewOf, 0, /*outside=*/ghost);
  }
 
  void gatherColour(std::vector<VofBlock>& blocks, SField cLocal) override {
    SField loc[1] = {cLocal};
    gatherImpl(blocks, 0, /*useExtended=*/false, 1, loc, &colourViewOf, ghostOf(blocks));
  }
 
  void scatterColourMax(std::vector<VofBlock>& blocks, SField cLocal) override {
    if (!deviceStaging) {
      scatterColourMaxHost(blocks, cLocal);
      return;
    }
    SField loc[1] = {cLocal};
    scatterImpl(blocks, 1, loc, &colourViewOf, ghostOf(blocks), /*op=*/0);
  }
 
  void scatterForceSum(std::vector<VofBlock>& blocks, SField fx, SField fy, SField fz) override {
    SField loc[3] = {fx, fy, fz};
    scatterImpl(blocks, 3, loc, &forceViewOf, ghostOf(blocks), /*op=*/1, /*gatherLocal=*/true);
  }
 
  void scatterColourSum(std::vector<VofBlock>& blocks, SField sLocal) override {
    SField loc[1] = {sLocal};
    scatterImpl(blocks, 1, loc, &colourViewOf, ghostOf(blocks), /*op=*/1);
  }
 
  void migrateColour(std::vector<VofBlock>& blocks, const std::vector<int>& oldMaster) override {
    mBytes_ = 0;
    mMsgs_ = 0;
    bufSeq_ = 0;
#ifdef PECLET_FLOW_MPI
    if (size_ <= 1)
      return;
    requireComm();
    std::vector<MPI_Request> reqs;
    std::vector<SField> dbuf(blocks.size());
    std::vector<HostView> hbuf(blocks.size());
    std::vector<std::vector<double>> hAux(blocks.size());
    for (std::size_t i = 0; i < blocks.size(); ++i) {
      const int from = oldMaster[i], to = blocks[i].master;
      if (from == to || (from != rank_ && to != rank_))
        continue;
      const I3 e = blocks[i].advector().extent();
      const long len = static_cast<long>(e.x) * e.y * e.z;
      dbuf[i] = buffer(len);
      hbuf[i] = Kokkos::create_mirror_view(dbuf[i]);
      hAux[i].resize(VofBlock::kAuxLen, 0.0);
      if (from == rank_) {
        Kokkos::deep_copy(dbuf[i], blocks[i].advector().colour());
        Kokkos::deep_copy(hbuf[i], dbuf[i]);
        blocks[i].serializeAux(hAux[i].data());
        reqs.push_back(MPI_REQUEST_NULL);
        MPI_Isend(hbuf[i].data(), static_cast<int>(len), MPI_DOUBLE, to, tagOf(blocks[i].id, 2),
                  comm_, &reqs.back());
        reqs.push_back(MPI_REQUEST_NULL);
        MPI_Isend(hAux[i].data(), VofBlock::kAuxLen, MPI_DOUBLE, to, tagOf(blocks[i].id, 3), comm_,
                  &reqs.back());
      } else {
        reqs.push_back(MPI_REQUEST_NULL);
        MPI_Irecv(hbuf[i].data(), static_cast<int>(len), MPI_DOUBLE, from, tagOf(blocks[i].id, 2),
                  comm_, &reqs.back());
        reqs.push_back(MPI_REQUEST_NULL);
        MPI_Irecv(hAux[i].data(), VofBlock::kAuxLen, MPI_DOUBLE, from, tagOf(blocks[i].id, 3),
                  comm_, &reqs.back());
      }
      mBytes_ += (len + VofBlock::kAuxLen) * 8;
      mMsgs_ += 2;
    }
    if (!reqs.empty())
      MPI_Waitall(static_cast<int>(reqs.size()), reqs.data(), MPI_STATUSES_IGNORE);
    for (std::size_t i = 0; i < blocks.size(); ++i)
      if (oldMaster[i] != blocks[i].master && blocks[i].master == rank_) {
        Kokkos::deep_copy(dbuf[i], hbuf[i]);
        Kokkos::deep_copy(blocks[i].advector().colour(), dbuf[i]);
        blocks[i].deserializeAux(hAux[i].data());
      }
    Kokkos::fence();
#else
    (void)blocks;
    (void)oldMaster;
#endif
  }
  long migrateBytes() const { return mBytes_; }
  long migrateMessages() const { return mMsgs_; }
 
  SField buffer(long n) {
    const std::size_t k = static_cast<std::size_t>(bufSeq_++);
    if (bufPool_.size() <= k)
      bufPool_.resize(k + 1);
    if (bufPool_[k].extent(0) != static_cast<std::size_t>(n))
      bufPool_[k] = SField(Kokkos::view_alloc("vof::block::buf", Kokkos::WithoutInitializing),
                           static_cast<std::size_t>(n));
    return bufPool_[k];
  }
 
  static int ghostOf(const std::vector<VofBlock>& blocks) {
    for (const auto& b : blocks)
      if (b.mine())
        return b.advector().ghost();
    return 3;
  }
 
  // ------------------------------------------------------- gather (rung W0, HOST-staged)
  //
  // Kept as the reference implementation and as the `deviceStaging = false` arm of the W1
  // measurement: it mirrors the whole local patch to the host, packs with host loops, and mirrors
  // the block arrays back. The device path below produces BIT-IDENTICAL buffers (every step is a
  // copy of a double), which is what the W1 ctest gates.
  void gatherFaceVelHost(std::vector<VofBlock>& blocks, int ghost) {
    gBytes_ = 0;
    gMsgs_ = 0;
    auto hu = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), lu_);
    auto hv = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), lv_);
    auto hw = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), lw_);
 
    std::vector<std::vector<double>> sendBuf(blocks.size());
    std::vector<HostView> hbu(blocks.size()), hbv(blocks.size()), hbw(blocks.size());
    struct Recv {
      std::size_t bi;
      int from;
      std::vector<double> buf;
    };
    std::vector<Recv> recvs;
#ifdef PECLET_FLOW_MPI
    std::vector<MPI_Request> reqs;
#endif
 
    // 1. post receives (master side) and pack+send (owner side); do the master's own cells locally.
    for (std::size_t bi = 0; bi < blocks.size(); ++bi) {
      VofBlock& b = blocks[bi];
      const VofBox eb = b.extended(ghost);
      vofBuildPieces(eb, gs_, per_, rankBox_, pieces_, /*outside=*/ghost);
      const bool master = (b.master == rank_);
      if (master) {
        hbu[bi] = Kokkos::create_mirror_view(b.advector().faceU());
        hbv[bi] = Kokkos::create_mirror_view(b.advector().faceV());
        hbw[bi] = Kokkos::create_mirror_view(b.advector().faceW());
        zero(hbu[bi]);
        zero(hbv[bi]);
        zero(hbw[bi]);
      }
      // per-rank aggregated counts, in the canonical piece order
      std::map<int, long> counts;
      for (const auto& p : pieces_)
        counts[p.rank] += p.cells();
      if (master) {
        for (const auto& kv : counts) {
          if (kv.first == rank_)
            continue;
          recvs.push_back(Recv{bi, kv.first, std::vector<double>(3 * kv.second)});
        }
      }
      if (counts.count(rank_)) {
        if (master) {
          // local: straight copy, no buffer
          long off = 0;
          for (const auto& p : pieces_) {
            if (p.rank != rank_)
              continue;
            copyPatchToBlock(p, hu, hv, hw, hbu[bi], hbv[bi], hbw[bi], b, ghost);
            off += p.cells();
          }
          (void)off;
        } else {
          std::vector<double>& sb = sendBuf[bi];
          sb.resize(3 * counts[rank_]);
          long off = 0;
          for (const auto& p : pieces_) {
            if (p.rank != rank_)
              continue;
            packPatch(p, hu, hv, hw, sb, off);
          }
        }
      }
    }
#ifdef PECLET_FLOW_MPI
    if (size_ > 1) {
      requireComm();
      for (auto& r : recvs) {
        reqs.push_back(MPI_REQUEST_NULL);
        MPI_Irecv(r.buf.data(), static_cast<int>(r.buf.size()), MPI_DOUBLE, r.from,
                  tagOf(blocks[r.bi].id, 0), comm_, &reqs.back());
        gBytes_ += static_cast<long>(r.buf.size()) * 8;
        ++gMsgs_;
      }
      for (std::size_t bi = 0; bi < blocks.size(); ++bi) {
        if (sendBuf[bi].empty())
          continue;
        reqs.push_back(MPI_REQUEST_NULL);
        MPI_Isend(sendBuf[bi].data(), static_cast<int>(sendBuf[bi].size()), MPI_DOUBLE,
                  blocks[bi].master, tagOf(blocks[bi].id, 0), comm_, &reqs.back());
        gBytes_ += static_cast<long>(sendBuf[bi].size()) * 8;
        ++gMsgs_;
      }
      if (!reqs.empty())
        MPI_Waitall(static_cast<int>(reqs.size()), reqs.data(), MPI_STATUSES_IGNORE);
    }
#endif
    // 2. unpack
    for (auto& r : recvs) {
      VofBlock& b = blocks[r.bi];
      vofBuildPieces(b.extended(ghost), gs_, per_, rankBox_, pieces_, /*outside=*/ghost);
      long off = 0;
      for (const auto& p : pieces_) {
        if (p.rank != r.from)
          continue;
        unpackBlock(p, r.buf, off, hbu[r.bi], hbv[r.bi], hbw[r.bi], b, ghost);
      }
    }
    for (std::size_t bi = 0; bi < blocks.size(); ++bi) {
      if (blocks[bi].master != rank_)
        continue;
      Kokkos::deep_copy(blocks[bi].advector().faceU(), hbu[bi]);
      Kokkos::deep_copy(blocks[bi].advector().faceV(), hbv[bi]);
      Kokkos::deep_copy(blocks[bi].advector().faceW(), hbw[bi]);
    }
  }
 
  // ------------------------------------------------------ scatter (rung W0, HOST-staged)
  void scatterColourMaxHost(std::vector<VofBlock>& blocks, SField cLocal) {
    sBytes_ = 0;
    sMsgs_ = 0;
    auto hc = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), cLocal);
    // UNPACK_MAX starts from an empty union: every inner cell not covered by a block is gas.
    const I3 e = patch_.e, n = patch_.n;
    const int g = patch_.g;
    for (int z = 0; z < n.z; ++z)
      for (int y = 0; y < n.y; ++y)
        for (int x = 0; x < n.x; ++x)
          hc(L3(x + g, y + g, z + g, e)) = 0.0;
 
    struct Recv {
      std::size_t bi;
      int from;
      std::vector<double> buf;
    };
    std::vector<Recv> recvs;
    std::vector<std::vector<double>> sendBuf;
    std::vector<int> sendTo;
    std::vector<long> sendTag;
#ifdef PECLET_FLOW_MPI
    std::vector<MPI_Request> reqs;
#endif
    std::vector<HostView> hbc(blocks.size());
 
    for (std::size_t bi = 0; bi < blocks.size(); ++bi) {
      VofBlock& b = blocks[bi];
      vofBuildPieces(b.box, gs_, per_, rankBox_, pieces_);
      const bool master = (b.master == rank_);
      if (master) {
        hbc[bi] = Kokkos::create_mirror_view(b.advector().colour());
        Kokkos::deep_copy(hbc[bi], b.advector().colour());
      }
      std::map<int, long> counts;
      for (const auto& p : pieces_)
        counts[p.rank] += p.cells();
      if (master) {
        for (const auto& kv : counts) {
          if (kv.first == rank_) {
            for (const auto& p : pieces_)
              if (p.rank == rank_)
                maxBlockIntoPatch(p, hbc[bi], b, hc);
            continue;
          }
          std::vector<double> sb(kv.second);
          long off = 0;
          for (const auto& p : pieces_)
            if (p.rank == kv.first)
              packBlockColour(p, hbc[bi], b, sb, off);
          sendBuf.push_back(std::move(sb));
          sendTo.push_back(kv.first);
          sendTag.push_back(b.id);
        }
      } else if (counts.count(rank_)) {
        recvs.push_back(Recv{bi, b.master, std::vector<double>(counts[rank_])});
      }
    }
#ifdef PECLET_FLOW_MPI
    if (size_ > 1) {
      requireComm();
      for (auto& r : recvs) {
        reqs.push_back(MPI_REQUEST_NULL);
        MPI_Irecv(r.buf.data(), static_cast<int>(r.buf.size()), MPI_DOUBLE, r.from,
                  tagOf(blocks[r.bi].id, 1), comm_, &reqs.back());
        sBytes_ += static_cast<long>(r.buf.size()) * 8;
        ++sMsgs_;
      }
      for (std::size_t k = 0; k < sendBuf.size(); ++k) {
        reqs.push_back(MPI_REQUEST_NULL);
        MPI_Isend(sendBuf[k].data(), static_cast<int>(sendBuf[k].size()), MPI_DOUBLE, sendTo[k],
                  tagOf(sendTag[k], 1), comm_, &reqs.back());
        sBytes_ += static_cast<long>(sendBuf[k].size()) * 8;
        ++sMsgs_;
      }
      if (!reqs.empty())
        MPI_Waitall(static_cast<int>(reqs.size()), reqs.data(), MPI_STATUSES_IGNORE);
    }
#endif
    for (auto& r : recvs) {
      VofBlock& b = blocks[r.bi];
      vofBuildPieces(b.box, gs_, per_, rankBox_, pieces_);
      long off = 0;
      for (const auto& p : pieces_)
        if (p.rank == rank_)
          maxBufIntoPatch(p, r.buf, off, hc);
    }
    Kokkos::deep_copy(cLocal, hc);
  }
 
  // -------------------------------------------------------------------------------- table
  void syncTable(std::vector<VofBlock>& blocks) override {
#ifdef PECLET_FLOW_MPI
    if (size_ <= 1)
      return;
    requireComm();
    const int nb = static_cast<int>(blocks.size());
    std::vector<int> mine(6 * nb, 0), all(6 * nb, 0);
    for (int i = 0; i < nb; ++i)
      if (blocks[i].master == rank_)
        for (int d = 0; d < 3; ++d) {
          mine[6 * i + d] = blocks[i].box.lo[d];
          mine[6 * i + 3 + d] = blocks[i].box.hi[d];
        }
    // exactly one rank contributes each block's row, so a SUM is a broadcast that needs no
    // per-block communicator and no ordering assumption.
    MPI_Allreduce(mine.data(), all.data(), 6 * nb, MPI_INT, MPI_SUM, comm_);
    for (int i = 0; i < nb; ++i)
      for (int d = 0; d < 3; ++d) {
        blocks[i].box.lo[d] = all[6 * i + d];
        blocks[i].box.hi[d] = all[6 * i + 3 + d];
      }
#else
    (void)blocks;
#endif
  }
 
 private:
  using HostView = decltype(Kokkos::create_mirror_view(SField()));
 
  static void zero(HostView v) {
    for (std::size_t i = 0; i < v.extent(0); ++i)
      v(i) = 0.0;
  }
  static int tagOf(long id, int phase) { return static_cast<int>((id % 4096) * 4 + phase) + 4096; }
 
#ifdef PECLET_FLOW_MPI
  void requireComm() const {
    if (!hasComm_)
      throw std::runtime_error("peclet::flow::vof::VofBlockExchange: no communicator set");
  }
#endif
 
  long patchIdx(int gx, int gy, int gz) const {
    return L3(gx - patch_.o.x + patch_.g, gy - patch_.o.y + patch_.g, gz - patch_.o.z + patch_.g,
              patch_.e);
  }
 
  template <class HV>
  void packPatch(const VofPiece& p, const HV& hu, const HV& hv, const HV& hw,
                 std::vector<double>& buf, long& off) const {
    for (int z = p.g.lo[2]; z < p.g.hi[2]; ++z)
      for (int y = p.g.lo[1]; y < p.g.hi[1]; ++y)
        for (int x = p.g.lo[0]; x < p.g.hi[0]; ++x) {
          const long i = patchIdx(x, y, z);
          buf[3 * off + 0] = hu(i);
          buf[3 * off + 1] = hv(i);
          buf[3 * off + 2] = hw(i);
          ++off;
        }
  }
 
  template <class HV, class HB>
  void copyPatchToBlock(const VofPiece& p, const HV& hu, const HV& hv, const HV& hw, HB bu, HB bv,
                        HB bw, const VofBlock& b, int ghost) const {
    const I3 be = b.advector().extent();
    for (int z = p.g.lo[2]; z < p.g.hi[2]; ++z)
      for (int y = p.g.lo[1]; y < p.g.hi[1]; ++y)
        for (int x = p.g.lo[0]; x < p.g.hi[0]; ++x) {
          const long i = patchIdx(x, y, z);
          // block-local EXTENDED index: `loc` is the index inside the extended box, whose cell
          // (0,0,0) sits at the advector's (0,0,0) — i.e. `ghost` cells before the inner region.
          const int lx = p.loc[0] + (x - p.g.lo[0]);
          const int ly = p.loc[1] + (y - p.g.lo[1]);
          const int lz = p.loc[2] + (z - p.g.lo[2]);
          const long j = L3(lx, ly, lz, be);
          bu(j) = hu(i);
          bv(j) = hv(i);
          bw(j) = hw(i);
        }
    (void)ghost;
  }
 
  template <class HB>
  void unpackBlock(const VofPiece& p, const std::vector<double>& buf, long& off, HB bu, HB bv,
                   HB bw, const VofBlock& b, int ghost) const {
    const I3 be = b.advector().extent();
    for (int z = p.g.lo[2]; z < p.g.hi[2]; ++z)
      for (int y = p.g.lo[1]; y < p.g.hi[1]; ++y)
        for (int x = p.g.lo[0]; x < p.g.hi[0]; ++x) {
          const int lx = p.loc[0] + (x - p.g.lo[0]);
          const int ly = p.loc[1] + (y - p.g.lo[1]);
          const int lz = p.loc[2] + (z - p.g.lo[2]);
          const long j = L3(lx, ly, lz, be);
          bu(j) = buf[3 * off + 0];
          bv(j) = buf[3 * off + 1];
          bw(j) = buf[3 * off + 2];
          ++off;
        }
    (void)ghost;
  }
 
  template <class HB>
  void packBlockColour(const VofPiece& p, const HB& bc, const VofBlock& b, std::vector<double>& buf,
                       long& off) const {
    const I3 be = b.advector().extent();
    const int g = b.advector().ghost();
    for (int z = p.g.lo[2]; z < p.g.hi[2]; ++z)
      for (int y = p.g.lo[1]; y < p.g.hi[1]; ++y)
        for (int x = p.g.lo[0]; x < p.g.hi[0]; ++x) {
          const int lx = p.loc[0] + (x - p.g.lo[0]);
          const int ly = p.loc[1] + (y - p.g.lo[1]);
          const int lz = p.loc[2] + (z - p.g.lo[2]);
          buf[off++] = bc(L3(lx + g, ly + g, lz + g, be));
        }
  }
 
  template <class HB, class HC>
  void maxBlockIntoPatch(const VofPiece& p, const HB& bc, const VofBlock& b, HC& hc) const {
    const I3 be = b.advector().extent();
    const int g = b.advector().ghost();
    for (int z = p.g.lo[2]; z < p.g.hi[2]; ++z)
      for (int y = p.g.lo[1]; y < p.g.hi[1]; ++y)
        for (int x = p.g.lo[0]; x < p.g.hi[0]; ++x) {
          const int lx = p.loc[0] + (x - p.g.lo[0]);
          const int ly = p.loc[1] + (y - p.g.lo[1]);
          const int lz = p.loc[2] + (z - p.g.lo[2]);
          const double v = bc(L3(lx + g, ly + g, lz + g, be));
          double& d = hc(patchIdx(x, y, z));
          if (v > d)
            d = v;
        }
  }
 
  template <class HC>
  void maxBufIntoPatch(const VofPiece& p, const std::vector<double>& buf, long& off, HC& hc) const {
    for (int z = p.g.lo[2]; z < p.g.hi[2]; ++z)
      for (int y = p.g.lo[1]; y < p.g.hi[1]; ++y)
        for (int x = p.g.lo[0]; x < p.g.hi[0]; ++x) {
          const double v = buf[off++];
          double& d = hc(patchIdx(x, y, z));
          if (v > d)
            d = v;
        }
  }
 
  I3 gs_{0, 0, 0};
  std::array<bool, 3> per_{true, true, true};
  std::vector<VofBox> rankBox_;
  int rank_ = 0, size_ = 1;
  Patch patch_;
  SField lu_, lv_, lw_;
  mutable std::vector<VofPiece> pieces_;
  std::vector<SField> bufPool_;
  long bufSeq_ = 0;
  long gBytes_ = 0, sBytes_ = 0, gMsgs_ = 0, sMsgs_ = 0, mBytes_ = 0, mMsgs_ = 0;
#ifdef PECLET_FLOW_MPI
  MPI_Comm comm_ = MPI_COMM_NULL;
  bool hasComm_ = false;
#endif
};
 
}  // namespace peclet::flow::vof
 
#endif  // PECLET_FLOW_VOF_BLOCK_EXCHANGE_HPP