core — VoF rung V5b (WO-S): the theta-consistent solid-band fill (static contact angle). More...

#include <Kokkos_Core.hpp>
#include <Kokkos_MathematicalFunctions.hpp>
#include "peclet/core/vof/curvature.hpp"
#include "peclet/core/vof/plic.hpp"

Namespaces
namespace	peclet

namespace	peclet::core

namespace	peclet::core::vof

Enumerations
enum	peclet::core::vof::VofWettingPivot : int { peclet::core::vof::kVofPivotVolume = 0 , peclet::core::vof::kVofPivotInterface = 1 , peclet::core::vof::kVofPivotWallNormal = 2 , peclet::core::vof::kVofPivotContactLine = 3 }
	How the theta-plane is anchored in the fluid cell. More...

enum	peclet::core::vof::VofWettingBranch : int { peclet::core::vof::kVofWetNone = 0 , peclet::core::vof::kVofWetTheta = 1 , peclet::core::vof::kVofWetNeighbour = 2 , peclet::core::vof::kVofWetPure = 3 , peclet::core::vof::kVofWetParallel = 4 , peclet::core::vof::kVofWetNeutral = 5 , peclet::core::vof::kVofWetCount = 6 }
	Which rule produced a band cell's colour, for `contact_angle_diagnostics()`. More...

Functions
bool	peclet::core::vof::youngsNormalFluidOnly (const double c[27], const unsigned char fluid[27], double m[3])
	Youngs' 27-point weighted gradient restricted to non-solid cells.

void	peclet::core::vof::vofPlicCentroid (const double m[3], double alpha, double p[3])
	The centroid of the PLIC polygon of cell `f` in the cell's own unit-cube frame, for the plane `m .`

int	peclet::core::vof::vofWettingPlane (const double mfIn[3], double cf, const double nwIn[3], double cosT, double sinT, double sdfF, int pivot, double tEps, double mth[3], double &alphaTh, double &cosApp, const VofMetric &g)
	Build the theta-plane in the FLUID cell's own unit-cube frame.

int	peclet::core::vof::vofWettingPlane (const double mfIn[3], double cf, const double nwIn[3], double cosT, double sinT, double sdfF, int pivot, double tEps, double mth[3], double &alphaTh, double &cosApp)
	Unit-metric overload (isotropic cells) — the pre-Phase-3 signature, unchanged arithmetic.

double	peclet::core::vof::vofWettingFraction (const double mth[3], double alphaTh, const int ds[3])
	The liquid fraction of the solid cell at integer offset `ds` from the fluid cell, under the theta-plane returned by `vofWettingPlane`.

Detailed Description

core — VoF rung V5b (WO-S): the theta-consistent solid-band fill (static contact angle).

Container-free (the plic.hpp rule: scalars and small arrays only, no View, no indexing), so every rule below is a pure function the device kernels in advect_wy.hpp call and a host oracle can call verbatim. It replaces ONLY pass 1 of WO-Q's three-pass solid-band fill; passes 2-3 (the mean of already-filled neighbours) and the shrinking depth budget are unchanged.

What the fill has to do

The MYC 3^3 stencil and the V3 height-function columns of a cell next to a wall reach INTO the solid. WO-Q's neutral rule (the mean of the fluid face neighbours) continues the colour into the wall with zero slope — the 90 deg limit of the Afkhami & Bussmann height-function boundary condition (IJNMF 57:453, 2008). To impose a contact angle theta we instead write into the band the fractions of the PLANE that continues the fluid-side interface into the solid at angle theta. The unmodified V3 curvature cascade then returns the curvature of an interface that meets the wall at theta, and the V4 balanced force does the rest — no force is added at the wall.

Sign conventions (all three are load-bearing; they are checked in <tt>test_vof_wetting.cpp</tt>)

n_w = grad(sdf)/|grad(sdf)| points from the SOLID INTO the fluid (sdf > 0 is fluid).
the PLIC normal m points from the LIQUID INTO the gas, liquid side m . x < alpha (plic.hpp).
theta is measured THROUGH THE LIQUID, and at the contact line **m . n_w = cos(theta)**. theta = 0 (complete wetting) gives m = n_w: the liquid side m . x < alpha is then the SOLID side, i.e. a liquid film wets the wall and the band fills with liquid. theta = 180 gives m = -n_w and an empty band. theta = 90 gives m perpendicular to the wall — the interface standing on the wall, which is the neutral case.

The construction

For a solid band cell s, walk along n_w(s) to the first FLUID cell f, then

