ExchangeLaw.hpp

Go to the documentation of this file.

 
#ifndef OPENSWMM_ENGINE_2D_EXCHANGE_LAW_HPP
#define OPENSWMM_ENGINE_2D_EXCHANGE_LAW_HPP
 
#include <cmath>
 
#if defined(KOKKOS_INLINE_FUNCTION)
#define OPENSWMM_EXCHANGE_HD KOKKOS_INLINE_FUNCTION
#else
#define OPENSWMM_EXCHANGE_HD inline
#endif
 
namespace openswmm::twoD::exchange {
 
constexpr double kGravity = 9.80665; 
constexpr double kSigmaBand = 0.1;   

struct Geometry {
    double z_crest = 0.0; 
    double l_weir = 0.0;  
    double a_orif = 0.0;  
};

struct Coefficients {
    double c_w = 0.540;  
    double c_sw = 0.360; 
    double c_oi = 0.167; 
    double c_oo = 0.160; 
    double h_eps = 0.02; 
};

OPENSWMM_EXCHANGE_HD double phi(double x, double eps) {
    if (!(x > 0.0))
        return 0.0;
    if (x >= eps)
        return std::sqrt(x);
    const double inv_sqrt_e = 1.0 / std::sqrt(eps);
    return (1.5 * inv_sqrt_e) * x - (0.5 * inv_sqrt_e / eps) * x * x;
}

OPENSWMM_EXCHANGE_HD double dphi(double x, double eps) {
    if (!(x > 0.0))
        return 1.5 / std::sqrt(eps);
    if (x >= eps)
        return 0.5 / std::sqrt(x);
    const double inv_sqrt_e = 1.0 / std::sqrt(eps);
    return 1.5 * inv_sqrt_e - (inv_sqrt_e / eps) * x;
}

OPENSWMM_EXCHANGE_HD double smoothstep(double t) {
    t = t < 0.0 ? 0.0 : (t > 1.0 ? 1.0 : t);
    return t * t * (3.0 - 2.0 * t);
}

OPENSWMM_EXCHANGE_HD double softmin(double a, double b, double& wa, double& wb) {
    wa = wb = 0.0;
    if (!(a > 0.0) || !(b > 0.0))
        return 0.0;
    const bool a_min = a <= b;
    const double lo = a_min ? a : b;
    const double r = lo / (a_min ? b : a);
    double r16 = r * r;
    r16 *= r16;
    r16 *= r16;
    r16 *= r16;
    const double root = std::sqrt(std::sqrt(std::sqrt(std::sqrt(1.0 + r16))));
    const double s = 1.0 / root; // Q / min
    const double q = lo * s;
    // ∂Q/∂x = (Q/x)^17 for each argument x.
    double s17 = s * s;
    s17 *= s17;
    s17 *= s17;
    s17 *= s17;
    s17 *= s;
    const double rs = r * s;
    double rs17 = rs * rs;
    rs17 *= rs17;
    rs17 *= rs17;
    rs17 *= rs17;
    rs17 *= rs;
    if (a_min) {
        wa = s17;
        wb = rs17;
    } else {
        wb = s17;
        wa = rs17;
    }
    return q;
}

struct Result {
    double q = 0.0;
    double g = 0.0;
};

OPENSWMM_EXCHANGE_HD Result exchange(const Geometry& geo, const Coefficients& c, double h2,
                                     double h1) {
    Result r;
    const double s2g = std::sqrt(2.0 * kGravity);
    const double z = geo.z_crest;
    const double eps = c.h_eps;
 
    if (h2 > h1) {
        const double hs = h2 - z;
        if (!(hs > 0.0))
            return r;
        const double dh = h2 - h1;
        const double S = (h1 > z) ? smoothstep((h1 - z) / (kSigmaBand * hs)) : 0.0;
 
        const double weir_free = (2.0 / 3.0) * c.c_w * geo.l_weir * s2g * hs * phi(hs, eps);
        const double weir_sub = c.c_sw * geo.l_weir * s2g * hs * phi(dh, eps);
        const double orif_free = c.c_oi * geo.a_orif * s2g * phi(hs, eps);
        const double orif_sub = c.c_oi * geo.a_orif * s2g * phi(dh, eps);
        const double qw = (1.0 - S) * weir_free + S * weir_sub;
        const double qo = (1.0 - S) * orif_free + S * orif_sub;
 
        double ww = 0.0, wo = 0.0;
        r.q = softmin(qw, qo, ww, wo);
        // Only the submerged forms depend on h1.
        const double dqw = S * c.c_sw * geo.l_weir * s2g * hs * dphi(dh, eps);
        const double dqo = S * c.c_oi * geo.a_orif * s2g * dphi(dh, eps);
        r.g = ww * dqw + wo * dqo;
        return r;
    }
    if (h1 > h2 && h1 > z) {
        const double head_dry = h1 - z;
        const double S = (h2 > z) ? smoothstep((h2 - z) / (kSigmaBand * head_dry)) : 0.0;
        const double k = c.c_oo * geo.a_orif * s2g;
        r.q = -k * ((1.0 - S) * phi(head_dry, eps) + S * phi(h1 - h2, eps));
        r.g = k * ((1.0 - S) * dphi(head_dry, eps) + S * dphi(h1 - h2, eps));
        return r;
    }
    return r;
}
 
} // namespace openswmm::twoD::exchange
 
#undef OPENSWMM_EXCHANGE_HD
 
#endif // OPENSWMM_ENGINE_2D_EXCHANGE_LAW_HPP