SimulationContext.hpp

Go to the documentation of this file.

 
#ifndef OPENSWMM_ENGINE_SIMULATION_CONTEXT_HPP
#define OPENSWMM_ENGINE_SIMULATION_CONTEXT_HPP
 
#include <cmath>
#include <functional>
#include "FilePathPair.hpp"
#include "../data/GageData.hpp"
#include "../data/LinkData.hpp"
#include "../data/NameIndex.hpp"
#include "../data/NodeData.hpp"
#include "../data/NodeSubtypes.hpp"
#include "../data/LinkSubtypes.hpp"
#include "../data/PollutantData.hpp"
#include "../data/SubcatchData.hpp"
#include "../data/TableData.hpp"
#include "SimulationOptions.hpp"
#include "SpatialFrame.hpp"
#include "UserFlags.hpp"
#include "../data/QualityData.hpp"
#include "../data/InflowData.hpp"
#include "../data/InfraData.hpp"
#include "../hydraulics/Transect.hpp"
#include "../data/HydrologyData.hpp"
#include "../data/ForcingData.hpp"
#include "../hydrology/Climate.hpp"
 
#include <algorithm>
#include <cstdint>
#include <string>
#include <vector>
 
namespace openswmm {
 
// ============================================================================
// Plugin specification (from [PLUGINS] section)
// ============================================================================

struct PluginSpec {
    std::string              path;       
    std::vector<std::string> init_args;  
};
 
// ============================================================================
// [FILES] section spec — secondary file references
// ============================================================================

enum class FileMode {
    NONE,
    SAVE,
    USE
};

struct HotstartSaveEntry {
    FilePathPair path;
    double       datetime = 0.0;
};
 
struct FilesSpec {
    FileMode     rainfall_mode = FileMode::NONE;
    FilePathPair rainfall_path;
 
    FileMode     runoff_mode   = FileMode::NONE;
    FilePathPair runoff_path;
 
    FileMode     rdii_mode     = FileMode::NONE;
    FilePathPair rdii_path;

    FilePathPair inflows_path;

    FilePathPair outflows_path;

    FilePathPair hotstart_use_path;

    std::vector<HotstartSaveEntry> hotstart_saves;

    [[nodiscard]] bool has_any() const noexcept {
        return rainfall_mode != FileMode::NONE
            || runoff_mode   != FileMode::NONE
            || rdii_mode     != FileMode::NONE
            || !inflows_path.empty()
            || !outflows_path.empty()
            || !hotstart_use_path.empty()
            || !hotstart_saves.empty();
    }
};
 
// ============================================================================
// Engine state enumeration
// ============================================================================

enum class EngineState : int32_t {
    CREATED      = 0,  
    OPENED       = 1,  
    INITIALIZED  = 2,  
    RUNNING      = 3,  
    PAUSED       = 4,  
    ENDED        = 5,  
    REPORTED     = 6,  
    CLOSED       = 7,  
    ERROR_STATE  = 8,  
    BUILDING     = 9   
};
 
// ============================================================================
// StateAccessors
// ============================================================================

struct StateAccessors {
    std::function<bool(int subcatch_index, int& out_model, double* out_infil)> get_infil_state;

    std::function<bool(int subcatch_index, int model, const double* infil)> set_infil_state;

    std::function<bool(int subcatch_index, double& out_theta, double& out_lower_depth)> get_gw_state;

    std::function<bool(int subcatch_index, double theta, double lower_depth)> set_gw_state;

    bool can_read() const noexcept {
        return static_cast<bool>(get_infil_state) || static_cast<bool>(get_gw_state);
    }

    bool can_write() const noexcept {
        return static_cast<bool>(set_infil_state) || static_cast<bool>(set_gw_state);
    }
};
 
// ============================================================================
// 2D model-IO bridge (forward declarations)
// ============================================================================
 
// Plain 2D data structs (all header-only and define-free; see
// src/engine/2d/data/). Forward-declared so SimulationContext can carry
// non-owning pointers without pulling the 2D module into every TU.
namespace twoD {
struct MeshData;
struct SolverOptions2D;
struct BoundaryData;
struct PendingBoundaryRow;
struct PendingEdgeConveyanceRow;
} // namespace twoD
 
// ============================================================================
// SimulationContext
// ============================================================================

struct SimulationContext {
 
    // =========================================================================
    // Engine state
    // =========================================================================

