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