openswmm::fv::kernels Namespace Reference

Classes
struct	FaceFlux
	Result of a face flux evaluation. More...
struct	FaceState
	One side of a face after hydrostatic reconstruction. More...

Functions
OPENSWMM_KERNEL_FN double	sectionArea (const FvGeometry &g, double h) noexcept
OPENSWMM_KERNEL_FN double	sectionWidth (const FvGeometry &g, double h) noexcept
OPENSWMM_KERNEL_FN double	sectionHydRad (const FvGeometry &g, double h) noexcept
OPENSWMM_KERNEL_FN double	slotRamp (double s) noexcept
OPENSWMM_KERNEL_FN double	slotRampIntegral (double s) noexcept
OPENSWMM_KERNEL_FN double	areaOfDepth (const FvGeometry &g, double h) noexcept
OPENSWMM_KERNEL_FN double	widthOfDepth (const FvGeometry &g, double h) noexcept
	Top width dA/dh at depth h, INCLUDING the tapered slot.
OPENSWMM_KERNEL_FN double	hydRadOfDepth (const FvGeometry &g, double h) noexcept
	Hydraulic radius at depth h.
OPENSWMM_KERNEL_FN double	i1OfDepth (const FvGeometry &g, double h, double area_at_h) noexcept
	Hydrostatic first moment I₁(h) = ∫₀ʰ A(η)dη.
OPENSWMM_KERNEL_FN double	depthOfAreaBracketed (const FvGeometry &g, double a) noexcept
	Invert A → h — the EXACT inverse of areaOfDepth above.
OPENSWMM_KERNEL_FN double	depthOfArea (const FvGeometry &g, double a) noexcept
	Invert A → h. Same root as depthOfAreaBracketed, found far faster.
OPENSWMM_KERNEL_FN double	celerity (double a, double t) noexcept
OPENSWMM_KERNEL_FN void	waveSpeeds (const FaceState &L, const FaceState &R, double &sl, double &sr) noexcept
OPENSWMM_KERNEL_FN void	physicalFlux (const FaceState &S, double &fa, double &fq) noexcept
	Physical flux of one state.
OPENSWMM_KERNEL_FN FaceFlux	riemannFlux (const FaceState &L, const FaceState &R) noexcept
	Flux for the conservative St. Venant system, plus the contact speed.
OPENSWMM_KERNEL_FN double	speciesFlux (const FaceState &L, const FaceState &R, const FaceFlux &f, double phi_l, double phi_r, bool hllc) noexcept
	Species flux for the A(phi) component.
OPENSWMM_KERNEL_FN double	frictionUpdate (double q, double u, double r, double dt, double rough_factor) noexcept
	Semi-implicit Manning friction — unconditionally stable, so friction imposes no time-step restriction.
OPENSWMM_KERNEL_FN double	localLossUpdate (double q, double u, double k, double dx, double dt) noexcept
	Semi-implicit entrance/exit loss at a node-coupling face.