    EngineState state = EngineState::CREATED;
 
    // =========================================================================
    // Project title / notes
    // =========================================================================

    std::vector<std::string> title_notes;
 
    // =========================================================================
    // Options & configuration
    // =========================================================================

    SimulationOptions options;

    climate::ClimateState climate_state;
 
    // =========================================================================
    // Simulation clock
    // =========================================================================

    double current_time = 0.0;

    double current_date = 0.0;

    double elapsed_ms = 0.0;

    double old_elapsed_ms = 0.0;

    double next_report_ms = 0.0;

    double dt_controls_remaining = 0.0;
 
    // =========================================================================
    // Object data stores (Structure-of-Arrays)
    // =========================================================================

    NodeData nodes;

    NodeSubtypes node_subtypes;

    LinkData links;

    LinkSubtypes link_subtypes;

    SubcatchData subcatches;

    GageData gages;

    PollutantData pollutants;

    TableData tables;
 
    // =========================================================================
    // Name-to-index lookup (O(1))
    // =========================================================================

    NameIndex node_names;

    NameIndex link_names;

    NameIndex subcatch_names;

    NameIndex gage_names;

    NameIndex pollutant_names;

    NameIndex table_names;
 
    // =========================================================================
    // Quality data (landuse, buildup, washoff, treatment)
    // =========================================================================
 
    NameIndex        landuse_names;
    LanduseData      landuses;
    BuildupData      buildup;
    WashoffData      washoff;
    TreatmentData    treatment;
 
    // =========================================================================
    // Inflow data (external, DWF, RDII, patterns)
    // =========================================================================
 
    ExtInflowData    ext_inflows;
    DwfData          dwf_inflows;
    RDIIAssignData   rdii_assigns;
    UnitHydData      unit_hyds;      
    RDIIDecayData    rdii_decay;     
    PatternData      patterns;
 
    // =========================================================================
    // Infrastructure data (transects, streets, inlets, controls)
    // =========================================================================
 
    TransectStore    transects;
    std::vector<transect::TransectData> transect_tables;
    std::uint64_t xsect_generation = 0;

    bool gpkg_units_internal = false;
 
    StreetStore      streets;
    InletStore       inlets;
    InletUsageStore  inlet_usages;
    ControlRuleStore control_rules;
 
    // =========================================================================
    // Events (from [EVENTS] section)
    // =========================================================================

    struct Event {
        double start = 0.0;  
        double end   = 0.0;  
    };
    std::vector<Event> events;
 
    // =========================================================================
    // Subcatchment adjustment patterns (from [ADJUSTMENTS] section)
    // =========================================================================

    double adjust_temp[12]   = {0,0,0,0,0,0,0,0,0,0,0,0};
    double adjust_evap[12]   = {1,1,1,1,1,1,1,1,1,1,1,1};
    double adjust_rain[12]   = {1,1,1,1,1,1,1,1,1,1,1,1};
    double adjust_hydcon[12] = {1,1,1,1,1,1,1,1,1,1,1,1};

    std::vector<int> subcatch_n_perv_pattern;
    std::vector<int> subcatch_d_store_pattern;
    std::vector<int> subcatch_infil_pattern;

    std::vector<double> base_n_perv;
    std::vector<double> base_ds_perv;
    bool has_subcatch_adj_patterns = false; 
 
    // =========================================================================
    // Hydrology data (snowpacks, aquifers, LID)
    // =========================================================================
 
    SnowpackStore    snowpacks;
    NameIndex        snowpack_names;
    AquiferStore     aquifers;
    NameIndex        aquifer_names;
    LidControlStore  lid_controls;
    NameIndex        lid_names;
    LidUsageStore    lid_usage;
 
    // =========================================================================
    // Spatial data
    // =========================================================================

    SpatialFrame spatial;
 
    // =========================================================================
    // User-defined flags
    // =========================================================================

    UserFlags user_flags;
 
    // =========================================================================
    // Object tags (from [TAGS] section)
    // =========================================================================
 
    // Tags from the [TAGS] section now live per-index on
    // NodeData::tags / LinkData::tags / SubcatchData::tags. Storing
    // them name-keyed here was a latent rename bug — a tagged node
    // would lose its tag the moment `swmm_node_rename` was called.
 
    // =========================================================================
    // Runtime forcing data
    // =========================================================================

    ForcingData forcing;
 
