InertialKernels.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_2D_INERTIAL_KERNELS_HPP
#define OPENSWMM_ENGINE_2D_INERTIAL_KERNELS_HPP
 
#include <algorithm>
#include <cmath>
 
#include "../data/MeshData.hpp"
#include "../data/SolverOptions2D.hpp"
#include "../mesh/VfrClosure.hpp"
#include "../mesh/QuadVfr.hpp"
 
// Portable kernel-function marker (P5 Kokkos port): host builds get plain
// `inline`; the GPU plugin defines OPENSWMM_KERNEL_FN to
// KOKKOS_INLINE_FUNCTION *before* including this header so the identical
// scalar bodies compile for the device (CUDA/HIP/SYCL) as well. The math uses
// std:: functions deliberately — CUDA/HIP provide device overloads and Kokkos
// enables relaxed-constexpr; keeping one spelling preserves single-source
// bit-identity with the serial marcher.
#ifndef OPENSWMM_KERNEL_FN
#define OPENSWMM_KERNEL_FN inline
#endif
 
namespace openswmm::twoD::inertial {
 
inline constexpr double kGravity = 9.80665; 

inline constexpr double kEtaDeadband = 1.0e-12;

OPENSWMM_KERNEL_FN double momentumBlendFraction(double theta, double dt,
                                                 double reference_step) noexcept {
    if (!(reference_step > 0.0))
        return 1.0 - theta;
    if (!(dt > 0.0) || theta == 1.0)
        return 0.0;
    return -std::expm1(std::log(theta) * (dt / reference_step));
}

OPENSWMM_KERNEL_FN double qMagnitude(double x, double y) noexcept {
    return std::sqrt(x * x + y * y);
}

OPENSWMM_KERNEL_FN void etaDepthScalar(double area, double cz, double z1, double z2, double z3,
                                       bool vfr, double vfr_min_wet_frac, double V, double& eta,
                                       double& depth) noexcept {
    const double v = (V > 0.0) ? V : 0.0;
    depth = (area > 1.0e-30) ? v / area : 0.0;
    if (vfr) {
        vfrSort3(z1, z2, z3);
        eta = vfrEtaFromMeanDepth(z1, z2, z3, depth, vfr_min_wet_frac);
    } else {
        eta = cz + depth;
    }
}

OPENSWMM_KERNEL_FN double volumeFromEtaScalar(double area, double cz, double z1, double z2,
                                              double z3, bool vfr, double vfr_min_wet_frac,
                                              double eta) noexcept {
    if (vfr) {
        vfrSort3(z1, z2, z3);
        return area * vfrMeanDepthFromEta(z1, z2, z3, eta, vfr_min_wet_frac);
    }
    const double d = eta - cz;
    return (d > 0.0) ? area * d : 0.0;
}

OPENSWMM_KERNEL_FN void etaDepthQuadScalar(double area, double cz, const double* zs, double A1,
                                           double A2, bool vfr, double vfr_min_wet_frac, double V,
                                           double& eta, double& depth) noexcept {
    const double v = (V > 0.0) ? V : 0.0;
    depth = (area > 1.0e-30) ? v / area : 0.0;
    if (vfr)
        eta = quadEtaFromMeanDepth(zs, A1, A2, depth, vfr_min_wet_frac);
    else
        eta = cz + depth;
}

OPENSWMM_KERNEL_FN double volumeFromEtaQuadScalar(double area, double cz, const double* zs,
                                                  double A1, double A2, bool vfr,
                                                  double vfr_min_wet_frac, double eta) noexcept {
    if (vfr)
        return area * quadMeanDepthFromEta(zs, A1, A2, eta, vfr_min_wet_frac);
    const double d = eta - cz;
    return (d > 0.0) ? area * d : 0.0;
}
 
#if defined(__GNUC__) || defined(__clang__)
#define OPENSWMM_2D_NOINLINE __attribute__((noinline))
#elif defined(_MSC_VER)
#define OPENSWMM_2D_NOINLINE __declspec(noinline)
#else
#define OPENSWMM_2D_NOINLINE
#endif

OPENSWMM_2D_NOINLINE inline void cellEtaDepthQuad(const MeshData& m, const SolverOptions2D& o,
                                                  int i, double V, double& eta,
                                                  double& depth) noexcept {
    etaDepthQuadScalar(
        m.tri_area[i], m.tri_cz[i], &m.quad_vfr_z[static_cast<std::size_t>(i) * kQuadVfrZ],
        m.quad_vfr_a[static_cast<std::size_t>(i) * 2],
        m.quad_vfr_a[static_cast<std::size_t>(i) * 2 + 1], true, o.vfr_min_wet_frac, V, eta, depth);
}

inline void cellEtaDepth(const MeshData& m, const SolverOptions2D& o, int i, double V, double& eta,
                         double& depth) noexcept {
    // Branch BEFORE the geometry loads: under the default FLAT closure the
    // three vertex-index loads and the three scattered vz gathers below are
    // dead, and this is the hottest call in the marcher (every cell, every
    // firing).
    if (o.cell_closure != CellClosure2D::VFR) {
        etaDepthScalar(m.tri_area[i], m.tri_cz[i], 0.0, 0.0, 0.0, false, o.vfr_min_wet_frac, V, eta,
                       depth);
        return;
    }
    if (m.cell_nv[static_cast<std::size_t>(i)] == 4) {
        cellEtaDepthQuad(m, o, i, V, eta, depth);
        return;
    }
    etaDepthScalar(m.tri_area[i], m.tri_cz[i], m.vz[m.cell_vertex(i, 0)], m.vz[m.cell_vertex(i, 1)],
                   m.vz[m.cell_vertex(i, 2)], true, o.vfr_min_wet_frac, V, eta, depth);
}

inline double cellWetArea(const MeshData& m, const SolverOptions2D& o, int i,
                          double volume) noexcept {
    if (!(volume > 0))
        return 0;
    if (o.cell_closure != CellClosure2D::VFR)
        return m.tri_area[i];
    double eta = 0, depth = 0;
    cellEtaDepth(m, o, i, volume, eta, depth);
    double fraction;
    if (m.cell_nv[std::size_t(i)] == 4)
        fraction = quadWetFraction(&m.quad_vfr_z[std::size_t(i) * kQuadVfrZ],
                                   m.quad_vfr_a[std::size_t(i) * 2],
                                   m.quad_vfr_a[std::size_t(i) * 2 + 1], eta);
    else {
        double a = m.vz[m.cell_vertex(i, 0)], b = m.vz[m.cell_vertex(i, 1)],
               c = m.vz[m.cell_vertex(i, 2)];
        vfrSort3(a, b, c);
        fraction = vfrWetFraction(a, b, c, eta);
    }
    return m.tri_area[i] * std::clamp(fraction, 0., 1.);
}

inline double cellVolumeFromEta(const MeshData& m, const SolverOptions2D& o, int i,
                                double eta) noexcept {
    if (m.cell_nv[static_cast<std::size_t>(i)] == 4) {
        return volumeFromEtaQuadScalar(
            m.tri_area[i], m.tri_cz[i], &m.quad_vfr_z[static_cast<std::size_t>(i) * kQuadVfrZ],
            m.quad_vfr_a[static_cast<std::size_t>(i) * 2],
            m.quad_vfr_a[static_cast<std::size_t>(i) * 2 + 1], o.cell_closure == CellClosure2D::VFR,
            o.vfr_min_wet_frac, eta);
    }
    return volumeFromEtaScalar(m.tri_area[i], m.tri_cz[i], m.vz[m.cell_vertex(i, 0)],
                               m.vz[m.cell_vertex(i, 1)], m.vz[m.cell_vertex(i, 2)],
                               o.cell_closure == CellClosure2D::VFR, o.vfr_min_wet_frac, eta);
}

OPENSWMM_KERNEL_FN double faceFlowDepth(double etaL, double etaR, double zface) noexcept {
    return std::max(etaL, etaR) - zface;
}

OPENSWMM_KERNEL_FN double faceDepthFromEta(double eta, double z_lo, double z_hi) noexcept {
    if (eta <= z_lo)
        return 0.0;
    const double dz = z_hi - z_lo;
    if (dz < 1.0e-9)
        return eta - z_lo; // level edge
    if (eta <= z_hi) {
        const double t = eta - z_lo;
        return t * t / (2.0 * dz);
    }
    return eta - 0.5 * (z_lo + z_hi);
}

OPENSWMM_KERNEL_FN double faceFlowDepthVfr(double etaL, double etaR, double z_lo,
                                           double z_hi) noexcept {
    return faceDepthFromEta(std::max(etaL, etaR), z_lo, z_hi);
}

OPENSWMM_KERNEL_FN double inertialFaceUpdate(double q, double qhat, double hf, double dt,
                                             double slope, double n2, double q_mag,
                                             double adv = 0.0) noexcept {
    const double h73 = hf * hf * std::cbrt(hf);
    const double num = qhat - dt * (kGravity * hf * slope + adv);
    const double den = 1.0 + kGravity * dt * n2 * q_mag / h73;
    return num / den;
}

OPENSWMM_KERNEL_FN double inertialAdvection(double q_f, double unL, double hL, double unR,
                                            double hR, double inv_dx_normal) noexcept {
    const double FL = unL * ((unL > 0.0) ? hL * unL : q_f);
    const double FR = unR * ((unR > 0.0) ? q_f : hR * unR);
    return (FR - FL) * inv_dx_normal;
}

OPENSWMM_KERNEL_FN double froudeCap(double q, double hf, double fr_max) noexcept {
    const double qcap = fr_max * hf * std::sqrt(kGravity * hf);
    return std::clamp(q, -qcap, qcap);
}

OPENSWMM_KERNEL_FN double cellCflDt(double alpha, double lchar, double h, double speed) noexcept {
    const double c = std::sqrt(kGravity * h) + speed;
    return (c > 1.0e-12) ? alpha * lchar / c : 1.0e30;
}

OPENSWMM_KERNEL_FN double positivityScale(double V, double outflow_m3s, double dt,
                                          double beta) noexcept {
    if (outflow_m3s <= 0.0 || dt <= 0.0)
        return 1.0;
    const double avail = beta * std::max(V, 0.0);
    const double take = outflow_m3s * dt;
    return (take <= avail) ? 1.0 : avail / take;
}
 
} // namespace openswmm::twoD::inertial
 
#endif // OPENSWMM_ENGINE_2D_INERTIAL_KERNELS_HPP