take the fluid-side normal m_f from a fluid-only Youngs gradient (youngsNormalFluidOnly) — the 27-point Youngs weights restricted to cells that are not solid and renormalized per half-plane. Restricting is what breaks the circularity: the ordinary MYC normal of f reads the very band cells this fill is writing, so using it would make the fill a fixed-point iteration on its own output. (Basilisk's contact.h breaks the same loop with a second reconstruction pass; that is the WO-S fallback if the gate fails.)
keep from m_f only its component IN the wall, t = m_f - (m_f . n_w) n_w, and build **m_theta = cos(theta) n_w + sin(theta) t_hat** — the unique unit vector in the (n_w, t_hat) plane making angle theta with n_w on the same side as m_f. When |t| is below tEps the interface is parallel to the wall and no rotation is defined, so m_f is kept.
anchor the plane, and
C_s = the fraction of cell s's unit cube on the liquid side of it.

Only the IN-WALL half of <tt>m_f</tt> is used, and that is deliberate

A fluid-only estimator cannot measure the WALL-NORMAL component of the interface normal: below the first fluid row there is no colour to difference against, and the one-sided substitute reads a saturating profile over half the distance. Measured on EXACT plane fractions next to a flat wall (gate G0e): the fluid-only normal's angle to the wall is wrong by up to 23 deg at theta = 30, against 2.3 deg for the same Youngs stencil with the solid rows present. The IN-WALL components are two-sided and complete, and their azimuth is accurate to **~1 deg** over the same sweep. Since the whole point of the rung is to OVERWRITE the apparent angle with the prescribed one, the construction above uses the azimuth (which is measured well) and discards the wall-normal component (which is not). This is what makes the rung's accuracy independent of the estimator defect the fill would otherwise inherit.

The anchor's own column is not enough (measured)

The construction as the work order states it walks from the band cell to the first fluid cell ALONG n_w and uses that cell's interface. Where that cell is PURE phase the rule degenerates to C_s = C_f, and that is wrong exactly where it matters: at a wetting angle the continued interface reaches into the solid BEYOND the column where the contact line crosses the first fluid row, so the columns just outside the contact circle have a pure-gas anchor and receive nothing, while the true continuation puts liquid there. MEASURED on an equilibrium theta = 60 cap (D/dx = 20, wall on a cell face): the band column one cell outside the contact column reads 0 where the continued plane gives 0.63 and 0.05 in the two rows below the wall, so the height function reads a wall slope closer to 90 deg than to 60 — and the equilibrium contact angle of the whole run comes out 5 deg biased towards 90.

The fix, shipped: when the anchor is pure phase, average the theta-plane fractions of the anchor's MIXED fluid neighbours (its 3^3 stencil, so at most one column away) evaluated in cell s; with no mixed neighbour the pure continuation stands. Every read stays inside the halo the exchange has already made owner-consistent, so the fill is still decomposition-independent.

The anchor: a measured deviation from the work order

WO-S anchors the plane at the contact point c = p_f - sdf(p_f) n_w, i.e. the fluid cell's PLIC centroid projected onto the wall ALONG n_w. That point is on the wall but it is not on the interface plane unless the interface is perpendicular to the wall: projecting along n_w changes the plane's offset by -sdf(p_f) cos(theta), an O(1) shift of the interface (up to ~1 cell at theta = 30 deg) with the wrong sign — it removes liquid from the band for a wetting angle. The consequence is that the fill is not idempotent: an interface already meeting the wall at exactly theta is not reproduced, theta is not a fixed point of the scheme, and the equilibrium angle is biased. Measured (gate G0a/G0b): the band fraction is off by up to 0.26 at theta = 60, against 1e-15 for the shipped rule.

The shipped default is kVofPivotVolume: the plane with normal m_theta whose liquid volume in cell f is exactly C_f (plicAlpha). It needs no pivot point at all, it is volume-consistent with the cell it is anchored on, it is exactly idempotent (if the true interface is a plane at angle theta then plicAlpha returns its offset exactly, so the band receives the exact fractions of the continued interface — that identity is gate G0a), and it is the only variant that does not read the poorly-measured wall-normal part of m_f at all. kVofPivotInterface (the plane through the PLIC centroid p_f reconstructed from m_f — the Afkhami-Bussmann / Basilisk contact.h rule), kVofPivotWallNormal (WO-S as written) and kVofPivotContactLine ship as measured ablations (set_contact_angle_pivot).

Known approximation, consistent with WO-Q. C_f is the liquid fraction of the cell's FLUID volume, but the reconstruction is done on the WHOLE unit cube (exactly as WO-Q's flux is), so in a CUT anchor cell the plane is displaced by roughly (1 - eps_f)/2 times the wall-normal component of m_theta — 0.14 cells at theta = 60 on a half-cut cell, zero at theta = 90 and zero in an uncut cell. Removing it needs the solid-clipped volume relation (Huang, JCP 2025/2026), the same refinement WO-Q's flux defers.

**Layer L1 (VOF_PLAN.md §11): promoted from flow/src/vof/ to core (WO-W0, 2026-09-02).** These are the container-free VoF math kernels — one copy shared by every container: flow's structured colour field, the per-bubble block container (Part III) and the future AMR path. flow/src/vof/wetting.hpp is now a thin include + a using-directive, so every existing peclet::flow::vof:: spelling keeps working unchanged.

Definition in file wetting.hpp.

Namespaces

Enumerations