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;
 
// ===========================================================================
// Section evaluation — the ONE place the cross-section machinery is called
// ===========================================================================
//
// These three go through the geometry's own evaluator (XSectKernels.hpp) — the
// SAME bodies XSection.cpp's public accessors and XSectBatch's kernels run.
// Only where the evaluator's tables live differs between the host solver and
// the device backend, so the §6.8 parity harness compares two instantiations of
// one implementation rather than two implementations.

OPENSWMM_KERNEL_FN double sectionArea(const FvGeometry& g, double h) noexcept {
    if (h <= 0.0) return 0.0;
    return g.barrel_scale * g.eval->getAofY(g.xs, (h < g.y_full) ? h : g.y_full);
}

OPENSWMM_KERNEL_FN double sectionWidth(const FvGeometry& g, double h) noexcept {
    if (h <= 0.0) return 0.0;
    return g.barrel_scale * g.eval->getWofY(g.xs, (h < g.y_full) ? h : g.y_full);
}

OPENSWMM_KERNEL_FN double sectionHydRad(const FvGeometry& g, double h) noexcept {
    if (h <= 0.0) return 0.0;
    if (h >= g.y_full) return g.r_full;
    return g.eval->getRofY(g.xs, h);
}
 
// ===========================================================================
// 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 {
    if (h <= 0.0) return 0.0;
    if (h >= g.y_full) return g.a_crown + g.t_slot * (h - g.y_full);
    const double band = g.y_full - g.y_crown;
    const double ax = sectionArea(g, h);
    if (band <= 0.0) return ax;                 // open section: no taper band
    const double s = (h - g.y_crown) / band;
    return ax + g.t_slot * band * slotRampIntegral(s);
}

OPENSWMM_KERNEL_FN double widthOfDepth(const FvGeometry& g, double h) noexcept {
    if (h <= 0.0) return 0.0;
    if (h >= g.y_full) return g.t_slot;
    const double band = g.y_full - g.y_crown;
    const double wx = sectionWidth(g, h);
    if (band <= 0.0) return wx;
    const double s = (h - g.y_crown) / band;
    return wx + g.t_slot * slotRamp(s);
}

OPENSWMM_KERNEL_FN double hydRadOfDepth(const FvGeometry& g, double h) noexcept {
    if (h <= 0.0) return 0.0;
    if (h >= g.y_full) return g.r_full;
    return sectionHydRad(g, h);
}

OPENSWMM_KERNEL_FN double i1OfDepth(const FvGeometry& g, double h,
                                    double area_at_h) noexcept {
    if (h <= 0.0) return 0.0;
    if (h >= g.y_full) {
        const double d = h - g.y_full;
        return g.i1_crown + g.a_crown * d + 0.5 * g.t_slot * d * d;
    }
    const int n = static_cast<int>(kI1Samples);
    const double dh = g.y_full / static_cast<double>(n - 1);
    int i = static_cast<int>(h / dh);
    if (i < 0) i = 0;
    if (i > n - 2) i = n - 2;
    const double h_i = static_cast<double>(i) * dh;
    // i1_tbl stores I₁ at the sample points; the companion area sample is the
    // second half of the same buffer (see buildI1Table).
    const double i1_i = g.i1_tbl[static_cast<std::size_t>(i)];
    const double a_i  = g.i1_tbl[static_cast<std::size_t>(n + i)];
    return i1_i + 0.5 * (a_i + area_at_h) * (h - h_i);
}

