Runoff.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_RUNOFF_HPP
#define OPENSWMM_RUNOFF_HPP
 
#include "../data/SubcatchData.hpp"
#include "surface/InfilBank.hpp"
#include <vector>
#include <array>
#include <functional>
#include <utility>
 
namespace openswmm::detail {
class CompletedSourceArchive;
}
 
namespace openswmm {
 
struct SimulationContext;
 
namespace runoff {

struct RunoffSourceForcing {
    int subcatch = -1;
    double rain = 0.0, pet = 0.0;
    std::array<double, 3> surface_return{}; // ft/s by IMPERV0/IMPERV1/PERV; never rain
    double raw_precipitation = -1;          // before subcatchment scales and phase partition
    double pet_imperv = -1, pet_perv = -1;  // optional exposed atmospheric footprints
};

struct CompletedSubareaWater {
    double area = 0, opening = 0, rain = 0, runon = 0, evaporation = 0;
    double infiltration = 0, spatial_infiltration = 0, outflow = 0, outlet = 0, closing = 0;
};
 
// ============================================================================
// Constants
// ============================================================================
 
constexpr double MEXP =
    1.6666667; 
constexpr double ODETOL = 0.0001; 
constexpr double PHI = 1.486;     
 
// ============================================================================
// Per-subcatchment subarea state (SoA for vectorization)
// ============================================================================
 
struct RunoffSoA {
    int n_subcatch = 0;
 
    // Per-subcatchment properties (set at init)
    std::vector<double> area;       
    std::vector<double> width;      
    std::vector<double> slope;      
    std::vector<double> imperv_pct; 
    std::vector<double> frac_imperv0; 
    std::vector<double> frac_imperv1; 
 
    // Per-subarea SoA: alpha = runoff coefficient
    std::vector<double> alpha_imperv; 
    std::vector<double> alpha_perv;   
 
    // Per-subarea SoA: depression storage (ft)
    std::vector<double> ds_imperv; 
    std::vector<double> ds_perv;   
 
    // Per-subarea SoA: Manning's n
    std::vector<double> n_imperv; 
    std::vector<double> n_perv;   
 
    // Per-subarea SoA: ponded depth (state, updated each step)
    std::vector<double> depth_imperv0; 
    std::vector<double> depth_imperv1; 
    std::vector<double> depth_perv;    
 
    // Per-subarea SoA: previous-step runoff rates (ft/sec, area-averaged)
    // Used for inter-subarea routing (legacy subcatch_getRunon)
    std::vector<double> old_runoff_imperv0; 
    std::vector<double> old_runoff_imperv1; 
    std::vector<double> old_runoff_perv;    
 
    // Per-subcatchment: computed runoff (output)
    std::vector<double> runoff;    
    std::vector<double> evap_loss; 
    std::vector<double>
        infil_loss; 
    std::vector<double>
        perv_evap_vol; 
    std::vector<double>
        actual_perv_evap_vol; 
    std::vector<double> spatial_infil_vol, native_infil_rate;
    std::vector<double> infil_vol; 
    std::vector<double> subarea_runoff_rate;
 
    // Per-subcatchment: per-subarea runoff CFS from non-LID area (Gap #23)
    // Used by SWMMEngine to compute LID unit inflow from impervious/pervious fractions.
    std::vector<double> imperv_runoff_cfs; 
    std::vector<double> perv_runoff_cfs;   
    std::vector<double> outflow_vol;
    std::vector<std::array<CompletedSubareaWater, 3>> completed;

    void resize(int n);
    void computeAlpha();
};
 
// ============================================================================
// Runoff solver
// ============================================================================

class RunoffSolver {
public:
    using InfiltrationBoundary = std::function<bool(int, double, double, double&)>;
    void init(SimulationContext& ctx,
              const std::vector<std::pair<int, double>>& spatial_areas = {});
    std::vector<std::pair<int, double>> reviewedAreas() const;
    void execute(SimulationContext& ctx, double dt, double evap_rate = 0.0,
                 double infil_factor = 1.0, double recovery_factor = 1.0, int month = -1,
                 const InfiltrationBoundary* boundary = nullptr,
                 const std::vector<RunoffSourceForcing>* source_forcing = nullptr,
                 const std::vector<double>* spatial_fractions = nullptr,
                 const std::vector<int>* source_owners = nullptr);

    double nativeInfilFullLid(SimulationContext& ctx, int i, double dt, double recovery_factor);

    double infilFactorUsed(int i) const { return infil_factor_used_[static_cast<std::size_t>(i)]; }

    const RunoffSoA& soa() const { return soa_; }
    RunoffSoA& soa() { return soa_; }
    void checkpoint(detail::CompletedSourceArchive& archive, bool model);
    void validateCheckpoint(const RunoffSolver& model) const;
    double pendingRoutingVolume(const SimulationContext& ctx, int source) const;
 
    // -----------------------------------------------------------------------
    // Hot start helpers — Gap #54
    // -----------------------------------------------------------------------

    void infil_get_state(int i, int& model, double state[6]) const noexcept;

    void infil_set_state(int i, int model, const double state[6]) noexcept;
 
private:
    RunoffSoA soa_;
    std::vector<double>
        spatial_full_area_ft2_; 
 
    // Infiltration state (one per subcatchment)
    surface::InfilBank infil_bank_;
    std::vector<double>
        infil_factor_used_; 
    std::vector<double> completed_step_seconds_; 

    double infilGetInfil(SimulationContext& ctx, int i, double precip, double runon, double depth,
                         double dt, double local_infil, double recovery_factor);
 
    // Working buffers (reused each step, sized to n_subcatch)
    std::vector<double> precip_;
    std::vector<double> evap_rate_;
    std::vector<double> infil_rate_;

    static void updatePondedDepth(double& depth, double inflow, double alpha, double dStore,
                                  double dt, double& t_runoff);
};
 
} // namespace runoff
} // namespace openswmm
 
#endif // OPENSWMM_RUNOFF_HPP