SubcatchData.hpp

Go to the documentation of this file.

 
#ifndef OPENSWMM_ENGINE_SUBCATCH_DATA_HPP
#define OPENSWMM_ENGINE_SUBCATCH_DATA_HPP
 
#include <vector>
#include <string>
#include <cstdint>
#include <algorithm>
#include <limits>
 
namespace openswmm {
 
// ============================================================================
// SubcatchData — SoA layout
// ============================================================================
 
struct SubcatchData {
 
    // -----------------------------------------------------------------------
    // Static properties — set at parse time
    // -----------------------------------------------------------------------
 
    std::vector<int>    outlet_node;
 
    std::vector<int>    outlet_subcatch;
 
    std::vector<std::string> outlet_name;
 
    std::vector<int>    gage;
 
    std::vector<double> area;
 
    std::vector<double> width;
 
    std::vector<double> slope;
 
    std::vector<double> curb_length;
 
    std::vector<double> frac_imperv;
 
    std::vector<double> frac_imperv_no_store;
 
    std::vector<double> n_imperv;
 
    std::vector<double> n_perv;
 
    std::vector<double> ds_imperv;
 
    std::vector<double> ds_perv;
 
    // -----------------------------------------------------------------------
    // Inter-subarea routing (from [SUBAREAS] section)
    // -----------------------------------------------------------------------
 
    std::vector<int>    subarea_routing;
 
    std::vector<double> pct_routed;
 
    // -----------------------------------------------------------------------
    // Infiltration parameters (per subcatchment)
    // -----------------------------------------------------------------------
 
    std::vector<int>    infil_model;
 
    std::vector<double> infil_p1;
    std::vector<double> infil_p2;
    std::vector<double> infil_p3;
    std::vector<double> infil_p4;
    std::vector<double> infil_p5;
 
    // -----------------------------------------------------------------------
    // State variables — updated each timestep
    // -----------------------------------------------------------------------
 
    std::vector<double> runoff;
 
    std::vector<double> rainfall;
 
    std::vector<double> evap_loss;
 
    std::vector<double> infil_loss;
 
    std::vector<double> ponded_depth;
 
    std::vector<double> gw_flow;
 
    // -----------------------------------------------------------------------
    // Previous-step state
    // -----------------------------------------------------------------------
 
    std::vector<double> old_runoff;
 
    std::vector<double> old_gw_flow;
 
    // -----------------------------------------------------------------------
    // Runon coupling — subcatch-to-subcatch routing
    // -----------------------------------------------------------------------
 
    std::vector<double> runon_inflow;
 
    std::vector<double> old_runon_inflow;
 
    std::vector<double> outfall_runon_vol;
 
    // -----------------------------------------------------------------------
    // Groundwater surface water head coupling
    // -----------------------------------------------------------------------
 
    std::vector<double> gw_sw_head;
 
    std::vector<double> gw_node_avail_flow;
 
    std::vector<double> gw_max_infil_vol;
 
    // -----------------------------------------------------------------------
    // Quality washoff output
    // -----------------------------------------------------------------------
 
    std::vector<double> washoff_load;
 
    // -----------------------------------------------------------------------
    // Per-subcatch quality state — flat 2D: [subcatch * n_pollutants + pollutant]
    // -----------------------------------------------------------------------
 
    std::vector<double> conc;
 
    std::vector<double> conc_old;
 
    std::vector<double> ponded_qual;
 
    int                 conc_n_pollutants = 0;
 
    // -----------------------------------------------------------------------
    // Per-object INP comment
    // -----------------------------------------------------------------------
 
    std::vector<std::string> comments;
 
    std::vector<std::string> tags;
 
    // -----------------------------------------------------------------------
    // Report flag — per-object output filter
    // -----------------------------------------------------------------------
 
    std::vector<char>       rpt_flag;
 
    // -----------------------------------------------------------------------
    // Cumulative statistics
    // -----------------------------------------------------------------------
 
