FvKernels.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: Apache-2.0
//
// Copyright 2026 Caleb Buahin
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

 
#ifndef OPENSWMM_ENGINE_FV_KERNELS_HPP
#define OPENSWMM_ENGINE_FV_KERNELS_HPP
 
#include <algorithm>
#include <cmath>
 
#include "FvOptions.hpp"
#include "NetworkMeshData.hpp"
 
// Portable kernel-function marker — identical convention to
// 2d/solver/InertialKernels.hpp:45. Host builds get plain `inline`; the GPU
// plugin defines OPENSWMM_KERNEL_FN to KOKKOS_INLINE_FUNCTION *before*
// including this header so the same scalar bodies compile for the device.
#ifndef OPENSWMM_KERNEL_FN
#define OPENSWMM_KERNEL_FN inline
#endif
 
namespace openswmm::fv::kernels {

inline constexpr double kGravity = 32.2;

inline constexpr double kDryDepth = 1.0e-7;

inline constexpr double kDryArea = 1.0e-12;

inline constexpr double kEtaDeadband = 1.0e-12;
 
// ===========================================================================
// Preissmann slot taper (plan §3.3.2)
// ===========================================================================

OPENSWMM_KERNEL_FN double slotRamp(double s) noexcept {
    if (s <= 0.0)
        return 0.0;
    if (s >= 1.0)
        return 1.0;
    return s * s * (3.0 - 2.0 * s);
}

OPENSWMM_KERNEL_FN double slotRampIntegral(double s) noexcept {
    if (s <= 0.0)
        return 0.0;
    if (s >= 1.0)
        return s - 0.5; // 1 − ½ from the taper, then linear
    const double s3 = s * s * s;
    return s3 - 0.5 * s3 * s; // s³ − s⁴/2
}
 
// ===========================================================================
// Closure — one continuous geometry from dry bed to full pressurization
// ===========================================================================
//
// This is the memoryless SLOT closure — every cell, every state, no history.
// TPA exception (issue #156): under FV_PRESSURE_CLOSURE TPA a cell whose
// regime flag is set evaluates the tpa* kernels at the end of this file
// instead. The flag is physical air-pathway history (cleared on cold start
// and hotstart restore), not numerical relaxation; unflagged cells — and the
// entire SLOT closure path — are untouched.

OPENSWMM_KERNEL_FN double areaOfDepth(const FvGeometry& g, double h) noexcept {
    return closureArea(g.closure_tbl, h);
}

OPENSWMM_KERNEL_FN double widthOfDepth(const FvGeometry& g, double h) noexcept {
    return closureWidth(g.closure_tbl, h);
}

OPENSWMM_KERNEL_FN double hydRadOfDepth(const FvGeometry& g, double h) noexcept {
    return closureHydRad(g.closure_tbl, h);
}

OPENSWMM_KERNEL_FN double i1OfDepth(const FvGeometry& g, double h, double area_at_h) noexcept {
    (void)area_at_h;
    return closureI1(g.closure_tbl, h);
}

OPENSWMM_KERNEL_FN double depthOfArea(const FvGeometry& g, double a) noexcept {
    return closureDepthOfArea(g.closure_tbl, a);
}

OPENSWMM_KERNEL_FN double celerity(double a, double t) noexcept {
    if (a <= kDryArea || t <= 0.0)
        return 0.0;
    return std::sqrt(kGravity * a / t);
}
 
// ===========================================================================
// Riemann solver
// ===========================================================================

struct FaceState {
    double a = 0.0;  
    double q = 0.0;  
    double u = 0.0;  
    double c = 0.0;  
    double i1 = 0.0; 
    uint8_t press = 0;
};

struct FaceFlux {
    double mass = 0.0;  
    double mom = 0.0;   
    double sstar = 0.0; 
    double sl = 0.0;    
    double sr = 0.0;    
};

OPENSWMM_KERNEL_FN void waveSpeeds(const FaceState& L, const FaceState& R, double& sl,
                                   double& sr) noexcept {
    const bool wetL = (L.a > kDryArea);
    const bool wetR = (R.a > kDryArea);
    if (wetL && wetR) {
        sl = std::min(L.u - L.c, R.u - R.c);
        sr = std::max(L.u + L.c, R.u + R.c);
        // Slot/free-surface interface (one side in the slot, the other not).
        // Davis's symmetric estimate carries the slot's acoustic celerity into
        // BOTH waves, but no signal enters the free-surface side at the slot
        // speed: the wave that crosses into it is the filling/emptying bore,
        // whose speed is the Rankine–Hugoniot jump (ΔQ/ΔA), bounded below by
        // that side's own u ± c and above by the Davis bound. With the slot
        // celerity on both waves the HLL diffusion term at a pressurized
        // entrance is O(c_slot·ΔA) — measured 52 cfs of spurious mass flux and
        // a −245 ft⁴/s² momentum sink on a 3 ft culvert whose junction then
        // flooded at 108 of 120 cfs however high its head rose.
        // A vented junction's ghost (press == 2) is excluded: bounding the
        // wave into a pressurized cell beside it cut the acoustic exchange at
        // every lateral junction of a pressurized tunnel (the node solve then
        // clamped at its rim and booked 63 000 m³ of flooding — Klaver
        // Example_02, 2026-09-13), and bounding the wave into a free cell
        // beside it drove a 2× over-conveying limit cycle at an overloaded
        // entrance. The bound acts at the INTERIOR slot/free faces, which is
        // where the entrance lock lived.
        if (L.press != R.press && L.press < 2 && R.press < 2) {
            if (L.press) {
                const double da = L.a - R.a;
                const double w = (da > 0.0) ? (L.q - R.q) / da : sr;
                sr = std::max(R.u + R.c, std::min(sr, w));
            } else {
                const double da = R.a - L.a;
                const double w = (da > 0.0) ? (R.q - L.q) / da : sl;
                sl = std::min(L.u - L.c, std::max(sl, w));
            }
        }
    } else if (wetR) { // dry left
        sl = R.u - 2.0 * R.c;
        sr = R.u + R.c;
    } else if (wetL) { // dry right
        sl = L.u - L.c;
        sr = L.u + 2.0 * L.c;
    } else {
        sl = 0.0;
        sr = 0.0;
    }
}

OPENSWMM_KERNEL_FN void physicalFlux(const FaceState& S, double& fa, double& fq) noexcept {
    fa = S.q;
    fq = S.q * S.u + kGravity * S.i1;
}

OPENSWMM_KERNEL_FN double riemannMassFlux(const FaceState& L, const FaceState& R) noexcept {
    if (L.a <= kDryArea && R.a <= kDryArea)
        return 0.0;
    double sl = 0.0, sr = 0.0;
    waveSpeeds(L, R, sl, sr);
    double fal = 0.0, fql = 0.0, far = 0.0, fqr = 0.0;
    physicalFlux(L, fal, fql);
    physicalFlux(R, far, fqr);
    if (sl >= 0.0)
        return fal;
    if (sr <= 0.0)
        return far;
    return (sr * fal - sl * far + sl * sr * (R.a - L.a)) / (sr - sl);
}

OPENSWMM_KERNEL_FN FaceFlux riemannFlux(const FaceState& L, const FaceState& R) noexcept {
    FaceFlux out;
    if (L.a <= kDryArea && R.a <= kDryArea)
        return out;
 
    double sl = 0.0, sr = 0.0;
    waveSpeeds(L, R, sl, sr);
    out.sl = sl;
    out.sr = sr;
 
    double fal = 0.0, fql = 0.0, far = 0.0, fqr = 0.0;
    physicalFlux(L, fal, fql);
    physicalFlux(R, far, fqr);
 
    if (sl >= 0.0) {
        out.mass = fal;
        out.mom = fql;
        out.sstar = sl;
        return out;
    }
    if (sr <= 0.0) {
        out.mass = far;
        out.mom = fqr;
        out.sstar = sr;
        return out;
    }
 
    const double dsr = sr - sl;
    out.mass = (sr * fal - sl * far + sl * sr * (R.a - L.a)) / dsr;
    out.mom = (sr * fql - sl * fqr + sl * sr * (R.q - L.q)) / dsr;
 
    // Contact speed (plan §3.2). Falls back to the HLL-averaged velocity when
    // the denominator degenerates, which happens only when both sides are
    // vanishing — where the species flux is zero anyway.
    const double den = R.a * (R.u - sr) - L.a * (L.u - sl);
    if (std::fabs(den) > 1.0e-14) {
        out.sstar = (sl * R.a * (R.u - sr) - sr * L.a * (L.u - sl)) / den;
    } else {
        const double a_hll = (sr * R.a - sl * L.a - (far - fal)) / dsr;
        out.sstar = (a_hll > kDryArea) ? out.mass / a_hll : 0.0;
    }
    return out;
}

OPENSWMM_KERNEL_FN double speciesFlux(const FaceState& L, const FaceState& R, const FaceFlux& f,
                                      double phi_l, double phi_r, bool hllc) noexcept {
    if (hllc)
        return f.mass * ((f.sstar >= 0.0) ? phi_l : phi_r);
 
    // HLL on the A(phi) component, with the same supersonic branches as the
    // hydrodynamic flux so the two can never disagree about which side is
    // upstream of the whole fan.
    if (f.sl >= 0.0)
        return L.q * phi_l;
    if (f.sr <= 0.0)
        return R.q * phi_r;
    const double dsr = f.sr - f.sl;
    return (f.sr * L.q * phi_l - f.sl * R.q * phi_r + f.sl * f.sr * (R.a * phi_r - L.a * phi_l)) /
           dsr;
}
 
// ===========================================================================
// Source terms
// ===========================================================================

OPENSWMM_KERNEL_FN double frictionUpdate(double q, double u, double r, double dt,
                                         double rough_factor) noexcept {
    if (r <= 0.0 || rough_factor <= 0.0)
        return q;
    const double r43 = r * std::cbrt(r);
    const double den = 1.0 + dt * rough_factor * std::fabs(u) / r43;
    return q / den;
}

OPENSWMM_KERNEL_FN double localLossUpdate(double q, double u, double k, double dx,
                                          double dt) noexcept {
    if (k <= 0.0 || dx <= 0.0)
        return q;
    return q / (1.0 + dt * k * std::fabs(u) / (2.0 * dx));
}

OPENSWMM_KERNEL_FN double ufUpdate(double q, double a, double u_old, double k3, double grad_term,
                                   double dt) noexcept {
    if (k3 <= 0.0 || a <= 0.0)
        return q;
    // Dead-band (measured, issue #156): the implicit fold acts as added
    // inertia — it resists velocity DECAY too, so applied to the ~mm/s
    // numerical ripple of a storage-coupled pool it sustained noise that
    // steady friction was correctly killing (0.004 → 0.011 cfs on the
    // at-rest fixture). UF correlations are calibrated for real transients
    // (paper velocities O(0.1–1 m/s)); below 0.01 ft/s the term is noise
    // amplification, not physics. Both the old and candidate velocities must
    // clear the floor, so a genuinely at-rest deck stays bit-identical.
    constexpr double kUfVelFloor = 0.01; // ft/s, internal units
    if (std::fabs(u_old) < kUfVelFloor && std::fabs(q / a) < kUfVelFloor)
        return q;
    double qn = (q + k3 * a * u_old) / (1.0 + k3);
    qn -= dt * k3 * a * grad_term;
    const double dq = qn - q;
    const double cap = 0.5 * std::fabs(q);
    if (std::fabs(dq) > cap)
        qn = q + ((dq > 0.0) ? cap : -cap);
    return qn;
}
 
// ---------------------------------------------------------------------------
// TPA — two-component pressure approach (issue #156 Phase 4)
// ---------------------------------------------------------------------------
// Vasconcelos, Wright & Roe (2006). For a cell whose regime FLAG is set
// (state.cell_tpa — physical air-pathway history, updated once per substep by
// the solver, NOT evaluated here), the closure is the slot line extended to
// BOTH signs of ΔA = A − a_crown:
//     h(A) = y_full + (A − a_crown)/t_slot        (hs = h − y_full, signed)
//     A(h) = a_crown + t_slot·(h − y_full)
//     I₁(h) = i1_crown + a_crown·(h − y_full)     (paper Eq. 12b: the slot's
//                                                  ½·t_slot·d² numerical
//                                                  storage pressure is dropped)
//     T = t_slot,  R = r_full,  c = √(g·A/t_slot) ≈ acoustic celerity a.
// Free-surface (unflagged) cells use the table closure above UNCHANGED,
// including the crown taper. The pair is continuous at A = a_crown, h = y_full.

OPENSWMM_KERNEL_FN double tpaDepthOfArea(const FvGeometry& g, double a) noexcept {
    return g.y_full + (a - g.a_crown) / g.t_slot;
}

OPENSWMM_KERNEL_FN double tpaAreaOfDepth(const FvGeometry& g, double h) noexcept {
    return g.a_crown + g.t_slot * (h - g.y_full);
}

OPENSWMM_KERNEL_FN double tpaI1OfDepth(const FvGeometry& g, double h) noexcept {
    return g.i1_crown + g.a_crown * (h - g.y_full);
}

OPENSWMM_KERNEL_FN double faceCflDt(double cfl, double dx, double u, double c) noexcept {
    const double s = std::fabs(u) + c;
    return (s > 1.0e-12) ? cfl * dx / s : 1.0e30;
}

OPENSWMM_KERNEL_FN double positivityScale(double vol, double outflow, double dt) noexcept {
    if (outflow <= 0.0 || dt <= 0.0)
        return 1.0;
    const double take = outflow * dt;
    const double avail = (vol > 0.0) ? vol : 0.0;
    return (take <= avail) ? 1.0 : avail / take;
}
 
} // namespace openswmm::fv::kernels
 
#endif // OPENSWMM_ENGINE_FV_KERNELS_HPP