flow 1.3.0
Kokkos cut-cell IBM incompressible Navier-Stokes solver + pnm pore extraction
Loading...
Searching...
No Matches
peclet::flow::vof Namespace Reference

Namespaces

namespace  bcdetail
 

Classes

struct  McVertex
 One dual-cube corner = one cell centre. More...
 
class  MomentumConsistentAdvector
 Momentum-consistent companion to WyAdvector. More...
 
struct  PcGfmRow
 WO-P3g item 2 — the Gibou–Fedkiw second-order one-sided (ghost-fluid) row. More...
 
struct  PcGradFit
 Accumulator for dT/dn fitted through the interface value: T - T_G = G * phi, so G = sum_j w_j phi_j (T_j - T_G) / sum_j w_j phi_j^2. More...
 
class  VofBlock
 One bubble's block. More...
 
class  VofBlockExchange
 The gather / scatter / table replication of the block container. More...
 
struct  VofBlockExchangeBase
 Gather / scatter / table replication. More...
 
struct  VofBlockJob
 One master block's share of a batched container launch. More...
 
class  VofBlockSet
 The container: a replicated table of blocks, the per-step orchestration (gather → advect → re-centre → scatter/union) and the per-bubble statistics. More...
 
struct  VofBlockStats
 Per-block Lagrangian census, computed on the master (VOF_PLAN §10 W1 outputs, available at W0). More...
 
struct  VofBlockTable
 The job table of one launch, by value. More...
 
struct  VofBox
 Half-open global cell box [lo, hi). More...
 
struct  VofCountTable
 The exact offset table of a host list launch: job k's entries are [o[k], o[k + 1]). More...
 
class  VofCurvature
 Curvature of the colour field on an extended (inner + ghost) block. More...
 
struct  VofCurvFallbackJob
 One block's share of a batched tier-3 launch (VofCurvature::fallbackBatch): its fields, its compacted interfacial list, and every argument its own fallbackPass would pass to curvFallbackCell. More...
 
struct  VofCurvFallbackTable
 The job table of one batched launch, passed to the kernel BY VALUE (namespace scope: nvcc's extended lambdas may not capture a function-local type). More...
 
struct  VofCurvJob
 One master block's cascade (VofCurvature, worklist mode) for the batched launches: its planes, curvature, branch and force fields, and every argument the per-cascade passes give the shared cell bodies (wyReconstructCell, curvHeightCell, curvFallbackCell, the clip, the census, csfFaceCurvature + csfFaceForce). More...
 
struct  VofCurvTable
 
struct  VofDebrisFlags
 Flags of the pass, uniform over the blocks. More...
 
struct  VofDynamicWetting
 The driver: produces WyAdvector::contactAngle() from the current colour, geometry, wall SDF, fluid-only normals and cell-centre velocity, once per band fill. More...
 
class  VofEnergyAdvector
 
class  VofInterfaceArea
 Interfacial area per cell from the V3 curvature cascade's own geometry. More...
 
class  VofMcArea
 Interfacial area per cell from the marching-tetrahedra surface of the C = 1/2 level set (or of the PLIC-reconstructed signed distance — set_phase_change_area modes 4-7). More...
 
struct  VofMoveJob
 One (piece, component) copy of VofBlockExchange::movePiece – patch cells [s, s + n) to block cells [q, q + n) – for the host's one-launch move (H-4e). More...
 
struct  VofMoveTable
 The job table of one move launch, by value. More...
 
struct  VofPiece
 One contiguous transfer: the global cells g of a block that rank r owns, and where they sit in the block's EXTENDED local index space. More...
 
struct  VofPvCached
 The cache-backed accessor: slot[j] >= 0 iff cell j of the grown region is interfacial, and cache[slot[j]] is its entry (written by the planes pass of the same cascade). More...
 
struct  VofPvOnTheFly
 The on-the-fly accessor of the dense mode (useWorklist = false, no list and no slot): the interfacial test is vofIsInterface(c(j)), and the entry is built from the planes. More...
 
struct  VofPvTeamScratch
 
struct  VofRawField
 A raw-pointer field accessor with the View<double*> call syntax, for kernels that walk SEVERAL blocks' fields in one launch (a device-side array of Views would carry reference counts). More...
 
struct  VofRawMask
 Read-only raw accessor of an unsigned-char field (the WO-R out-of-domain mask). More...
 
struct  VofRowTable
 The row offsets of a host region launch: job k's rows are [off[k], off[k + 1]). More...
 
struct  VofRun
 One contiguous run of GLOBAL indices covered by a block's index range on one axis, with the block-local index the run starts at. More...
 
struct  VofStatsPrev
 What measure() reads besides the colour: each job's previous centroid and whether it is valid (the host's prevCentroid_ / hasPrev_, current when the launch is issued). More...
 
struct  VofSumN
 N double sums in one team reduction (Kokkos::Sum<VofSumN<N>>). More...
 
class  WyAdvector
 Weymouth-Yue split advection of a colour field on an extended (inner + ghost) block. More...
 

Typedefs

using UCField = Kokkos::View< unsigned char *, SMem >
 
using LField = Kokkos::View< long *, SMem >
 
using VofPvTerms = Kokkos::View< PvTerm *, SExec::scratch_memory_space, Kokkos::MemoryUnmanaged >
 Team scratch of the device tier-3 body: the compacted stencil's canonical offset indices, their terms, and the 27 accumulated entries (about 10 KB at the 5^3 stencil).
 
using VofPvInts = Kokkos::View< int *, SExec::scratch_memory_space, Kokkos::MemoryUnmanaged >
 
using VofPvDbls = Kokkos::View< double *, SExec::scratch_memory_space, Kokkos::MemoryUnmanaged >
 

Enumerations

enum  VofDebrisSlot : int {
  kDbSumD = 0 , kDbSumR , kDbW , kDbVA ,
  kDbDV , kDbActed , kDbND , kDbNR ,
  kDbPending , kDbSlots
}
 The per-job slots of the debris packet (out(stride * (base + k) + slot)). More...
 
enum class  VofMasterAssign { RoundRobin = 0 , Lpt = 1 , WeightedOrb = 2 }
 How a block's MASTER rank is chosen. More...
 
enum  InterfaceAreaMode : int {
  kAreaPlic = 0 , kAreaMetric = 1 , kAreaNormal = 2 , kAreaFootprint = 3 ,
  kAreaMcColour = 4 , kAreaMcColourSplit = 5 , kAreaMcPlic = 6 , kAreaMcPlicSplit = 7
}
 Which construction turns the cascade's normal into an area. More...
 
enum  McSource : int { kMcSrcColour = 0 , kMcSrcPlic = 1 }
 What the edge crossing interpolates. More...
 
enum  McDeposit : int { kMcDepositCentroid = 0 , kMcDepositSplit = 1 }
 Which cell a triangle's area is booked to. More...
 
enum  VofDynamicState : int {
  kVofDynNone = 0 , kVofDynStatic = 1 , kVofDynPinned = 2 , kVofDynAdvancing = 3 ,
  kVofDynReceding = 4 , kVofDynStateCount = 5
}
 Which rule produced a band cell's imposed angle (vof_dynamic_field(4)). More...
 

Functions

bool wyIsMixed (double c)
 The mixed-cell predicate.
 
bool wyIsMixed (double c, double eps)
 WO-R2 item 4 — the same predicate with a WISP TOLERANCE.
 
bool wyColourJump (double a, double b, double eps)
 The interface-BAND jump test of maxCourantInterface: two neighbouring cells carry a DIFFERENT colour.
 
template<class SF >
void wyReconstructCell (const SF &c, long i, long sy, long sz, SF mx, SF my, SF mz, SF alpha)
 PLIC reconstruction of one cell: MYC normal from the 3^3 colour stencil + the analytic plane offset.
 
template<class SF >
double wyFaceFlux (double a, long p, long sd, int dir, const SF &c, const SF &mx, const SF &my, const SF &mz, const SF &alpha, double eps=0.0)
 Signed Eulerian donor-cell flux through the dir-face between cell p and cell p + sd, as a fraction of a cell volume, positive along +dir.
 
template<class SF , class MF >
double wyFaceFluxBc (double a, long p, long sd, int dir, const SF &c, const SF &mx, const SF &my, const SF &mz, const SF &alpha, const MF &outside, double eps=0.0)
 WO-R sibling of wyFaceFlux for a face whose DONOR lies outside the global domain.
 
int vofJobOf (const long *off, long t)
 Job of flat index t (off[0] <= t < off[nj]): a walk over at most 16 offsets.
 
template<class F >
long vofSetOffsets (VofBlockTable &T, F size)
 Fill T.off from a per-job range size.
 
long vofJobLen (const VofBlockJob &J)
 
long vofJobInner (const VofBlockJob &J)
 
long vofJobGrown1 (const VofBlockJob &J)
 The inner region grown by one cell per side: WyAdvector::reconstructImpl's range.
 
template<class Exec , class Cnt , class F >
void vofHostListFor (const char *name, int nj, Cnt cnt, F body)
 H-4(a): a host list kernel over the EXACT entry counts, cnt(k) (read directly: the count Views are HostSpace), with a dynamic schedule in chunks of 16 – so the active entries spread over every thread instead of sitting at the head of one block's static share.
 
template<class Tab >
auto vofListCount (const Tab &T, LField start, LField end)
 The [start, end) count of job k of a list launch.
 
template<class Exec , class Rows , class F >
void vofHostRowFor (const char *name, int nj, Rows rows, F row)
 H-4(b): a host region kernel over ROWS – rows(k) x-runs for job k, row(k, r) handling job k's local row r with the x loop inside – so the job search and the index math run once per row, not once per cell.
 
template<class Exec , class Rows , class RowOf , class Hit >
void vofHostRowScan (const char *name, int nj, int base, Rows rows, RowOf rowOf, Hit hit, LField list, LField start, LField end)
 H-4(d): a host compaction scan over ROWS.
 
void vofBatchCfl (const VofBlockTable &T, SField out)
 One team per job: WyAdvector::maxCourant or maxCourantInterface (per the job's interfaceCfl), the same predicates and the same fmax chain.
 
template<class Exec = SExec>
void vofBatchFreeze (const VofBlockTable &T)
 WyAdvector::freezeDilationFlag over every job's extended block.
 
template<class Exec = SExec>
void vofBatchWorklist (const VofBlockTable &T, LField list, LField start, LField end)
 WyAdvector::reconstructImpl's worklist scan, over the concatenated grown-by-one regions (T.off = vofJobGrown1).
 
template<class Exec = SExec>
void vofBatchPlic (const VofBlockTable &T, LField list, LField start, LField end)
 wyReconstructCell over each job's compacted list: launched over the upper bound (the job's region), a thread past its job's count exits.
 
long vofJobFaces (const VofBlockJob &J, int d)
 A sweep's face range of one job: the inner region with one extra layer on the low side of d.
 
void vofBatchFluxCell (const VofBlockJob &J, int d, long p)
 The flux of d-face p of job J (the per-cell body of vofBatchFlux, both launch forms).
 
template<class Exec = SExec>
void vofBatchFlux (const VofBlockTable &T, int d)
 WyAdvector::computeFluxesImpl (uncut): one flux per d-face, wyFaceFlux (or, for a job carrying the WO-R mask, wyFaceFluxBc).
 
void vofBatchBcLedger (const VofBlockTable &T, int d, SField out, int stride, int slot0)
 The WO-R boundary ledger (WyAdvector::accumulateBcFaceVolume) of the jobs carrying the mask: one thread per (job, side), summing its plane of d-fluxes in (p0 fastest, p1) row order into out(stride (base + k) + slot0 + face) – a fixed order, which is a C2 order change against the per-block MDRange2 reduction (recorded; doc/vof_step_performance_design.md §5.9).
 
void vofBatchUpdateCell (const VofBlockJob &J, int d, long i)
 The sweep update of inner cell i of job J (the per-cell body of vofBatchUpdate).
 
template<class Exec = SExec>
void vofBatchUpdate (const VofBlockTable &T, int d)
 WyAdvector::applySweepImpl (uncut) over the inner regions (T.off = vofJobInner); the expression verbatim.
 
template<class Exec = SExec>
void vofBatchGhostZero (const VofBlockTable &T)
 Pass 1: every ghost cell of every job -> 0.
 
template<class Exec = SExec>
void vofBatchPeriodic (const VofBlockTable &T, int a)
 Pass 2 on axis a (colour_field.hpp::periodicFill, one axis): only the jobs whose block spans the periodic axis a are in T; T.off = the transverse plane sizes dims[b] * dims[c].
 
void vofBatchClampCell (const VofBlockJob &J, int x, int y, int z, int q0, int q1, int q2, bool p0, bool p1, bool p2)
 The clamp of extended cell (x, y, z) of job J (the per-cell body of vofBatchClamp).
 
template<class Exec = SExec>
void vofBatchClamp (const VofBlockTable &T, I3 gs, bool px, bool py, bool pz)
 Pass 3 (colour_field.hpp::clampFill, skip = 0): the part of each extended block outside a NON-periodic domain takes the globally clamped value.
 
template<class SF >
bool vofDebrisCell (const SF &c, long i, long sy, long sz, double ieps, double cfull)
 The debris predicate of doc/vof_overlap_design.md §11 on cell i: interfacial, and NO cell of the marker at or above cfull within two cells (the 5^3 fit stencil).
 
template<class SF >
bool vofAttachedCell (const SF &c, long i, long sy, long sz, double ieps, double cfull)
 The ATTACHED predicate (§5.3 step 2): interfacial with a cell at or above cfull within two.
 
bool vofResidueCell (double ci, double weps)
 The residue predicate (§13): 0 < |C| <= wispEps.
 
template<class Exec = SExec>
void vofBatchMark (const VofBlockTable &T, int kind, double ieps, double cfull, double weps, LField list, LField start, LField end)
 One mark scan over the concatenated inner regions (T.off = vofJobInner), in each block's inner-box index order – the per-block scan's order, so every list is the per-block list.
 
