advect_wy.hpp File Reference

flow — VoF rung V1: Weymouth & Yue (2010) directionally-split, exactly volume-conserving geometric advection of a colour field on a structured block. More...

#include <chrono>
#include <cstdio>
#include <functional>
#include <Kokkos_Core.hpp>
#include <stdexcept>
#include "mac_stencils.hpp"
#include "policy.hpp"
#include "vof/cutcell.hpp"
#include "vof/plic.hpp"
#include "vof/wetting.hpp"

Go to the source code of this file.

Classes
class	peclet::flow::vof::WyAdvector
	Weymouth-Yue split advection of a colour field on an extended (inner + ghost) block. More...
struct	peclet::flow::vof::WyAdvector::Diagnostics
	Per-step census. More...
struct	peclet::flow::vof::WyAdvector::Timing

Namespaces
namespace	peclet
namespace	peclet::flow
namespace	peclet::flow::vof

Typedefs
using	peclet::flow::vof::UCField = Kokkos::View< unsigned char *, SMem >
using	peclet::flow::vof::LField = Kokkos::View< long *, SMem >

Functions
bool	peclet::flow::vof::wyIsMixed (double c)
	The mixed-cell predicate.
bool	peclet::flow::vof::wyIsMixed (double c, double eps)
	WO-R2 item 4 — the same predicate with a WISP TOLERANCE.
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.
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.
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.
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.

Variables
constexpr int	peclet::flow::vof::kWySweepPerm [6][3]
	The six permutations of (x, y, z), cycled by step index so no direction is systematically favoured (perm[step % 6]).

Detailed Description

flow — VoF rung V1: Weymouth & Yue (2010) directionally-split, exactly volume-conserving geometric advection of a colour field on a structured block.

Standalone. This header owns its own extended field block with its own ghost width (g = 3, suite/docs/VOF_PLAN.md §3 design rule 1) and knows nothing about IbmSolver, the solver's G = 2 machinery, MPI, or cut cells. The ghost refresh is a host-side callback (exchange), so the same header drives the single-block periodic case in tests/kokkos and the multi-rank core::halo::GridHaloTopology case in tests/kokkos_mpi with no code change.

The algorithm (Weymouth & Yue, JCP 229:2853, 2010)

Primary source used here: G. D. Weymouth, Physics and learning based computational models for breaking bow waves based on new boundary immersion approaches, PhD thesis, MIT 2008 (dspace 1721.1/44754) — its §2.2.2 and Appendix A are the JCP paper's own derivation and the full boundedness proof. Cross-checked against the restatement in Arrufat et al., Computers & Fluids 215:104785 (2021), arXiv:1811.12327, §3.3.4.

Split the colour equation dC/dt + div(u C) = C div(u) into three one-dimensional sweeps

C^{l+1}_i = C^l_i - (F_{i+1/2} - F_{i-1/2}) + c_i (a_{i+1/2} - a_{i-1/2}) ,
a_{i+1/2} = uf_{i+1/2} dt / h  (the face Courant number),

where F is the Eulerian donor-cell geometric flux: the fluid volume of the upwind cell's PLIC polyhedron inside the un-stretched slab of width |a| adjacent to the face (thesis Fig. 2-4; plicSlabVolume). The interface is re-reconstructed before every sweep from the current C (thesis: "to avoid the possibility of fluxing the same fluid into two different cells the surface must be reconstructed after each sweep").

The dilation coefficient is the frozen cell-centre indicator

c_i = H(C^n_i - 1/2)     (thesis eq. A.28: 1 if C^n_i > 1/2, else 0)

evaluated ONCE per step from C^n and used unchanged by all three sweeps. This is the whole method. Because c_i is a constant of the step and direction-independent, summing the three sweeps gives

C^{n+1}_i - C^n_i = -sum_faces F + c_i (dt/h) sum_d (uf_{d+} - uf_{d-}) ,

whose flux part telescopes to zero over the domain and whose second part is c_i times the discrete divergence — zero for a discretely solenoidal face field. Hence exact conservation, to the accuracy with which the discrete face divergence vanishes (thesis §2.2.2 requirements 1–3). Recomputing c_i between sweeps destroys exactly this cancellation, and the damage is quiet: conservation degrades from ~1e-15 to ~1e-10 rather than failing outright.

Two further consequences of the same structure, both relied on here:

• a cell whose whole 1D neighbourhood is full (C = 1) is exactly stationary in floating point: its fluxes are the algebraic 1 * a_{+/-} and the dilation term is the negation of the flux difference, and IEEE subtraction is antisymmetric — so (a_- - a_+) + (a_+ - a_-) is an exact zero. Only interface cells accumulate round-off, which is why the conservation floor scales with the interface area and not with the domain volume. To keep this, the flux and the dilation term must scale the SAME uf by the SAME dt/h (see applySweep).
• empty cells (C = 0) stay exactly empty.

Boundedness needs a CFL cap. Thesis Appendix A bounds the flux by max(0, a - C) <= F <= min(a, C) and shows that c = H(C - 1/2) is the only quadrature that prevents both over-filling and over-emptying, provided |a| <= 1/(2(N-1)) for N-dimensional flow (thesis eq. A.33 / eq. 2.23) — 1/2 in 2D, 1/4 in 3D. The widely-quoted CFL < 0.5 (including in this rung's work order) is the 2D value; quoting it for a 3D solver is a transcription error the WO inherited. cflLimit therefore defaults to the proven 3D bound 0.25 and advect() aborts at or above it; raise it deliberately (cflLimit = 0.5) for 2D work or to probe the gap, never as a way to take bigger steps in 3D. Note that conservation is independent of boundedness: the telescoping above holds whatever C does, so an over-CFL run loses 0 <= C <= 1, not volume — which is exactly why the failure is easy to miss (volume still closes to round-off while C leaves [0,1]). Empirically this margin is not tight — a sweep to CFL 0.48 on the LeVeque field never left [0,1] — but the target regime (capillary-limited dt, VOF_PLAN.md §4 V4) sits far below both bounds, so there is nothing to buy by defaulting past what is proven.

No clipping is applied at this rung (VOF_PLAN.md §4 V1): conservation must close to round-off with nothing hiding the error. diagnostics() reports the wisp census instead.

Definition in file advect_wy.hpp.