    // =========================================================================
    // Plugin specifications (from [PLUGINS] section)
    // =========================================================================

    std::vector<PluginSpec> plugin_specs;

    FilesSpec files;

    StateAccessors state_accessors;

    struct TwoDModelIO {
        twoD::MeshData*                              mesh       = nullptr;
        twoD::SolverOptions2D*                       options    = nullptr;
        twoD::BoundaryData*                          boundary   = nullptr;
        std::vector<twoD::PendingBoundaryRow>*       pending_bc = nullptr;
        std::vector<twoD::PendingEdgeConveyanceRow>* pending_ec = nullptr;
    } twod_io;
 
    // =========================================================================
    // Error / warning tracking
    // =========================================================================

    int error_code = 0;

    int warning_code = 0;

    std::string error_message;

    std::vector<std::string> warnings;

    std::vector<std::string> errors;

    std::vector<std::uint8_t> coupled_node;
 
    // =========================================================================
    // Mass balance accumulators (SoA — vectorisable batch updates)
    // =========================================================================

    struct MassBalance {
        // Runoff totals
        double runoff_rainfall   = 0.0;  
        double runoff_evap       = 0.0;  
        double runoff_infil      = 0.0;  
        double runoff_runoff     = 0.0;  
        double runoff_snowremov  = 0.0;  
        double runoff_init_store = 0.0;  
        double runoff_final_store= 0.0;  
 
        // Routing totals
        double routing_dry_weather   = 0.0;
        double routing_wet_weather   = 0.0;
        double routing_gw_inflow     = 0.0;
        double routing_rdii          = 0.0;
        double routing_external      = 0.0;
        double routing_flooding      = 0.0;
        double routing_outflow       = 0.0;
        double routing_evap_loss     = 0.0;
        double routing_seep_loss     = 0.0;
        double routing_init_storage  = 0.0;
        double routing_final_storage = 0.0;
 
        // User-forced volumes (diagnostic — subset of routing_external)
        double routing_forcing_inflow = 0.0; 
 
        // User-forced quality mass (diagnostic — cumulative mass injected via user_conc_mass_flux)
        std::vector<double> routing_forcing_qual_inflow; 
 
        // Groundwater mass balance totals (all in feet — depth per unit area)
        // Matches legacy TGwaterTotals
        double gw_infil         = 0.0; 
        double gw_upper_evap    = 0.0; 
        double gw_lower_evap    = 0.0; 
        double gw_lower_perc    = 0.0; 
        double gw_lateral_flow  = 0.0; 
        double gw_init_storage  = 0.0; 
        double gw_final_storage = 0.0; 
 
        // Per-step accumulators (reset each step for reporting)
        double step_flooding     = 0.0;
        double step_outflow      = 0.0;
        double step_dw_inflow    = 0.0;
        double step_gw_inflow    = 0.0;
        double step_rdii_inflow  = 0.0;
        double step_ext_inflow   = 0.0;
 
        // Quality mass balance (per-pollutant, in mass units)
        std::vector<double> qual_init_buildup;   
        std::vector<double> qual_final_buildup;  
        std::vector<double> qual_surface_buildup; 
        std::vector<double> qual_wet_deposition; 
        std::vector<double> qual_sweeping;       
        std::vector<double> qual_bmp_removal;    
        std::vector<double> qual_infil_loss;     
        std::vector<double> qual_runoff_load;    
        std::vector<double> qual_routing_wet;    
        std::vector<double> qual_routing_outflow;
        std::vector<double> qual_routing_flood;  
        std::vector<double> qual_routing_init;   
        std::vector<double> qual_routing_final;  
        std::vector<double> qual_routing_reacted;
        std::vector<double> qual_routing_ii_in;  
        std::vector<double> qual_routing_dw_in;  
        std::vector<double> qual_routing_gw_in;  
        std::vector<double> qual_routing_ex_in;  
        std::vector<double> qual_routing_seep;   
        std::vector<double> qual_routing_evap;   
 
