PondModel.hpp

Go to the documentation of this file.

 
#ifndef OPENSWMM_ENGINE_2D_POND_MODEL_HPP
#define OPENSWMM_ENGINE_2D_POND_MODEL_HPP
 
#include <cmath>
 
#if defined(KOKKOS_INLINE_FUNCTION)
#define OPENSWMM_POND_HD KOKKOS_INLINE_FUNCTION
#else
#define OPENSWMM_POND_HD inline
#endif
 
namespace openswmm::twoD::pond {

OPENSWMM_POND_HD double area(double x, double a_eq, double band) {
    if (!(x > 0.0))
        return 0.0;
    if (!(band > 0.0) || x >= band)
        return a_eq;
    const double t = x / band;
    return a_eq * t * t * (3.0 - 2.0 * t);
}

OPENSWMM_POND_HD double volume(double x, double a_eq, double band) {
    if (!(x > 0.0))
        return 0.0;
    if (!(band > 0.0))
        return a_eq * x;
    if (x >= band)
        return a_eq * (x - 0.5 * band);
    const double t = x / band;
    return a_eq * band * t * t * t * (1.0 - 0.5 * t);
}

OPENSWMM_POND_HD double depthForVolume(double v, double a_eq, double band) {
    if (!(v > 0.0) || !(a_eq > 0.0))
        return 0.0;
    if (!(band > 0.0))
        return v / a_eq;
    const double full = 0.5 * a_eq * band; // volume at x = band
    if (v >= full)
        return v / a_eq + 0.5 * band;
    // t³ − t⁴/2 = r on [0, 1], monotone: safeguarded Newton from t ≈ ∛r.
    const double r = v / (a_eq * band);
    double lo = 0.0, hi = 1.0, t = std::cbrt(r);
    if (t > 1.0)
        t = 1.0;
    for (int it = 0; it < 60; ++it) {
        const double f = t * t * t * (1.0 - 0.5 * t) - r;
        if (f > 0.0)
            hi = t;
        else
            lo = t;
        const double df = t * t * (3.0 - 2.0 * t);
        double tn = (df > 0.0) ? t - f / df : 0.5 * (lo + hi);
        if (!(tn > lo && tn < hi))
            tn = 0.5 * (lo + hi);
        if (std::fabs(tn - t) <= 1e-15) {
            t = tn;
            break;
        }
        t = tn;
    }
    return t * band;
}

OPENSWMM_POND_HD double equalizeCap(double dh, double window, double a1) {
    if (!(window > 0.0) || !(a1 > 0.0))
        return HUGE_VAL;
    return std::fabs(dh) * a1 / window;
}
 
} // namespace openswmm::twoD::pond
 
#undef OPENSWMM_POND_HD
 
#endif // OPENSWMM_ENGINE_2D_POND_MODEL_HPP