    std::vector<double> stat_precip_vol;
    std::vector<double> stat_evap_vol;     
    std::vector<double> stat_infil_vol;    
    std::vector<double> stat_imperv_vol;   
    std::vector<double> stat_perv_vol;     
 
    std::vector<double> stat_runoff_vol;
 
    std::vector<double> stat_max_runoff;
 
    // Gap #63: Per-subcatch groundwater statistics (accumulated each GW step)
    std::vector<double> stat_gw_infil_vol;      
    std::vector<double> stat_gw_upper_evap_vol; 
    std::vector<double> stat_gw_lower_evap_vol; 
    std::vector<double> stat_gw_deep_perc_vol;  
    std::vector<double> stat_gw_flow_vol;       
    std::vector<double> stat_gw_max_flow;       
    std::vector<double> stat_gw_sum_theta;      
    std::vector<double> stat_gw_sum_depth;      
    std::vector<double> stat_gw_final_theta;    
    std::vector<double> stat_gw_final_depth;    
    std::vector<long>   stat_gw_steps;          
 
    // -----------------------------------------------------------------------
    // Groundwater parameters (from [GROUNDWATER] section)
    // -----------------------------------------------------------------------
 
    std::vector<int>    gw_aquifer;
 
    std::vector<int>    gw_node;
 
    std::vector<double> gw_surf_elev;
 
    std::vector<double> gw_a1;
 
    std::vector<double> gw_b1;
 
    std::vector<double> gw_a2;
 
    std::vector<double> gw_b2;
 
    std::vector<double> gw_a3;
 
    std::vector<double> gw_tw;
 
    std::vector<double> gw_hstar;
 
    // -----------------------------------------------------------------------
    // Snowpack assignment (index into SnowpackStore, -1 = none)
    // -----------------------------------------------------------------------
 
    std::vector<int>    snowpack;
 
    std::vector<double> snow_net_imperv;
 
    std::vector<double> snow_net_perv;
 
    std::vector<double> total_lid_area_ft2;
 
    std::vector<double> lid_return_to_perv_cfs;
 
    std::vector<double> lid_drain_runon_cfs;
 
    // -----------------------------------------------------------------------
    // Cumulative pollutant washoff loads — flat 2D: [subcatch * n_pollutants + pollutant]
    // -----------------------------------------------------------------------
 
    std::vector<double> total_load;
 
    int total_load_n_pollutants = 0;
 
    void resize_total_load(int n_sc, int n_poll) {
        total_load_n_pollutants = n_poll;
        total_load.assign(
            static_cast<std::size_t>(n_sc) * static_cast<std::size_t>(n_poll), 0.0);
    }
 
    // -----------------------------------------------------------------------
    // Land-use coverage — flat 2D: [subcatch * n_landuses + landuse]
    // -----------------------------------------------------------------------
 
    std::vector<double> coverage;
 
    int coverage_n_landuses = 0;
 
    std::vector<double> sweep_last_swept;
 
    void resize_coverage(int n_sc, int n_lu) {
        coverage_n_landuses = n_lu;
        auto sz = static_cast<std::size_t>(n_sc) * static_cast<std::size_t>(n_lu);
        coverage.assign(sz, 0.0);
        sweep_last_swept.assign(sz, 0.0);
    }
 
    // -----------------------------------------------------------------------
    // Capacity management
    // -----------------------------------------------------------------------
 
    int count() const noexcept { return static_cast<int>(area.size()); }
 