        void resize_quality(int n_pollutants) {
            auto np = static_cast<std::size_t>(n_pollutants);
            qual_init_buildup.assign(np, 0.0);
            qual_final_buildup.assign(np, 0.0);
            qual_surface_buildup.assign(np, 0.0);
            qual_wet_deposition.assign(np, 0.0);
            qual_sweeping.assign(np, 0.0);
            qual_bmp_removal.assign(np, 0.0);
            qual_infil_loss.assign(np, 0.0);
            qual_runoff_load.assign(np, 0.0);
            qual_routing_wet.assign(np, 0.0);
            qual_routing_outflow.assign(np, 0.0);
            qual_routing_flood.assign(np, 0.0);
            qual_routing_init.assign(np, 0.0);
            qual_routing_final.assign(np, 0.0);
            qual_routing_reacted.assign(np, 0.0);
            qual_routing_ii_in.assign(np, 0.0);
            qual_routing_dw_in.assign(np, 0.0);
            qual_routing_gw_in.assign(np, 0.0);
            qual_routing_ex_in.assign(np, 0.0);
            qual_routing_seep.assign(np, 0.0);
            qual_routing_evap.assign(np, 0.0);
            routing_forcing_qual_inflow.assign(np, 0.0);
        }
 
        void reset() {
            // Save quality vectors, reset scalars, restore vectors
            auto qi = std::move(qual_init_buildup);
            auto qf = std::move(qual_final_buildup);
            auto qsb = std::move(qual_surface_buildup);
            auto qwd = std::move(qual_wet_deposition);
            auto qsw = std::move(qual_sweeping);
            auto qbmp = std::move(qual_bmp_removal);
            auto qil = std::move(qual_infil_loss);
            auto qrl = std::move(qual_runoff_load);
            auto qrw = std::move(qual_routing_wet);
            auto qro = std::move(qual_routing_outflow);
            auto qrf = std::move(qual_routing_flood);
            auto qri = std::move(qual_routing_init);
            auto qrfi = std::move(qual_routing_final);
            auto qrr = std::move(qual_routing_reacted);
            auto qrii = std::move(qual_routing_ii_in);
            auto qrdw = std::move(qual_routing_dw_in);
            auto qrgw = std::move(qual_routing_gw_in);
            auto qrex = std::move(qual_routing_ex_in);
            auto qrseep = std::move(qual_routing_seep);
            auto qrevap = std::move(qual_routing_evap);
            auto qrfqi = std::move(routing_forcing_qual_inflow);
            *this = MassBalance{};
            qual_init_buildup = std::move(qi);
            qual_final_buildup = std::move(qf);
            qual_surface_buildup = std::move(qsb);
            qual_wet_deposition = std::move(qwd);
            qual_sweeping = std::move(qsw);
            qual_bmp_removal = std::move(qbmp);
            qual_infil_loss = std::move(qil);
            qual_runoff_load = std::move(qrl);
            qual_routing_wet = std::move(qrw);
            qual_routing_outflow = std::move(qro);
            qual_routing_flood = std::move(qrf);
            qual_routing_init = std::move(qri);
            qual_routing_final = std::move(qrfi);
            qual_routing_reacted = std::move(qrr);
            qual_routing_ii_in = std::move(qrii);
            qual_routing_dw_in = std::move(qrdw);
            qual_routing_gw_in = std::move(qrgw);
            qual_routing_ex_in = std::move(qrex);
            qual_routing_seep = std::move(qrseep);
            qual_routing_evap = std::move(qrevap);
            routing_forcing_qual_inflow = std::move(qrfqi);
            // Zero out the quality vectors (except init_buildup which is
            // computed once during initQuality and must survive reset)
            for (auto* v : {&qual_final_buildup,
                            &qual_surface_buildup, &qual_wet_deposition,
                            &qual_sweeping, &qual_bmp_removal,
                            &qual_infil_loss, &qual_runoff_load,
                            &qual_routing_wet, &qual_routing_outflow,
                            &qual_routing_flood, &qual_routing_init,
                            &qual_routing_final, &qual_routing_reacted,
                            &qual_routing_ii_in,
                            &qual_routing_dw_in,
                            &qual_routing_gw_in,
                            &qual_routing_ex_in,
                            &qual_routing_seep,
                            &qual_routing_evap,
                            &routing_forcing_qual_inflow}) {
                std::fill(v->begin(), v->end(), 0.0);
            }
        }

        double runoff_error() const {
            double total_in = runoff_rainfall + runoff_init_store;
            double total_out = runoff_evap + runoff_infil + runoff_runoff + runoff_final_store;
            return (total_in > 0.0) ? (total_in - total_out) / total_in : 0.0;
        }

