SweKernels.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_SWE_KERNELS_HPP
#define OPENSWMM_ENGINE_2D_SWE_KERNELS_HPP
 
#include <algorithm>
#include <cmath>
 
#ifndef OPENSWMM_KERNEL_FN
#define OPENSWMM_KERNEL_FN inline
#endif
 
namespace openswmm::twoD::swe {
 
inline constexpr double kGravity = 9.80665;

struct FaceSide {
    double h = 0.0;
    double un = 0.0;
    double ut = 0.0;
};

struct FaceFlux {
    double mass = 0.0;
    double mx = 0.0;
    double my = 0.0;
    double sstar = 0.0;
};

OPENSWMM_KERNEL_FN void physicalFlux(const FaceSide& s, double& fh, double& fn,
                                     double& ft) noexcept {
    fh = s.h * s.un;
    fn = s.h * s.un * s.un + 0.5 * kGravity * s.h * s.h;
    ft = s.h * s.un * s.ut;
}

OPENSWMM_KERNEL_FN void waveSpeeds(const FaceSide& L, const FaceSide& R, double& sl,
                                   double& sr) noexcept {
    const double cl = std::sqrt(kGravity * L.h);
    const double cr = std::sqrt(kGravity * R.h);
    if (L.h > 0.0 && R.h > 0.0) {
        const double ustar = 0.5 * (L.un + R.un) + cl - cr;
        const double cstar = 0.5 * (cl + cr) + 0.25 * (L.un - R.un);
        sl = std::min(L.un - cl, ustar - cstar);
        sr = std::max(R.un + cr, ustar + cstar);
    } else if (L.h > 0.0) { // dry right
        sl = L.un - cl;
        sr = L.un + 2.0 * cl;
    } else if (R.h > 0.0) { // dry left
        sl = R.un - 2.0 * cr;
        sr = R.un + cr;
    } else {
        sl = sr = 0.0;
    }
}

OPENSWMM_KERNEL_FN void hllcFlux(const FaceSide& L, const FaceSide& R, double& fh, double& fn,
                                 double& ft, double& sstar) noexcept {
    if (L.h <= 0.0 && R.h <= 0.0) {
        fh = fn = ft = 0.0;
        sstar = 0.0;
        return;
    }
    double sl, sr;
    waveSpeeds(L, R, sl, sr);
    double fhl, fnl, ftl, fhr, fnr, ftr;
    // Contact speed (Toro Eq. 10.70) — falls back to the mean normal
    // velocity when the denominator degenerates (both sides at rest).
    const double dl = L.h * (L.un - sl), dr = R.h * (R.un - sr);
    const double den = dr - dl;
    sstar = (std::fabs(den) > 1.0e-14) ? (sl * dr - sr * dl) / den : 0.5 * (L.un + R.un);
 
    if (sl >= 0.0) {
        physicalFlux(L, fhl, fnl, ftl);
        fh = fhl;
        fn = fnl;
        ft = ftl;
    } else if (sr <= 0.0) {
        physicalFlux(R, fhr, fnr, ftr);
        fh = fhr;
        fn = fnr;
        ft = ftr;
    } else {
        // HLL middle state for mass and normal momentum.
        physicalFlux(L, fhl, fnl, ftl);
        physicalFlux(R, fhr, fnr, ftr);
        const double inv = 1.0 / (sr - sl);
        fh = (sr * fhl - sl * fhr + sl * sr * (R.h - L.h)) * inv;
        fn = (sr * fnl - sl * fnr + sl * sr * (R.h * R.un - L.h * L.un)) * inv;
        // Tangential momentum: passive scalar upwinded on the contact.
        ft = fh * ((sstar >= 0.0) ? L.ut : R.ut);
    }
}

OPENSWMM_KERNEL_FN bool faceFlux(double etaL, double hL, double qxL, double qyL, double etaR,
                                 double hR, double qxR, double qyR, double nx, double ny,
                                 double h_dry, FaceFlux& out, double& corrL_x, double& corrL_y,
                                 double& corrR_x, double& corrR_y) noexcept {
    // Beds consistent with the storage: z = η − h̄ (FLAT closure: the cell
    // mean bed; VFR: the flat-equivalent bed, which is what keeps a uniform η
    // an exact rest state under either closure).
    const double zL = etaL - hL, zR = etaR - hR;
    const double zf = (zL > zR) ? zL : zR;
    FaceSide L, R;
    L.h = (etaL - zf > 0.0) ? etaL - zf : 0.0;
    R.h = (etaR - zf > 0.0) ? etaR - zf : 0.0;
    // A side whose reconstructed depth is negligible carries no velocity
    // (avoids q/h blow-ups at fronts); the hydrostatic reconstruction alone
    // already removes the "dry cell standing higher" flux.
    const double tx = -ny, ty = nx;
    if (hL > h_dry && L.h > 0.0) {
        const double ux = qxL / hL, uy = qyL / hL;
        L.un = ux * nx + uy * ny;
        L.ut = ux * tx + uy * ty;
    }
    if (hR > h_dry && R.h > 0.0) {
        const double ux = qxR / hR, uy = qyR / hR;
        R.un = ux * nx + uy * ny;
        R.ut = ux * tx + uy * ty;
    }
    // Hydrostatic pressure remains on the wet side of a blocked face.
    // Dropping it when BOTH reconstructed depths vanish accelerates a lake
    // away from its shoreline even though no water can cross this face.
    const double cL = 0.5 * kGravity * (L.h * L.h - hL * hL);
    const double cR = 0.5 * kGravity * (R.h * R.h - hR * hR);
    corrL_x = cL * nx;
    corrL_y = cL * ny;
    corrR_x = -cR * nx;
    corrR_y = -cR * ny;
    if (L.h <= 0.0 && R.h <= 0.0) {
        out = FaceFlux{};
        return cL != 0.0 || cR != 0.0;
    }
    double fh, fn, ft, sstar;
    hllcFlux(L, R, fh, fn, ft, sstar);
    out.mass = fh;
    out.mx = fn * nx + ft * tx;
    out.my = fn * ny + ft * ty;
    out.sstar = sstar;
    return true;
}

OPENSWMM_KERNEL_FN bool faceFluxReconBed(double etaLf, double uxLf, double uyLf, double zL,
                                         double hL_cell, double etaRf, double uxRf, double uyRf,
                                         double zR, double hR_cell, double nx, double ny,
                                         double h_dry, FaceFlux& out, double& corrL_x,
                                         double& corrL_y, double& corrR_x, double& corrR_y,
                                         double zLf, double zRf) noexcept {
    const double zf = (zLf > zRf) ? zLf : zRf;
    FaceSide L, R;
    L.h = (etaLf - zf > 0.0) ? etaLf - zf : 0.0;
    R.h = (etaRf - zf > 0.0) ? etaRf - zf : 0.0;
    const double tx = -ny, ty = nx;
    if (hL_cell > h_dry && L.h > 0.0) {
        L.un = uxLf * nx + uyLf * ny;
        L.ut = uxLf * tx + uyLf * ty;
    }
    if (hR_cell > h_dry && R.h > 0.0) {
        R.un = uxRf * nx + uyRf * ny;
        R.ut = uxRf * tx + uyRf * ty;
    }
    const double hLf = std::max(0.0, etaLf - zLf);
    const double hRf = std::max(0.0, etaRf - zRf);
    // Hydrostatic face correction plus a centered within-cell bed source.
    const double cL = 0.5 * kGravity * (L.h * L.h - hLf * hLf - (hL_cell + hLf) * (zLf - zL));
    const double cR = 0.5 * kGravity * (R.h * R.h - hRf * hRf - (hR_cell + hRf) * (zRf - zR));
    corrL_x = cL * nx;
    corrL_y = cL * ny;
    corrR_x = -cR * nx;
    corrR_y = -cR * ny;
    if (L.h <= 0.0 && R.h <= 0.0) {
        out = FaceFlux{};
        return cL != 0.0 || cR != 0.0;
    }
    double fh, fn, ft, sstar;
    hllcFlux(L, R, fh, fn, ft, sstar);
    out.mass = fh;
    out.mx = fn * nx + ft * tx;
    out.my = fn * ny + ft * ty;
    out.sstar = sstar;
    return true;
}

OPENSWMM_KERNEL_FN bool faceFluxRecon(double etaLf, double uxLf, double uyLf, double zL,
                                      double hL_cell, double etaRf, double uxRf, double uyRf,
                                      double zR, double hR_cell, double nx, double ny, double h_dry,
                                      FaceFlux& out, double& corrL_x, double& corrL_y,
                                      double& corrR_x, double& corrR_y) noexcept {
    return faceFluxReconBed(etaLf, uxLf, uyLf, zL, hL_cell, etaRf, uxRf, uyRf, zR, hR_cell, nx, ny,
                            h_dry, out, corrL_x, corrL_y, corrR_x, corrR_y, zL, zR);
}

OPENSWMM_KERNEL_FN double bjLimiter(double w, double wf, double wmin, double wmax) noexcept {
    const double d = wf - w;
    if (d > 1.0e-14)
        return std::min(1.0, (wmax - w) / d);
    else if (d < -1.0e-14)
        return std::min(1.0, (wmin - w) / d);
    return 1.0;
}

OPENSWMM_KERNEL_FN void frictionUpdate(double& qx, double& qy, double h, double n,
                                       double dt) noexcept {
    if (!(h > 0.0)) {
        qx = qy = 0.0;
        return;
    }
    const double qm = std::sqrt(qx * qx + qy * qy);
    if (qm == 0.0)
        return;
    const double h73 = h * h * std::cbrt(h);
    const double den = 1.0 + kGravity * dt * n * n * qm / h73;
    qx /= den;
    qy /= den;
}

OPENSWMM_KERNEL_FN double criticalDepth(double q) noexcept {
    return std::cbrt(q * q / kGravity);
}

OPENSWMM_KERNEL_FN double depthFromInvariantAndDischarge(double r, double q,
                                                         double h_interior) noexcept {
    const double rg = std::sqrt(kGravity);
    double s = std::sqrt(h_interior > 0.0 ? h_interior : 1.0e-6);
    for (int it = 0; it < 40; ++it) {
        const double f = 2.0 * rg * s * s * s - r * s * s + q;
        const double df = 6.0 * rg * s * s - 2.0 * r * s;
        if (!(std::fabs(df) > 1.0e-12))
            break;
        const double s1 = s - f / df;
        if (!(s1 > 0.0) || !std::isfinite(s1))
            break;
        const double step = std::fabs(s1 - s);
        s = s1;
        if (step < 1.0e-12 * (1.0 + s)) {
            const double h = s * s;
            return (h > 0.0 && std::isfinite(h)) ? h : h_interior;
        }
    }
    const double h = s * s;
    // Accept only a converged, physical root; otherwise keep the interior
    // depth (transmissive), which is what the pre-invariant code did.
    return (h > 0.0 && std::isfinite(h) &&
            std::fabs(2.0 * rg * s * s * s - r * s * s + q) < 1.0e-6 * (1.0 + std::fabs(q)))
               ? h
               : h_interior;
}
 
} // namespace openswmm::twoD::swe
 
#endif // OPENSWMM_ENGINE_2D_SWE_KERNELS_HPP