    void resize(int n) {
        const auto un = static_cast<std::size_t>(n);
 
        outlet_node.assign(un, -1);
        outlet_subcatch.assign(un, -1);
        outlet_name.resize(un);
        gage.assign(un, -1);
        area.assign(un, 0.0);
        width.assign(un, 0.0);
        slope.assign(un, 0.0);
        curb_length.assign(un, 0.0);
        frac_imperv.assign(un, 0.0);
        frac_imperv_no_store.assign(un, 0.0);
        n_imperv.assign(un, 0.013);
        n_perv.assign(un, 0.1);
        ds_imperv.assign(un, 0.0);
        ds_perv.assign(un, 0.0);
        subarea_routing.assign(un, 0);  // TO_OUTLET
        pct_routed.assign(un, 0.0);
 
        infil_model.assign(un, 0);
        infil_p1.assign(un, 0.0);
        infil_p2.assign(un, 0.0);
        infil_p3.assign(un, 0.0);
        infil_p4.assign(un, 0.0);
        infil_p5.assign(un, 0.0);
 
        runoff.assign(un, 0.0);
        rainfall.assign(un, 0.0);
        evap_loss.assign(un, 0.0);
        infil_loss.assign(un, 0.0);
        ponded_depth.assign(un, 0.0);
        gw_flow.assign(un, 0.0);
        old_runoff.assign(un, 0.0);
        old_gw_flow.assign(un, 0.0);
        runon_inflow.assign(un, 0.0);
        old_runon_inflow.assign(un, 0.0);
        outfall_runon_vol.assign(un, 0.0);
        gw_sw_head.assign(un, 0.0);
        gw_node_avail_flow.assign(un, 0.0);
        gw_max_infil_vol.assign(un, std::numeric_limits<double>::max());
 
        comments.assign(un, std::string{});
        tags.assign(un, std::string{});
 
        rpt_flag.assign(un, 0);
 
        stat_precip_vol.assign(un, 0.0);
        stat_evap_vol.assign(un, 0.0);
        stat_infil_vol.assign(un, 0.0);
        stat_imperv_vol.assign(un, 0.0);
        stat_perv_vol.assign(un, 0.0);
        stat_runoff_vol.assign(un, 0.0);
        stat_max_runoff.assign(un, 0.0);
        stat_gw_infil_vol.assign(un, 0.0);
        stat_gw_upper_evap_vol.assign(un, 0.0);
        stat_gw_lower_evap_vol.assign(un, 0.0);
        stat_gw_deep_perc_vol.assign(un, 0.0);
        stat_gw_flow_vol.assign(un, 0.0);
        stat_gw_max_flow.assign(un, 0.0);
        stat_gw_sum_theta.assign(un, 0.0);
        stat_gw_sum_depth.assign(un, 0.0);
        stat_gw_final_theta.assign(un, 0.0);
        stat_gw_final_depth.assign(un, 0.0);
        stat_gw_steps.assign(un, 0L);
 
        gw_aquifer.assign(un, -1);
        gw_node.assign(un, -1);
        gw_surf_elev.assign(un, 0.0);
        gw_a1.assign(un, 0.0);
        gw_b1.assign(un, 0.0);
        gw_a2.assign(un, 0.0);
        gw_b2.assign(un, 0.0);
        gw_a3.assign(un, 0.0);
        gw_tw.assign(un, 0.0);
        gw_hstar.assign(un, 0.0);
        snowpack.assign(un, -1);
        snow_net_imperv.assign(un, -1.0);
        snow_net_perv.assign(un, -1.0);
        total_lid_area_ft2.assign(un, 0.0);
        lid_return_to_perv_cfs.assign(un, 0.0);
        lid_drain_runon_cfs.assign(un, 0.0);
    }
 