void vofBatchDebrisSums (const VofBlockTable &T, VofDebrisFlags F, LField listD, LField listR, LField listA, LField sD0, LField eD0, LField sR0, LField eR0, LField sA0, LField eA0, SField out, int stride, int slot0)
 Step 2 of the pass: one thread per job, every sum in list order (the per-block kernel's single thread and order, so the same bits).
 
long vofDebrisActCount (const VofDebrisFlags &F, LField sA0, LField eA0, SField out, long ob, int j)
 The number of act entries of job j (its D removed, then R, then A), 0 when it did not act: the per-block write kernel's range.
 
void vofDebrisActEntry (const VofBlockTable &T, const VofDebrisFlags &F, int k, long t, LField listD, LField listR, LField listA, LField sD0, LField sR0, LField sA0, LField eA0, SField out, int stride, int slot0, SField ex, Kokkos::View< int *, SMem > capf)
 Act entry t of job k (the per-block write kernel's expressions verbatim); t past the job's count, or a job that did not act, is a no-op.
 
template<class Exec = SExec>
void vofBatchDebrisWrite (const VofBlockTable &T, VofDebrisFlags F, LField listD, LField listR, LField listA, LField sD0, LField sR0, LField sA0, LField eA0, SField out, int stride, int slot0, SField ex, Kokkos::View< int *, SMem > capf)
 Step 4, the act: removed cells -> 0, attached C += dV C(1-C)/W capped at 1 (the per-block write kernel's expressions verbatim).
 
void vofBatchDebrisLost (const VofBlockTable &T, LField sA0, LField eA0, SField out, int stride, int slot0, SField ex, Kokkos::View< int *, SMem > capf)
 lost: the capped excess of an acted job summed in list order on one thread (only where a cap fired, as the per-block path; 0 otherwise).
 
template<class Exec = SExec>
void vofBatchBox (const VofBlockTable &T, double eps, SField out, int stride, int slot0)
 One team per job: the tight box of |C| > eps over the inner region, in the block's LOCAL inner frame (bubbleBox's six Min/Max reductions; order-free).
 
template<class F >
long vofSetOffsets (VofCurvTable &T, F size)
 
long vofCurvInner (const VofCurvJob &J)
 
long vofCurvGrown (const VofCurvJob &J)
 The inner region grown by kPvHalf: every cell VofCurvature::reconstructPlanes writes.
 
long vofCurvCell (const VofCurvJob &J, long r, bool grown)
 Cell r of job J's grown (grown) or inner region, in the region's x-fastest order.
 
long vofCurvRows (const VofCurvJob &J, bool grown)
 Rows of job J's grown (grown) or inner region, for the host row forms (H-4b).
 
long vofCurvRow (const VofCurvJob &J, long r, bool grown, int &x0, int &x1)
 The first cell of row r of job J's grown or inner region, and the region's x range.
 
template<class Exec = SExec>
void vofCurvCompact (const VofCurvTable &T, bool grown, LField list, LField start, LField end)
 VofCurvature::compact: the interfacial cells of the grown (grown) or inner region, in the per-cascade scan's order; each job's range lands in start/end (absolute list positions).
 
template<class Exec = SExec>
void vofCurvPlanesZero (const VofCurvTable &T)
 reconstructPlanes (worklist mode), part 1: zero the four plane fields and clear the slot map (design G §5.2) over the grown region.
 
void vofCurvPlaneEntry (const VofCurvTable &T, int k, long q, LField list)
 Entry q (an absolute grown-list position) of job k of the planes pass (below).
 
template<class Exec = SExec>
void vofCurvPlanesList (const VofCurvTable &T, LField list, LField start, LField end)
 Part 2: wyReconstructCell over the grown interfacial list (upper bound: the grown region), and each cell's moment cache entry from the plane it just wrote, at its list position q, which becomes its slot (design G §5.2).
 
template<class Exec = SExec>
void vofCurvHfReset (const VofCurvTable &T)
 heightPass (worklist mode), part 1: every inner cell to kappa = 0, branch = kCurvNone.
 
void vofCurvFallbackTeams (const VofCurvTable &T, LField list, LField start, LField end)
 
void vofCurvListEntry (const VofCurvTable &T, int pass, int k, long q, LField list, LField cnt)
 Entry q (an absolute list position) of job k of a cascade list pass (see below): the body shared by the device and host launch forms.
 
template<class Exec = SExec>
void vofCurvListPass (const VofCurvTable &T, int pass, LField list, LField start, LField end, LField cnt)
 The per-list-cell passes of the cascade over the inner interfacial list (upper bound: the inner region): pass 0 = tiers 1-2 (curvHeightCell), 1 = tier 3 (curvFallbackCell), 2 = the admissibility clip (its count into cnt(8 (base + k) + 7)), 3 = the branch census (its seven counters into cnt(8 (base + k) + 0..6), in VofCurvature::census's order).
 
template<class Exec = SExec>
void vofCurvListPasses (const VofCurvTable &T, bool anyClip, LField list, LField start, LField end, LField cnt)
 The whole cascade list stage of a chunk: tiers 1-2, tier 3, the clip (only if some job clips), the census.
 
void vofCsfForceCell (const VofCurvJob &J, long i, double sig, double w0, double w1, double w2)
 The three CSF face-force components at the low faces of inner cell i of job J (the per-cell body of vofCsfForceBatch).
 
template<class Exec = SExec>
void vofCsfForceBatch (const VofCurvTable &T, double sig, double w0, double w1, double w2)
 VofBlockSet::buildCsfForce for every job, the three components per cell: the V4 balanced-force CSF at the LOW face of each inner cell, the per-block kernel's expressions verbatim.
 
void vofBatchStats (const VofBlockTable &T, const VofStatsPrev &P, double hh, double aeps, double dt, SField out, long off)
 VofBlockSet::measure for every job in one launch, one team per job: the first moments, then the central second moments and the PLIC interface area (VofBlockSet::interfaceArea's body) together, each as ONE team reduction over the inner box, with the per-cell expressions of the per-block reductions verbatim.
 
int vofAxisRuns (int a, int b, int gs, bool periodic, VofRun out[5], int outside=0)
 Decompose [a, b) on one axis into global runs.
 
VofBox vofClampBox (VofBox b, I3 gs, const std::array< bool, 3 > &per)
 Clamp a block box so it is representable: at most the whole axis on a periodic axis, inside the domain on a non-periodic one.
 
void vofAssignLpt (const std::vector< long > &w, int size, std::vector< int > &master)
 Longest-processing-time greedy (Graham 1969): heaviest block first onto the currently least-loaded rank, ties to the LOWEST rank id.
 
void vofAssignOrb (const std::vector< long > &w, int size, std::vector< int > &master)
 Weighted ORB over a 1-D block space: the blocks are laid out as a 1-D "grid" of N cells in block-id order carrying their weights, and core::decomp::BlockDecomposer<1> recursively bisects it into size CONTIGUOUS runs of balanced weight.
 
bool vofMoveOverlap (const VofMoveJob &a, const VofMoveJob &b)
 Debug check of H-4(e)'s premise: two jobs writing the same block array overlap.
 
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)
 Build the pieces of one block box against every rank's owned box, in canonical order.
 
void buildOutsideMask (UCField m, I3 e, int g, I3 o, I3 gs, bool px, bool py, bool pz)
 1 for every ghost cell of the block that lies OUTSIDE the global domain on a non-periodic axis, 0 everywhere else (inner cells, and ghosts that are interior halo on a periodic axis or across a rank boundary).
 
void bcColourConst (SField f, I3 e, int g, int a, int s, double value)
 Overwrite ALL g ghost layers of domain face (axis a, side s) with a constant colour.
 
void bcColourProfile (SField f, I3 e, int g, int a, int s, SField prof, int nc)
 Same, with a per-position colour on the face's (b, c) plane.
 
void bcColourBackflow (SField f, I3 e, int g, int a, int s, SField un, double backflow)
 inletOutlet on an outflow face: where the boundary face velocity points INTO the domain, the ghost band takes backflow; where it points out, the band is left as clampFill wrote it (zero-gradient, so what leaves is what is inside).
 
void copyFaceVelocity (SField dst, I3 de, int dg, SField src, I3 se, int sg, int dir)
 Embed the ENTIRE source extended block (inner + its sg ghost layers) into a LARGER destination extended block, converting the staggered face index of velocity component dir from the solver's low-face convention to the advector's high-face convention:
 
void copyAdvectedVelocity (SField dst, I3 de, int dg, SField src, I3 se, int sg, int dir, I3 n)
 The INVERSE of copyFaceVelocity, restricted to the inner region: bring a face-indexed field back from the advector's g = 3 block to the solver's G = 2 block, undoing the high-face -> low-face shift.
 
void periodicFill (SField f, I3 e, int g, bool px, bool py, bool pz)
 Periodic wrap of a single-block extended field at ghost width g (the single-rank ghost path; under MPI core::halo::GridHalo does this and the cross-rank exchange in one call).
 
void clampFill (SField f, I3 e, int g, I3 o, I3 gs, bool px, bool py, bool pz, int skip=0)
 Zero-gradient (Neumann) fill of the ghost cells that lie OUTSIDE the global domain on a non-periodic axis: each such cell takes the value at its globally CLAMPED index.
 
bool vofIsInterface (double c, double eps)
 wyIsMixed with a wisp threshold: the cell carries an interface only while eps < C < 1 - eps.
 
template<class SF >
void curvHeightCell (long i, SF c, SF mx, SF my, SF mz, SF al, SF kap, SF br, long s0, long s1, long s2, double mtol, double ptW, double ieps, bool forceFb, bool oneDir, bool useFit, VofMetric g, double peps)
 WO-V9 — the per-cell bodies of the two cascade tiers, lifted VERBATIM out of the kernels they used to be written inline in, so the dense (whole-region) and the compacted (interfacial-cell list) kernels below can share ONE body instead of two copies that could drift apart.
 
template<class SF >
void curvFallbackStore (long i, const PvFit &fit, SF kap, SF br)
 Tier 3's last step: solve the accumulated PV system and write the target's kappa and branch.
 
PECLET_VOF_DEVICE_NOINLINE void curvMomentsBuild (double mx, double my, double mz, double alpha, const VofMetric &g, PvMoments &P)
 pvMomentsBuild, the one call every planes pass and the dense accessor make (out of line on a device, see PECLET_VOF_DEVICE_NOINLINE).
 
template<class SF , class Acc >
void curvFallbackCell (long i, const Acc &acc, SF kap, SF br, long sy, long sz, int gr, double dW, double cmin, VofMetric g)
 Tier 3 for one target, one thread (host): the frame from the target's entry (pvFrameMoments: the physical normal, its tangents, the polygon's area centroid as the origin), the 5^3 stencil in canonical order (pvTermMoments), lower-triangle accumulation, the solve.
 
template<class Member , class SF , class Acc >
PECLET_VOF_DEVICE_NOINLINE void curvFallbackTeam (const Member &tm, const VofPvTeamScratch &sc, long i, const Acc &acc, SF kap, SF br, long sy, long sz, int gr, double dW, double cmin, VofMetric g)
 curvFallbackCell for one target on a whole team (design G §5.4; the team shape of C3, doc/vof_step_performance_design.md §4.7, §5.11).
 
double hfAreaElement (const double h[9])
 The height function's area element sqrt(1 + h_x^2 + h_y^2), with the same central differences hfPatchKappa uses — so the metric and the curvature are two derivatives of one interpolant, not two independent estimates.
 
double hfAreaElement (const double h[9], const VofMetric &g, int d)
 The same on an ANISOTROPIC cell: the graph is the PHYSICAL f = h_d * h_idx over the physical transverse coordinates, so the two slopes carry the ratios of hfPatchKappa(h, g, d) (V2.3).
 
double hfFootprintArea (const double h[9])
 Variant C — the height function's own footprint times its own metric (kAreaFootprint).
 
double hfFootprintArea (const double h[9], const VofMetric &g, int d)
 The same on an ANISOTROPIC cell.
 
void hfSurfaceNormal (const double h[9], int d, double m[3])
 The unit surface normal implied by the SAME 3x3 patch of heights, in world components.
 
void hfSurfaceNormal (const double h[9], int d, const VofMetric &g, double m[3])
 The same on an ANISOTROPIC cell — the PHYSICAL unit normal of the graph, i.e.
 
void pvSurfaceNormal (const double a[6], const double t1[3], const double t2[3], const double nn[3], double m[3])
 The unit normal of the fitted paraboloid z = a0 + a1 x + a2 y + ... AT THE ORIGIN of its own frame (the target cell's PLIC centroid), in world components: (-a1, -a2, 1) rotated by (t1, t2, nn) and normalized.
 
double interfaceAreaFromNormal (int mode, double mx, double my, double mz, double alpha, double c, const double ns[3], const VofMetric &g)
 The projection axis of kAreaMetric is the one ns is most aligned with — the column direction on tiers 1/2 by construction, and the same choice on tier 3 — because the footprint factor |n.e_d| is the one that stays farthest from zero there.
 
double interfaceAreaFromNormal (int mode, double mx, double my, double mz, double alpha, double c, const double ns[3])
 Unit-metric overload (isotropic cells) — the pre-Phase-3 signature, unchanged arithmetic.
 
void mcCornerPos (int k, double p[3])
 Local coordinates of dual-cube corner k (bit 0 = x, bit 1 = y, bit 2 = z).
 
double mcEdgeT (const McVertex &a, const McVertex &b, const double dv[3], int src)
 The parameter t in [0,1] at which the level set cuts the edge from corner a to corner b.
 
double mcTwiceArea (const double p0[3], const double p1[3], const double p2[3], const VofMetric &g)
 Twice the area of the triangle (p0, p1, p2) — i.e. |(p1-p0) x (p2-p0)|.
 
double mcPolygonArea (const double p[8][3], int n, const VofMetric &g)
 Area of a PLANAR polygon of n vertices, 1/2 |sum p_i x p_{i+1}|.
 
int mcClipHalf (double p[8][3], int n, int axis, double c, bool keepAbove)
 One Sutherland-Hodgman clip of a convex polygon against the half-space keepAbove ? x_axis >= c : x_axis <= c.
 
double mcClipToCorner (const double p0[3], const double p1[3], const double p2[3], int k, const VofMetric &g)
 The area of the part of one triangle inside corner k's octant of the dual cube.
 
int mcNearestMasked (const double p[3], int mask)
 The corner in mask nearest the point p; the lowest index breaks a tie.
 
int mcRetarget (int k, int mask)
 The retarget rule.
 
double mcTriangleToCorner (const double p0[3], const double p1[3], const double p2[3], int lc, int deposit, int mask, const VofMetric &g)
 The part of one triangle that belongs to local corner lc of the dual cube.
 
double mcCubeCornerArea (const McVertex v[8], int lc, int src, int deposit, const VofMetric &g)
 The interfacial area inside one dual cube that belongs to its local corner lc.
 
double mcCubeCornerArea (const McVertex v[8], int lc, int src, int deposit)
 Unit-metric overload (isotropic cells) — the pre-Phase-3 signature, unchanged arithmetic.
 
double vofPhaseRho (double rhoG, double rhoL, double c)
 Two-phase density from a fraction, written so that C = 1 returns rho_l and C = 0 returns rho_g EXACTLY (a plain rho_g + (rho_l - rho_g) C does not).
 
double vofMinmod (double a, double b)
 
double vofCellBox (const SField &c, const SField &mx, const SField &my, const SField &mz, const SField &al, long q, const double lo[3], const double hi[3], double eps=0.0)
 Liquid volume of the PLIC polyhedron of cell q inside a box given as per-axis intervals, as a fraction of a whole cell.
 
double plicArea (double mx, double my, double mz, double alpha)
 Area of the PLIC polygon {x in [0,1]^3 : m.x = alpha}, as a multiple of h^2.
 
double plicAreaMetric (double mx, double my, double mz, double alpha, const VofMetric &g)
 The PHYSICAL area of the PLIC polygon of an ANISOTROPIC cell, in hRef^2 (flow/doc/anisotropic_vof.md §2, §7 V5.1).
 
double pcUnitNormal (double mx, double my, double mz, double n[3])
 Unit (L2) interface normal from a PLIC normal of any scale.
 
double pcCentreDistance (double mx, double my, double mz, double alpha)
 Signed distance (in cells) from the CELL CENTRE to the PLIC plane, measured along n and POSITIVE ON THE GAS SIDE: phi(x) = (m.x - alpha)/|m|_2, evaluated at x = (1/2,1/2,1/2).
 
double pcOffsetDistance (double phic, const double n[3], double dx, double dy, double dz)
 Signed distance from the plane of the neighbour at INTEGER cell offset d, given the centre distance phic of the reference cell: phi_j = phi_c + d.n (cells, positive in the gas).
 
double pcGradWeight (double dx, double dy, double dz, const double n[3])
 Malan's collinearity weight for a sample at cell offset d: w = xi / |d| with xi = (d.n)^2/|d|^2 — i.e.
 
void pcGradAdd (PcGradFit &f, double w, double phi, double T, double Tg)
 Add one pure-phase sample at normal distance phi with temperature T.
 
double pcGradSolve (const PcGradFit &f)
 The fitted normal derivative dT/dn (0 if the side carried no usable sample).
 
double pcCurvedDistance (double phi, double dx, double dy, double dz, const double n[3], double kappa)
 WO-P3f — the sample's distance to the CURVED interface, not to the tangent plane.
 
double pcGradSolve2 (const PcGradFit &f)
 The QUADRATIC one-sided fit (WO-P23): T - T_G = G phi + Q phi^2, returning G.
 
double pcMassFlux (double kg, double gradG, double kl, double gradL, double hlv)
 mdot = (k_g dT_g/dn - k_l dT_l/dn)/h_lv with n the PLIC normal (liquid -> gas).
 
double pcRegressVolume (double mdot, double area, double dt, double rhoL)
 Liquid volume (fraction of a cell) removed from an interfacial cell in one step: dV = mdot A_G dt / rho_l.
 
double pcDivSource (double mdot, double area, double rhoG, double rhoL)
 The volumetric divergence source of one interfacial cell (1/s, cell volume = 1): S = mdot A_G (1/rho_g - 1/rho_l).
 
double pcInterfaceTemperature (double Tsat, double mdot, double Rint)
 Interfacial temperature under the Schrage-derived interfacial heat-transfer resistance (Bureš & Sato): T_G = T_sat + mdot R_int.
 
double pcCarriedValue (double Tg, double grad, double phic)
 The plane-anchored (ghost-fluid) interfacial Dirichlet condition (WO-P23, the first item).
 
double pcGfmTheta (double phiC, double nd, double s, double thetaMin, double thetaMax)
 thetaMin keeps the coefficient finite when the plane grazes the cell centre (the standard GFM clamp); thetaMax bounds it when the plane sits deep inside the far cell (|phi_c| <= sqrt(3)/2, so theta <= 1 + sqrt(3)/2) AND when the interface is nearly PARALLEL to this face's axis, where |phi_i|/|n_d| diverges and the honest statement is "the interface is far along this line".
 
double pcFaceConductivity (double ki, double kj, bool mi, bool mj)
 Face conductivity of the variable-coefficient energy operator, k(C) arithmetic in the colour with ONE exception: where one of the two cells carries the per-cell Dirichlet condition (an interfacial cell pinned at the plane-anchored T_G) the face takes the OTHER cell's k.
 
double pcPhaseMix (double vg, double vl, double c)
 Linear-in-C phase property (k, rho c_p), written so that C = 0 and C = 1 return the endpoint EXACTLY — the same construction as vofPhaseRho in momentum_advect.hpp.
 
bool pcPushWeights (const double n[3], double sgn, const bool avail[3], int step[3], double w[3])
 The clip-and-redistribute allocation: which face neighbours a clipped cell's residue is pushed into, and in what fractions.
 
double pcGfmThetaK (double phiC, double nd, double s, double kappa, double thetaMin, double thetaMax)
 WO-P3g item 3 — the curvature-consistent GFM distance.
 
double pcGfmThetaKAniso (double phiC, const double n[3], int d, double s, double kappa, double thetaMin, double thetaMax, const VofMetric &g)
 The same on an ANISOTROPIC cell (Phase 3, flow/doc/anisotropic_vof.md §7, V5.3).
 
PcGfmRow pcGfmRow (double theta, bool behind, int order)
 
double pcOperatorMassFlux (double q, double coefSum, double area, double hlv, double Rint)
 WO-P3g item 1 — mdot from the ENERGY OPERATOR's own interfacial flux.
 
bool pcIsInterfacial (double c, double eps)
 Is this a cell the regression acts on? The PLIC contract of advect_wy.hpp (wyIsMixed) with an explicit wisp guard: a cell whose colour is round-off residue has no meaningful plane, and giving it an area would inject an unbounded mdot into the census.
 
bool csfKappaDefined (double branch)
 Is the curvature of a cell carrying branch code branch usable?
 
bool csfFaceCurvature (double k0, double b0, double k1, double b1, double &kf)
 Face curvature for the balanced-force CSF, from the two cells sharing the face.
 
double csfFaceForce (double sigma, double kf, double dC, double hGrad)
 The balanced-force CSF contribution to the momentum RHS at one staggered velocity unknown, in the same units and at the same place as the incremental scheme's -(P(i) - P(i - s_c)).
 
double capillaryDt (double rhoSum, double h, double sigma)
 Brackbill capillary time-step limit sqrt((ρ₁+ρ₂) h³ / (4πσ)).
 
double coxVoinovAngle (double thetaE, double ca, double logRatio, double thetaMin, double thetaMax)
 Cox-Voinov / Afkhami-Zaleski-Bussmann grid-scale apparent angle, in RADIANS.
 
int vofHysteresisBase (double thetaApp, double thetaA, double thetaR, double &base)
 The hysteresis selector: which base angle the Cox-Voinov correction is applied to, and whether the line is pinned.
 
int vofDynamicContactAngle (double thetaApp, double uCl, double muL, double sigma, double logRatio, bool hyst, double thetaA, double thetaR, double thetaE, double thetaMin, double thetaMax, double &thetaOut, double &caOut)
 The whole rung in one pure function: hysteresis selector, then Cox-Voinov on the selected base.
 
bool vofWallTangent (const double mfIn[3], const double nwIn[3], double tEps, double that[3], double &cosApp, const VofMetric &g)
 The in-wall (tangential) unit direction of the interface and its apparent cosine.
 
bool vofWallTangent (const double mfIn[3], const double nwIn[3], double tEps, double that[3], double &cosApp)
 Unit-metric overload (isotropic cells) — the pre-Phase-3 signature, unchanged arithmetic.
 
double vofContactLineSpeed (const double u[3], const double that[3])
 U_cl = +u .
 

Variables

constexpr int kWySweepPerm [6][3]
 The six permutations of (x, y, z), cycled by step index so no direction is systematically favoured (perm[step % 6]).
 
constexpr int kVofBlockBatch = 16
 Jobs per launch; a container with more master blocks runs each stage in chunks.
 
template<class Exec >
constexpr bool kVofHostExec = std::is_same_v<typename Exec::memory_space, Kokkos::HostSpace>
 The host branches select on this (§14.2).
 
constexpr int kVofTier3MaxTeams = 1 << 16
 Tier 3 of the batched cascade on a DEVICE (C3, doc/vof_step_performance_design.md §4.7, §5.11; launch shape: coordinator decision 2026-10-02, option C).
 
constexpr int kVofStats = 14
 Doubles per job in the statistics packet: volume, centroid[3], velocity[3], area, moment[6].
 
constexpr int kVofWarp = 32
 The team size of a device tier-3 launch: one warp / wavefront.
 

Typedef Documentation

◆ UCField

Definition at line 89 of file advect_wy.hpp.

◆ LField

Definition at line 90 of file advect_wy.hpp.

◆ VofPvTerms

using peclet::flow::vof::VofPvTerms = typedef Kokkos::View<PvTerm*, SExec::scratch_memory_space, Kokkos::MemoryUnmanaged>

Team scratch of the device tier-3 body: the compacted stencil's canonical offset indices, their terms, and the 27 accumulated entries (about 10 KB at the 5^3 stencil).

Definition at line 277 of file curvature_field.hpp.

◆ VofPvInts

using peclet::flow::vof::VofPvInts = typedef Kokkos::View<int*, SExec::scratch_memory_space, Kokkos::MemoryUnmanaged>

Definition at line 278 of file curvature_field.hpp.

◆ VofPvDbls

using peclet::flow::vof::VofPvDbls = typedef Kokkos::View<double*, SExec::scratch_memory_space, Kokkos::MemoryUnmanaged>

Definition at line 279 of file curvature_field.hpp.

Enumeration Type Documentation

◆ VofDebrisSlot

The per-job slots of the debris packet (out(stride * (base + k) + slot)).

Enumerator
kDbSumD 

sum_D C in list order (the removed / census debris volume)

kDbSumR 

sum_R C in list order (signed residue)

kDbW 

sum_A C(1 - C)

kDbVA 

sum_A C

kDbDV 

the volume to return (D if removing, then R); replaced by lost once acted

kDbActed 

1 when the return ran (pending and VA >= kDebrisMinAttached and W > 0)

kDbND 

debris count

kDbNR 

residue count

kDbPending 

1 when there was something to remove (removeD or nR > 0)

kDbSlots 

Definition at line 817 of file block_batch.hpp.

◆ VofMasterAssign

How a block's MASTER rank is chosen.

The assignment is deliberately independent of where the bubble's cells live (that is the whole point of the container: bubble work is load-balanced separately from the flow decomposition), so it is a pure function of the REPLICATED block table and needs no communication in any of the three modes.

Enumerator
RoundRobin 

id % size — rung W0's assignment; blind to block SIZE

Lpt 

longest-processing-time greedy on the block cell counts (rung W1 default)

WeightedOrb 

weighted ORB over a 1-D "block space" (core::decomp::BlockDecomposer<1>)

Definition at line 207 of file block_container.hpp.

◆ InterfaceAreaMode

Which construction turns the cascade's normal into an area.

Exposed because the choice is a MEASUREMENT (see the sphere probe in the WO-P3c findings), not a taste.

Enumerator
kAreaPlic 

rung P0/P1: plicArea on the MYC normal. The old numbers.

kAreaMetric 

variant A: the PLIC footprint x the cascade's slope. DEFAULT.

kAreaNormal 

variant B: plicArea rebuilt on the cascade's normal.

kAreaFootprint 

variant C: the height function's OWN footprint x its own metric.

kAreaMcColour 

the C = 1/2 level set, whole triangles to the centroid's cell

kAreaMcColourSplit 

the same sheet, triangles CLIPPED to each cell's cube

kAreaMcPlic 

the PLIC signed-distance level set, centroid deposit

kAreaMcPlicSplit 

the PLIC signed-distance level set, clipped deposit

Definition at line 65 of file interface_area.hpp.

◆ McSource

What the edge crossing interpolates.

Enumerator
kMcSrcColour 

the C = 1/2 level set (the WO's design)

kMcSrcPlic 

the zero of the PLIC-reconstructed signed distance (the smooth variant)

Definition at line 96 of file marching_cubes.hpp.

◆ McDeposit

Which cell a triangle's area is booked to.

Enumerator
kMcDepositCentroid 

the whole triangle to the cell holding its centroid

kMcDepositSplit 

clipped to each cell's cube, each cell gets its own piece

Definition at line 102 of file marching_cubes.hpp.

◆ VofDynamicState

Which rule produced a band cell's imposed angle (vof_dynamic_field(4)).

Enumerator
kVofDynNone 

not a contact cell (no anchor / no usable normal): theta = the base

kVofDynStatic 

Cox-Voinov applied to the static base (no hysteresis configured)

kVofDynPinned 

theta_r <= theta_app <= theta_a: the apparent angle is imposed

kVofDynAdvancing 

theta_app > theta_a: theta_a (+ Cox-Voinov)

kVofDynReceding 

theta_app < theta_r: theta_r (+ Cox-Voinov)

kVofDynStateCount 

Definition at line 84 of file wetting_dynamic.hpp.

Function Documentation

◆ wyIsMixed() [1/2]

bool peclet::flow::vof::wyIsMixed ( double  c)

The mixed-cell predicate.

Reconstruction and fluxing MUST agree on it: the flux reads (m, alpha) exactly for the cells the reconstruction pass wrote, so the two tests are the same function by construction (that is also what makes the worklist a pure optimization).

Definition at line 95 of file advect_wy.hpp.

Referenced by peclet::flow::vof::WyAdvector::computeFluxesCut(), peclet::flow::vof::WyAdvector::diagnostics(), peclet::flow::vof::WyAdvector::maxCourantInterface(), peclet::flow::vof::WyAdvector::maxCourantInterfaceCut(), peclet::flow::vof::WyAdvector::reconstructImpl(), vofBatchCfl(), vofBatchWorklist(), vofCellBox(), wyFaceFlux(), and wyFaceFluxBc().

◆ wyIsMixed() [2/2]

bool peclet::flow::vof::wyIsMixed ( double  c,
double  eps 
)

WO-R2 item 4 — the same predicate with a WISP TOLERANCE.

eps == 0 is wyIsMixed(c) bit for bit (1.0 - 0.0 == 1.0), which is what makes WyAdvector::wispEps = 0 the V1 scheme verbatim.

Why it exists. An open-boundary domain can DRAIN completely, and no earlier VoF configuration could. When it does, the colour field is nothing but Weymouth-Yue round-off residue — values of order 1e-18 — and 0 < c < 1 still calls those cells MIXED. The MYC normal of a stencil of +-1e-18 values is degenerate and plicAlpha divides by it; WO-R measured the consequence on a drained slug run (nvidia-cuda, np = 1): sum C -1.5e-17 at step 185, -inf at 186, NaN at 187. Above the threshold a cell is treated as PURE for reconstruction and for fluxing, and its flux is the ALGEBRAIC C_donor * a of its ACTUAL colour — so nothing is thrown away and the exact telescoping conservation is untouched. 1e-8 is the same threshold the V3 curvature predicate already uses under surface tension (VofCurvature::interfaceEps).

Definition at line 111 of file advect_wy.hpp.

◆ wyColourJump()

bool peclet::flow::vof::wyColourJump ( double  a,
double  b,
double  eps 
)

The interface-BAND jump test of maxCourantInterface: two neighbouring cells carry a DIFFERENT colour.

WO-R2 item 4b — with an exact != the Weymouth-Yue round-off residue (min C ~ -3.8e-17 in the wake of a passing interface) keeps the whole wake inside the band and the interface-local limiter degenerates towards the global max: measured on Zalesak, vof_last_courant() reads 0.2545 after step 1 and 0.3110 by step 1000 on a case whose interface never exceeds 0.255, which forces the caller to raise set_vof_cfl_limit. eps == 0 is a != b bit for bit (|a-b| > 0 is false exactly when a == b, NaN aside, and a NaN colour has already failed every other gate).

Definition at line 123 of file advect_wy.hpp.

Referenced by peclet::flow::vof::WyAdvector::maxCourantInterface(), peclet::flow::vof::WyAdvector::maxCourantInterfaceCut(), and vofBatchCfl().

◆ wyReconstructCell()

template<class SF >
void peclet::flow::vof::wyReconstructCell ( const SF &  c,
long  i,
long  sy,
long  sz,
SF  mx,
SF  my,
SF  mz,
SF  alpha 
)

PLIC reconstruction of one cell: MYC normal from the 3^3 colour stencil + the analytic plane offset.

Container-free apart from the flat view read, so it stays a thin wrapper over plic.hpp. A template on the field type (anything with double& operator()(long): a View, or the raw-pointer VofRawField a batched multi-block launch walks) – the body is one expression sequence whatever the accessor.

Definition at line 139 of file advect_wy.hpp.

Referenced by peclet::flow::vof::WyAdvector::reconstructImpl(), peclet::flow::vof::VofCurvature::reconstructPlanes(), peclet::flow::vof::VofInterfaceArea::reconstructPlanes(), vofBatchPlic(), and vofCurvPlaneEntry().

◆ wyFaceFlux()

template<class SF >
double peclet::flow::vof::wyFaceFlux ( double  a,
long  p,
long  sd,
int  dir,
const SF &  c,
const SF &  mx,
const SF &  my,
const SF &  mz,
const SF &  alpha,
double  eps = 0.0 
)

Signed Eulerian donor-cell flux through the dir-face between cell p and cell p + sd, as a fraction of a cell volume, positive along +dir.

a is the face Courant number.

Definition at line 157 of file advect_wy.hpp.

References wyIsMixed().

Referenced by peclet::flow::vof::WyAdvector::computeFluxesCut(), peclet::flow::vof::WyAdvector::computeFluxesImpl(), and vofBatchFluxCell().

◆ wyFaceFluxBc()

template<class SF , class MF >
double peclet::flow::vof::wyFaceFluxBc ( double  a,
long  p,
long  sd,
int  dir,
const SF &  c,
const SF &  mx,
const SF &  my,
const SF &  mz,
const SF &  alpha,
const MF &  outside,
double  eps = 0.0 
)

WO-R sibling of wyFaceFlux for a face whose DONOR lies outside the global domain.

Identical in every branch except one: when the donor is marked outside, the flux is the ALGEBRAIC C_donor * a instead of the PLIC slab volume of a reconstructed plane.

Why the donor's reconstruction must not be used there (see vof/colour_bc.hpp for the long version): the ghost band on an inflow face is a uniform, prescribed DATUM — it has no MYC normal worth the name, and a fractional inflow colour states what fraction of the incoming FLUX is liquid, not where a sub-cell interface sits. C_donor * a is exactly that statement. At a wall the two rules agree trivially (the normal face velocity is zero, so a == 0), and at an outflow face the donor is the inner cell unless the flow reverses, so this branch is the backflow rule and nothing else.

With outside(donor) == 0 this returns wyFaceFlux's value BIT FOR BIT — the same expressions in the same order — which is what makes the mask branch inert (gate G5).

Definition at line 188 of file advect_wy.hpp.

References wyIsMixed().

Referenced by peclet::flow::vof::WyAdvector::computeFluxesCut(), peclet::flow::vof::WyAdvector::computeFluxesImpl(), and vofBatchFluxCell().

◆ vofJobOf()

int peclet::flow::vof::vofJobOf ( const long *  off,
long  t 
)

Job of flat index t (off[0] <= t < off[nj]): a walk over at most 16 offsets.

Definition at line 69 of file block_batch.hpp.

Referenced by vofBatchClamp(), vofBatchDebrisWrite(), vofBatchFlux(), vofBatchFreeze(), vofBatchPlic(), vofBatchUpdate(), vofBatchWorklist(), vofCsfForceBatch(), vofCurvCompact(), vofCurvListPass(), vofCurvPlanesList(), and vofHostListFor().

◆ vofSetOffsets() [1/2]

template<class F >
long peclet::flow::vof::vofSetOffsets ( VofBlockTable &  T,
F  size 
)
inline

◆ vofJobLen()

long peclet::flow::vof::vofJobLen ( const VofBlockJob &  J)
inline

◆ vofJobInner()

long peclet::flow::vof::vofJobInner ( const VofBlockJob &  J)
inline

◆ vofJobGrown1()

long peclet::flow::vof::vofJobGrown1 ( const VofBlockJob &  J)
inline

The inner region grown by one cell per side: WyAdvector::reconstructImpl's range.

Definition at line 94 of file block_batch.hpp.

Referenced by peclet::flow::vof::VofBlockSet::advectBatched().

◆ vofHostListFor()

template<class Exec , class Cnt , class F >
void peclet::flow::vof::vofHostListFor ( const char *  name,
int  nj,
Cnt  cnt,
F  body 
)
inline

H-4(a): a host list kernel over the EXACT entry counts, cnt(k) (read directly: the count Views are HostSpace), with a dynamic schedule in chunks of 16 – so the active entries spread over every thread instead of sitting at the head of one block's static share.

body(k, t) is entry t (0-based) of job k: the device form's entry q = start + t, which exits for t >= cnt.

Definition at line 124 of file block_batch.hpp.

References kVofBlockBatch, and vofJobOf().

◆ vofListCount()

template<class Tab >
auto peclet::flow::vof::vofListCount ( const Tab &  T,
LField  start,
LField  end 
)
inline

The [start, end) count of job k of a list launch.

Definition at line 144 of file block_batch.hpp.

Referenced by vofBatchPlic(), vofCurvListPass(), and vofCurvPlanesList().

◆ vofHostRowFor()

template<class Exec , class Rows , class F >
void peclet::flow::vof::vofHostRowFor ( const char *  name,
int  nj,
Rows  rows,
F  row 
)
inline

H-4(b): a host region kernel over ROWS – rows(k) x-runs for job k, row(k, r) handling job k's local row r with the x loop inside – so the job search and the index math run once per row, not once per cell.

Static schedule: the rows are uniform, as in ccFor3.

Definition at line 169 of file block_batch.hpp.

References kVofBlockBatch.

◆ vofHostRowScan()

template<class Exec , class Rows , class RowOf , class Hit >
void peclet::flow::vof::vofHostRowScan ( const char *  name,
int  nj,
int  base,
Rows  rows,
RowOf  rowOf,
Hit  hit,
LField  list,
LField  start,
LField  end 
)
inline

H-4(d): a host compaction scan over ROWS.

rowOf(k, r, x0, x1) returns the index of x = 0 on job k's local row r and sets its x range; a row's value is its count of hit(k, i) cells, and the final pass writes the row's hits at [prefix, prefix + count) in x order – so the list, and each job's start(base + k) / end(base + k), are the device scan's exactly.

Definition at line 187 of file block_batch.hpp.

References kVofBlockBatch.

◆ vofBatchCfl()

void peclet::flow::vof::vofBatchCfl ( const VofBlockTable &  T,
SField  out 
)
inline

One team per job: WyAdvector::maxCourant or maxCourantInterface (per the job's interfaceCfl), the same predicates and the same fmax chain.

Max is order-free, so the team reduction returns the per-block reduction's value exactly. out(base + k) = the job's Courant number (not yet all-reduced: a block advector has no globalMax).

Definition at line 225 of file block_batch.hpp.

References peclet::flow::L3(), vofBatchCfl(), wyColourJump(), and wyIsMixed().

Referenced by peclet::flow::vof::VofBlockSet::advectBatched(), and vofBatchCfl().

◆ vofBatchFreeze()

template<class Exec = SExec>
void peclet::flow::vof::vofBatchFreeze ( const VofBlockTable &  T)
inline

WyAdvector::freezeDilationFlag over every job's extended block.

T.off = extended lengths. Host: the extended block's rows (H-4b).

Definition at line 282 of file block_batch.hpp.

References vofBatchFreeze(), and vofJobOf().

Referenced by peclet::flow::vof::VofBlockSet::advectBatched(), and vofBatchFreeze().

◆ vofBatchWorklist()

template<class Exec = SExec>
void peclet::flow::vof::vofBatchWorklist ( const VofBlockTable &  T,
LField  list,
LField  start,
LField  end 
)
inline

WyAdvector::reconstructImpl's worklist scan, over the concatenated grown-by-one regions (T.off = vofJobGrown1).

The list of job k is written at list(listBase + start(k) ...) in the per-block scan's order; start(base + k) / end(base + k) receive the job's list range (absolute positions), so its count never leaves the device.

Definition at line 315 of file block_batch.hpp.

References peclet::flow::L3(), vofBatchWorklist(), vofJobOf(), and wyIsMixed().

Referenced by peclet::flow::vof::VofBlockSet::advectBatched(), and vofBatchWorklist().

◆ vofBatchPlic()

template<class Exec = SExec>
void peclet::flow::vof::vofBatchPlic ( const VofBlockTable &  T,
LField  list,
LField  start,
LField  end 
)
inline

wyReconstructCell over each job's compacted list: launched over the upper bound (the job's region), a thread past its job's count exits.

Host: the exact counts (H-4a).

Definition at line 360 of file block_batch.hpp.

References vofBatchPlic(), vofJobOf(), vofListCount(), and wyReconstructCell().

Referenced by peclet::flow::vof::VofBlockSet::advectBatched(), and vofBatchPlic().

◆ vofJobFaces()

long peclet::flow::vof::vofJobFaces ( const VofBlockJob &  J,
int  d 
)
inline

A sweep's face range of one job: the inner region with one extra layer on the low side of d.

Definition at line 391 of file block_batch.hpp.

References vofJobFaces().

Referenced by peclet::flow::vof::VofBlockSet::advectBatched(), and vofJobFaces().

◆ vofBatchFluxCell()

void peclet::flow::vof::vofBatchFluxCell ( const VofBlockJob &  J,
int  d,
long  p 
)

The flux of d-face p of job J (the per-cell body of vofBatchFlux, both launch forms).

Definition at line 396 of file block_batch.hpp.

References vofBatchFluxCell(), wyFaceFlux(), and wyFaceFluxBc().

Referenced by vofBatchFlux(), and vofBatchFluxCell().

◆ vofBatchFlux()

template<class Exec = SExec>
void peclet::flow::vof::vofBatchFlux ( const VofBlockTable &  T,
int  d 
)
inline

WyAdvector::computeFluxesImpl (uncut): one flux per d-face, wyFaceFlux (or, for a job carrying the WO-R mask, wyFaceFluxBc).

T.off = vofJobFaces(d). Host: the face range's rows (H-4b).

Definition at line 412 of file block_batch.hpp.

References peclet::flow::L3(), vofBatchFlux(), vofBatchFluxCell(), and vofJobOf().

Referenced by peclet::flow::vof::VofBlockSet::advectBatched(), and vofBatchFlux().

◆ vofBatchBcLedger()

void peclet::flow::vof::vofBatchBcLedger ( const VofBlockTable &  T,
int  d,
SField  out,
int  stride,
int  slot0 
)
inline

The WO-R boundary ledger (WyAdvector::accumulateBcFaceVolume) of the jobs carrying the mask: one thread per (job, side), summing its plane of d-fluxes in (p0 fastest, p1) row order into out(stride (base + k) + slot0 + face) – a fixed order, which is a C2 order change against the per-block MDRange2 reduction (recorded; doc/vof_step_performance_design.md §5.9).

The block container installs no mask today, so the path is inert on every shipped configuration.

Definition at line 453 of file block_batch.hpp.

References vofBatchBcLedger().

Referenced by peclet::flow::vof::VofBlockSet::advectBatched(), and vofBatchBcLedger().

◆ vofBatchUpdateCell()

void peclet::flow::vof::vofBatchUpdateCell ( const VofBlockJob &  J,
int  d,
long  i 
)

The sweep update of inner cell i of job J (the per-cell body of vofBatchUpdate).

Definition at line 479 of file block_batch.hpp.

References vofBatchUpdateCell().

Referenced by vofBatchUpdate(), and vofBatchUpdateCell().

◆ vofBatchUpdate()

template<class Exec = SExec>
void peclet::flow::vof::vofBatchUpdate ( const VofBlockTable &  T,
int  d 
)
inline

WyAdvector::applySweepImpl (uncut) over the inner regions (T.off = vofJobInner); the expression verbatim.

Host: the inner rows (H-4b).

Definition at line 495 of file block_batch.hpp.

References peclet::flow::L3(), vofBatchUpdate(), vofBatchUpdateCell(), and vofJobOf().

Referenced by peclet::flow::vof::VofBlockSet::advectBatched(), and vofBatchUpdate().

◆ vofBatchGhostZero()

template<class Exec = SExec>
void peclet::flow::vof::vofBatchGhostZero ( const VofBlockTable &  T)
inline

Pass 1: every ghost cell of every job -> 0.

T.off = extended lengths. Host: only the ghost shell's rows, in CutcellMG::fillWrap's three-slab split (H-4b) – the z-ghost planes (every y), the y-ghost rows of the inner z range, then the two x-ghost runs of each inner row.

Definition at line 532 of file block_batch.hpp.

References peclet::flow::L3(), and vofBatchGhostZero().

Referenced by peclet::flow::vof::VofBlockSet::batchGhostFill(), and vofBatchGhostZero().

◆ vofBatchPeriodic()

template<class Exec = SExec>
void peclet::flow::vof::vofBatchPeriodic ( const VofBlockTable &  T,
int  a 
)
inline

Pass 2 on axis a (colour_field.hpp::periodicFill, one axis): only the jobs whose block spans the periodic axis a are in T; T.off = the transverse plane sizes dims[b] * dims[c].

Host: only the rows of the two a-ghost slabs (H-4b) – for a = x the g cells at each end of every extended row, for a = y or z whole x-rows. Each ghost cell is written once, from a cell inside along a, as in the device form.

Definition at line 601 of file block_batch.hpp.

References peclet::flow::L3(), N, PECLET_FLOW_OMP_SIMD, and vofBatchPeriodic().

Referenced by peclet::flow::vof::VofBlockSet::batchGhostFill(), and vofBatchPeriodic().

◆ vofBatchClampCell()

void peclet::flow::vof::vofBatchClampCell ( const VofBlockJob &  J,
int  x,
int  y,
int  z,
int  q0,
int  q1,
int  q2,
bool  p0,
bool  p1,
bool  p2 
)

The clamp of extended cell (x, y, z) of job J (the per-cell body of vofBatchClamp).

Definition at line 670 of file block_batch.hpp.

References peclet::flow::L3(), and vofBatchClampCell().

Referenced by vofBatchClamp(), and vofBatchClampCell().

◆ vofBatchClamp()

template<class Exec = SExec>
void peclet::flow::vof::vofBatchClamp ( const VofBlockTable &  T,
I3  gs,
bool  px,
bool  py,
bool  pz 
)
inline

Pass 3 (colour_field.hpp::clampFill, skip = 0): the part of each extended block outside a NON-periodic domain takes the globally clamped value.

T.off = extended lengths. Host: the extended rows, skipping a row inside the domain on every non-periodic axis – no cell of it writes (H-4b). A source cell is inside on every non-periodic axis, so it is never written.

Definition at line 687 of file block_batch.hpp.

References vofBatchClamp(), vofBatchClampCell(), and vofJobOf().

Referenced by peclet::flow::vof::VofBlockSet::batchGhostFill(), and vofBatchClamp().

◆ vofDebrisCell()

template<class SF >
bool peclet::flow::vof::vofDebrisCell ( const SF &  c,
long  i,
long  sy,
long  sz,
double  ieps,
double  cfull 
)

The debris predicate of doc/vof_overlap_design.md §11 on cell i: interfacial, and NO cell of the marker at or above cfull within two cells (the 5^3 fit stencil).

Shared by the per-block scan and the batched one.

Definition at line 728 of file block_batch.hpp.

References vofDebrisCell(), and vofIsInterface().

Referenced by peclet::flow::vof::VofBlockSet::debrisPassBatchImpl(), vofBatchMark(), and vofDebrisCell().

◆ vofAttachedCell()

template<class SF >
bool peclet::flow::vof::vofAttachedCell ( const SF &  c,
long  i,
long  sy,
long  sz,
double  ieps,
double  cfull 
)

The ATTACHED predicate (§5.3 step 2): interfacial with a cell at or above cfull within two.

Definition at line 742 of file block_batch.hpp.

References vofAttachedCell(), and vofIsInterface().

Referenced by peclet::flow::vof::VofBlockSet::debrisPassBatchImpl(), vofAttachedCell(), and vofBatchMark().

◆ vofResidueCell()

bool peclet::flow::vof::vofResidueCell ( double  ci,
double  weps 
)

The residue predicate (§13): 0 < |C| <= wispEps.

Definition at line 755 of file block_batch.hpp.

References vofResidueCell().

Referenced by peclet::flow::vof::VofBlockSet::debrisPassBatchImpl(), vofBatchMark(), and vofResidueCell().

◆ vofBatchMark()

template<class Exec = SExec>
void peclet::flow::vof::vofBatchMark ( const VofBlockTable &  T,
int  kind,
double  ieps,
double  cfull,
double  weps,
LField  list,
LField  start,
LField  end 
)
inline

One mark scan over the concatenated inner regions (T.off = vofJobInner), in each block's inner-box index order – the per-block scan's order, so every list is the per-block list.

kind: 0 debris, 1 residue, 2 attached. Each job's range lands in start/end (absolute).

Definition at line 763 of file block_batch.hpp.

References vofAttachedCell(), vofBatchMark(), vofDebrisCell(), and vofResidueCell().

Referenced by peclet::flow::vof::VofBlockSet::debrisBatched(), and vofBatchMark().

◆ vofBatchDebrisSums()

void peclet::flow::vof::vofBatchDebrisSums ( const VofBlockTable &  T,
VofDebrisFlags  F,
LField  listD,
LField  listR,
LField  listA,
LField  sD0,
LField  eD0,
LField  sR0,
LField  eR0,
LField  sA0,
LField  eA0,
SField  out,
int  stride,
int  slot0 
)
inline

Step 2 of the pass: one thread per job, every sum in list order (the per-block kernel's single thread and order, so the same bits).

A job with nothing to remove computes only its census sum.

Definition at line 838 of file block_batch.hpp.

References F, kDbActed, kDbDV, kDbND, kDbNR, kDbPending, kDbSumD, kDbSumR, kDbVA, kDbW, and vofBatchDebrisSums().

Referenced by peclet::flow::vof::VofBlockSet::debrisBatched(), and vofBatchDebrisSums().

◆ vofDebrisActCount()

long peclet::flow::vof::vofDebrisActCount ( const VofDebrisFlags &  F,
LField  sA0,
LField  eA0,
SField  out,
long  ob,
int  j 
)

The number of act entries of job j (its D removed, then R, then A), 0 when it did not act: the per-block write kernel's range.

Definition at line 887 of file block_batch.hpp.

References F, kDbActed, kDbND, kDbNR, and vofDebrisActCount().

Referenced by vofDebrisActCount().

◆ vofDebrisActEntry()

void peclet::flow::vof::vofDebrisActEntry ( const VofBlockTable &  T,
const VofDebrisFlags &  F,
int  k,
long  t,
LField  listD,
LField  listR,
LField  listA,
LField  sD0,
LField  sR0,
LField  sA0,
LField  eA0,
SField  out,
int  stride,
int  slot0,
SField  ex,
Kokkos::View< int *, SMem >  capf 
)

Act entry t of job k (the per-block write kernel's expressions verbatim); t past the job's count, or a job that did not act, is a no-op.

Definition at line 899 of file block_batch.hpp.

References F, kDbActed, kDbDV, kDbND, kDbNR, kDbW, and vofDebrisActEntry().

Referenced by vofBatchDebrisWrite(), and vofDebrisActEntry().

◆ vofBatchDebrisWrite()

template<class Exec = SExec>
void peclet::flow::vof::vofBatchDebrisWrite ( const VofBlockTable &  T,
VofDebrisFlags  F,
LField  listD,
LField  listR,
LField  listA,
LField  sD0,
LField  sR0,
LField  sA0,
LField  eA0,
SField  out,
int  stride,
int  slot0,
SField  ex,
Kokkos::View< int *, SMem >  capf 
)
inline

Step 4, the act: removed cells -> 0, attached C += dV C(1-C)/W capped at 1 (the per-block write kernel's expressions verbatim).

Over the upper bound (the inner region: D, R and A are disjoint subsets of it); a job that did not act, or a thread past its count, exits. capf(j) is raised when some cell of job j was capped. Host: the exact counts (H-4a).

Definition at line 944 of file block_batch.hpp.

References F, vofBatchDebrisWrite(), vofDebrisActEntry(), and vofJobOf().

Referenced by peclet::flow::vof::VofBlockSet::debrisBatched(), and vofBatchDebrisWrite().

◆ vofBatchDebrisLost()

void peclet::flow::vof::vofBatchDebrisLost ( const VofBlockTable &  T,
LField  sA0,
LField  eA0,
SField  out,
int  stride,
int  slot0,
SField  ex,
Kokkos::View< int *, SMem >  capf 
)
inline

lost: the capped excess of an acted job summed in list order on one thread (only where a cap fired, as the per-block path; 0 otherwise).

Written over the job's kDbDV slot.

Definition at line 972 of file block_batch.hpp.

References kDbActed, kDbDV, and vofBatchDebrisLost().

Referenced by peclet::flow::vof::VofBlockSet::debrisBatched(), and vofBatchDebrisLost().

◆ vofBatchBox()

template<class Exec = SExec>
void peclet::flow::vof::vofBatchBox ( const VofBlockTable &  T,
double  eps,
SField  out,
int  stride,
int  slot0 
)
inline

One team per job: the tight box of |C| > eps over the inner region, in the block's LOCAL inner frame (bubbleBox's six Min/Max reductions; order-free).

Written as doubles (exact integers) to out(stride (base + k) + slot0 + {0..5}) = lo x, y, z, hi x, y, z; an empty block reads hi < lo (the reducers' identities), exactly as bubbleBox tests it.

Host (H-4c): TeamPolicy(nj, AUTO) is one thread per block there. Instead one row-parallel pass over every job's inner rows: a row takes its six extrema over x, then integer atomic min/max into the job's slots, which start at the same identities. Order-free, so bitwise.

Definition at line 1004 of file block_batch.hpp.

References kVofBlockBatch, and vofBatchBox().

Referenced by peclet::flow::vof::VofBlockSet::recentreBatched(), and vofBatchBox().

◆ vofSetOffsets() [2/2]

template<class F >
long peclet::flow::vof::vofSetOffsets ( VofCurvTable &  T,
F  size 
)
inline

Definition at line 1112 of file block_batch.hpp.

References kVofBlockBatch, and vofSetOffsets().

◆ vofCurvInner()

long peclet::flow::vof::vofCurvInner ( const VofCurvJob &  J)
inline

Definition at line 1120 of file block_batch.hpp.

References vofCurvInner().

Referenced by peclet::flow::vof::VofBlockSet::computeCsfBatched(), and vofCurvInner().

◆ vofCurvGrown()

long peclet::flow::vof::vofCurvGrown ( const VofCurvJob &  J)
inline

The inner region grown by kPvHalf: every cell VofCurvature::reconstructPlanes writes.

Definition at line 1124 of file block_batch.hpp.

References vofCurvGrown().

Referenced by peclet::flow::vof::VofBlockSet::computeCsfBatched(), and vofCurvGrown().

◆ vofCurvCell()

long peclet::flow::vof::vofCurvCell ( const VofCurvJob &  J,
long  r,
bool  grown 
)

Cell r of job J's grown (grown) or inner region, in the region's x-fastest order.

Definition at line 1129 of file block_batch.hpp.

References peclet::flow::L3(), and vofCurvCell().

Referenced by vofCsfForceBatch(), vofCurvCell(), and vofCurvCompact().

◆ vofCurvRows()

long peclet::flow::vof::vofCurvRows ( const VofCurvJob &  J,
bool  grown 
)
inline

Rows of job J's grown (grown) or inner region, for the host row forms (H-4b).

Definition at line 1139 of file block_batch.hpp.

References vofCurvRows().

Referenced by vofCurvRows().

◆ vofCurvRow()

long peclet::flow::vof::vofCurvRow ( const VofCurvJob &  J,
long  r,
bool  grown,
int &  x0,
int &  x1 
)
inline

The first cell of row r of job J's grown or inner region, and the region's x range.

Definition at line 1145 of file block_batch.hpp.

References peclet::flow::L3(), and vofCurvRow().

Referenced by vofCurvRow().

◆ vofCurvCompact()

template<class Exec = SExec>
void peclet::flow::vof::vofCurvCompact ( const VofCurvTable &  T,
bool  grown,
LField  list,
LField  start,
LField  end 
)
inline

VofCurvature::compact: the interfacial cells of the grown (grown) or inner region, in the per-cascade scan's order; each job's range lands in start/end (absolute list positions).

Definition at line 1156 of file block_batch.hpp.

References vofCurvCell(), vofCurvCompact(), vofIsInterface(), and vofJobOf().

Referenced by peclet::flow::vof::VofBlockSet::computeCsfBatched(), and vofCurvCompact().

◆ vofCurvPlanesZero()

template<class Exec = SExec>
void peclet::flow::vof::vofCurvPlanesZero ( const VofCurvTable &  T)
inline

reconstructPlanes (worklist mode), part 1: zero the four plane fields and clear the slot map (design G §5.2) over the grown region.

Host: the grown rows (H-4b).

Definition at line 1193 of file block_batch.hpp.

References vofCurvPlanesZero().

Referenced by peclet::flow::vof::VofBlockSet::computeCsfBatched(), and vofCurvPlanesZero().

◆ vofCurvPlaneEntry()

void peclet::flow::vof::vofCurvPlaneEntry ( const VofCurvTable &  T,
int  k,
long  q,
LField  list 
)

Entry q (an absolute grown-list position) of job k of the planes pass (below).

Definition at line 1228 of file block_batch.hpp.

References curvMomentsBuild(), vofCurvPlaneEntry(), and wyReconstructCell().

Referenced by vofCurvPlaneEntry(), and vofCurvPlanesList().

◆ vofCurvPlanesList()

template<class Exec = SExec>
void peclet::flow::vof::vofCurvPlanesList ( const VofCurvTable &  T,
LField  list,
LField  start,
LField  end 
)
inline

Part 2: wyReconstructCell over the grown interfacial list (upper bound: the grown region), and each cell's moment cache entry from the plane it just wrote, at its list position q, which becomes its slot (design G §5.2).

Host: the exact counts (H-4a).

Definition at line 1242 of file block_batch.hpp.

References vofCurvPlaneEntry(), vofCurvPlanesList(), vofJobOf(), and vofListCount().

Referenced by peclet::flow::vof::VofBlockSet::computeCsfBatched(), and vofCurvPlanesList().

◆ vofCurvHfReset()

template<class Exec = SExec>
void peclet::flow::vof::vofCurvHfReset ( const VofCurvTable &  T)
inline

heightPass (worklist mode), part 1: every inner cell to kappa = 0, branch = kCurvNone.

Host: the inner rows (H-4b).

Definition at line 1263 of file block_batch.hpp.

References vofCurvHfReset().

Referenced by peclet::flow::vof::VofBlockSet::computeCsfBatched(), and vofCurvHfReset().

◆ vofCurvFallbackTeams()

void peclet::flow::vof::vofCurvFallbackTeams ( const VofCurvTable &  T,
LField  list,
LField  start,
LField  end 
)
inline

◆ vofCurvListEntry()

void peclet::flow::vof::vofCurvListEntry ( const VofCurvTable &  T,
int  pass,
int  k,
long  q,
LField  list,
LField  cnt 
)

Entry q (an absolute list position) of job k of a cascade list pass (see below): the body shared by the device and host launch forms.

pass 0 = tiers 1-2, 1 = tier 3, 2 = the clip, 3 = the census, 4 = all four in that order (the host's fused pass, design G §5.5).

Definition at line 1336 of file block_batch.hpp.

References csfKappaDefined(), curvFallbackCell(), curvHeightCell(), and vofCurvListEntry().

Referenced by vofCurvListEntry(), and vofCurvListPass().

◆ vofCurvListPass()

template<class Exec = SExec>
void peclet::flow::vof::vofCurvListPass ( const VofCurvTable &  T,
int  pass,
LField  list,
LField  start,
LField  end,
LField  cnt 
)
inline

The per-list-cell passes of the cascade over the inner interfacial list (upper bound: the inner region): pass 0 = tiers 1-2 (curvHeightCell), 1 = tier 3 (curvFallbackCell), 2 = the admissibility clip (its count into cnt(8 (base + k) + 7)), 3 = the branch census (its seven counters into cnt(8 (base + k) + 0..6), in VofCurvature::census's order).

The counts are integer atomics: exact whatever the order. Host: the exact counts (H-4a).

Definition at line 1391 of file block_batch.hpp.

References vofCurvFallbackTeams(), vofCurvListEntry(), vofCurvListPass(), vofJobOf(), and vofListCount().

Referenced by vofCurvListPass().

◆ vofCurvListPasses()

template<class Exec = SExec>
void peclet::flow::vof::vofCurvListPasses ( const VofCurvTable &  T,
bool  anyClip,
LField  list,
LField  start,
LField  end,
LField  cnt 
)
inline

The whole cascade list stage of a chunk: tiers 1-2, tier 3, the clip (only if some job clips), the census.

A device runs the four passes (tier 3 as teams); a HOST runs ONE fused list pass that calls them in that order per entry (design G §5.5) — bitwise, because each entry's tier 3 reads only its own branch (written by its own tiers 1-2) plus planes, colour and cache written before the launch, the clip touches only its own kappa, and the census counts are integer atomics.

Definition at line 1424 of file block_batch.hpp.

References vofCurvListPasses().

Referenced by peclet::flow::vof::VofBlockSet::computeCsfBatched(), and vofCurvListPasses().

◆ vofCsfForceCell()

void peclet::flow::vof::vofCsfForceCell ( const VofCurvJob &  J,
long  i,
double  sig,
double  w0,
double  w1,
double  w2 
)

The three CSF face-force components at the low faces of inner cell i of job J (the per-cell body of vofCsfForceBatch).

Definition at line 1439 of file block_batch.hpp.

References vofCsfForceCell().

Referenced by vofCsfForceBatch(), and vofCsfForceCell().

◆ vofCsfForceBatch()

template<class Exec = SExec>
void peclet::flow::vof::vofCsfForceBatch ( const VofCurvTable &  T,
double  sig,
double  w0,
double  w1,
double  w2 
)
inline

VofBlockSet::buildCsfForce for every job, the three components per cell: the V4 balanced-force CSF at the LOW face of each inner cell, the per-block kernel's expressions verbatim.

Host: the inner rows (H-4b).

Definition at line 1461 of file block_batch.hpp.

References vofCsfForceBatch(), vofCsfForceCell(), vofCurvCell(), and vofJobOf().

Referenced by peclet::flow::vof::VofBlockSet::computeCsfBatched(), and vofCsfForceBatch().

◆ vofBatchStats()

void peclet::flow::vof::vofBatchStats ( const VofBlockTable &  T,
const VofStatsPrev &  P,
double  hh,
double  aeps,
double  dt,
SField  out,
long  off 
)
inline

VofBlockSet::measure for every job in one launch, one team per job: the first moments, then the central second moments and the PLIC interface area (VofBlockSet::interfaceArea's body) together, each as ONE team reduction over the inner box, with the per-cell expressions of the per-block reductions verbatim.

The summation ORDER is the team reduction's, not the per-block MDRange's: a recorded change (C2) of these diagnostics, which feed no state. Job k's values land at out(off + kVofStats * (T.base + k) + ...); the host applies them later (deferred read).

Definition at line 1530 of file block_batch.hpp.

References peclet::flow::vof::VofStatsPrev::c, kVofStats, peclet::flow::L3(), and peclet::flow::vof::VofStatsPrev::valid.

Referenced by peclet::flow::vof::VofBlockSet::measureBatched().

◆ vofAxisRuns()

int peclet::flow::vof::vofAxisRuns ( int  a,
int  b,
int  gs,
bool  periodic,
VofRun  out[5],
int  outside = 0 
)
inline

Decompose [a, b) on one axis into global runs.

Returns the count.

The runs partition the range in BLOCK-LOCAL index, not in global index: a range longer than the periodic axis (which is exactly what the EXTENDED box of a block spanning the whole axis is) revisits the same global cells, and that is correct — the extra local cells are the block's periodic ghosts and must carry the wrapped owner's value.

Definition at line 150 of file block_container.hpp.

Referenced by vofBuildPieces().

◆ vofClampBox()

VofBox peclet::flow::vof::vofClampBox ( VofBox  b,
I3  gs,
const std::array< bool, 3 > &  per 
)
inline

Clamp a block box so it is representable: at most the whole axis on a periodic axis, inside the domain on a non-periodic one.

Definition at line 177 of file block_container.hpp.

Referenced by peclet::flow::vof::VofBlockSet::recentreTo(), peclet::flow::vof::VofBlockSet::seedBox(), peclet::flow::vof::VofBlockSet::seedBoxWithColour(), and peclet::flow::vof::VofBlockSet::seedSphere().

◆ vofAssignLpt()

void peclet::flow::vof::vofAssignLpt ( const std::vector< long > &  w,
int  size,
std::vector< int > &  master 
)
inline

Longest-processing-time greedy (Graham 1969): heaviest block first onto the currently least-loaded rank, ties to the LOWEST rank id.

4/3 − 1/(3p) of optimal in the worst case, and — the property the bitwise gates need — a deterministic pure function of the replicated weights, identical on every rank without an exchange.

Definition at line 217 of file block_container.hpp.

Referenced by peclet::flow::vof::VofBlockSet::plannedMasters().

◆ vofAssignOrb()

void peclet::flow::vof::vofAssignOrb ( const std::vector< long > &  w,
int  size,
std::vector< int > &  master 
)
inline

Weighted ORB over a 1-D block space: the blocks are laid out as a 1-D "grid" of N cells in block-id order carrying their weights, and core::decomp::BlockDecomposer<1> recursively bisects it into size CONTIGUOUS runs of balanced weight.

This is the work order's BlockDecomposer::init(…, weights) route, using the same partitioner the flow decomposition uses. Its constraint (and its weakness here) is that a rank's blocks must be contiguous in id, which LPT is free of.

Definition at line 246 of file block_container.hpp.

Referenced by peclet::flow::vof::VofBlockSet::plannedMasters().

◆ vofMoveOverlap()

bool peclet::flow::vof::vofMoveOverlap ( const VofMoveJob &  a,
const VofMoveJob &  b 
)
inline

Debug check of H-4(e)'s premise: two jobs writing the same block array overlap.

Definition at line 88 of file block_exchange.hpp.

Referenced by peclet::flow::vof::VofBlockExchange::moveLocalPieces().

◆ vofBuildPieces()

void peclet::flow::vof::vofBuildPieces ( const VofBox &  box,
I3  gs,
const std::array< bool, 3 > &  per,
const std::vector< VofBox > &  rankBox,
std::vector< VofPiece > &  out,
int  outside = 0 
)
inline

Build the pieces of one block box against every rank's owned box, in canonical order.

outside > 0 (the face-velocity gather only) also covers up to outside cells beyond each face of a non-periodic domain, owned by the rank whose box touches that face and read from its patch ghost layers (so outside must not exceed the patch ghost width). The rank boxes tile the domain, so the grown boxes tile the grown domain and the pieces still partition the block box.

Definition at line 101 of file block_exchange.hpp.

References peclet::flow::vof::VofBox::intersect(), and vofAxisRuns().

Referenced by peclet::flow::vof::VofBlockExchange::gatherFaceVelHost(), peclet::flow::vof::VofBlockExchange::gatherImpl(), peclet::flow::vof::VofBlockExchange::scatterColourMaxHost(), and peclet::flow::vof::VofBlockExchange::scatterImpl().

◆ buildOutsideMask()

void peclet::flow::vof::buildOutsideMask ( UCField  m,
I3  e,
int  g,
I3  o,
I3  gs,
bool  px,
bool  py,
bool  pz 
)
inline

1 for every ghost cell of the block that lies OUTSIDE the global domain on a non-periodic axis, 0 everywhere else (inner cells, and ghosts that are interior halo on a periodic axis or across a rank boundary).

o is the GLOBAL index of the block's inner cell (0,0,0) and gs the global grid size — the same pair clampFill uses, and for the same reason: the test must be a GLOBAL one, so that a rank in the middle of a decomposition marks nothing and the mask is decomposition-independent by construction.

Definition at line 84 of file colour_bc.hpp.

References peclet::flow::L3().

Referenced by peclet::flow::Solver< Grid >::vofRebuildBcBlock().

◆ bcColourConst()

void peclet::flow::vof::bcColourConst ( SField  f,
I3  e,
int  g,
int  a,
int  s,
double  value 
)
inline

Overwrite ALL g ghost layers of domain face (axis a, side s) with a constant colour.

Call AFTER the exchange and clampFill; the caller is responsible for the ownership test (touchesGlobalFace) so a rank whose block does not touch that global face writes nothing.

Definition at line 101 of file colour_bc.hpp.

References peclet::flow::vof::bcdetail::axisDims().

Referenced by peclet::flow::Solver< Grid >::vofApplyColourBc(), and peclet::flow::Solver< Grid >::vofBcPropGhosts().

◆ bcColourProfile()

void peclet::flow::vof::bcColourProfile ( SField  f,
I3  e,
int  g,
int  a,
int  s,
SField  prof,
int  nc 
)
inline

Same, with a per-position colour on the face's (b, c) plane.

prof is indexed prof(p0 * nc + p1) over the block's FULL ghost-inclusive transverse extents — i.e. it has already been resampled with the clamp rule IbmSolver::setDomainBcProfile uses, so the kernel indexes it directly by face position and needs no bounds arithmetic.

Definition at line 123 of file colour_bc.hpp.

References peclet::flow::vof::bcdetail::axisDims().

Referenced by peclet::flow::Solver< Grid >::vofApplyColourBc(), and peclet::flow::Solver< Grid >::vofBcPropGhosts().

◆ bcColourBackflow()

void peclet::flow::vof::bcColourBackflow ( SField  f,
I3  e,
int  g,
int  a,
int  s,
SField  un,
double  backflow 
)
inline

inletOutlet on an outflow face: where the boundary face velocity points INTO the domain, the ghost band takes backflow; where it points out, the band is left as clampFill wrote it (zero-gradient, so what leaves is what is inside).

un is the advector's own face-velocity component a (HIGH-face convention, un(i) is the +a face of cell i), so the boundary face of side 0 is the +a face of the last ghost cell (a-index g-1) and the boundary face of side 1 is the +a face of the last inner cell (a-index na - g - 1). Inward is + on side 0 and − on side 1. The test is on the velocity the advector is about to flux with, which is why this must run AFTER bridgeVelocityToVof — see IbmSolver::vofApplyColourBc.

Definition at line 152 of file colour_bc.hpp.

References peclet::flow::vof::bcdetail::axisDims().

Referenced by peclet::flow::Solver< Grid >::vofApplyColourBc().

◆ copyFaceVelocity()

void peclet::flow::vof::copyFaceVelocity ( SField  dst,
I3  de,
int  dg,
SField  src,
I3  se,
int  sg,
int  dir 
)
inline

Embed the ENTIRE source extended block (inner + its sg ghost layers) into a LARGER destination extended block, converting the staggered face index of velocity component dir from the solver's low-face convention to the advector's high-face convention:

dst(x + s - [dir==0], y + s - [dir==1], z + s - [dir==2]) = src(x, y, z),   s = dg - sg

i.e. a plain concentric embed on the two transverse axes and one cell less on the component's own axis (u_adv(i) = u_solver(i + s_dir) — see the file header; omitting this shift is a silent, 35 %-scale conservation defect). Requires equal inner extents and dg > sg, which the only caller satisfies (g = 3 colour block from the G = 2 velocity block).

Definition at line 65 of file colour_field.hpp.

References peclet::flow::L3().

Referenced by peclet::flow::Solver< Grid >::bridgeVelocityToVof(), and peclet::flow::Solver< Grid >::buildVofGeometry().

◆ copyAdvectedVelocity()

void peclet::flow::vof::copyAdvectedVelocity ( SField  dst,
I3  de,
int  dg,
SField  src,
I3  se,
int  sg,
int  dir,
I3  n 
)
inline

The INVERSE of copyFaceVelocity, restricted to the inner region: bring a face-indexed field back from the advector's g = 3 block to the solver's G = 2 block, undoing the high-face -> low-face shift.

dst(x + dg, y + dg, z + dg) = src(x + sg - [dir==0], y + sg - [dir==1], z + sg - [dir==2])

for x,y,z over the inner region n. Used at rung V2b (WO-K) to hand the momentum-consistent advected velocity back to the momentum RHS. The source index runs one cell BELOW the advector's inner region along dir — that is exactly the momentum control volume the solver's first velocity unknown on that axis belongs to, and MomentumConsistentAdvector::cvRange computes there for this reason (the two index sets are in one-to-one correspondence, so nothing is dropped or duplicated across a rank boundary).

Definition at line 87 of file colour_field.hpp.

References peclet::flow::L3().

◆ periodicFill()

void peclet::flow::vof::periodicFill ( SField  f,
I3  e,
int  g,
bool  px,
bool  py,
bool  pz 
)
inline

Periodic wrap of a single-block extended field at ghost width g (the single-rank ghost path; under MPI core::halo::GridHalo does this and the cross-rank exchange in one call).

Sequential axis passes, so edges and corners come out right.

Definition at line 101 of file colour_field.hpp.

References N.

Referenced by peclet::flow::vof::VofBlockSet::fillBlockGhosts(), peclet::flow::Solver< Grid >::vofExchangeRaw(), and peclet::flow::Solver< Grid >::vofExchangeScalar().

◆ clampFill()

void peclet::flow::vof::clampFill ( SField  f,
I3  e,
int  g,
I3  o,
I3  gs,
bool  px,
bool  py,
bool  pz,
int  skip = 0 
)
inline

Zero-gradient (Neumann) fill of the ghost cells that lie OUTSIDE the global domain on a non-periodic axis: each such cell takes the value at its globally CLAMPED index.

o is the GLOBAL index of this block's inner cell (0,0,0) and gs the global grid size. The clamp is per axis and global, which is what makes it decomposition-independent (WO-E finding 2: sequential zero-gradient axis passes over the block are NOT — a corner cell outside the domain in y whose x-neighbourhood is an interior halo gets a different value on a different decomposition, and that difference reaches C through the MYC stencil of the inner corner cell). The clamped source is provably inside the block's extended range and is either an inner cell or a ghost the exchange has already filled, so call this AFTER the halo exchange / periodic fill. skip (ISSUES sweep item 3) is a per-FACE bitmask (bit 2a+s = face -x,+x,-y,+y,-z,+z): a ghost cell that lies outside the global domain across a SKIPPED face is left ALONE instead of receiving the zero-gradient copy. That band is owned by something else — the theta-consistent solid-band fill, when a contact angle is imposed on a DOMAIN wall, whose write set is exactly those cells and which would otherwise be overwritten by the fill's own closing exchange. skip = 0 (the default) is the shipped rule bit for bit.

Definition at line 139 of file colour_field.hpp.

References peclet::flow::L3().

Referenced by peclet::flow::vof::VofBlockSet::fillBlockGhosts(), peclet::flow::Solver< Grid >::vofExchangeRaw(), and peclet::flow::Solver< Grid >::vofExchangeScalar().

◆ vofIsInterface()

◆ curvHeightCell()

template<class SF >
void peclet::flow::vof::curvHeightCell ( long  i,
SF  c,
SF  mx,
SF  my,
SF  mz,
SF  al,
SF  kap,
SF  br,
long  s0,
long  s1,
long  s2,
double  mtol,
double  ptW,
double  ieps,
bool  forceFb,
bool  oneDir,
bool  useFit,
VofMetric  g,
double  peps 
)

WO-V9 — the per-cell bodies of the two cascade tiers, lifted VERBATIM out of the kernels they used to be written inline in, so the dense (whole-region) and the compacted (interfacial-cell list) kernels below can share ONE body instead of two copies that could drift apart.

The only substitution is the cell index: i is a parameter instead of L3(x, y, z, e), and tier 3's neighbour index L3(x+ox, y+oy, z+oz, e) becomes the identical i + ox + oy*sy + oz*sz (the same integers, by the definition of L3). Nothing else changed, which is what lets the ctest gate the compaction as bitwise.

Definition at line 59 of file curvature_field.hpp.

References vofIsInterface().

Referenced by peclet::flow::vof::VofCurvature::heightPass(), and vofCurvListEntry().

◆ curvFallbackStore()

template<class SF >
void peclet::flow::vof::curvFallbackStore ( long  i,
const PvFit &  fit,
SF  kap,
SF  br 
)

Tier 3's last step: solve the accumulated PV system and write the target's kappa and branch.

Definition at line 172 of file curvature_field.hpp.

Referenced by curvFallbackCell(), and curvFallbackTeam().

◆ curvMomentsBuild()

PECLET_VOF_DEVICE_NOINLINE void peclet::flow::vof::curvMomentsBuild ( double  mx,
double  my,
double  mz,
double  alpha,
const VofMetric &  g,
PvMoments &  P 
)

pvMomentsBuild, the one call every planes pass and the dense accessor make (out of line on a device, see PECLET_VOF_DEVICE_NOINLINE).

Definition at line 199 of file curvature_field.hpp.

Referenced by peclet::flow::vof::VofPvOnTheFly< SF >::entry(), peclet::flow::vof::VofCurvature::reconstructPlanes(), and vofCurvPlaneEntry().

◆ curvFallbackCell()

template<class SF , class Acc >
void peclet::flow::vof::curvFallbackCell ( long  i,
const Acc &  acc,
SF  kap,
SF  br,
long  sy,
long  sz,
int  gr,
double  dW,
double  cmin,
VofMetric  g 
)

Tier 3 for one target, one thread (host): the frame from the target's entry (pvFrameMoments: the physical normal, its tangents, the polygon's area centroid as the origin), the 5^3 stencil in canonical order (pvTermMoments), lower-triangle accumulation, the solve.

acc is VofPvCached or VofPvOnTheFly.

Definition at line 246 of file curvature_field.hpp.

References curvFallbackStore().

Referenced by peclet::flow::vof::VofCurvature::fallbackBatch(), peclet::flow::vof::VofCurvature::fallbackPass(), and vofCurvListEntry().

◆ curvFallbackTeam()

template<class Member , class SF , class Acc >
PECLET_VOF_DEVICE_NOINLINE void peclet::flow::vof::curvFallbackTeam ( const Member &  tm,
const VofPvTeamScratch &  sc,
long  i,
const Acc &  acc,
SF  kap,
SF  br,
long  sy,
long  sz,
int  gr,
double  dW,
double  cmin,
VofMetric  g 
)

curvFallbackCell for one target on a whole team (design G §5.4; the team shape of C3, doc/vof_step_performance_design.md §4.7, §5.11).

Every lane computes the frame (deterministic, so no broadcast); then (1) a team scan COMPACTS the (2gr+1)^3 stencil to its interfacial cells, keeping the canonical order k = ((oz+gr)(2gr+1) + (oy+gr))(2gr+1) + (ox+gr); (2) the lanes map those n cells to PvTerms (pvTermMoments); (3) ENTRY-PARALLEL accumulation: lane e sums entry e of the 27 (lower triangle + b) over the accepted terms in canonical order, from 0.0, with pvFitAccum's expression; and one lane assembles the fit and solves. Each entry is the one-thread body's sum, term for term in the same order; under nvcc's FMA contraction the per-lane expression rounds differently from the single-lane fold, so on CUDA this is part of WO-3's RECORDED round-off change (with --fmad=false the two are bitwise). A warp reduction is NOT used (register: not bitwise). Every lane must call it (it holds team barriers). Out of line on a device (PECLET_VOF_DEVICE_NOINLINE): the batched and the per-block kernels run one compiled body.

Definition at line 311 of file curvature_field.hpp.

References curvFallbackStore().

Referenced by peclet::flow::vof::VofCurvature::fallbackBatch(), and vofCurvFallbackTeams().

◆ hfAreaElement() [1/2]

double peclet::flow::vof::hfAreaElement ( const double  h[9])

The height function's area element sqrt(1 + h_x^2 + h_y^2), with the same central differences hfPatchKappa uses — so the metric and the curvature are two derivatives of one interpolant, not two independent estimates.

h is the 3x3 patch of heights indexed p + 3 q exactly as hfPatchKappa takes it (p along the first transverse axis d1 = (d+1)%3, q along d2 = (d+2)%3). The heights are SIGNED (they grow with liquid whichever end the liquid is at) so the interface position is orient * h; the area element is even in that sign and needs no orientation.

Definition at line 87 of file interface_area.hpp.

References hfAreaElement().

Referenced by hfAreaElement(), hfAreaElement(), and hfFootprintArea().

◆ hfAreaElement() [2/2]

double peclet::flow::vof::hfAreaElement ( const double  h[9],
const VofMetric &  g,
int  d 
)

The same on an ANISOTROPIC cell: the graph is the PHYSICAL f = h_d * h_idx over the physical transverse coordinates, so the two slopes carry the ratios of hfPatchKappa(h, g, d) (V2.3).

Unit metric == the expression above, bit for bit.

Definition at line 96 of file interface_area.hpp.

References hfAreaElement().

◆ hfFootprintArea() [1/2]

double peclet::flow::vof::hfFootprintArea ( const double  h[9])

Variant C — the height function's own footprint times its own metric (kAreaFootprint).

This is the one construction whose cell pieces TILE, and the measurement says that is the whole question. Write the interface as the graph x_d = f(x_d1, x_d2) linearized on the patch, f = h_0 + h_x x + h_y y in cell-centred coordinates. The piece of it that belongs to THIS cell projects on the transverse square to R = {(x,y) in [-1/2,1/2]^2 : |f| <= 1/2}, and the pieces of the cells of one column partition that square EXACTLY — which is what makes the column sums telescope and the total converge:

A = |R| sqrt(1 + h_x^2 + h_y^2) ,   |R| = F(+1/2) - F(-1/2) ,
F(t) = fraction of the square with f <= t = plicVolume(h_x, h_y, 0, t - h_0 + (h_x+h_y)/2)

(the 2-D PLIC volume fraction, which plicVolume returns for a zero third component).

Why the other two variants cannot converge, measured. They both take the FOOTPRINT from the cell's own PLIC polygon, and on a CURVED interface the neighbouring cells' planes are different planes, so their projections overlap or leave a gap and the column's footprints do not sum to one. The resulting deficit is FIRST order in h/R and it is not a normal error: an exact- fraction circle, with the EXACT radial normal and an analytic chord, sums to -2.86 / -1.47 / -1.42 / -0.84 % of 2 pi R at R = 8/12/20/28 (the code's cylinder probe reads -1.42 % at R = 20 on the same footing). Replacing the slope changes that by 0.02 pp; replacing the FOOTPRINT is what removes it.

On a plane it is exact, and for an axis-aligned plane it is exact BITWISE: h_x = h_y = 0 makes plicVolume take its degenerate-normal branch (1 and 0), so |R| = 1 and the metric is 1.

Definition at line 128 of file interface_area.hpp.

Referenced by peclet::flow::vof::VofInterfaceArea::heightPass().

◆ hfFootprintArea() [2/2]

double peclet::flow::vof::hfFootprintArea ( const double  h[9],
const VofMetric &  g,
int  d 
)

The same on an ANISOTROPIC cell.

The FOOTPRINT FRACTION f — which part of the transverse unit square this cell's piece of the graph covers — is cell-crossing bookkeeping in INDEX units and is unchanged; the PHYSICAL area is that fraction times the physical footprint h_d1 h_d2 times the physical slope factor. Unit metric: f * 1.0 * sqrt(1+hx^2+hy^2).

Definition at line 141 of file interface_area.hpp.

References hfAreaElement().

◆ hfSurfaceNormal() [1/2]

void peclet::flow::vof::hfSurfaceNormal ( const double  h[9],
int  d,
double  m[3] 
)

The unit surface normal implied by the SAME 3x3 patch of heights, in world components.

x_d = orient h(x_d1, x_d2) has normal (-h_x, -h_y, 1) in the frame (d1, d2, d), up to a global sign that no consumer here cares about (plicArea is invariant under m -> -m with alpha re-solved from the same colour: the two planes are point reflections of each other through the cell centre, an isometry of the cube — the identity V(a) + V(1-a) = 1 of plicArea's derivation).

Definition at line 157 of file interface_area.hpp.

References hfSurfaceNormal().

Referenced by peclet::flow::vof::VofInterfaceArea::heightPass(), hfSurfaceNormal(), and hfSurfaceNormal().

◆ hfSurfaceNormal() [2/2]

void peclet::flow::vof::hfSurfaceNormal ( const double  h[9],
int  d,
const VofMetric &  g,
double  m[3] 
)

The same on an ANISOTROPIC cell — the PHYSICAL unit normal of the graph, i.e.

the same expression with the physical slopes. Unit metric == the overload above, bit for bit.

Definition at line 169 of file interface_area.hpp.

References hfSurfaceNormal().

◆ pvSurfaceNormal()

void peclet::flow::vof::pvSurfaceNormal ( const double  a[6],
const double  t1[3],
const double  t2[3],
const double  nn[3],
double  m[3] 
)

The unit normal of the fitted paraboloid z = a0 + a1 x + a2 y + ... AT THE ORIGIN of its own frame (the target cell's PLIC centroid), in world components: (-a1, -a2, 1) rotated by (t1, t2, nn) and normalized.

This is the tier-3 twin of hfSurfaceNormal and it is the same derivative of the same fit paraboloidKappa differentiates twice.

Definition at line 184 of file interface_area.hpp.

References pvSurfaceNormal().

Referenced by peclet::flow::vof::VofInterfaceArea::fallbackPass(), and pvSurfaceNormal().

◆ interfaceAreaFromNormal() [1/2]

double peclet::flow::vof::interfaceAreaFromNormal ( int  mode,
double  mx,
double  my,
double  mz,
double  alpha,
double  c,
const double  ns[3],
const VofMetric &  g 
)

The projection axis of kAreaMetric is the one ns is most aligned with — the column direction on tiers 1/2 by construction, and the same choice on tier 3 — because the footprint factor |n.e_d| is the one that stays farthest from zero there.

kAreaPlic returns plicArea(m, alpha) unchanged (times the metric), so the old behaviour is one branch and not a special case.

Turn an accurate PHYSICAL unit normal ns into the cell's PHYSICAL interfacial area (in hRef^2), given the cell's MYC plane (m, alpha) in INDEX space and its colour c. mode selects the construction (see the file header).

**Phase 3 (flow/doc/anisotropic_vof.md §7, V5.1).** A unit-cube polygon area A_xi is the physical area A_xi * det(H) * s(m) (plicAreaMetric); ns is a physical direction, so the footprint ratio of kAreaMetric is taken between the two PHYSICAL normals. Unit metric: det = 1, s = q/q = 1, and the physical normal of m is m/|m| — today's arithmetic.

Definition at line 205 of file interface_area.hpp.

References interfaceAreaFromNormal(), kAreaNormal, kAreaPlic, and plicAreaMetric().

Referenced by peclet::flow::vof::VofInterfaceArea::fallbackPass(), peclet::flow::vof::VofInterfaceArea::heightPass(), interfaceAreaFromNormal(), and interfaceAreaFromNormal().

◆ interfaceAreaFromNormal() [2/2]

double peclet::flow::vof::interfaceAreaFromNormal ( int  mode,
double  mx,
double  my,
double  mz,
double  alpha,
double  c,
const double  ns[3] 
)

Unit-metric overload (isotropic cells) — the pre-Phase-3 signature, unchanged arithmetic.

Definition at line 238 of file interface_area.hpp.

References interfaceAreaFromNormal().

◆ mcCornerPos()

void peclet::flow::vof::mcCornerPos ( int  k,
double  p[3] 
)

Local coordinates of dual-cube corner k (bit 0 = x, bit 1 = y, bit 2 = z).

Definition at line 125 of file marching_cubes.hpp.

Referenced by mcCubeCornerArea(), mcNearestMasked(), and mcRetarget().

◆ mcEdgeT()

double peclet::flow::vof::mcEdgeT ( const McVertex &  a,
const McVertex &  b,
const double  dv[3],
int  src 
)

The parameter t in [0,1] at which the level set cuts the edge from corner a to corner b.

kMcSrcColour is the plain linear interpolation of psi. kMcSrcPlic solves each endpoint cell's own plane along the edge — phi_a(t) = d_a + t (n_a . dv) and phi_b(t) = d_b + (t-1)(n_b . dv) — and averages the roots the two cells supply, which is what makes a plane exact from both sides; it falls back to the psi interpolation when neither endpoint carries a usable plane (two pure cells of opposite phase in contact, i.e. an interface with no mixed cell at all).

Definition at line 139 of file marching_cubes.hpp.

References kMcSrcPlic.

Referenced by mcCubeCornerArea().

◆ mcTwiceArea()

double peclet::flow::vof::mcTwiceArea ( const double  p0[3],
const double  p1[3],
const double  p2[3],
const VofMetric &  g 
)

Twice the area of the triangle (p0, p1, p2) — i.e. |(p1-p0) x (p2-p0)|.

Definition at line 205 of file marching_cubes.hpp.

Referenced by mcTriangleToCorner().

◆ mcPolygonArea()

double peclet::flow::vof::mcPolygonArea ( const double  p[8][3],
int  n,
const VofMetric &  g 
)

Area of a PLANAR polygon of n vertices, 1/2 |sum p_i x p_{i+1}|.

Definition at line 219 of file marching_cubes.hpp.

Referenced by mcClipToCorner().

◆ mcClipHalf()

int peclet::flow::vof::mcClipHalf ( double  p[8][3],
int  n,
int  axis,
double  c,
bool  keepAbove 
)

One Sutherland-Hodgman clip of a convex polygon against the half-space keepAbove ? x_axis >= c : x_axis <= c.

Returns the new vertex count.

Definition at line 241 of file marching_cubes.hpp.

Referenced by mcClipToCorner().

◆ mcClipToCorner()

double peclet::flow::vof::mcClipToCorner ( const double  p0[3],
const double  p1[3],
const double  p2[3],
int  k,
const VofMetric &  g 
)

The area of the part of one triangle inside corner k's octant of the dual cube.

Definition at line 273 of file marching_cubes.hpp.

References mcClipHalf(), and mcPolygonArea().

Referenced by mcTriangleToCorner().

◆ mcNearestMasked()

int peclet::flow::vof::mcNearestMasked ( const double  p[3],
int  mask 
)

The corner in mask nearest the point p; the lowest index breaks a tie.

mask == 0 returns -1.

Definition at line 289 of file marching_cubes.hpp.

References mcCornerPos().

Referenced by mcRetarget(), and mcTriangleToCorner().

◆ mcRetarget()

int peclet::flow::vof::mcRetarget ( int  k,
int  mask 
)

The retarget rule.

A triangle (or a clipped piece of one) can land in a cell the wisp predicate calls PURE — measured: exactly the 12 axis-tangent "pole" cells of a sphere, where the cell just outside the interface is exactly C = 1 while the C = 1/2 sheet still passes through its cube (1.5 h^2, i.e. 0.19 % of the area at R = 8). The phase-change consumer reads A only on interfacial cells, so that area would be SILENTLY DROPPED from int mdot dA — and, worse, the drop partly cancels the sheet's own error, which is exactly the kind of cancellation this campaign has twice been caught by. So a piece whose octant's cell carries no interface is booked to the nearest cell that does. mask is the 8-bit interfacial census of the cube's corners; mask == 0 (an interface with no mixed cell at all) keeps the geometric target.

Definition at line 317 of file marching_cubes.hpp.

References mcCornerPos(), and mcNearestMasked().

Referenced by mcTriangleToCorner().

◆ mcTriangleToCorner()

double peclet::flow::vof::mcTriangleToCorner ( const double  p0[3],
const double  p1[3],
const double  p2[3],
int  lc,
int  deposit,
int  mask,
const VofMetric &  g 
)

The part of one triangle that belongs to local corner lc of the dual cube.

kMcDepositCentroid books the whole triangle to the corner nearest its centroid (a cube's 8 corner OCTANTS partition it, so exactly one corner claims it); kMcDepositSplit clips the triangle to that corner's octant [0,1/2] x ... and returns the clipped area. Summed over the 8 corners both rules return the triangle's own area exactly — including the retargeted pieces, which are moved between corners and never lost.

Definition at line 333 of file marching_cubes.hpp.

References kMcDepositCentroid, mcClipToCorner(), mcNearestMasked(), mcRetarget(), and mcTwiceArea().

Referenced by mcCubeCornerArea().

◆ mcCubeCornerArea() [1/2]

double peclet::flow::vof::mcCubeCornerArea ( const McVertex  v[8],
int  lc,
int  src,
int  deposit,
const VofMetric &  g 
)

The interfacial area inside one dual cube that belongs to its local corner lc.

Kuhn's 6 tetrahedra of the cube, each cut by the level set into nothing, one triangle (a 1-vs-3 corner split) or a quad (2-vs-2, emitted as two triangles). Every vertex sits on a tet edge and therefore on a cube edge, a face diagonal or the body diagonal — all of which are shared with the neighbouring cubes in the same orientation, so the sheet is closed.

Definition at line 367 of file marching_cubes.hpp.

References mcCornerPos(), mcCubeCornerArea(), mcEdgeT(), and mcTriangleToCorner().

Referenced by peclet::flow::vof::VofMcArea::areaPass(), mcCubeCornerArea(), and mcCubeCornerArea().

◆ mcCubeCornerArea() [2/2]

double peclet::flow::vof::mcCubeCornerArea ( const McVertex  v[8],
int  lc,
int  src,
int  deposit 
)

Unit-metric overload (isotropic cells) — the pre-Phase-3 signature, unchanged arithmetic.

Definition at line 445 of file marching_cubes.hpp.

References mcCubeCornerArea().

◆ vofPhaseRho()

double peclet::flow::vof::vofPhaseRho ( double  rhoG,
double  rhoL,
double  c 
)

Two-phase density from a fraction, written so that C = 1 returns rho_l and C = 0 returns rho_g EXACTLY (a plain rho_g + (rho_l - rho_g) C does not).

That exactness is what makes the single-phase reduction density-ratio-independent to round-off.

Definition at line 139 of file momentum_advect.hpp.

Referenced by peclet::flow::vof::MomentumConsistentAdvector::diagnostics(), peclet::flow::vof::MomentumConsistentAdvector::freezeShiftedState(), peclet::flow::vof::MomentumConsistentAdvector::momentumUpdate(), and peclet::flow::vof::MomentumConsistentAdvector::momentumUpdateCut().

◆ vofMinmod()

◆ vofCellBox()

double peclet::flow::vof::vofCellBox ( const SField &  c,
const SField &  mx,
const SField &  my,
const SField &  mz,
const SField &  al,
long  q,
const double  lo[3],
const double  hi[3],
double  eps = 0.0 
)

Liquid volume of the PLIC polyhedron of cell q inside a box given as per-axis intervals, as a fraction of a whole cell.

Non-mixed cells take the algebraic branch C * volume — the same contract wyFaceFlux uses, which is what keeps full/empty cells exactly stationary.

Definition at line 152 of file momentum_advect.hpp.

References wyIsMixed().

Referenced by peclet::flow::vof::MomentumConsistentAdvector::buildShiftedColour(), peclet::flow::vof::MomentumConsistentAdvector::buildShiftedColourCut(), peclet::flow::vof::MomentumConsistentAdvector::momentumFluxes(), and peclet::flow::vof::MomentumConsistentAdvector::momentumFluxesCut().

◆ plicArea()

double peclet::flow::vof::plicArea ( double  mx,
double  my,
double  mz,
double  alpha 
)

Area of the PLIC polygon {x in [0,1]^3 : m.x = alpha}, as a multiple of h^2.

The plane<->volume relation gives it for free: the volume between the planes at alpha and alpha + d(alpha) is a slab of thickness d(alpha)/|m|_2 on the polygon, so

A = |m|_2 * dV/d(alpha) .

V(alpha) is the Scardovelli–Zaleski piecewise cubic, so dV/d(alpha) is a piecewise quadratic and this is ANALYTIC, not a finite difference. With the components mirrored into the positive octant (an isometry: the area is unchanged), L1-normalized (n1+n2+n3 = 1) and sorted ascending, and with w = min(a, 1-a) (the point reflection x -> 1-x gives V(a)+V(1-a) = 1, hence ‘V’(a) = V'(1-a)‘):

6 n1 n2 n3 V(w) = w^3 - sum_i <w-n_i>^3 + sum_i <w-(1-n_i)>^3        (<y> = max(y,0))
2 n1 n2 n3 V&rsquo;(w) = w^2 - sum_i <w-n_i>^2 + sum_i <w-(1-n_i)>^2 .

Two rearrangements make it cancellation-free as n1 -> 0 (a nearly axis-aligned plane, which is EXACTLY the planar configuration rungs P0/P1 run on):

  • w^2 - <w-n1>^2 is evaluated as n1(2w-n1) where w > n1 — a product, not a difference of two O(1) numbers;
  • the terms <w-(1-n_2)>^2 = <w-n1-n3>^2 and <w-n3>^2 are PAIRED (they differ by O(n1)) and their difference is evaluated in closed form. Likewise for the n2 pair.
  • <w-(1-n1)>^2 = <w-n2-n3>^2 is identically zero: n1 <= 1/3 (it is the minimum of three numbers summing to 1) so n2+n3 >= 2/3 > 1/2 >= w. The degenerate branches (n1 = 0: the polygon is a rectangle; n2 = 0: an axis-aligned face) are taken exactly, so a grid-aligned interface returns A = 1 to the last bit.

Definition at line 88 of file phase_change.hpp.

References plicArea().

Referenced by plicArea(), and plicAreaMetric().

◆ plicAreaMetric()

double peclet::flow::vof::plicAreaMetric ( double  mx,
double  my,
double  mz,
double  alpha,
const VofMetric &  g 
)

The PHYSICAL area of the PLIC polygon of an ANISOTROPIC cell, in hRef^2 (flow/doc/anisotropic_vof.md §2, §7 V5.1).

The unit-cube polygon area A_xi and the physical one are related by the cofactor transformation of the area vector, A n = det(H) H^-T (A_xi m_hat), i.e.

A_phys = A_xi * det(H) * s(m),      s(m) = |H^-1 m| / |m|

which is exactly plicArea times two factors that are both EXACTLY 1.0 at equal spacings (det = 1*1*1, and s is one L2 norm divided by the same L2 norm) — so the isotropic call is plicArea(...) * 1.0 * 1.0, bit for bit today's number.

Definition at line 137 of file phase_change.hpp.

References plicArea(), and plicAreaMetric().

Referenced by peclet::flow::vof::VofInterfaceArea::fallbackPass(), interfaceAreaFromNormal(), peclet::flow::Solver< Grid >::pcBuildInterface(), plicAreaMetric(), and peclet::flow::Solver< Grid >::vofInterfaceArea().

◆ pcUnitNormal()

double peclet::flow::vof::pcUnitNormal ( double  mx,
double  my,
double  mz,
double  n[3] 
)

Unit (L2) interface normal from a PLIC normal of any scale.

Points INTO THE GAS (liquid -> gas), the plic.hpp convention. Returns the L2 norm of the input (0 => n untouched).

Definition at line 154 of file phase_change.hpp.

Referenced by peclet::flow::vof::VofMcArea::planePass().

◆ pcCentreDistance()

double peclet::flow::vof::pcCentreDistance ( double  mx,
double  my,
double  mz,
double  alpha 
)

Signed distance (in cells) from the CELL CENTRE to the PLIC plane, measured along n and POSITIVE ON THE GAS SIDE: phi(x) = (m.x - alpha)/|m|_2, evaluated at x = (1/2,1/2,1/2).

A cell that is mostly liquid has phi_c < 0.

Definition at line 167 of file phase_change.hpp.

References pcCentreDistance().

Referenced by peclet::flow::Solver< Grid >::pcBuildInterface(), pcCentreDistance(), peclet::flow::Solver< Grid >::pcUpdateThermalMask(), and peclet::flow::vof::VofMcArea::planePass().

◆ pcOffsetDistance()

double peclet::flow::vof::pcOffsetDistance ( double  phic,
const double  n[3],
double  dx,
double  dy,
double  dz 
)

Signed distance from the plane of the neighbour at INTEGER cell offset d, given the centre distance phic of the reference cell: phi_j = phi_c + d.n (cells, positive in the gas).

Definition at line 176 of file phase_change.hpp.

References pcOffsetDistance().

Referenced by peclet::flow::Solver< Grid >::pcBuildInterface(), pcOffsetDistance(), and peclet::flow::Solver< Grid >::pcUpdateThermalMask().

◆ pcGradWeight()

double peclet::flow::vof::pcGradWeight ( double  dx,
double  dy,
double  dz,
const double  n[3] 
)

Malan's collinearity weight for a sample at cell offset d: w = xi / |d| with xi = (d.n)^2/|d|^2 — i.e.

w = (d.n)^2 / |d|^3. Samples lying along the normal dominate; samples in the tangent plane (which carry no normal-gradient information) get zero weight.

Definition at line 201 of file phase_change.hpp.

References pcGradWeight().

Referenced by peclet::flow::Solver< Grid >::pcBuildInterface(), pcGradWeight(), and peclet::flow::Solver< Grid >::pcUpdateThermalMask().

◆ pcGradAdd()

void peclet::flow::vof::pcGradAdd ( PcGradFit &  f,
double  w,
double  phi,
double  T,
double  Tg 
)

Add one pure-phase sample at normal distance phi with temperature T.

Definition at line 210 of file phase_change.hpp.

References pcGradAdd().

Referenced by peclet::flow::Solver< Grid >::pcBuildInterface(), pcGradAdd(), and peclet::flow::Solver< Grid >::pcUpdateThermalMask().

◆ pcGradSolve()

double peclet::flow::vof::pcGradSolve ( const PcGradFit &  f)

The fitted normal derivative dT/dn (0 if the side carried no usable sample).

Definition at line 221 of file phase_change.hpp.

References pcGradSolve().

Referenced by peclet::flow::Solver< Grid >::pcBuildInterface(), pcGradSolve(), pcGradSolve2(), and peclet::flow::Solver< Grid >::pcUpdateThermalMask().

◆ pcCurvedDistance()

double peclet::flow::vof::pcCurvedDistance ( double  phi,
double  dx,
double  dy,
double  dz,
const double  n[3],
double  kappa 
)

WO-P3f — the sample's distance to the CURVED interface, not to the tangent plane.

pcOffsetDistance returns the distance to the interfacial cell's PLIC PLANE. The one-sided fit then models T as a function of that, which is exact for a plane and FIRST ORDER in h/R for anything curved: a sample at lateral offset rho from the normal line sits rho^2/(2R) further from a sphere than from its tangent plane, and since T grows away from the interface every off-axis sample is HOTTER than the plane model expects, so the fitted dT/dn — and with it mdot — comes out systematically HIGH. Measured a priori on an exact sphere carrying an exactly linear profile (vof_scriven.py --mdot-probe --mdot-prof linear), 128^3, sub = 16: +19.2 / +12.1 / +8.8 / +6.2 % at R = 6 / 10 / 14 / 20 with observed order 0.91 / 0.93 / 0.98 in h/R — a clean first-order curvature bias, and the largest single error in the P3 rung.

The correction. With kappa = div(n) (the mean curvature of the level sets of the gas-positive distance function; -2/R for a spherical gas bubble, whose n points inward), the signed distance to the SURFACE of a point at plane distance phi and lateral offset rho is phi + kappa rho^2/4 to second order. Derivation, for a gas sphere of radius R with the foot point at R z_hat: a point displaced zeta along n = -z_hat and rho tangentially sits at radius sqrt((R-zeta)^2 + rho^2) ~ R - zeta + rho^2/(2R), so the gas-positive distance R - r is zeta - rho^2/(2R) where the plane says zeta; and kappa/4 = -1/(2R).

Definition at line 245 of file phase_change.hpp.

Referenced by peclet::flow::Solver< Grid >::pcBuildInterface(), and peclet::flow::Solver< Grid >::pcUpdateThermalMask().

◆ pcGradSolve2()

double peclet::flow::vof::pcGradSolve2 ( const PcGradFit &  f)

The QUADRATIC one-sided fit (WO-P23): T - T_G = G phi + Q phi^2, returning G.

The linear fit above is a straight line through the interface value fitted to samples that start ~1 cell away and reach ~2.5, so a profile with curvature ‘T’'biases it by O(T'' h)— a clean FIRST-ORDER error inmdot`, and (once the plane-anchored Dirichlet has removed the cell-centre mismatch) the leading error of the whole rung. This is the Aslam-quadratic lever VOF_PLAN §9 item 1 names, in its least-squares rather than PDE form: the same samples, one more basis function, no extra stencil reach and no sweeps.

Normal equations of the 2-parameter weighted fit: [S2 S3; S3 S4] [G; Q] = [b1; b2] , S_k = sum w phi^k , b_k = sum w phi^k (T - T_G). Falls back to the linear solve when the system is ill-conditioned (fewer than two distinct normal distances — the determinant is then zero to round-off), so a side that carries a single usable sample still returns something rather than a division by noise.

Definition at line 266 of file phase_change.hpp.

References pcGradSolve().

Referenced by peclet::flow::Solver< Grid >::pcBuildInterface(), and peclet::flow::Solver< Grid >::pcUpdateThermalMask().

◆ pcMassFlux()

double peclet::flow::vof::pcMassFlux ( double  kg,
double  gradG,
double  kl,
double  gradL,
double  hlv 
)

mdot = (k_g dT_g/dn - k_l dT_l/dn)/h_lv with n the PLIC normal (liquid -> gas).

Positive = evaporation (liquid consumed). See the file header for the sign derivation.

Definition at line 279 of file phase_change.hpp.

Referenced by peclet::flow::Solver< Grid >::pcBuildInterface().

◆ pcRegressVolume()

double peclet::flow::vof::pcRegressVolume ( double  mdot,
double  area,
double  dt,
double  rhoL 
)

Liquid volume (fraction of a cell) removed from an interfacial cell in one step: dV = mdot A_G dt / rho_l.

Positive = evaporation.

Definition at line 286 of file phase_change.hpp.

Referenced by peclet::flow::Solver< Grid >::pcRegress().

◆ pcDivSource()

double peclet::flow::vof::pcDivSource ( double  mdot,
double  area,
double  rhoG,
double  rhoL 
)

The volumetric divergence source of one interfacial cell (1/s, cell volume = 1): S = mdot A_G (1/rho_g - 1/rho_l).

Definition at line 292 of file phase_change.hpp.

Referenced by peclet::flow::Solver< Grid >::pcBuildInterface().

◆ pcInterfaceTemperature()

double peclet::flow::vof::pcInterfaceTemperature ( double  Tsat,
double  mdot,
double  Rint 
)

Interfacial temperature under the Schrage-derived interfacial heat-transfer resistance (Bureš & Sato): T_G = T_sat + mdot R_int.

R_int = 0 is the hard Dirichlet.

Definition at line 298 of file phase_change.hpp.

Referenced by peclet::flow::Solver< Grid >::pcBuildInterface(), and peclet::flow::Solver< Grid >::pcUpdateThermalMask().

◆ pcCarriedValue()

double peclet::flow::vof::pcCarriedValue ( double  Tg,
double  grad,
double  phic 
)

The plane-anchored (ghost-fluid) interfacial Dirichlet condition (WO-P23, the first item).

Rungs P0/P1 pinned the whole interfacial CELL at T_G, so the numerical thermal boundary sat at the CELL CENTRE while the mass-flux gradient is fitted from the PLIC PLANE. The mismatch is up to half a cell and it CHANGES SIGN as the interface sweeps through a cell — the first-order, oscillating error component WO-P01 measured (P1 orders 1.07 on 64->128 and 1.50 on 128->256).

The repair is a PER-FACE condition, not a per-cell value. For a PURE cell i whose neighbour j across the face in direction d is interfacial, the plane of cell j sits at signed normal distance phi_i = phi_c(j) - s n_d(j) from cell i's centre (s the face step from i to j), so along the grid line it is theta = |phi_i| / |n_d| cells away. That face's row therefore becomes the DIRICHLET row k open (T_i - T_G)/theta instead of the interior coupling k open (T_i - T_j), and the interfacial cell's own value is never read.

Why the per-CELL value the work order describes does NOT work, measured. Giving the interfacial cell the value the one-sided profile takes at its centre is right for the neighbour on the side the fit was taken from and WRONG for the neighbour on the other side: on the P1 Stefan ladder it heats the saturated LIQUID through the interfacial cell's liquid-side face, which gives the liquid a spurious gradient, which feeds straight back into mdot = (k_g G_g - k_l G_l)/h_lv. The error saturates instead of converging — measured +6.20 / +5.62 / +5.42 % at N = 64/128/256, observed order 0.10, against +1.31 / +0.59 / +0.20 % and order 1.37 for the cell-centre pinning it was meant to improve on. One cell-centred value cannot serve two sides; the per-face form below can, and does (the liquid-side face then reads T_G at its own theta, i.e. exactly zero flux for a saturated liquid).

The value to CARRY in an interfacial cell (WO-P23). With the plane-anchored rows above, no neighbour reads it — the solve's only use for it is that the cell becomes PURE a few steps later as the interface sweeps past, and then it IS read, both by the energy operator and by the one-sided gradient fit. Leaving it at T_G makes every newly exposed cell start at the saturation temperature, i.e. a cold spot that the diffusion then has to remove: measured on the P1 Stefan ladder it holds the interface back by a clean first order (-1.31 / -0.72 / -0.36 % at N = 64/128/256, order 0.93). The value the one-sided linear profile takes at the cell centre, T_G + G phi_c, removes it. (This is the per-cell value the work order describes; it is right HERE, where nothing reads it across a face, and wrong as a boundary condition — see above.)

Definition at line 336 of file phase_change.hpp.

Referenced by peclet::flow::Solver< Grid >::pcUpdateThermalMask().

◆ pcGfmTheta()

double peclet::flow::vof::pcGfmTheta ( double  phiC,
double  nd,
double  s,
double  thetaMin,
double  thetaMax 
)

thetaMin keeps the coefficient finite when the plane grazes the cell centre (the standard GFM clamp); thetaMax bounds it when the plane sits deep inside the far cell (|phi_c| <= sqrt(3)/2, so theta <= 1 + sqrt(3)/2) AND when the interface is nearly PARALLEL to this face's axis, where |phi_i|/|n_d| diverges and the honest statement is "the interface is far along this line".

The n_d -> 0 limit has to be thetaMax, and it has to be reached CONTINUOUSLY. The first version fell back to theta = 1 below |n_d| = 1e-6, i.e. a JUMP of 1 -> 1.9 in that face's coefficient across a threshold that a transverse MYC normal component crosses on ROUND-OFF: in a quasi-1D scene the energy solve's red-black parity asymmetry gives symmetric columns ~1e-16 differences in T and hence ~1e-8 transverse normal components (WO-P01 finding 3), so the transverse faces of every interfacial cell sit on that threshold. Measured on the P2 sucking-interface MPI gate at np = 1 — where the two solvers differ only by round-off — max|C_dist - C_ref| went 3.3e-16 at 3 steps -> 2.4e-4 at 12 -> 4.5e-4 at 55 (a plateau) with the jump, against 4.1e-14 / 7.4e-14 for the rung P0/P1 treatment on the same scene. A discontinuous coefficient is a round-off amplifier even when the jump itself is bounded.

Definition at line 356 of file phase_change.hpp.

Referenced by pcGfmThetaK().

◆ pcFaceConductivity()

double peclet::flow::vof::pcFaceConductivity ( double  ki,
double  kj,
bool  mi,
bool  mj 
)

Face conductivity of the variable-coefficient energy operator, k(C) arithmetic in the colour with ONE exception: where one of the two cells carries the per-cell Dirichlet condition (an interfacial cell pinned at the plane-anchored T_G) the face takes the OTHER cell's k.

The exception is not cosmetic. A Dirichlet row is an identity row, so the only thing this face coefficient does is set the conductance with which the PURE neighbour draws heat to the interface. That conductance must be the pure neighbour's own phase conductivity — an arithmetic mean with k(C_iface) puts a slab of the WRONG phase between the pure cell and a boundary condition that already sits (plane-anchored) at the interface, and at the water/steam ratio k_l/k_g ~ 26 that is a 50 % error in the wall heat flux of every interfacial face.

Definition at line 375 of file phase_change.hpp.

Referenced by peclet::flow::scalarBuildDiffusionVarK().

◆ pcPhaseMix()

double peclet::flow::vof::pcPhaseMix ( double  vg,
double  vl,
double  c 
)

Linear-in-C phase property (k, rho c_p), written so that C = 0 and C = 1 return the endpoint EXACTLY — the same construction as vofPhaseRho in momentum_advect.hpp.

Definition at line 385 of file phase_change.hpp.

Referenced by peclet::flow::vof::VofEnergyAdvector::energyUpdate(), peclet::flow::Solver< Grid >::pcCarryDeposit(), and peclet::flow::Solver< Grid >::pcUpdateEnergyProps().

◆ pcPushWeights()

bool peclet::flow::vof::pcPushWeights ( const double  n[3],
double  sgn,
const bool  avail[3],
int  step[3],
double  w[3] 
)

The clip-and-redistribute allocation: which face neighbours a clipped cell's residue is pushed into, and in what fractions.

sgn = -1 pushes along -n (a LIQUID deficit goes to the liquid side), sgn = +1 along +n (a condensation excess goes to the gas side). step[d] is the integer face step (-1, 0 or +1) and w[d] the fraction; both are zero for n_d == 0.

The raw weight is n_d^2 (the direction cosines squared, so a plane normal to an axis puts everything on that axis), restricted to the neighbours avail[d] says can absorb it — for a liquid deficit, a neighbour that still HOLDS liquid. That restriction is not cosmetic: a slightly-tilted plane otherwise sends a n_t^2 share of the deficit into an in-plane neighbour that is already empty, which leaves a permanent negative colour wisp there (measured -2.5e-6 on the P1 Stefan ladder before the restriction; the transverse tilt came from the red-black smoother's parity asymmetry in the energy solve, so it is not avoidable upstream).

Returns false when NO neighbour can absorb the residue; the unrestricted weights are then returned, so the redistribution stays exactly conservative and the caller counts the event.

Definition at line 404 of file phase_change.hpp.

Referenced by peclet::flow::Solver< Grid >::pcRegress().

◆ pcGfmThetaK()

double peclet::flow::vof::pcGfmThetaK ( double  phiC,
double  nd,
double  s,
double  kappa,
double  thetaMin,
double  thetaMax 
)

WO-P3g item 3 — the curvature-consistent GFM distance.

pcGfmTheta measures cell i's distance to the interfacial neighbour j's PLIC PLANE. On a curved interface that is first order in h/R, exactly as it is for the one-sided fit (pcCurvedDistance), and WO-P3f measured the consequence on the row: the GFM flux carries an interface-curvature bias about TWICE the fit's (+12.4 % against +6.2 % at R = 20).

The correction is the same one, and it is exact to second order for a sphere. Cell i sits at integer offset d = -s e_d from the interfacial cell j, so its plane distance is phi_i = phi_c(j) - s n_d and its LATERAL offset from the plane's foot point is rho^2 = |d|^2 - (d.n)^2 = 1 - n_d^2 (the foot point differs from j's centre by a purely NORMAL displacement, so the tangential part of the offset is unchanged). The signed distance to the surface is then phi_i + kappa rho^2 / 4 with kappa = div(n) (-2/R for a gas sphere), and the distance ALONG the grid line is that divided by |n_d| to the same order.

kappa == 0 returns pcGfmTheta bitwise — that is the shipped default and the inertness contract of this rung.

Definition at line 446 of file phase_change.hpp.

References pcGfmTheta().

◆ pcGfmThetaKAniso()

double peclet::flow::vof::pcGfmThetaKAniso ( double  phiC,
const double  n[3],
int  d,
double  s,
double  kappa,
double  thetaMin,
double  thetaMax,
const VofMetric &  g 
)

The same on an ANISOTROPIC cell (Phase 3, flow/doc/anisotropic_vof.md §7, V5.3).

theta is the distance to the interface ALONG a grid line, measured in CELLS of that axis — an axis ratio, hence metric-free in form. What changes is that the caller now carries the PHYSICAL centre distance phiC and the PHYSICAL unit normal n, so both are pulled back to the axis before the ratio is formed:

index centre distance along `d`   :  phiC / h_d      (phiC is a length in hRef)
index normal component along `d`  :  n_d / h_d  ...  renormalized? NO — see below.

The clean way to see it: the plane's distance from cell i's centre, measured along the grid line, is |phi_i| / |n . e_d^phys| where phi_i is the PHYSICAL plane distance of that centre and e_d^phys is the UNIT physical direction of the axis — which is just e_d. So the numerator is |phiC - s h_d n_d| (the physical distance of the neighbour's centre, one cell of size h_d away) and the denominator is |n_d| h_d (the physical step per cell along d projected on the normal). The curvature correction's lateral offset is likewise physical: rho^2 = h_d^2 (1 - n_d^2).

At g.h = {1,1,1} every factor is 1.0 and the expression is pcGfmThetaK term for term.

Definition at line 476 of file phase_change.hpp.

Referenced by peclet::flow::Solver< Grid >::pcBudgetPost(), and peclet::flow::Solver< Grid >::pcBuildInterface().

◆ pcGfmRow()

PcGfmRow peclet::flow::vof::pcGfmRow ( double  theta,
bool  behind,
int  order 
)

◆ pcOperatorMassFlux()

double peclet::flow::vof::pcOperatorMassFlux ( double  q,
double  coefSum,
double  area,
double  hlv,
double  Rint 
)

WO-P3g item 1 — mdot from the ENERGY OPERATOR's own interfacial flux.

q is the heat the operator's Dirichlet rows actually draw INTO the interface from its pure face neighbours, q = sum_p a_Gamma(p) k_p o_p (T_p - T_Gamma) (positive = the surrounding phases are superheated and feed evaporation), coefSum = sum_p a_Gamma(p) k_p o_p, area the interfacial area of the cell and Rint the Schrage/IHTR resistance of T_Gamma = T_sat + mdot R_int. The interfacial energy balance mdot h_lv A = q then closes in ONE step:

mdot (h_lv A + R_int coefSum) = sum_p a_Gamma(p) k_p o_p (T_p - T_sat) .

Why this and not the least-squares fit: the heat the energy equation LOSES and the mass the regression PRODUCES become the same discrete number, so the flux consistency -q_gfm / (mdot h_lv A) is 1 by construction instead of the measured 0.95…1.04 (WO-P3f), and the mdot A products that drive the regression and the divergence source no longer contain the interfacial AREA at all — it cancels, dV = q dt/(h_lv rho_l).

Definition at line 546 of file phase_change.hpp.

Referenced by peclet::flow::Solver< Grid >::pcBuildInterface().

◆ pcIsInterfacial()

bool peclet::flow::vof::pcIsInterfacial ( double  c,
double  eps 
)

Is this a cell the regression acts on? The PLIC contract of advect_wy.hpp (wyIsMixed) with an explicit wisp guard: a cell whose colour is round-off residue has no meaningful plane, and giving it an area would inject an unbounded mdot into the census.

Definition at line 555 of file phase_change.hpp.

Referenced by peclet::flow::Solver< Grid >::pcBandDivergence(), peclet::flow::Solver< Grid >::pcBuildInterface(), peclet::flow::Solver< Grid >::pcUpdateThermalMask(), and peclet::flow::Solver< Grid >::vofInterfaceArea().

◆ csfKappaDefined()

bool peclet::flow::vof::csfKappaDefined ( double  branch)

Is the curvature of a cell carrying branch code branch usable?

True for branches 1..5 (the four shipped cascade tiers). FALSE for kCurvNone (0 — not an interfacial cell, where kappa is legitimately 0 and must not be read as a curvature) and for kCurvNoEstimate (6 — the cascade failed, where kappa is 0 for want of anything better). The branch field rides the ordinary field registry as a double, hence the rounded compare.

Definition at line 76 of file surface_tension.hpp.

Referenced by peclet::flow::Solver< Grid >::addCsfRhsCellInterp(), peclet::flow::vof::VofCurvature::clipPass(), csfFaceCurvature(), and vofCurvListEntry().

◆ csfFaceCurvature()

bool peclet::flow::vof::csfFaceCurvature ( double  k0,
double  b0,
double  k1,
double  b1,
double &  kf 
)

Face curvature for the balanced-force CSF, from the two cells sharing the face.

k0/b0 are the LOW cell (i - s_c) and k1/b1 the HIGH cell (i) — the two cells the face difference C(i) - C(i - s_c) reads, in that order. Returns false (and kf = 0) when neither cell carries a curvature; a caller that sees this on a face with ΔC ≠ 0 has an orphan face and should say so.

Definition at line 87 of file surface_tension.hpp.

References csfKappaDefined().

Referenced by peclet::flow::Solver< Grid >::addCsfRhs(), peclet::flow::addFaceAccelCsf(), peclet::flow::vof::VofBlockSet::buildCsfForce(), and peclet::flow::Solver< Grid >::csfDiagnostics().

◆ csfFaceForce()

double peclet::flow::vof::csfFaceForce ( double  sigma,
double  kf,
double  dC,
double  hGrad 
)

The balanced-force CSF contribution to the momentum RHS at one staggered velocity unknown, in the same units and at the same place as the incremental scheme's -(P(i) - P(i - s_c)).

dC = C(i) - C(i - s_c) is formed by the caller from the SAME colour field the density closure reads, with the projection's own face difference.

**hGrad is the PRESSURE-GRADIENT WEIGHT of this component's axis, and that is the whole point** (Phase 3, flow/doc/anisotropic_vof.md §5, decision V3.1). Substituting the Phase 1 conversions into the physical force sigma kappa dC/h_a — ‘sigma’ = sigma tRef^2/(rhoRef hRef^3),kappa' = kappa hRef,F_a' = F_a tRef^2/(rhoRef h_a)‘ — gives

F_a&rsquo;(i) = sigma' * kappa_f'(i) * (C(i) - C(i - s_a)) / h_a'^2,     h_a' = h_a/hRef

i.e. **exactly the weight ‘w_a = 1/h_a’^2the momentum RHS applies to-(P(i) - P(i - s_a))** (Phase 2,doc/anisotropic_metric.md§1.3). That identity is what keeps the balanced-force property on a stretched grid: with a constantkappa'the force is again the discrete (axis-weighted) gradient ofsigma' kappa' C`, so it lies in the range of the operator the projection inverts, the projection annihilates it, and a static droplet stays at machine zero. A different weight here than in the pressure gradient would put a floor under the spurious currents no curvature accuracy could remove.

hGrad is 1.0 in cell units and on every isotropic run (x / 1.0 == x), which is why the name change is not a numerical change: the argument used to be the uniform spacing, always 1.0.

Definition at line 129 of file surface_tension.hpp.

Referenced by peclet::flow::Solver< Grid >::addCsfRhs(), peclet::flow::addFaceAccelCsf(), peclet::flow::vof::VofBlockSet::buildCsfForce(), and peclet::flow::Solver< Grid >::csfDiagnostics().

◆ capillaryDt()

double peclet::flow::vof::capillaryDt ( double  rhoSum,
double  h,
double  sigma 
)

Brackbill capillary time-step limit sqrt((ρ₁+ρ₂) h³ / (4πσ)).

Returns +inf for σ <= 0 (no constraint). rhoSum is the SUM of the two phase densities — see the file header for why it is the sum and not a mean.

**h is the SMALLEST cell size** (Phase 3, decision V3.2): the constraint is the period of the shortest RESOLVABLE capillary wave, whose wavelength is 2 h_min — the finest axis sets it, and dt_sigma ~ h^{3/2} makes that the binding one. Basilisk takes the same min. On an isotropic grid min_a h_a is the one spacing, so nothing moves.

Definition at line 142 of file surface_tension.hpp.

Referenced by peclet::flow::Solver< Grid >::capillaryDtInternal().

◆ coxVoinovAngle()

double peclet::flow::vof::coxVoinovAngle ( double  thetaE,
double  ca,
double  logRatio,
double  thetaMin,
double  thetaMax 
)

Cox-Voinov / Afkhami-Zaleski-Bussmann grid-scale apparent angle, in RADIANS.

theta^3 = thetaE^3 + 9 Ca logRatio,      logRatio = ln(Delta/lambda)

Clamped into [thetaMin, thetaMax]; a non-positive cube (the film-entrainment branch of a fast receding line) returns thetaMin.

Definition at line 99 of file wetting_dynamic.hpp.

Referenced by vofDynamicContactAngle().

◆ vofHysteresisBase()

int peclet::flow::vof::vofHysteresisBase ( double  thetaApp,
double  thetaA,
double  thetaR,
double &  base 
)

The hysteresis selector: which base angle the Cox-Voinov correction is applied to, and whether the line is pinned.

All angles in RADIANS.

Definition at line 110 of file wetting_dynamic.hpp.

References kVofDynAdvancing, kVofDynPinned, and kVofDynReceding.

Referenced by vofDynamicContactAngle().

◆ vofDynamicContactAngle()

int peclet::flow::vof::vofDynamicContactAngle ( double  thetaApp,
double  uCl,
double  muL,
double  sigma,
double  logRatio,
bool  hyst,
double  thetaA,
double  thetaR,
double  thetaE,
double  thetaMin,
double  thetaMax,
double &  thetaOut,
double &  caOut 
)

The whole rung in one pure function: hysteresis selector, then Cox-Voinov on the selected base.

Parameters
thetaAppthe measured apparent angle (rad)
uClthe contact-line speed, + advancing
muL,sigmathe liquid viscosity and the surface tension (Ca_cl = muL uCl / sigma)
logRatioln(Delta/lambda); pass 0 to switch the Cox-Voinov correction off
hystwhether thetaA/thetaR are configured
thetaEthe static base angle (rad) used when hyst is false
thetaOutout: the angle to impose (rad)
caOutout: Ca_cl
Returns
the VofDynamicState

Definition at line 135 of file wetting_dynamic.hpp.

References coxVoinovAngle(), kVofDynPinned, kVofDynStatic, and vofHysteresisBase().

Referenced by peclet::flow::vof::VofDynamicWetting::impose().

◆ vofWallTangent() [1/2]

bool peclet::flow::vof::vofWallTangent ( const double  mfIn[3],
const double  nwIn[3],
double  tEps,
double  that[3],
double &  cosApp,
const VofMetric &  g 
)

The in-wall (tangential) unit direction of the interface and its apparent cosine.

mfIn is the fluid-side normal (into the gas), nwIn the wall normal (solid -> fluid); both any scale. Returns false when the interface is parallel to the wall (|t| < tEps), in which case that is zeroed and no contact-line direction is defined.

Definition at line 160 of file wetting_dynamic.hpp.

Referenced by peclet::flow::vof::VofDynamicWetting::measure(), and vofWallTangent().

◆ vofWallTangent() [2/2]

bool peclet::flow::vof::vofWallTangent ( const double  mfIn[3],
const double  nwIn[3],
double  tEps,
double  that[3],
double &  cosApp 
)

Unit-metric overload (isotropic cells) — the pre-Phase-3 signature, unchanged arithmetic.

Definition at line 192 of file wetting_dynamic.hpp.

References vofWallTangent().

◆ vofContactLineSpeed()

double peclet::flow::vof::vofContactLineSpeed ( const double  u[3],
const double  that[3] 
)

U_cl = +u .

t_hat: positive when the LIQUID ADVANCES.

This sign is a CORRECTION of the work order, which writes "sign positive when the liquid advances (velocity along `-t_hat`, since `m` points into the gas)". m points into the gas, so its in-wall part t_hat points from the LIQUID side towards the DRY side along the wall; the liquid advancing over dry wall is therefore motion along +t_hat, not -t_hat. Concretely, a vertical interface with liquid at x < x0 has m = +x_hat and t_hat = +x_hat, and a flow u = +U x_hat pushes the liquid onto the dry wall — advancing, Ca_cl > 0, which Cox-Voinov must turn into theta_D > theta_e (viscous bending RETARDS spreading; the Young force cos theta_e - cos theta_D then opposes the motion). With the work order's sign the same flow would report a receding line and Cox-Voinov would ACCELERATE the spreading — an unstable feedback, which is what the gate measures.

Definition at line 209 of file wetting_dynamic.hpp.

Referenced by peclet::flow::vof::VofDynamicWetting::measure().

Variable Documentation

◆ kWySweepPerm

constexpr int peclet::flow::vof::kWySweepPerm[6][3]
inlineconstexpr
Initial value:
= {{0, 1, 2}, {1, 2, 0}, {2, 0, 1},
{0, 2, 1}, {2, 1, 0}, {1, 0, 2}}

The six permutations of (x, y, z), cycled by step index so no direction is systematically favoured (perm[step % 6]).

Definition at line 130 of file advect_wy.hpp.

Referenced by peclet::flow::vof::WyAdvector::advect(), peclet::flow::vof::VofEnergyAdvector::advect(), peclet::flow::vof::MomentumConsistentAdvector::advect(), and peclet::flow::vof::VofBlockSet::advectBatched().

◆ kVofBlockBatch

◆ kVofHostExec

template<class Exec >
constexpr bool peclet::flow::vof::kVofHostExec = std::is_same_v<typename Exec::memory_space, Kokkos::HostSpace>
inlineconstexpr

The host branches select on this (§14.2).

Definition at line 112 of file block_batch.hpp.

◆ kVofTier3MaxTeams

constexpr int peclet::flow::vof::kVofTier3MaxTeams = 1 << 16
inlineconstexpr

Tier 3 of the batched cascade on a DEVICE (C3, doc/vof_step_performance_design.md §4.7, §5.11; launch shape: coordinator decision 2026-10-02, option C).

The interfacial counts stay on the device (§4.6), so the league cannot be the number of entries; instead a FIXED league of persistent warp-teams strides over the concatenated device-counted entries g = league_rank, league_rank + league_size, ..., each entry handled by curvFallbackTeam (the same body as VofCurvature::fallbackBatch's team path: neighbour compaction + canonical-order accumulation, design G §5.4, bit for bit the one-thread curvFallbackCell). No host read, no empty teams.

League size: one team per hardware warp slot, concurrency() / kVofWarp (on CUDA = SMs x resident threads per SM / 32), capped by kVofTier3MaxTeams. Teams beyond the resident number simply start later; the stride covers every entry whatever the league.

Definition at line 1303 of file block_batch.hpp.

Referenced by vofCurvFallbackTeams().

◆ kVofStats

constexpr int peclet::flow::vof::kVofStats = 14
inlineconstexpr

Doubles per job in the statistics packet: volume, centroid[3], velocity[3], area, moment[6].

Definition at line 1522 of file block_batch.hpp.

Referenced by peclet::flow::vof::VofBlockSet::advectBatched(), peclet::flow::vof::VofBlockSet::applyStats(), peclet::flow::vof::VofBlockSet::flushStats(), peclet::flow::vof::VofBlockSet::measureBatched(), and vofBatchStats().

◆ kVofWarp

constexpr int peclet::flow::vof::kVofWarp = 32
inlineconstexpr

The team size of a device tier-3 launch: one warp / wavefront.

Definition at line 376 of file curvature_field.hpp.

Referenced by peclet::flow::vof::VofCurvature::fallbackBatch(), and vofCurvFallbackTeams().