        double routing_error() const {
            double total_in = routing_dry_weather + routing_wet_weather +
                              routing_gw_inflow + routing_rdii + routing_external +
                              routing_init_storage;
            double total_out = routing_flooding + routing_outflow +
                               routing_evap_loss + routing_seep_loss +
                               routing_final_storage;
            return (total_in > 0.0) ? (total_in - total_out) / total_in : 0.0;
        }

        double gw_error() const {
            double total_in  = gw_infil + gw_init_storage;
            double total_out = gw_upper_evap + gw_lower_evap + gw_lower_perc +
                               gw_lateral_flow + gw_final_storage;
            return (total_in > 0.0) ? (total_in - total_out) / total_in : 0.0;
        }
    } mass_balance;

    struct MassBalance2D {
        double init_storage          = 0.0;  
        double final_storage         = 0.0;  
        double rainfall_in           = 0.0;  
        double coupling_1d_to_2d_in  = 0.0;  
        double coupling_2d_to_1d_out = 0.0;  
        double outfall_in            = 0.0;  
        double outfall_out           = 0.0;  
        double boundary_in           = 0.0;  
        double boundary_out          = 0.0;  
        double evap_out              = 0.0;  
        bool   active                = false;

        double error() const {
            double total_in  = rainfall_in + coupling_1d_to_2d_in + outfall_in
                               + boundary_in + init_storage;
            double total_out = coupling_2d_to_1d_out + outfall_out + boundary_out
                               + evap_out + final_storage;
            return (total_in > 0.0) ? (total_in - total_out) / total_in : 0.0;
        }
    } mass_balance_2d;
 
    // =========================================================================
    // Routing time-step statistics
    // =========================================================================

    struct MaxStats {
        int    obj_type = -1;  
        int    index    = -1;  
        double value    = 0.0; 
    };
 
    static constexpr int MAX_STATS = 5;

    MaxStats max_courant_crit[MAX_STATS] = {};
    MaxStats max_flow_turns[MAX_STATS] = {};
    MaxStats max_non_converged[MAX_STATS] = {};
    MaxStats max_mass_bal_errs[MAX_STATS] = {};
 
    struct RoutingStepStats {
        double min_step  = 1.0e30; 
        double max_step  = 0.0;    
        double sum_step  = 0.0;    
        long   n_steps   = 0;      
        double steady_pct = 0.0;   

        static constexpr int N_TIME_BINS = 5;
        long   step_counts[N_TIME_BINS + 1] = {};   
        double step_intervals[N_TIME_BINS + 1] = {}; 

        long   n_non_converged = 0;   
        double sum_iterations  = 0.0; 
        double max_courant     = 0.0; 
 
        void update(double dt) {
            min_step = std::min(min_step, dt);
            max_step = std::max(max_step, dt);
            sum_step += dt;
            ++n_steps;
        }

        void update_iterations(int iters, bool converged) {
            sum_iterations += iters;
            if (!converged) ++n_non_converged;
        }

        void init_histogram(double route_step, double min_route_step) {
            if (route_step <= 0.0) return;
            if (min_route_step <= 0.0) min_route_step = route_step;
            double log_hi = std::log10(route_step);
            double log_lo = std::log10(min_route_step);
            double delta = (log_hi - log_lo) / static_cast<double>(N_TIME_BINS);
            step_intervals[0] = route_step;
            for (int i = 1; i <= N_TIME_BINS; ++i)
                step_intervals[i] = std::pow(10.0, log_hi - i * delta);
            step_intervals[N_TIME_BINS] = min_route_step;
        }

        void build_histogram() {
            if (n_steps == 0 || max_step <= 0.0) return;
            double hi = max_step;
            double lo = (min_step < 1.0e30) ? min_step : 0.0;
            if (lo <= 0.0) lo = hi;
            init_histogram(hi, lo);
        }

        void record_step_bin(double dt) {
            for (int i = 0; i < N_TIME_BINS; ++i) {
                if (dt >= step_intervals[i+1]) {
                    step_counts[i]++;
                    return;
                }
            }
            step_counts[N_TIME_BINS - 1]++;
        }
 
        double avg_step() const {
            return (n_steps > 0) ? sum_step / static_cast<double>(n_steps) : 0.0;
        }
 
        double pct_non_converged() const {
            return (n_steps > 0) ? 100.0 * static_cast<double>(n_non_converged) / static_cast<double>(n_steps) : 0.0;
        }
 