    void grow_to(int n) {
        if (n <= count()) return;
        const auto un = static_cast<std::size_t>(n);
        auto g = [&](auto& vec, auto def) { vec.resize(un, def); };
        g(outlet_node, -1); g(outlet_subcatch, -1);
        outlet_name.resize(un); g(gage, -1);
        g(area, 0.0); g(width, 0.0); g(slope, 0.0); g(curb_length, 0.0);
        g(frac_imperv, 0.0); g(frac_imperv_no_store, 0.0);
        g(n_imperv, 0.013); g(n_perv, 0.1);
        g(ds_imperv, 0.0); g(ds_perv, 0.0);
        g(subarea_routing, 0); g(pct_routed, 0.0);
        g(infil_model, 0);
        g(infil_p1, 0.0); g(infil_p2, 0.0); g(infil_p3, 0.0);
        g(infil_p4, 0.0); g(infil_p5, 0.0);
        g(runoff, 0.0); g(rainfall, 0.0);
        g(evap_loss, 0.0); g(infil_loss, 0.0);
        g(ponded_depth, 0.0); g(gw_flow, 0.0);
        g(old_runoff, 0.0); g(old_gw_flow, 0.0);
        g(runon_inflow, 0.0); g(old_runon_inflow, 0.0);
        g(gw_sw_head, 0.0); g(gw_node_avail_flow, 0.0);
        g(gw_max_infil_vol, std::numeric_limits<double>::max());
        g(outfall_runon_vol, 0.0);
        comments.resize(un, std::string{});
        tags.resize(un, std::string{});
 
        g(rpt_flag, static_cast<char>(0));
        g(stat_precip_vol, 0.0); g(stat_evap_vol, 0.0);
        g(stat_infil_vol, 0.0); g(stat_imperv_vol, 0.0);
        g(stat_perv_vol, 0.0); g(stat_runoff_vol, 0.0);
        g(stat_max_runoff, 0.0);
        g(stat_gw_infil_vol, 0.0); g(stat_gw_upper_evap_vol, 0.0);
        g(stat_gw_lower_evap_vol, 0.0); g(stat_gw_deep_perc_vol, 0.0);
        g(stat_gw_flow_vol, 0.0); g(stat_gw_max_flow, 0.0);
        g(stat_gw_sum_theta, 0.0); g(stat_gw_sum_depth, 0.0);
        g(stat_gw_final_theta, 0.0); g(stat_gw_final_depth, 0.0);
        stat_gw_steps.resize(un, 0L);
        g(gw_aquifer, -1); g(gw_node, -1); g(gw_surf_elev, 0.0);
        g(gw_a1, 0.0); g(gw_b1, 0.0); g(gw_a2, 0.0); g(gw_b2, 0.0);
        g(gw_a3, 0.0); g(gw_tw, 0.0); g(gw_hstar, 0.0);
        g(snowpack, -1);
        g(snow_net_imperv, -1.0);
        g(snow_net_perv,   -1.0);
        g(total_lid_area_ft2, 0.0);
        g(lid_return_to_perv_cfs, 0.0);
        g(lid_drain_runon_cfs, 0.0);
        // Note: conc, conc_old, ponded_qual, washoff_load handled by resize_quality()
        // Note: coverage, total_load handled separately
    }
 
