flow — VoF rung V2b (WO-K): momentum-consistent transport of rho^c u_c on the half-shifted (MAC) momentum control volumes, driven by the SAME geometric fluxes as the colour advection of the same step. More...

#include <Kokkos_Core.hpp>
#include <stdexcept>
#include "policy.hpp"
#include "vof/advect_wy.hpp"
#include "vof/cutcell.hpp"

Classes
class	peclet::flow::vof::MomentumConsistentAdvector
	Momentum-consistent companion to `WyAdvector`. More...

struct	peclet::flow::vof::MomentumConsistentAdvector::Diagnostics

Namespaces
namespace	peclet

namespace	peclet::flow

namespace	peclet::flow::vof

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

double	peclet::flow::vof::vofMinmod (double a, double b)

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.

Detailed Description

flow — VoF rung V2b (WO-K): momentum-consistent transport of rho^c u_c on the half-shifted (MAC) momentum control volumes, driven by the SAME geometric fluxes as the colour advection of the same step.

Why

With mass and momentum advected by different fluxes, a mixed cell multiplies the gas acceleration by the liquid density — a spurious interfacial momentum source of order d(rho). Rudman (1998) and Arrufat et al. (Computers & Fluids 215:104785, 2021) show the naive scheme breaks down around ratio 1000 unless the resolution is absurd; with consistency a raindrop at ratio 831.8 is accurate within 15 % at 15 cells/diameter instead of ~200.

Indexing

The momentum control volume of component e is the pressure cell shifted half a cell in e. Everything here is indexed in **flow's own low-face convention**: CV_e(i) is the control volume centred on the -e face of cell i, spanning [x_{i-1}, x_i] (the +e half of cell i-1 plus the -e half of cell i), which is exactly where the solver's velocity unknown u_e(i) lives. Three things follow and all three are load-bearing:

the owned control volumes are exactly the block's INNER region, so C^e and rho^e u_e are ordinary cell fields on the advector's g = 3 block and WyAdvector::exchange carries their ghosts unchanged (a +e-shifted ownership would have put an owned value inside the halo's ghost band, where the exchange would overwrite it);
the bridge back to the solver is a plain copyInner, with no shift to get wrong;
WyAdvector's own face velocity uf_e(i) sits on the +e face of cell i (colour_field.hpp), so the velocity AT CV_e(i) is uf_e(i - s_e).

The construction

**C^e by clipping, never by interpolation.** C^e(i) = slab(cell i-s_e, e, [1/2, 1]) + slab(cell i, e, [0, 1/2]), each term the PLIC polyhedron of that pressure cell truncated to the half cell (plicSlabVolume — a rescale, not a clipper). Interpolating C destroys sharpness and is the classic error in this rung.
The same fluxes. For a sweep in direction d the flux slab of the shifted CV is
- d == e: a slab of the SINGLE pressure cell whose centre the CV face sits on (|a| < 1/2 keeps it inside that cell's half);
- d != e: the union of two HALF cells — [1/2,1] of cell p-s_e and [0,1/2] of cell p, each |a| deep in d. Two axes at once, which is why plicBoxVolume exists. Both read the planes WyAdvector reconstructed for THIS sweep, between its computeFluxes and its applySweep — the planes are overwritten every sweep, so the momentum sweeps are interleaved with the colour sweeps rather than run after advect() returns. The face Courant number of the shifted CV is the e-average of the two staggered faces; summed over the three sweeps that is exactly 1/2 (div_{p-s_e} + div_p), so the dilation term vanishes with the projection's own divergence residual — the colour's conservation floor.
Weymouth's flux bounds, on the shifted volume's OWN colour. The geometric flux above is bounded by the liquid content of the donor CV as the current cell planes see it, which is not the same as the ADVECTED C^e once the second sweep re-reconstructs from an updated C. The gap is O(a^2) and does not accumulate — measured on a tilted plane at CFL 0.2/0.1/0.05/ 0.025: C^e reached -2.6e-2 / -6.0e-3 / -1.0e-3 / -3.1e-5 — but at density ratio 10^4 a -2.6e-2 undershoot drives rho^e to -255, and the recovery divides by it. The flux is therefore clamped into Weymouth's own admissible interval max(0, |a| - (1 - C^e_don)) <= |F| <= min(|a|, C^e_don) (thesis Appendix A), which is exactly the hypothesis of his boundedness proof and restores 0 <= C^e <= 1. The clamp is inactive wherever the geometry is self-consistent (clampedFluxes() counts it), it is applied to the ONE face value both neighbours share so conservation still telescopes bit-exactly, and it is applied BEFORE the momentum flux is formed so the two updates keep sharing one F — the consistency identity is untouched. For C^e_don exactly 0 or 1 the interval collapses to the algebraic value, so full and empty control volumes stay exactly stationary.
One frozen flag for the pair. flag^e = H(C^{e,n} - 1/2), frozen once per step from the freshly built C^e and used unchanged by all three sweeps in BOTH the C^e update and the rho^e u_e update. That shared flag is what makes the two updates telescope identically. (The work order says "the same frozen dilation flag as the C advection". The flag that has to be shared is the one of the pair that must telescope, and both members of that pair live on the SHIFTED CV, so the structural analogue H(C^e - 1/2) is what is shared — WyAdvector::dilationFlag(), the pressure-cell flag, belongs to the pressure-cell pair. useCellDilationFlag switches to the literal reading as a measured ablation.)
Recovery u_e = (rho^e u_e) / max(rho^e, floor); see rhoFloor.

The consistency identity (why the uniform-velocity gate is exact)

Per sweep, with Psi = rho_g a + (rho_l - rho_g) F the mass flux and rho^ the frozen-flag density, the conservative pair is

C^e_i   += dC = (F_- - F_+) + flag_i (a_+ - a_-)
M_i     += (Psi_- u^_- - Psi_+ u^_+) + rho^_i u_i (a_+ - a_-)

with rho^e = rho_g (1 - C^e) + rho_l C^e, so d(rho^e) = drho dC — the SAME dC. If u_e == U everywhere then every u^ and every u_i is U and the two updates differ by exactly the factor U: dM = U d(rho^e), so M = U rho^e is preserved and u = M/rho^e = U at ANY density ratio, with no tuning knob. An inconsistent scheme is wrong at O(d rho).

The identity holds whatever the fluxes are, provided the two updates use the same ones. That is the whole point of driving both from one set of planes, and it is why the flux clamp of point 3 is free: it changes F, not the fact that there is one F. momentumUpdate evolves the algebraically identical u_new = u_old + dev/rho_new form so that the identity survives in FLOATING POINT too, not only in exact arithmetic — the note there explains why that matters.

The transported velocity in the flux

Cut cells (rung V5a, WO-Q item 8)

With an immersed solid the half-shifted control volume gets a fluid volume and open face areas of its own, by the same half-cell averaging that defines it:

fluid volume eps^e(i) = 1/2 (eps(i - s_e) + eps(i));
transverse face (d != e) openness 1/2 (o_d(i - s_e) + o_d(i)) — the average is what makes the dilation sum reproduce 1/2 (div_open(i - s_e) + div_open(i)), i.e. the same half-and-half of the projection's own constraint the uncut construction relies on;
AXIAL face (d == e) openness eps(p), the CELL fraction of the pressure cell whose CENTRE the face sits on — the axial CV face IS that cell-centre plane, and its fluid area fraction is the cell's fluid volume fraction to the same (uniform-solid-in-cell) accuracy as everything else here. The liquid content of the CV is eps(i - s_e) slabHi(i - s_e) + eps(i) slabLo(i) and C^e = that / eps^e: each half cell's PLIC slab fraction weighted by ITS OWN fluid volume, which is the same "spread the solid uniformly over the cell" approximation the colour flux makes, and which returns C^e = 1 exactly for a uniformly liquid field whatever the geometry. Both updates are then done in fluid-volume units (divide by eps^e), exactly as the colour is.

The consistency identity survives all of it verbatim, and that is the point: every term of dev is a DIFFERENCE OF VELOCITIES, so a uniform field gives dev = 0 bit for bit whatever eps^e and the openness are, and u_new = u_old exactly — in cut cells too, not only away from them. Measured on a sphere packing at ratio 1e3: 0.0.

u^ is the CURRENT transported velocity of the donor CV (plain donor-cell upwind by default; momentumMuscl adds a MinMod-limited slope to the flux slab's centroid, and the note on that flag says why it is NOT the default at high density ratio). With a uniform field u^ = U bit for bit either way, so the gate above is untouched. The velocity is exchanged between sweeps alongside C^e — fluxing the FROZEN u^n instead would be one exchange cheaper but leaves a term with gain rho_l/rho_g on any control volume a sweep empties (see the derivation in momentumUpdate), which is not something to carry into a high-ratio rung to save a halo exchange.

Definition in file momentum_advect.hpp.

Classes

Namespaces