OPENSWMM_KERNEL_FN double depthOfAreaBracketed(const FvGeometry& g,
                                               double a) noexcept {
    if (a <= 0.0) return 0.0;
    if (a >= g.a_crown) return g.y_full + (a - g.a_crown) / g.t_slot;
 
    const int n = static_cast<int>(kI1Samples);
    const double dh = g.y_full / static_cast<double>(n - 1);
 
    int lo = 0, hi = n - 1;
    while (hi - lo > 1) {
        const int mid = (lo + hi) / 2;
        if (g.i1_tbl[static_cast<std::size_t>(n + mid)] <= a) lo = mid;
        else                                                  hi = mid;
    }
 
    double xa = static_cast<double>(lo) * dh;
    double xb = static_cast<double>(hi) * dh;
    double fa = g.i1_tbl[static_cast<std::size_t>(n + lo)] - a;
    double fb = g.i1_tbl[static_cast<std::size_t>(n + hi)] - a;
    if (fa >= 0.0) return xa;
    if (fb <= 0.0) return xb;
 
    for (int it = 0; it < 40; ++it) {
        double x = xb - fb * (xb - xa) / (fb - fa);
        if (!(x > xa && x < xb)) x = 0.5 * (xa + xb);
        const double f = areaOfDepth(g, x) - a;
        if (f == 0.0) return x;
        if (f < 0.0) { xa = x; fa = f; fb *= 0.5; }   // Illinois down-weighting
        else         { xb = x; fb = f; fa *= 0.5; }
        if (xb - xa <= 1.0e-15 * g.y_full) break;
    }
    return 0.5 * (xa + xb);
}

OPENSWMM_KERNEL_FN double depthOfArea(const FvGeometry& g, double a) noexcept {
    if (a <= 0.0) return 0.0;
    if (a >= g.a_crown) return g.y_full + (a - g.a_crown) / g.t_slot;
 
    const int n = static_cast<int>(kI1Samples);
    const double da = g.a_crown / static_cast<double>(n - 1);
    if (!(da > 0.0)) return depthOfAreaBracketed(g, a);
 
    int j = static_cast<int>(a / da);
    if (j < 1) j = 1;
    if (j > n - 3) j = n - 3;
 
    double xa = g.h_tbl[static_cast<std::size_t>(j - 1)];
    double xb = g.h_tbl[static_cast<std::size_t>(j + 2)];
    if (!(xb > xa)) return depthOfAreaBracketed(g, a);
 
    double fa = areaOfDepth(g, xa) - a;
    double fb = areaOfDepth(g, xb) - a;
    if (fa > 0.0 || fb < 0.0) return depthOfAreaBracketed(g, a);
    if (fa == 0.0) return xa;
    if (fb == 0.0) return xb;
 
    const double tol = 1.0e-15 * g.y_full;
    double c = xa, fc = fa, d = xb - xa, e = d;
    for (int it = 0; it < 60; ++it) {
        if (fb * fc > 0.0) { c = xa; fc = fa; d = xb - xa; e = d; }
        if (std::fabs(fc) < std::fabs(fb)) {
            xa = xb; xb = c; c = xa;
            fa = fb; fb = fc; fc = fa;
        }
        const double m = 0.5 * (c - xb);
        if (std::fabs(m) <= tol || fb == 0.0) return xb;
 
        if (std::fabs(e) < tol || std::fabs(fa) <= std::fabs(fb)) {
            d = m; e = m;                                   // bisect
        } else {
            const double sfb = fb / fa;
            double p, q;
            if (xa == c) {                                  // secant
                p = 2.0 * m * sfb;
                q = 1.0 - sfb;
            } else {                                        // inverse quadratic
                const double qq = fa / fc;
                const double r  = fb / fc;
                p = sfb * (2.0 * m * qq * (qq - r) - (xb - xa) * (r - 1.0));
                q = (qq - 1.0) * (r - 1.0) * (sfb - 1.0);
            }
            if (p > 0.0) q = -q; else p = -p;
            if (2.0 * p < ((3.0 * m * q - std::fabs(tol * q) < std::fabs(e * q))
                               ? 3.0 * m * q - std::fabs(tol * q)
                               : std::fabs(e * q))) {
                e = d; d = p / q;
            } else {
                d = m; e = m;
            }
        }
        xa = xb; fa = fb;
        xb += (std::fabs(d) > tol) ? d : ((m > 0.0) ? tol : -tol);
        fb = areaOfDepth(g, xb) - a;
    }
    return xb;
}

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;  
};

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);
    } 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 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