    void erase_at(int idx) {
        const auto ui = static_cast<std::size_t>(idx);
        auto e = [&](auto& v) { if (ui < v.size()) v.erase(v.begin() + static_cast<std::ptrdiff_t>(idx)); };
 
        e(outlet_node); e(outlet_subcatch); e(outlet_name); e(gage);
        e(area); e(width); e(slope); e(curb_length);
        e(frac_imperv); e(frac_imperv_no_store); e(n_imperv); e(n_perv);
        e(ds_imperv); e(ds_perv); e(subarea_routing); e(pct_routed);
 
        e(infil_model); e(infil_p1); e(infil_p2); e(infil_p3); e(infil_p4); e(infil_p5);
 
        e(runoff); e(rainfall); e(evap_loss); e(infil_loss); e(ponded_depth);
        e(gw_flow); e(old_runoff); e(old_gw_flow);
        e(runon_inflow); e(old_runon_inflow); e(outfall_runon_vol);
        e(gw_sw_head); e(gw_node_avail_flow); e(gw_max_infil_vol);
        e(comments); e(tags); e(rpt_flag);
 
        e(stat_precip_vol); e(stat_evap_vol); e(stat_infil_vol);
        e(stat_imperv_vol); e(stat_perv_vol); e(stat_runoff_vol); e(stat_max_runoff);
        e(stat_gw_infil_vol); e(stat_gw_upper_evap_vol); e(stat_gw_lower_evap_vol);
        e(stat_gw_deep_perc_vol); e(stat_gw_flow_vol); e(stat_gw_max_flow);
        e(stat_gw_sum_theta); e(stat_gw_sum_depth); e(stat_gw_final_theta);
        e(stat_gw_final_depth); e(stat_gw_steps);
 
        e(gw_aquifer); e(gw_node); e(gw_surf_elev);
        e(gw_a1); e(gw_b1); e(gw_a2); e(gw_b2); e(gw_a3); e(gw_tw); e(gw_hstar);
        e(snowpack); e(snow_net_imperv); e(snow_net_perv);
        e(total_lid_area_ft2); e(lid_return_to_perv_cfs); e(lid_drain_runon_cfs);
 
        // Flat 2D quality arrays: [sc * np + p]
        if (conc_n_pollutants > 0) {
            const auto np = static_cast<std::size_t>(conc_n_pollutants);
            const auto base = ui * np;
            auto erase2d = [&](auto& v) {
                if (base + np <= v.size())
                    v.erase(v.begin() + static_cast<std::ptrdiff_t>(base),
                            v.begin() + static_cast<std::ptrdiff_t>(base + np));
            };
            erase2d(conc); erase2d(conc_old); erase2d(ponded_qual); erase2d(washoff_load);
        }
 
        // Flat 2D total load: [sc * np + p]
        if (total_load_n_pollutants > 0) {
            const auto np = static_cast<std::size_t>(total_load_n_pollutants);
            const auto base = ui * np;
            if (base + np <= total_load.size())
                total_load.erase(
                    total_load.begin() + static_cast<std::ptrdiff_t>(base),
                    total_load.begin() + static_cast<std::ptrdiff_t>(base + np));
        }
 
        // Flat 2D coverage/sweep matrices: [sc * n_lu + lu]
        if (coverage_n_landuses > 0) {
            const auto nlu = static_cast<std::size_t>(coverage_n_landuses);
            const auto base = ui * nlu;
            auto erase2d = [&](auto& v) {
                if (base + nlu <= v.size())
                    v.erase(v.begin() + static_cast<std::ptrdiff_t>(base),
                            v.begin() + static_cast<std::ptrdiff_t>(base + nlu));
            };
            erase2d(coverage); erase2d(sweep_last_swept);
        }
    }
 
    void resize_quality(int n_pollutants) {
        conc_n_pollutants = n_pollutants;
        if (n_pollutants > 0) {
            auto total = static_cast<std::size_t>(count()) *
                         static_cast<std::size_t>(n_pollutants);
            conc.assign(total, 0.0);
            conc_old.assign(total, 0.0);
            ponded_qual.assign(total, 0.0);
            washoff_load.assign(total, 0.0);
        }
    }
 
    void resize_washoff_load(int n_pollutants) {
        if (n_pollutants > 0) {
            washoff_load.assign(
                static_cast<std::size_t>(count()) *
                static_cast<std::size_t>(n_pollutants), 0.0);
        }
    }
 
