VfrClosure.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_VFR_CLOSURE_HPP
#define OPENSWMM_ENGINE_2D_VFR_CLOSURE_HPP
 
#include <cmath>
 
// Portable kernel-function marker (P5 Kokkos port) — see InertialKernels.hpp.
// Host builds: plain inline. The GPU plugin defines OPENSWMM_KERNEL_FN to
// KOKKOS_INLINE_FUNCTION before including so the identical closure bodies are
// device-callable.
#ifndef OPENSWMM_KERNEL_FN
#define OPENSWMM_KERNEL_FN inline
#endif
 
namespace openswmm::twoD {

inline constexpr double kVfrFlatRelief = 1.0e-9;

OPENSWMM_KERNEL_FN void vfrSort3(double& z1, double& z2, double& z3) noexcept {
    if (z1 > z2) {
        const double t = z1;
        z1 = z2;
        z2 = t;
    }
    if (z2 > z3) {
        const double t = z2;
        z2 = z3;
        z3 = t;
    }
    if (z1 > z2) {
        const double t = z1;
        z1 = z2;
        z2 = t;
    }
}

OPENSWMM_KERNEL_FN double vfrWetFraction(double z1, double z2, double z3, double eta) noexcept {
    const double relief = z3 - z1;
    if (relief < kVfrFlatRelief)
        return (eta > z1) ? 1.0 : 0.0;
    if (eta <= z1)
        return 0.0;
    if (eta >= z3)
        return 1.0;
    if (eta <= z2) {
        const double d21 = z2 - z1;
        if (d21 < kVfrFlatRelief)
            return 0.0; // measure-zero band z1==z2
        const double t = eta - z1;
        return t * t / (d21 * relief);
    }
    // z2 < eta < z3 (z3 > z2 strictly here)
    const double d32 = z3 - z2;
    const double t = z3 - eta;
    return 1.0 - t * t / (relief * d32);
}

OPENSWMM_KERNEL_FN double vfrMeanDepthFromEtaExact(double z1, double z2, double z3,
                                                   double eta) noexcept {
    const double zbar = (z1 + z2 + z3) / 3.0;
    const double relief = z3 - z1;
    if (relief < kVfrFlatRelief)
        return (eta > zbar) ? (eta - zbar) : 0.0;
    if (eta <= z1)
        return 0.0;
    if (eta >= z3)
        return eta - zbar;
    if (eta <= z2) {
        const double d21 = z2 - z1;
        if (d21 < kVfrFlatRelief)
            return 0.0; // measure-zero band z1==z2
        const double t = eta - z1;
        return t * t * t / (3.0 * d21 * relief);
    }
    const double d32 = z3 - z2;
    const double t = z3 - eta;
    return (eta - zbar) + t * t * t / (3.0 * relief * d32);
}

OPENSWMM_KERNEL_FN double vfrStageAtWetFraction(double z1, double z2, double z3,
                                                double eps) noexcept {
    const double relief = z3 - z1;
    const double d21 = z2 - z1;
    // Lower branch reaches wet fraction (z2−z1)/(z3−z1) at η = z2.
    if (eps <= d21 / relief)
        return z1 + std::sqrt(eps * d21 * relief);
    const double d32 = z3 - z2;
    return z3 - std::sqrt((1.0 - eps) * relief * d32);
}

OPENSWMM_KERNEL_FN double vfrMeanDepthFromEta(double z1, double z2, double z3, double eta,
                                              double eps) noexcept {
    const double relief = z3 - z1;
    if (eps <= 0.0 || relief < kVfrFlatRelief)
        return vfrMeanDepthFromEtaExact(z1, z2, z3, eta);
    const double eta_s = vfrStageAtWetFraction(z1, z2, z3, eps);
    if (eta >= eta_s)
        return vfrMeanDepthFromEtaExact(z1, z2, z3, eta);
    const double h_s = vfrMeanDepthFromEtaExact(z1, z2, z3, eta_s);
    const double h = h_s - eps * (eta_s - eta);
    return (h > 0.0) ? h : 0.0;
}

OPENSWMM_KERNEL_FN double vfrEtaFromMeanDepth(double z1, double z2, double z3, double mean_depth,
                                              double eps) noexcept {
    const double zbar = (z1 + z2 + z3) / 3.0;
    const double relief = z3 - z1;
 
    // Flat (or degenerate) cell: the flat closure is exact.
    if (relief < kVfrFlatRelief)
        return zbar + ((mean_depth > 0.0) ? mean_depth : 0.0);
 
    // Fully wet: flat closure is exact. Checked FIRST — it is the dominant case
    // in deep water (e.g. flooded urban meshes) and skips the ε-tail's sqrt/cubic
    // switch-point evaluation below, which every cell would otherwise pay each
    // RHS. The fully-wet threshold z3−z̄ always exceeds the ε-tail depth h_s, so
    // reordering is exact (the two branches never overlap).
    if (mean_depth >= z3 - zbar)
        return zbar + mean_depth;
 
    // Regularized tail: below the switch depth h_s the closure is linear.
    double eta_s = z1, h_s = 0.0;
    if (eps > 0.0) {
        eta_s = vfrStageAtWetFraction(z1, z2, z3, eps);
        h_s = vfrMeanDepthFromEtaExact(z1, z2, z3, eta_s);
        const double h = (mean_depth > 0.0) ? mean_depth : 0.0;
        if (h <= h_s)
            return eta_s - (h_s - h) / eps;
    } else if (!(mean_depth > 0.0)) {
        return z1; // exact relation, dry limit
    }
 
    // Lower branch (z1 < η ≤ z2): closed-form cube root.
    const double d21 = z2 - z1;
    const double h_at_z2 = d21 * d21 / (3.0 * relief);
    if (mean_depth <= h_at_z2)
        return z1 + std::cbrt(3.0 * mean_depth * d21 * relief);
 
    // z2 == z3 (to rounding): the upper branch is an empty interval — the
    // fully-wet check above and h_at_z2 meet at R/3 — but 1-ulp slivers must
    // not reach the 1/(relief·d32) divisions below.
    if (z3 - z2 < kVfrFlatRelief)
        return zbar + mean_depth;
 
    // Upper branch (z2 < η < z3): safeguarded Newton on the bracket [z2, z3].
    // h̄ is strictly increasing (dh̄/dη = A_wet/A > 0 here) so this converges
    // unconditionally. Mirrors the pre-existing render-path iteration.
    const double d32 = z3 - z2;
    const double denom = 3.0 * relief * d32;
    double lo = z2, hi = z3;
    double eta = zbar + mean_depth; // flat-closure initial guess
    if (eta <= lo || eta >= hi)
        eta = 0.5 * (lo + hi);
    for (int it = 0; it < 64; ++it) {
        const double dz3 = z3 - eta;
        const double f = (eta - zbar) + dz3 * dz3 * dz3 / denom - mean_depth;
        if (f > 0.0)
            hi = eta;
        else
            lo = eta;
        const double df = 1.0 - dz3 * dz3 / (relief * d32); // A_wet/A
        double next = (df > 1.0e-12) ? eta - f / df : 0.5 * (lo + hi);
        if (next <= lo || next >= hi)
            next = 0.5 * (lo + hi); // safeguard
        if (std::abs(next - eta) < 1.0e-12 * (1.0 + relief))
            return next;
        eta = next;
    }
    return eta;
}

OPENSWMM_KERNEL_FN double vfrDryEta(double z1, double z2, double z3, double eps) noexcept {
    return vfrEtaFromMeanDepth(z1, z2, z3, 0.0, eps);
}

OPENSWMM_KERNEL_FN double vfrDEtaDMeanDepth(double z1, double z2, double z3, double eta,
                                            double eps) noexcept {
    double w = vfrWetFraction(z1, z2, z3, eta);
    if (w < eps)
        w = eps;
    if (w < 1.0e-12)
        w = 1.0e-12; // eps == 0 safety (exact relation, dry cell)
    return 1.0 / w;
}
 
} // namespace openswmm::twoD
 
#endif // OPENSWMM_ENGINE_2D_VFR_CLOSURE_HPP