        double computed_avg_iterations() const {
            return (n_steps > 0) ? sum_iterations / static_cast<double>(n_steps) : 0.0;
        }
    } routing_stats;
 
    // =========================================================================
    // Control action log — Gap #67
    // Populated by ControlEngine::applyPendingActions() when rpt_controls is on.
    // =========================================================================

    struct ControlLogEntry {
        int         link_idx;    
        std::string rule_name;   
        double      new_setting; 
        double      date;        
    };

    std::vector<ControlLogEntry> control_log;
 
    // =========================================================================
    // Input file path (for model write / hot start)
    // =========================================================================

    std::string inp_file_path;
 
    // =========================================================================
    // Context-level operations
    // =========================================================================

    static void updateMaxStats(MaxStats arr[], int obj_type, int idx, double value) {
        MaxStats candidate;
        candidate.obj_type = obj_type;
        candidate.index    = idx;
        candidate.value    = value;
        for (int k = 0; k < MAX_STATS; ++k) {
            if (std::fabs(candidate.value) > std::fabs(arr[k].value)) {
                MaxStats tmp = arr[k];
                arr[k] = candidate;
                candidate = tmp;
            }
        }
    }

    void finalize_max_stats() {
        long step_count = routing_stats.n_steps;
        if (step_count <= 0) return;
        double inv_steps = 1.0 / static_cast<double>(step_count);
 
        // CFL-critical elements: percentage of steps each element was critical
        for (int j = 0; j < n_nodes(); ++j) {
            double x = nodes.stat_time_courant_critical[static_cast<std::size_t>(j)] * inv_steps;
            updateMaxStats(max_courant_crit, 0, j, 100.0 * x);
        }
        for (int j = 0; j < n_links(); ++j) {
            double x = links.stat_time_courant_critical[static_cast<std::size_t>(j)] * inv_steps;
            updateMaxStats(max_courant_crit, 1, j, 100.0 * x);
        }
 
        // Flow instability index (matching legacy normalization)
        long rpt_steps = routing_stats.n_steps;
        if (rpt_steps > 2) {
            double z = 100.0 / (2.0 / 3.0 * static_cast<double>(rpt_steps - 2));
            for (int j = 0; j < n_links(); ++j) {
                double x = static_cast<double>(links.stat_flow_turns[static_cast<std::size_t>(j)]) * z;
                updateMaxStats(max_flow_turns, 1, j, x);
            }
        }
 
        // Non-convergence: fraction of total steps each node failed to converge
        for (int j = 0; j < n_nodes(); ++j) {
            double x = static_cast<double>(nodes.stat_non_converged_count[static_cast<std::size_t>(j)]) * inv_steps;
            updateMaxStats(max_non_converged, 0, j, x);
        }
    }

    void reset() {
        state = EngineState::CREATED;
        control_log.clear();
        current_time = 0.0;
        current_date = 0.0;
        dt_controls_remaining = 0.0;
        elapsed_ms = 0.0;
        old_elapsed_ms = 0.0;
        next_report_ms = 0.0;
        error_code = 0;
        warning_code = 0;
        error_message.clear();
        warnings.clear();
        errors.clear();
        title_notes.clear();
 
        // Clear SoA stores
        nodes      = NodeData{};
        links      = LinkData{};
        subcatches = SubcatchData{};
        gages      = GageData{};
        pollutants = PollutantData{};
        tables     = TableData{};
 
        // Clear name indices
        node_names.clear();
        link_names.clear();
        subcatch_names.clear();
        gage_names.clear();
        pollutant_names.clear();
        table_names.clear();
 
        // Clear inflow-related stores that aren't reset by their owning solvers
        rdii_decay = RDIIDecayData{};
 
        // Clear daily climate state (re-initialized by SWMMEngine on next run)
        climate_state = climate::ClimateState{};
 
        // Clear spatial, flags, events, and forcing. Per-object tags
        // are owned by NodeData/LinkData/SubcatchData and cleared when
        // those SoAs are resized/cleared by the wider reset path.
        spatial    = SpatialFrame{};
        user_flags.clear();
        events.clear();
        std::fill(std::begin(adjust_temp), std::end(adjust_temp), 0.0);
        std::fill(std::begin(adjust_evap), std::end(adjust_evap), 1.0);
        std::fill(std::begin(adjust_rain), std::end(adjust_rain), 1.0);
        std::fill(std::begin(adjust_hydcon), std::end(adjust_hydcon), 1.0);
        subcatch_n_perv_pattern.clear();
        subcatch_d_store_pattern.clear();
        subcatch_infil_pattern.clear();
        base_n_perv.clear();
        base_ds_perv.clear();
        has_subcatch_adj_patterns = false;
        forcing    = ForcingData{};
    }