    void shrink_to_fit() {
        outlet_node.shrink_to_fit();
        outlet_subcatch.shrink_to_fit();
        outlet_name.shrink_to_fit();
        gage.shrink_to_fit();
        area.shrink_to_fit();
        width.shrink_to_fit();
        slope.shrink_to_fit();
        curb_length.shrink_to_fit();
        frac_imperv.shrink_to_fit();
        frac_imperv_no_store.shrink_to_fit();
        n_imperv.shrink_to_fit();
        n_perv.shrink_to_fit();
        ds_imperv.shrink_to_fit();
        ds_perv.shrink_to_fit();
        subarea_routing.shrink_to_fit();
        pct_routed.shrink_to_fit();
 
        infil_model.shrink_to_fit();
        infil_p1.shrink_to_fit();
        infil_p2.shrink_to_fit();
        infil_p3.shrink_to_fit();
        infil_p4.shrink_to_fit();
        infil_p5.shrink_to_fit();
 
        runoff.shrink_to_fit();
        rainfall.shrink_to_fit();
        evap_loss.shrink_to_fit();
        infil_loss.shrink_to_fit();
        ponded_depth.shrink_to_fit();
        gw_flow.shrink_to_fit();
        old_runoff.shrink_to_fit();
        old_gw_flow.shrink_to_fit();
        runon_inflow.shrink_to_fit();
        old_runon_inflow.shrink_to_fit();
        outfall_runon_vol.shrink_to_fit();
        gw_sw_head.shrink_to_fit();
        gw_node_avail_flow.shrink_to_fit();
        gw_max_infil_vol.shrink_to_fit();
 
        comments.shrink_to_fit();
        tags.shrink_to_fit();
 
        rpt_flag.shrink_to_fit();
 
        stat_precip_vol.shrink_to_fit();
        stat_evap_vol.shrink_to_fit();
        stat_infil_vol.shrink_to_fit();
        stat_imperv_vol.shrink_to_fit();
        stat_perv_vol.shrink_to_fit();
        stat_runoff_vol.shrink_to_fit();
        stat_max_runoff.shrink_to_fit();
 
        gw_aquifer.shrink_to_fit();
        gw_node.shrink_to_fit();
        gw_surf_elev.shrink_to_fit();
        gw_a1.shrink_to_fit();
        gw_b1.shrink_to_fit();
        gw_a2.shrink_to_fit();
        gw_b2.shrink_to_fit();
        gw_a3.shrink_to_fit();
        gw_tw.shrink_to_fit();
        gw_hstar.shrink_to_fit();
        snowpack.shrink_to_fit();
        snow_net_imperv.shrink_to_fit();
        snow_net_perv.shrink_to_fit();
 
        conc.shrink_to_fit();
        conc_old.shrink_to_fit();
        ponded_qual.shrink_to_fit();
        washoff_load.shrink_to_fit();
        total_load.shrink_to_fit();
        coverage.shrink_to_fit();
        sweep_last_swept.shrink_to_fit();
    }
 
    void save_state() noexcept {
        std::copy(runoff.begin(),        runoff.end(),        old_runoff.begin());
        std::copy(runon_inflow.begin(),  runon_inflow.end(),  old_runon_inflow.begin());
        std::copy(conc.begin(),          conc.end(),          conc_old.begin());
    }
 
    void reset_state() noexcept {
        std::fill(runoff.begin(),       runoff.end(),       0.0);
        std::fill(rainfall.begin(),     rainfall.end(),     0.0);
        std::fill(snow_net_imperv.begin(), snow_net_imperv.end(), -1.0);
        std::fill(snow_net_perv.begin(),   snow_net_perv.end(),   -1.0);
        std::fill(evap_loss.begin(),    evap_loss.end(),    0.0);
        std::fill(infil_loss.begin(),   infil_loss.end(),   0.0);
        std::fill(ponded_depth.begin(), ponded_depth.end(), 0.0);
        std::fill(old_runoff.begin(),   old_runoff.end(),   0.0);
        std::fill(old_gw_flow.begin(),  old_gw_flow.end(),  0.0);
        std::fill(runon_inflow.begin(), runon_inflow.end(), 0.0);
        std::fill(old_runon_inflow.begin(), old_runon_inflow.end(), 0.0);
        std::fill(gw_sw_head.begin(),   gw_sw_head.end(),   0.0);
        std::fill(gw_node_avail_flow.begin(), gw_node_avail_flow.end(), 0.0);
        std::fill(washoff_load.begin(), washoff_load.end(), 0.0);
        std::fill(conc.begin(), conc.end(), 0.0);
        std::fill(conc_old.begin(), conc_old.end(), 0.0);
    }
};
 
} /* namespace openswmm */
 
#endif /* OPENSWMM_ENGINE_SUBCATCH_DATA_HPP */