flow_ibm_mpi.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_FLOW_IBM_MPI_HPP
#define PECLET_FLOW_FLOW_IBM_MPI_HPP
 
namespace peclet::flow {
 
#ifdef PECLET_FLOW_MPI
template <class Grid>
void Solver<Grid>::initMpi(int gnx, int gny, int gnz, MPI_Comm comm) {
  int size = 1;
  MPI_Comm_size(comm, &size);
  weightedDec_ = false;  // the equal-weight partition; rebalanceByWeights sets it
  // Build the shared decomposition so the pressure MG can derive nested coarse levels
  // (CutcellMG::coarsened) — aligned ORB by default, or coarse-first when `set_decomposition`
  // asks for it. Depends only on the global grid and that setting, so it matches mpi_block()'s
  // sizing exactly when the same values are passed there.
  initMpi(CutcellMG::decomposition(static_cast<std::size_t>(size), gnx, gny, gnz, decompLevels_,
                                   decompMaxImbalance_),
          comm);
}
 
template <class Grid>
void Solver<Grid>::initMpi(const peclet::core::decomp::BlockDecomposer<3>& dec, MPI_Comm comm) {
  // The FIRST initMpi must precede the geometry. setSolid builds the pressure MG, the velocity MG,
  // the SDF ghost ring and the cut-cell overlays for the partition in force when it runs; before
  // initMpi that is the single-rank one, so each rank's hierarchy treats its own block as a
  // periodic domain of its own. Nothing below rebuilds it, and the step then runs with a
  // distributed momentum solve but P independent pressure solves: every PCG converges, and the
  // velocity leaves each rank carrying the divergence of its block-local mean (measured, one
  // step: a Gaussian cell force at np = 4, 2x2x1, max|du| 4.0e-2 of max|u| 0.36; a pure-gradient
  // force at np = 2, max|div u| 7.5e-2 where single-rank gives 2e-15). np = 1 is exact by
  // accident -- its block IS the domain -- so the order is enforced at every rank count.
  // redistribute() re-enters here with distributed_ already set and rebuilds the geometry itself.
  if (geometryBuilt_ && !distributed_)
    throw std::runtime_error(
        "init_mpi: call BEFORE the geometry (set_solid / set_pressure_geometry / "
        "set_solid_from_scene) -- the geometry builds the pressure and velocity multigrids for the "
        "partition in force when it runs, so a geometry built first gives every rank an "
        "independent, block-periodic pressure solve (a wrong answer at np > 1, with no error)");
  pcInDomain_ = CCField();  // WO-P3g: which ghosts carry a row depends on the decomposition
  distributed_ = true;
  comm_ = comm;
  const auto& gs = dec.globalSize();
  gnx_ = (int)gs[0];
  gny_ = (int)gs[1];
  gnz_ = (int)gs[2];
  int rank = 0;
  MPI_Comm_rank(comm, &rank);
  std::array<bool, 3> per{true, true, true};
  velHalo_ = std::make_shared<GridHaloTopology<3>>();
  velHalo_->buildTopology(dec, rank, G, per, comm);
  velDev_ = std::make_shared<GridHalo<double>>();
  velDev_->setLabel("velocity g2");
  velDev_->init(*velHalo_);
  // Communication-avoiding momentum sweeps (see smoothComp): the velocity block is g=2 already,
  // so the CA pair needs no new topology — only this float exchange for the stencil ring
  // (operator coefficients are float). Eligible when every rank's block is >= 4 on every axis
  // (rank-uniform: the decomposition is replicated), gated by `set_comm_avoiding`.
  velDevF_ = std::make_shared<GridHalo<MReal>>();
  velDevF_->init(*velHalo_);
  long minExt = std::numeric_limits<long>::max();
  for (const auto& s : dec.sizes())
    for (int k = 0; k < 3; ++k)
      minExt = std::min(minExt, (long)s[k]);
  caMomentum_ = (caMode_ & kCaMomentum) && minExt >= 4;
  for (bool& d : momStencilDirty_)
    d = true;
  dec_ = std::make_shared<peclet::core::decomp::BlockDecomposer<3>>(dec);  // remember the partition
  const auto oig = velHalo_->indexer().originInclGhost();
  og_ = {(int)oig[0] + G, (int)oig[1] + G, (int)oig[2] + G};  // block inner origin -> global parity
  // The colour field's OWN g=3 topology follows the same partition (VOF_PLAN §3 rule 1). Rebuilt
  // here rather than only in enableVof() so the order enable_vof/init_mpi does not matter — and
  // because a redistribute() re-inits through this path. The advector's block carries no state
  // between steps (C's canonical storage is the G=2 registry field "C"), so reallocating it is
  // free of consequence.
  if (vofEnabled_)
    buildVofBlock();
}
 
template <class Grid>
void Solver<Grid>::redistribute(const peclet::core::decomp::BlockDecomposer<3>& newDec) {
  if (!distributed_ || !dec_)
    return;
  int rank = 0;
  MPI_Comm_rank(comm_, &rank);
  const auto ob = dec_->block(rank), nb = newDec.block(rank);
  const int oex = (int)ob.size[0] + 2 * G, oey = (int)ob.size[1] + 2 * G,
            oez = (int)ob.size[2] + 2 * G;
  const int nex = (int)nb.size[0] + 2 * G, ney = (int)nb.size[1] + 2 * G,
            nez = (int)nb.size[2] + 2 * G;
 
  // 1. gather the surviving registered fields to host padded buffers on the OLD block.
  const auto names = fields_.names();
  std::vector<std::vector<double>> oldHost(names.size()), newHost(names.size());
  for (std::size_t k = 0; k < names.size(); ++k) {
    CCField f = fields_.at(names[k]).data;
    auto h = Kokkos::create_mirror_view(f);
    Kokkos::deep_copy(h, f);
    oldHost[k].assign(h.data(), h.data() + (std::size_t)oex * oey * oez);
    newHost[k].assign((std::size_t)nex * ney * nez, 0.0);
  }
  // 2. redistribute each field OLD -> NEW (host, bit-exact pure data movement).
  std::vector<const double*> op(names.size());
  std::vector<double*> np(names.size());
  for (std::size_t k = 0; k < names.size(); ++k) {
    op[k] = oldHost[k].data();
    np[k] = newHost[k].data();
  }
  // 2b. The Chebyshev driver's kept power iterates (D3, doc/vof_step_performance_design.md §5.13)
  //    are cross-step state of the pressure MG, held on its g = 1 level-0 block: carried along
  //    in the same exchange (inner cells re-padded to G), handed back after the MG rebuild (5b).
  //    Dropping them would make the next re-estimate cold, and an np = 1 rebalance would stop
  //    being bit-exact. Collective decision: every rank carries them or none does.
  int warmLocal = mg_.eigenWarmReady() ? 1 : 0, warmAll = 0;
  MPI_Allreduce(&warmLocal, &warmAll, 1, MPI_INT, MPI_MIN, comm_);
  std::vector<double> eigOld[2], eigNew[2];
  if (warmAll) {
    const CCField src[2] = {mg_.eigenWarmMax(), mg_.eigenWarmMin()};
    const int o1x = (int)ob.size[0] + 2, o1y = (int)ob.size[1] + 2;
    for (int v = 0; v < 2; ++v) {
      auto h = Kokkos::create_mirror_view(src[v]);
      Kokkos::deep_copy(h, src[v]);
      eigOld[v].assign((std::size_t)oex * oey * oez, 0.0);
      eigNew[v].assign((std::size_t)nex * ney * nez, 0.0);
      for (int z = 0; z < (int)ob.size[2]; ++z)
        for (int y = 0; y < (int)ob.size[1]; ++y)
          for (int x = 0; x < (int)ob.size[0]; ++x)
            eigOld[v][(std::size_t)(x + G) + (std::size_t)(y + G) * oex +
                      (std::size_t)(z + G) * oex * oey] =
                h((std::size_t)(x + 1) + (std::size_t)(y + 1) * o1x +
                  (std::size_t)(z + 1) * o1x * o1y);
      op.push_back(eigOld[v].data());
      np.push_back(eigNew[v].data());
    }
  }
  peclet::core::decomp::redistributeGridFields<double>(*dec_, newDec, rank, G, op, np, comm_);
 
  // 3. reallocate every buffer to the new block; re-init the halo + MG on the new partition.
  //    `resizeForBlock` must run BETWEEN the two: `allocateBlock` only re-creates the buffers it
  //    owns by name, and `initMpi` -> `buildVofBlock` already writes through `cField_`, which is
  //    a registry alias. Without it the scatter below memcpy's the new padded extent into the
  //    old allocation (ASan heap-buffer-overflow -> `free(): invalid pointer`).
  allocateBlock((int)nb.size[0], (int)nb.size[1], (int)nb.size[2]);
  resizeForBlock();
  initMpi(newDec, comm_);
  // scatter a padded host buffer into a registered field's device buffer.
  auto scatterPadded = [&](const std::string& name, const std::vector<double>& src) {
    CCField f = fields_.at(name).data;
    auto h = Kokkos::create_mirror_view(f);
    std::memcpy(h.data(), src.data(), sizeof(double) * (std::size_t)nex * ney * nez);
    Kokkos::deep_copy(f, h);
  };
  // 4. scatter all migrated fields into the fresh (new-block) buffers.
  for (std::size_t k = 0; k < names.size(); ++k)
    scatterPadded(names[k], newHost[k]);
  // 5. rebuild geometry-derived state (openness/IBM/stencils/MG) from the migrated SDF. setSolid
  //    zeroes the velocity + pressure (it is the initial-geometry setup), so re-instate every
  //    non-SDF field afterward from the migrated data.
  {
    // A repartition is not an operator change: the balanced-force pre-projection keeps its
    // Chebyshev bounds across it, as the main solve keeps its own (setSolid's MG rebuild would
    // otherwise force a re-estimate, and an np = 1 rebalance would stop being bit-exact).
    const bool keepBfpBounds = bfpChebSet_;
    setSolid(gatherInner(sdf_), cutcellPressure_);
    bfpChebSet_ = keepBfpBounds;
  }
  for (std::size_t k = 0; k < names.size(); ++k)
    if (names[k] != "sdf")
      scatterPadded(names[k], newHost[k]);
  // 5b. hand the migrated Chebyshev power iterates back to the rebuilt MG (see 2b): inner cells
  //    into a fresh g = 1 level-0 vector (ghosts zero; no consumer reads them before a refill).
  if (warmAll) {
    const int n1x = (int)nb.size[0] + 2, n1y = (int)nb.size[1] + 2, n1z = (int)nb.size[2] + 2;
    CCField dst[2];
    for (int v = 0; v < 2; ++v) {
      dst[v] = CCField(v == 0 ? "ev_vmax" : "ev_vmin", (std::size_t)n1x * n1y * n1z);
      auto h = Kokkos::create_mirror_view(dst[v]);
      Kokkos::deep_copy(h, 0.0);
      for (int z = 0; z < (int)nb.size[2]; ++z)
        for (int y = 0; y < (int)nb.size[1]; ++y)
          for (int x = 0; x < (int)nb.size[0]; ++x)
            h((std::size_t)(x + 1) + (std::size_t)(y + 1) * n1x +
              (std::size_t)(z + 1) * n1x * n1y) =
                eigNew[v][(std::size_t)(x + G) + (std::size_t)(y + G) * nex +
                          (std::size_t)(z + G) * nex * ney];
      Kokkos::deep_copy(dst[v], h);
    }
    mg_.setEigenWarm(dst[0], dst[1]);
  }
  // 6. REFILL THE GHOSTS. `redistributeGridFields` moves the INNER cells only and says so
  //    ("the caller refills ghosts with a halo exchange afterwards"); the padded destination
  //    buffers were zero-initialised in step 1, so without this every migrated field enters the
  //    next step with a zero halo on the new partition. Measured at np = 4 on the VoF gate: one
  //    step after the move, du = 1.01e-01 against the single-rank reference (a run started
  //    directly on the SAME weighted partition reproduces it to 6.5e-17), and 0.00e+00 with the
  //    refill. np = 1/2 hid it because the new neighbour set happened to reproduce the old one.
  for (std::size_t k = 0; k < names.size(); ++k)
    fillGhosts(fields_.at(names[k]).data);
  for (auto& sc : scalars_)
    applyScalarBc(sc);  // a scalar's own domain BCs override the halo/periodic base
  // 6b. RE-SEED eps^n. It is STATE, not scratch: a re-partition is not a time step, so the
  //    porosity must not appear to have jumped across it -- eps^n := the migrated eps^{n+1}, i.e.
  //    d(eps)/dt = 0 over the redistribute, the only choice that leaves the projection RHS
  //    unchanged. It must happen HERE, after step 4 scattered the migrated eps into `epsField_`:
  //    it used to happen in `resizeForBlock` (step 3), where `epsField_` is the freshly allocated
  //    ZERO buffer whenever the block changes size, so eps^n = 0 and the first porous projection
  //    after a rebalance saw d(eps)/dt = eps/dt. Measured on coupling's test_mpi_rebalance (a
  //    32^3 heap, np = 2 and 4): pressure off by 2.2e+02 and velocity by 2.4 (|u| <= 0.07) one
  //    step after the move; a rebalance that keeps every block's size hid it (no reallocation).
  if (epsPrev_.extent(0) == n_ && epsField_.extent(0) == n_)
    Kokkos::deep_copy(epsPrev_, epsField_);
  // 7. RE-SEED THE COLLOCATED FACE FIELD. `uf_/vf_/wf_` are block scratch, not registry fields:
  //    step 3's `allocateBlock` handed us fresh ZERO buffers and nothing above migrates them, so
  //    the projected face field of the old partition is simply gone. Everything that rides it --
  //    momentum advection (ufAdvVelocity()), the colour transport, the scalars -- would read
  //    zeros on the step after a rebalance, and `advect_vof`'s divergence guard would wave them
  //    through (a zero field is perfectly solenoidal). Seeding from the migrated cell velocity
  //    with the same `centerToFace` map `project()` uses restores a consistent (not yet
  //    divergence-free) field, exactly as `set_state` does; the next `project()` makes it
  //    solenoidal again. A no-op on the staggered grid.
  seedFaceFieldFromCells();
}
 
template <class Grid>
int Solver<Grid>::rebalanceByWeights(const std::vector<peclet::core::Real>& w) {
  if (!distributed_)
    return 1;
  int size = 1, rank = 0;
  MPI_Comm_size(comm_, &size);
  MPI_Comm_rank(comm_, &rank);
  // The ALIGNED weighted ORB (amr/docs/amr_mg_core_boundary.md §6, §9.2, §11.4): split planes on
  // multiples of 2^a, a the largest alignment whose weight imbalance stays within the 1.05 budget
  // (a = 0: the plain weighted ORB, bit for bit). Every split is then even for a lifts, so the
  // pressure MG coarsens in place for a levels before any telescope stage can fire, and whatever
  // stage fires does so on a level 8^a smaller. A pure function of (ranks, grid, weights): every
  // rank computes the same partition, and so must any code that co-decomposes from `w` -- which is
  // why the chosen alignment is RETURNED: dem rebuilds this partition as the aligned weighted ORB
  // `init(size, G, w, {2^a,...})`, and at a = 0 as the plain weighted ORB (core guarantees the
  // two coincide bit for bit at align 1).
  const auto c = peclet::core::decomp::chooseAlignedWeighted(
      (std::size_t)size, peclet::core::IVec<3>{gnx_, gny_, gnz_}, w);
  if (mgDebugLevel() && rank == 0) {
    std::printf(
        "[mg] rebalanceByWeights: aligned weighted ORB a = %d (align %d), weight imbalance "
        "%.4f\n",
        c.a, 1 << c.a, (double)c.imbalance);
    std::fflush(stdout);
  }
  weightedDec_ = true;  // before redistribute: its setSolid rebuilds the pressure MG
  redistribute(c.dec);
  return 1 << c.a;
}
 
template <class Grid>
std::vector<CutcellMG::PlanRow> Solver<Grid>::pressureMgPlan() const {
  return mg_.builtPlan();
}
#endif
 
template <class Grid>
void Solver<Grid>::resizeForBlock() {
  reallocOwnedFields();
  rebindFieldAliases();
  resizeBlockScratch();
  for (int f = 0; f < 6; ++f)
    resampleBcProfile(f);  // the resampled inlet profile is indexed by LOCAL face position
}
 
template <class Grid>
void Solver<Grid>::reallocOwnedFields() {
  for (const auto& nm : fields_.names()) {
    auto& rec = fields_.at(nm);
    if (!rec.ownStorage || rec.data.extent(0) == n_)
      continue;
    rec.data = CCField(nm, n_);
  }
}
 
template <class Grid>
void Solver<Grid>::rebindFieldAliases() {
  auto bind = [this](CCField& f, const std::string& nm) {
    if (fields_.has(nm))
      f = fields_.at(nm).data;
  };
  bind(cField_, "C");          // VoF: the canonical G=2 colour mirror
  bind(kappaField_, "kappa");  // VoF: the curvature mirrors
  bind(kappaBranch_, "kappa_branch");
  bind(rhoField_, "rho");   // property closures
  bind(Pb_, "p_balanced");  // balanced-force projection split (doc/collocated_varrho_forces.md)
  bind(muField_, "mu");
  bind(epsField_, "eps");  // CFD-DEM: porosity + drag
  bind(dragBeta_, "drag_beta");
  bind(pcMdot_, "mdot");  // phase change
  bind(pcSrc_, "pc_source");
  bind(pcUser_, "div_source");
  bind(pcSink_, "div_sink");
  static const char* fn[3] = {"force_x", "force_y", "force_z"};
  for (int c = 0; c < 3; ++c)
    bind(cellForce_[c], fn[c]);
  for (auto& sc : scalars_) {
    bind(sc.c, sc.name);
    if (sc.cut && sc.cut->solid.extent(0) > 0)  // WO-7: a conjugate scalar's solid field
      bind(sc.cut->solid, sc.name + "_solid");
  }
  for (auto& cl : closures_) {  // property_closures.hpp keeps the registry keys for exactly this
    if (fields_.has(cl.outName))
      cl.out = fields_.at(cl.outName).data;
    if (fields_.has(cl.in0Name))
      cl.in0 = CCConst(fields_.at(cl.in0Name).data);
    if (!cl.in1Name.empty() && fields_.has(cl.in1Name))
      cl.in1 = CCConst(fields_.at(cl.in1Name).data);
  }
}
 
template <class Grid>
void Solver<Grid>::resizeIfAllocated(CCField& f, const char* label, std::size_t n) {
  if (f.extent(0) != 0 && f.extent(0) != n)
    f = CCField(label, n);
}
 
template <class Grid>
void Solver<Grid>::resizeBlockScratch() {
  const std::size_t n = n_, n1 = n1_;
  // --- velocity / pressure scratch --------------------------------------------------------
  resizeIfAllocated(velRes_, "velRes", n);
  resizeIfAllocated(zp1_, "zp1", n1);
  resizeIfAllocated(vmgTheta_, "vmgTheta", n);
  resizeIfAllocated(vmgClean_, "vmgClean", n);
  for (int c = 0; c < 3; ++c) {
    resizeIfAllocated(uStar_[c], "uStar", n);
    resizeIfAllocated(advRhs_[c], "advRhs", n);
    resizeIfAllocated(uBc_[c], "uBc", n);
    resizeIfAllocated(uwCell_[c], "uwCell", n);
    resizeIfAllocated(uwAdv_[c], "uwAdv", n);
    resizeIfAllocated(faceAcc_[c], "faceAcc", n);
  }
  resizeIfAllocated(pb1_, "pb1",
                    n1);  // balanced-force projection scratch (P_b itself is registered)
  haveUStar_ = false;     // u* belonged to the old block's last momentum solve
  haveAdvRhs_ = false;
  // --- ghost projection (its overlay is rebuilt by setSolid) ------------------------------
  resizeIfAllocated(oxb_, "oxb", n);
  resizeIfAllocated(oyb_, "oyb", n);
  resizeIfAllocated(ozb_, "ozb", n);
  resizeIfAllocated(sdfGp_, "peclet::flow::sdfGp", n);
  resizeIfAllocated(gpRh_, "gpRh", n1);
  resizeIfAllocated(gpT_, "gpT", n1);
  resizeIfAllocated(gpZ2_, "gpZ2", n1);
  resizeIfAllocated(gpX2_, "gpXg2", n);
  // --- VoF (the g=3 block itself is rebuilt by buildVofBlock) -----------------------------
  resizeIfAllocated(vofCs_, "vofCs", n);
  resizeIfAllocated(vofSolidG2_, "vofSolidG2", n);
  for (int c = 0; c < 3; ++c) {
    resizeIfAllocated(uAdv_[c], "uAdv", n);
    resizeIfAllocated(vofDynVel_[c], "vofDynVel", n);
    resizeIfAllocated(csfBlkF_[c], "vof::blockcsf", n);
  }
  // --- phase change -----------------------------------------------------------------------
  resizeIfAllocated(pcArea_, "pc_area", n);
  for (int c = 0; c < 3; ++c) {
    resizeIfAllocated(pcNrm_[c], "pc_nrm", n);
    resizeIfAllocated(pcGn_[c], "pc_gn", n);
  }
  resizeIfAllocated(pcDep_, "pc_dep", n);
  resizeIfAllocated(pcTgt_, "pc_tgt", n);
  resizeIfAllocated(pcCnew_, "pc_cnew", n);
  resizeIfAllocated(pcDefic_, "pc_defic", n);
  resizeIfAllocated(pcKcell_, "pc_k", n);
  resizeIfAllocated(pcRcp_, "pc_rcp", n);
  resizeIfAllocated(pcClsPrev_, "pc_cls_prev", n);
  resizeIfAllocated(pcCarrySrc_, "pc_carry_src", n);
  resizeIfAllocated(pcMdotFit_, "pc_mdot_fit", n);
  resizeIfAllocated(pcKappa_, "pc_kappa", n);
  resizeIfAllocated(pcAreaCg2_, "pc_area_cascade", n);
  resizeIfAllocated(pcTgam_, "pc_tgam", n);
  resizeIfAllocated(pcGphi_, "pc_gphi", n);
  pcInDomain_ = CCField();  // which ghosts carry a row depends on the decomposition (WO-P3g)
  pcMaskFresh_ = false;
  // --- transported scalars (their operator is rebuilt; the solution rides the registry) -----
  for (auto& sc : scalars_) {
    resizeIfAllocated(sc.cOld, "scalar_old", n);
    resizeIfAllocated(sc.b, "scalar_b", n);
    resizeIfAllocated(sc.AC, "scalar_AC", n);
    resizeIfAllocated(sc.AW, "scalar_AW", n);
    resizeIfAllocated(sc.AE, "scalar_AE", n);
    resizeIfAllocated(sc.AS, "scalar_AS", n);
    resizeIfAllocated(sc.AN, "scalar_AN", n);
    resizeIfAllocated(sc.AB, "scalar_AB", n);
    resizeIfAllocated(sc.AT, "scalar_AT", n);
    resizeIfAllocated(sc.dmask, "scalar_dmask", n);
    resizeIfAllocated(sc.dval, "scalar_dval", n);
    resizeIfAllocated(sc.gfmB, "scalar_gfmb", n);
    sc.stencilBuilt = false;  // the operator is a property of the block
  }
  // `ScalarField::kcell`/`rcp` are NOT the scalar's own storage: they ALIAS `pcKcell_`/`pcRcp_`
  // (set at `set_phase_change_energy`). Resizing them independently would give the energy
  // operator a different buffer from the one `pcUpdateEnergyProps` fills every step — measured
  // as dP 3.09e-04 against the never-rebalanced control at np = 4. Re-alias instead.
  if (pcEnergy_ && !pcTName_.empty() && hasScalar(pcTName_)) {
    ScalarField& sc = scalarField(pcTName_);
    sc.kcell = pcKcell_;
    sc.rcp = pcRcp_;
  }
  // --- variable density / porous (CFD-DEM) -------------------------------------------------
  resizeIfAllocated(rho1_, "rho1", n1);
  resizeIfAllocated(cx1_, "cx1", n1);
  resizeIfAllocated(cy1_, "cy1", n1);
  resizeIfAllocated(cz1_, "cz1", n1);
  resizeIfAllocated(eps1_, "eps1", n1);
  resizeIfAllocated(beta1_, "beta1", n1);
  resizeIfAllocated(depsdt_, "depsdt", n);
  resizeIfAllocated(divAdv_, "divAdv", n);
  resizeIfAllocated(epsRho_, "epsRho", n);
  resizeIfAllocated(epsPrev_, "epsPrev", n);  // eps^n is STATE: `redistribute` re-seeds it (6b)
  for (int c = 0; c < 3; ++c)
    for (int k = 0; k < 3; ++k)
      if (tEx_[c][k].extent(0) != 0 && tEx_[c][k].extent(0) != (std::size_t)nx_ * ny_ * nz_) {
        tEx_[c][k] = CCField();  // exact crossings are single-rank only; drop the stale set
        hasExactCross_ = false;
      }
  // `cutOwner_` is deliberately NOT resized: every consumer guards on
  // `extent(0) == nx_*ny_*nz_` and skips, so a stale-sized owner map is inert, while a
  // freshly-zeroed one would silently name instance 0 as the owner of every cut cell.
}
 
}  // namespace peclet::flow
 
#endif  // PECLET_FLOW_FLOW_IBM_MPI_HPP