    void save_state() noexcept {
        nodes.save_state();
        links.save_state();
        // NOTE: subcatches.save_state() is intentionally NOT called here.
        // Subcatchment old-state (old_runoff/old_runon/conc_old) is the
        // runoff-step snapshot used to linearly interpolate lateral inflow
        // between runoff evaluations (legacy subcatch_setOldState, called
        // ONLY inside runoff_execute — i.e. per WET/DRY step, not per routing
        // step). Saving it here, every routing step, clobbered old_runoff with
        // the current value so old==new and the lateral inflow jumped to the
        // full new runoff instantly instead of ramping. It is now saved in the
        // runoff-advance loop in SWMMEngine::stepRunoff().
    }

    void reset_state() noexcept {
        nodes.reset_state();
        links.reset_state();
        subcatches.reset_state();
        gages.reset_state();
        tables.reset_cursors();
        current_time = 0.0;
        current_date = options.start_date;
        dt_controls_remaining = 0.0;
        elapsed_ms = 0.0;
        old_elapsed_ms = 0.0;
        // Legacy swmm5.c:721 — ReportTime = 1000 * ReportStep (ms from SIM start)
        next_report_ms = 1000.0 * options.report_step;
    }

    void allocate_objects() {
        nodes.resize(node_names.size());
        links.resize(link_names.size());
        subcatches.resize(subcatch_names.size());
        gages.resize(gage_names.size());
        pollutants.resize_pollutants(pollutant_names.size());
        int np = static_cast<int>(pollutant_names.size());
        nodes.resize_quality(np);
        links.resize_quality(np);
        subcatches.resize_quality(np);
 
        // Resize spatial coordinate arrays
        spatial.node_x.assign(static_cast<std::size_t>(node_names.size()), 0.0);
        spatial.node_y.assign(static_cast<std::size_t>(node_names.size()), 0.0);
        spatial.link_x.assign(static_cast<std::size_t>(link_names.size()), 0.0);
        spatial.link_y.assign(static_cast<std::size_t>(link_names.size()), 0.0);
        spatial.subcatch_x.assign(static_cast<std::size_t>(subcatch_names.size()), 0.0);
        spatial.subcatch_y.assign(static_cast<std::size_t>(subcatch_names.size()), 0.0);
 
        // Resize link vertex and subcatchment polygon arrays
        spatial.link_vertices_x.resize(static_cast<std::size_t>(link_names.size()));
        spatial.link_vertices_y.resize(static_cast<std::size_t>(link_names.size()));
        spatial.subcatch_polygon_x.resize(static_cast<std::size_t>(subcatch_names.size()));
        spatial.subcatch_polygon_y.resize(static_cast<std::size_t>(subcatch_names.size()));
 
        // Resize gage coordinates
        spatial.gage_x.assign(static_cast<std::size_t>(gage_names.size()), 0.0);
        spatial.gage_y.assign(static_cast<std::size_t>(gage_names.size()), 0.0);
    }

    void shrink_all_to_fit() {
        nodes.shrink_to_fit();
        links.shrink_to_fit();
        subcatches.shrink_to_fit();
        gages.shrink_to_fit();
        pollutants.shrink_to_fit();
        landuses.shrink_to_fit();
        buildup.shrink_to_fit();
        washoff.shrink_to_fit();
        treatment.shrink_to_fit();
        spatial.shrink_to_fit();
    }
 
    // =========================================================================
    // Convenience accessors
    // =========================================================================

    int n_nodes()      const noexcept { return node_names.size(); }

    int n_links()      const noexcept { return link_names.size(); }

    int n_subcatches() const noexcept { return subcatch_names.size(); }

    int n_gages()      const noexcept { return gage_names.size(); }

    int n_pollutants() const noexcept { return pollutant_names.size(); }
    int n_landuses()   const noexcept { return landuse_names.size(); }

    int n_tables()     const noexcept { return table_names.size(); }
};
 
} /* namespace openswmm */
 
#endif /* OPENSWMM_ENGINE_SIMULATION_CONTEXT_HPP */