XSectLookup.hpp

Go to the documentation of this file.

 
#ifndef OPENSWMM_XSECT_LOOKUP_HPP
#define OPENSWMM_XSECT_LOOKUP_HPP
 
namespace openswmm::xsect {

inline double lookup_exact(double x, const double* t, int n) noexcept {
    const double delta = 1.0 / static_cast<double>(n - 1);   // same rounding as legacy
    int i = static_cast<int>(x / delta);                     // DIVIDE (legacy), not x*inv_delta
    if (i >= n - 1) return t[n - 1];
 
    double x0 = i * delta;
    double y  = t[i] + (x - x0) * (t[i + 1] - t[i]) / delta; // grouping matches legacy
 
    if (i < 2) {                                             // quadratic refinement for low x
        double x1 = (static_cast<double>(i) + 1.0) * delta;
        double y2 = y + (x - x0) * (x - x1) / (delta * delta) *
                        (t[i] / 2.0 - t[i + 1] + t[i + 2] / 2.0);
        if (y2 > 0.0) y = y2;
    }
    if (y < 0.0) y = 0.0;
    return y;
}
 
// ============================================================================
// Fast mode (opt-in, NOT bit-exact) — plan §6.
//
// Enabled per-build with -DSWMM_XSECT_FAST_LOOKUP. Substitutes the baseline
// reciprocal-multiply form: normalize by a precomputed 1/y_full and index /
// interpolate with `* inv_delta` (inv_delta = N-1 is an EXACT integer) instead
// of `/ delta`. This removes ~3 IEEE divisions per element (measured ~25% faster
// in the fused circular kernel on arm64) at the cost of ULP-level deviation from
// legacy (measured well inside the abs 1e-8 / rel 1e-7 tolerance gate — see
// tests/unit/engine/test_xsect_parity.cpp XSectFastMode). Default OFF: the
// bit-exact `lookup_exact` path above is the shipped default.
//
// `lookup_fast` is always compiled (so the tolerance test can assert its bound
// regardless of the active default), but is only *wired into the kernels* when
// SWMM_XSECT_FAST_LOOKUP is defined.
// ============================================================================

inline double lookup_fast(double x, const double* t, int n) noexcept {
    const double inv_delta = static_cast<double>(n - 1);   // exact integer
    const double delta     = 1.0 / inv_delta;
    if (x < 0.0) x = 0.0; else if (x > 1.0) x = 1.0;       // fast mode pre-clamps
    int i = static_cast<int>(x * inv_delta);
    if (i >= n - 1) return t[n - 1];
    double x0 = i * delta;
    double y  = t[i] + (x - x0) * (t[i + 1] - t[i]) * inv_delta;
    if (i < 2) {
        double x1 = (static_cast<double>(i) + 1.0) * delta;
        double y2 = y + (x - x0) * (x - x1) * (inv_delta * inv_delta) *
                        (t[i] / 2.0 - t[i + 1] + t[i + 2] / 2.0);
        if (y2 > 0.0) y = y2;
    }
    if (y < 0.0) y = 0.0;
    return y;
}
 
// Mode-aware normalization. In the default (bit-exact) build `param` is y_full
// and x = y/yFull via IEEE divide; under SWMM_XSECT_FAST_LOOKUP `param` is the
// precomputed inv_y_full and x = depth*inv (clamped). The kernels feed the
// matching per-link array (see XSectBatch.cpp norm_param()).
inline double norm_x(double depth, double param) noexcept {
#ifdef SWMM_XSECT_FAST_LOOKUP
    double x = depth * param;                              // param == inv_y_full
    if (x < 0.0) x = 0.0; else if (x > 1.0) x = 1.0;
    return x;
#else
    return (param > 0.0 && depth > 0.0) ? (depth / param) : 0.0;  // param == y_full
#endif
}

inline double norm_lookup(double depth, double param, const double* t, int n) noexcept {
#ifdef SWMM_XSECT_FAST_LOOKUP
    return lookup_fast(norm_x(depth, param), t, n);
#else
    return lookup_exact(norm_x(depth, param), t, n);
#endif
}
 
} // namespace openswmm::xsect
 
#endif // OPENSWMM_XSECT_LOOKUP_HPP