RDII.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_RDII_HPP
#define OPENSWMM_RDII_HPP
 
#include <array>
#include <vector>
#include <string>
#include <unordered_map>
 
#include "../core/StringCase.hpp"
 
namespace openswmm {
 
struct SimulationContext;
 
namespace rdii {
 
struct UnitHydParams {
    double r[12][3]     = {};  
    double tPeak[12][3] = {};  
    double tBase[12][3] = {};  
    double iaMax[12][3] = {};  
    double iaRecov[12][3] = {};
    double iaInit[12][3] = {}; 
};

struct ExpDecayParams {
    bool   active    = false;
    double k_dep     = 0.0;   
    double k_0       = 0.0;   
    double k_T       = 0.0;   
    double T_ref     = 10.0;  
    double theta_rec = 0.0;   
    double T_freeze  = 0.0;   
 
    // Optional degree-day snow model. When `snow_on`, precipitation at
    // T <= snow_T accumulates as SWE with no liquid input; at T > snow_T,
    // melt = min(SWE, snow_ddf*(T - snow_T)*dt_days) is added to rainfall
    // (rain-on-snow) before the IA depletion step.
    bool   snow_on   = false; 
    double snow_T    = 1.0;   
    double snow_ddf  = 0.0;   
};

struct UHResponseData {
    std::vector<double> past_rain;      
    std::vector<int>    past_month;     
    int    period      = 0;             
    int    max_periods = 0;             
    int    has_past_rain = 0;           
    double ia_used     = 0.0;           
    double swe         = 0.0;           
    long   dry_seconds = 0;             
 
    void allocate(int n) {
        max_periods = n;
        past_rain.assign(static_cast<std::size_t>(n), 0.0);
        past_month.assign(static_cast<std::size_t>(n), 0);
        period = 0;
        has_past_rain = 0;
        ia_used = 0.0;
        swe = 0.0;
        dry_seconds = static_cast<long>(n) * 300 + 1; // start dry
    }
};

double getRecoveryRate(const ExpDecayParams& dp, double T_celsius);

double updateIA_exp(const UnitHydParams& uh, UHResponseData& rd,
                    const ExpDecayParams& dp, int month, int response,
                    double rainDepth, double dt_sec,
                    const SimulationContext& ctx);
 
struct RDIIGroupSoA {
    int count = 0;
    std::vector<int>    node_idx;       
    std::vector<int>    uh_idx;         
    std::vector<int>    gage_idx;       
    std::vector<double> area;           

    std::vector<UHResponseData> uh_data;
 
    std::vector<int>    rain_interval;  
 
    // Strict-grid driver state (legacy createRdiiFile embedded at runtime).
    // Chunk starts advance in rain_interval steps from simulation start;
    // rates are recorded per chunk while the runoff substep containing the
    // chunk's start is current (the substep never crosses a gage entry
    // boundary, so the current gage rate IS the rate at the chunk start).
    std::vector<double> gage_elapsed;         
    std::vector<double> rate_recorded_until;  
    std::vector<std::vector<double>> pending_rates;  
 
    void resize(int n);
};
 
class RDIISolver {
public:
    void init(SimulationContext& ctx);

    int addUnitHydParams(const std::string& name, const UnitHydParams& params);

    int findUnitHyd(const std::string& name) const;

    using UhNameMap = std::unordered_map<std::string, int, CiHash, CiEqual>;

    const UhNameMap& uhNameIndex() const { return uh_name_to_idx_; }

    void advance(SimulationContext& ctx, double new_elapsed_sec);

    void applyRdiiInflows(SimulationContext& ctx, double elapsed_sec) const;
 
    std::vector<UnitHydParams> uh_params;

    std::vector<std::array<ExpDecayParams, 3>> decay_params;

    const std::vector<double>& nodeFlows() const { return node_rdii_flow_; }

    std::vector<int> rdiiNodeList() const;

    static int getRainInterval(const UnitHydParams& uh, double wet_step);

    static int getMaxPeriods(const UnitHydParams& uh, int response, int rainInterval);
 
private:
    RDIIGroupSoA groups_;
    UhNameMap uh_name_to_idx_;  
    std::vector<double> node_rdii_flow_;  
 
    // Strict grid (legacy createRdiiFile records, kept in memory).
    double grid_step_ = 0.0;        
    double next_tick_ = 0.0;        
    std::vector<int>   grid_node_;  
    std::vector<float> grid_flows_; 

    void emitTick(SimulationContext& ctx, double T);

    double uhOrdinate(const UnitHydParams& uh, int month, int response, double t) const;

    void validateExpDecay(SimulationContext& ctx) const;
};
 
} // namespace rdii
} // namespace openswmm
 
#endif // OPENSWMM_RDII_HPP