OPENSWMM_KERNEL_FN double	ufUpdate (double q, double a, double u_old, double k3, double grad_term, double dt) noexcept
	Unsteady-friction momentum update (issue #156).
OPENSWMM_KERNEL_FN double	tpaDepthOfArea (const FvGeometry &g, double a) noexcept
OPENSWMM_KERNEL_FN double	tpaAreaOfDepth (const FvGeometry &g, double h) noexcept
OPENSWMM_KERNEL_FN double	tpaI1OfDepth (const FvGeometry &g, double h) noexcept
OPENSWMM_KERNEL_FN double	faceCflDt (double cfl, double dx, double u, double c) noexcept
	CFL-limited step for one face: α·Δx/(|u| + c).
OPENSWMM_KERNEL_FN double	positivityScale (double vol, double outflow, double dt) noexcept
	Positivity scale for a cell about to export more volume than it holds.

Variables
constexpr double	kGravity = 32.2
constexpr double	kDryDepth = 1.0e-7
constexpr double	kDryArea = 1.0e-12
	Area floor paired with kDryDepth for the u = Q/A division.
constexpr double	kEtaDeadband = 1.0e-12

Function Documentation

areaOfDepth()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::areaOfDepth ( const FvGeometry & g, double h )

noexcept

Flow area at depth h, INCLUDING the tapered slot. Monotone in h for any section, which is what makes the depth inversion well posed.

Here is the call graph for this function:

Here is the caller graph for this function:

celerity()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::celerity ( double a, double t )

noexcept

Gravity-wave celerity √(g·A/T). Guarded so a vanishing top width (dry, or a section evaluated exactly at a cusp) cannot produce a non-finite dt bound.

depthOfArea()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::depthOfArea ( const FvGeometry & g, double a )

noexcept

Invert A → h. Same root as depthOfAreaBracketed, found far faster.

This is the solver's hottest kernel by a wide margin — profiling a Δx = 20 ft run put it and the closure evaluations it drives at 87 % of total time — so how it converges matters more than anywhere else in the scheme.

Why the obvious approach is slow. Illinois regula-falsi on the depth-uniform bracket does not converge superlinearly here: measured 16 closure evaluations per call on a circular pipe and 35 on a trapezoid, with the iteration count falling only linearly as the tolerance is relaxed — the signature of bisection. Seeding it better changes nothing, because the exit test is the BRACKET collapsing and Illinois replaces only one end per step.

Why not Newton. The width is dA/dh analytically, and Newton with it needs 3–5 evaluations. But for tabulated shapes W and A are INDEPENDENT legacy tabulations rather than an exact derivative pair (§7A.2), and the error that introduces is not small: measured round-trip error 4.7e-4 ft on a 3 ft circular pipe, five orders worse than the scheme needs and enough to break lake-at-rest.

Brent. Inverse quadratic interpolation with a secant fallback and a bisection safeguard — superlinear using function values ONLY, so the width inconsistency cannot mislead it. 5.9 / 3.1 / 6.3 evaluations on circular / rectangular / trapezoidal at full round-trip accuracy (≤ 2.7e-15 ft).

The bracket comes from the area-uniform inverse table widened by one panel each side and is then VERIFIED by evaluating both ends. The two evaluations that costs are why a rectangular pipe — which the old path nailed in 1.5 — now takes 3.1; that trade is worth it, and the guard falls back to the bracketed inverse if the table's bracket somehow fails to contain the root.

Here is the call graph for this function:

depthOfAreaBracketed()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::depthOfAreaBracketed ( const FvGeometry & g, double a )

noexcept

Invert A → h — the EXACT inverse of areaOfDepth above.

This must be a true inverse, not merely an accurate one. The solver carries A and derives the free surface as η = z_b + depthOfArea(A); if the composition depthOfArea∘areaOfDepth were not the identity, cells sitting at the same η but different bed elevations would reconstruct DIFFERENT surfaces, and lake-at-rest — the property the whole well-balanced construction exists to deliver — would fail on every partly-full pipe.

xsect::getYofA cannot be used for this. The legacy geometry tables are INDEPENDENT tabulations of the same shape — A_Circ gives area from depth, Y_Circ gives depth from area — and they round-trip only to table resolution (measured: 0.016 ft on a 3 ft circular pipe near the crown, ~0.5 % of the diameter). That is fine for the legacy solver, which never composes them; it is fatal here.

So the inverse is built from the forward closure itself: bracket on the A samples stored alongside the I₁ table (both were produced by areaOfDepth, so the bracket is exact), then Illinois-regula-falsi inside the panel. A is strictly increasing — a monotone section plus a strictly increasing slot term — so the bracket always holds and convergence is unconditional; within one panel A is nearly linear, so it typically takes four or five evaluations.

Above the crown A is exactly linear, so that branch is closed-form.

Here is the call graph for this function:

Here is the caller graph for this function:

faceCflDt()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::faceCflDt ( double cfl, double dx, double u, double c )

noexcept

CFL-limited step for one face: α·Δx/(|u| + c).

frictionUpdate()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::frictionUpdate ( double q, double u, double r, double dt, double rough_factor )

noexcept

Semi-implicit Manning friction — unconditionally stable, so friction imposes no time-step restriction.

Q^{n+1} = Q* / (1 + Δt·g·(n/φ)²·|u|/R^{4/3}), with g·(n/φ)² supplied as rough_factor (ConduitData::rough_factor, already adjusted for Courant lengthening). R^{4/3} is written as R·cbrt(R) to avoid a libm pow().

hydRadOfDepth()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::hydRadOfDepth ( const FvGeometry & g, double h )

noexcept

Hydraulic radius at depth h.

Here is the call graph for this function:

i1OfDepth()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::i1OfDepth ( const FvGeometry & g, double h, double area_at_h )

noexcept

Hydrostatic first moment I₁(h) = ∫₀ʰ A(η)dη.

Below the crown this reads the per-geometry table built at init (uniform on [0, y_full]) and refines with one trapezoid step using the exact A(h) the caller needs anyway — second-order accurate and, crucially, a single-valued function of h, which is all the well-balanced property requires.

Above the crown A is exactly linear (A = a_crown + t_slot·(h − y_full)), so the extension is analytic — deep surcharge stays exact with a small table.

localLossUpdate()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::localLossUpdate ( double q, double u, double k, double dx, double dt )

noexcept

Semi-implicit entrance/exit loss at a node-coupling face.

The local head loss K·u²/2g is spread over the end cell as a momentum sink g·A·S_loss with S_loss·Δx = K·u²/2g, then integrated implicitly in |u| for the same stability reason as friction. K = 0 leaves Q untouched, so calibrated models that never set losses are bit-unaffected.

physicalFlux()

OPENSWMM_KERNEL_FN void openswmm::fv::kernels::physicalFlux ( const FaceState & S, double & fa, double & fq )

noexcept

Physical flux of one state.

Here is the caller graph for this function:

positivityScale()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::positivityScale ( double vol, double outflow, double dt )

noexcept

Positivity scale for a cell about to export more volume than it holds.

Outgoing fluxes may take at most the cell's available volume over Δt; faces leaving the cell are scaled by λ and the identical scaled flux updates both incident cells, so conservation is untouched. Same contract as the 2D marcher's positivityScale.

riemannFlux()

OPENSWMM_KERNEL_FN FaceFlux openswmm::fv::kernels::riemannFlux ( const FaceState & L, const FaceState & R )

noexcept

Flux for the conservative St. Venant system, plus the contact speed.

The hydrodynamic flux is HLL, and that is not a shortcut — it is what HLLC reduces to here. The system U = [A, Q] is 2x2 with two genuinely nonlinear fields and NO middle wave: the contact appears only once a third component is carried, U = [A, Q, A(phi)], whose eigenvalues are u-c, u, u+c. Applying the Euler-style HLLC star-state construction to the 2x2 subsystem is inconsistent — the resulting momentum flux violates the HLL consistency condition, and on a pressurized/part-full interface (a surcharged manhole feeding a half-full pipe — the case this solver exists for) it disagrees with HLL by more than an order of magnitude and drives the flow backwards.

So S* is computed here for its actual job: selecting the upwind side for the ADVECTED SCALAR (see speciesFlux). That is exactly the role plan §3.2 assigns it — "λ = v IS the contact discontinuity that carries the species" — and it is why FV_RIEMANN changes the transport answer while leaving the hydraulics identical.

Here is the call graph for this function:

sectionArea()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::sectionArea ( const FvGeometry & g, double h )

noexcept

Section area at depth h, EXCLUDING the slot (h clamped to the crown). Scaled by the barrel count: a cell is the aggregate of the parallel barrels at a shared depth (see FvGeometry::barrel_scale).

Here is the caller graph for this function:

sectionHydRad()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::sectionHydRad ( const FvGeometry & g, double h )

noexcept

Section hydraulic radius at depth h. Frozen at r_full above the crown — the slot is a numerical device and must not contribute wetted perimeter (same convention as legacy dwflow.c::getHydRad).

Here is the caller graph for this function:

sectionWidth()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::sectionWidth ( const FvGeometry & g, double h )

noexcept

Section top width at depth h, EXCLUDING the slot. RECT_CLOSED returns 0 at exactly y_full (the crown is a point); the slot term added by widthOfDepth is what keeps the total width — and hence the celerity — finite there.

Here is the caller graph for this function:

slotRamp()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::slotRamp ( double s )

noexcept

Smoothstep ramp opening the slot over [y_crown, y_full]. C¹ at both ends, so dA/dh is continuous through the crown and the Riemann solver's wave-speed estimates do not see a jump. Above the crown the ramp saturates at 1.

Here is the caller graph for this function:

slotRampIntegral()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::slotRampIntegral ( double s )

noexcept

∫₀ˢ slotRamp — the normalized slot area accumulated up to s. Continuous and exactly matched to slotRamp so A(h) = ∫T(h)dh holds through the taper.

Here is the caller graph for this function:

speciesFlux()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::speciesFlux ( const FaceState & L, const FaceState & R, const FaceFlux & f, double phi_l, double phi_r, bool hllc )

noexcept

Species flux for the A(phi) component.

HLLC upwinds on the sign of the contact speed and multiplies by the SAME mass flux the water used. Flux consistency is a hard requirement, not a detail: computing this from a separately-evaluated velocity — the usual trap when transport is bolted onto a hydraulic solver — decouples solute mass from water mass and produces spurious extrema and non-conservation. Reusing the mass flux guarantees exact solute conservation, a discrete maximum principle, and that a uniform field stays uniform under any flow (plan §3.2).

HLL instead averages the whole fan into one intermediate state, smearing precisely the wave that carries the species. It still preserves a uniform field (the average is linear in phi when phi is constant), so it is a legitimate baseline — just a diffusive one. The §6.11(d) front-width comparison separates this layer's contribution from the reconstruction's.

tpaAreaOfDepth()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::tpaAreaOfDepth ( const FvGeometry & g, double h )

noexcept

tpaDepthOfArea()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::tpaDepthOfArea ( const FvGeometry & g, double a )

noexcept

tpaI1OfDepth()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::tpaI1OfDepth ( const FvGeometry & g, double h )

noexcept

ufUpdate()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::ufUpdate ( double q, double a, double u_old, double k3, double grad_term, double dt )

noexcept

Unsteady-friction momentum update (issue #156).

Pinto/Vasconcelos/Soares (2025) source term, split by stiffness: S_fu = (k3/g)·(∂V/∂t + c·sgn(V)·|∂V/∂x|) The local-acceleration half integrates IMPLICITLY (Δt cancels: ΔQ = −k3·A·(V^{n+1} − Vⁿ) ⇒ Q^{n+1} = (Q* + k3·A·Vⁿ)/(1 + k3)), so it is unconditionally stable like frictionUpdate. The convective half enters explicitly through grad_term = c·sgn(Vⁿ)·|∂V/∂x|ⁿ, precomputed by the solver from a consistent old-state snapshot (never from mid-update neighbors). The combined change is clamped to half the incoming momentum per substep — the paper reports instability for k3 ≳ 0.02 with no such guard; the clamp also makes UF exactly inert at rest (q = 0 ⇒ cap = 0, and u_old = grad = 0 anyway), preserving the well-balanced property. k3 = 0 returns q bit-unchanged.

waveSpeeds()

OPENSWMM_KERNEL_FN void openswmm::fv::kernels::waveSpeeds ( const FaceState & L, const FaceState & R, double & sl, double & sr )

noexcept

Davis/Einfeldt wave-speed estimates with dry-state handling. A dry side has no celerity of its own, so the wet side's rarefaction-tail speed (u ∓ 2c) is used — the standard dry-bed estimate that keeps the front speed correct.

Here is the caller graph for this function:

widthOfDepth()

OPENSWMM_KERNEL_FN double openswmm::fv::kernels::widthOfDepth ( const FvGeometry & g, double h )

noexcept

Top width dA/dh at depth h, INCLUDING the tapered slot.

Here is the call graph for this function:

Variable Documentation

kDryArea

double openswmm::fv::kernels::kDryArea = 1.0e-12

inlineconstexpr

Area floor paired with kDryDepth for the u = Q/A division.

kDryDepth

double openswmm::fv::kernels::kDryDepth = 1.0e-7

inlineconstexpr

Depth below which a cell is treated as dry: no velocity, no flux, no celerity contribution. Small enough to be hydraulically irrelevant (3e-5 ft ≈ 0.01 mm) and large enough to keep Q/A finite.

kEtaDeadband

double openswmm::fv::kernels::kEtaDeadband = 1.0e-12

inlineconstexpr

Free-surface difference below which a face is treated as exactly level. Without it the momentum update integrates 1-ulp closure round-trip noise; 1e-12 ft is far below any physical head.

kGravity

double openswmm::fv::kernels::kGravity = 32.2

inlineconstexpr

Gravity in internal units. Matches constants::GRAVITY — spelled here so the header stays dependency-free for the plugin build.