SimulationContext.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_ENGINE_SIMULATION_CONTEXT_HPP
#define OPENSWMM_ENGINE_SIMULATION_CONTEXT_HPP
 
#include <cmath>
#include <ctime>
#include <functional>
#include "FilePathPair.hpp"
#include "../data/ArdConfigData.hpp"
#include "../data/BedZoneData.hpp"
#include "../data/GageData.hpp"
#include "../data/HeatData.hpp"
#include "../data/WaterAgeData.hpp"
#include "../data/LidLayerSpeciesData.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/ReactionData.hpp"
#include "../data/SpeciesRegistry.hpp"
#include "../data/SubcatchData.hpp"
#include "../data/TableData.hpp"
#include "SimulationOptions.hpp"
#include "SpatialFrame.hpp"
#include "UserFlags.hpp"
#include "../data/QualityData.hpp"
#include "../data/InitialQualityData.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 <utility>
#include <vector>
 
namespace openswmm {
 
// ============================================================================
// Plugin specification (from [PLUGINS] section)
// ============================================================================

struct PluginSpec {
    std::string              path;       
    std::vector<std::string> init_args;  
};
 
// ============================================================================
// [PROCESS_COMPONENTS] section spec — process-component registrations
// ============================================================================

struct ProcessComponentSpec {
    std::string id;           
    std::string config_path;  
    std::vector<std::pair<std::string, std::string>> args;  

    std::string resolved_config_path;
};
 
// ============================================================================
// [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;
class  Infil2D;
struct PendingBoundaryRow;
struct PendingEdgeConveyanceRow;
struct PendingInitialQualityRow;
struct PendingBoundaryQualityRow;
struct GwTransportData;   // U4 (2026-09-07)
struct SubsurfaceConfig;  // G1: the [2D_AQUIFER*] rows
struct SubsurfaceState;   // G1: the running two-zone kernel state
} // namespace twoD
 
// ============================================================================
// SimulationContext
// ============================================================================

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

    EngineState state = EngineState::CREATED;

    std::time_t wall_start = 0;
 
    // =========================================================================
    // 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;

    SpeciesRegistry species_registry;

    ReactionData reactions;

    ArdConfigData ard_config;

    WaterAgeConfigData water_age_config;
    WaterAgeState      water_age_state;

    LidLayerSpeciesState lid_layer_state;

    HeatConfigData heat_config;
    HeatState      heat_state;

    BedZoneState bed_state;

    std::vector<std::string> reported_species_names;

    int n_reported_species() const noexcept {
        return static_cast<int>(reported_species_names.size());
    }

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

    NameIndex node_names;

    NameIndex link_names;

    NameIndex subcatch_names;

    NameIndex gage_names;

    NameIndex pollutant_names;
 
    // NOTE: there is deliberately NO shared name registry for tables.
    // Legacy keeps TSERIES and CURVE in separate hash tables, so a curve and
    // a timeseries may legally share one name; ctx.tables is the authority
    // (each row stores id + type) and lookups are kind-aware via
    // find_timeseries() / find_curve() below.
 
    // =========================================================================
    // Quality data (landuse, buildup, washoff, treatment)
    // =========================================================================
 
    NameIndex        landuse_names;
    LanduseData      landuses;
    BuildupData      buildup;
    WashoffData      washoff;
    TreatmentData    treatment;

    InitialQualityData initial_quality;
 
    // =========================================================================
    // 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;
 
    // =========================================================================
    // Transient deferred-resolution capture lists (input loading only)
    // =========================================================================
    // Sections whose referenced objects may be defined later in the .inp store
    // the raw name here so PostParseResolver can re-resolve after every section
    // has been parsed. Cleared by the resolver; empty outside input loading, so
    // they can never desync with the editing APIs.

    std::vector<std::pair<int, std::string>> pending_gw_nodes;

    std::vector<std::pair<int, std::pair<std::string, std::string>>>
        pending_link_nodes;
 
    // =========================================================================
    // 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;

    std::vector<ProcessComponentSpec> process_component_specs;

    std::vector<std::pair<std::string, std::vector<std::string>>>
        embedded_component_sections;

    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;
        std::vector<twoD::PendingInitialQualityRow>*  pending_iq = nullptr;
        std::vector<twoD::PendingBoundaryQualityRow>* pending_bq = nullptr;
        twoD::Infil2D*                               infil      = nullptr;
        twoD::GwTransportData*                       gw         = nullptr;
        twoD::SubsurfaceConfig*                      aquifer    = nullptr;
        std::vector<std::string>*                    aquifer_nodes = nullptr;
        twoD::SubsurfaceState*                       aquifer_state = nullptr;
    } twod_io;

    std::vector<std::pair<std::string, std::string>> deferred_section_rows;
 
    // =========================================================================
    // 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;

    double coupling_delivery_remaining = 0.0;
 
    // =========================================================================
    // Mass balance accumulators (SoA — vectorisable batch updates)
    // =========================================================================

    struct MassBalance {
        // Runoff totals
        double runoff_rainfall   = 0.0;  
        double runoff_runon      = 0.0;  
        double runoff_evap       = 0.0;  
        double runoff_infil      = 0.0;  
        double runoff_runoff     = 0.0;  
        double runoff_lid_drain  = 0.0;  
        double runoff_snowremov  = 0.0;  
        double runoff_init_store = 0.0;  
        double runoff_final_store= 0.0;  
        double runoff_init_snow  = 0.0;
        double runoff_final_snow = 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_coupling_out  = 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; 
        double gw_infil_2d_recharge = 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 {
            // Snow is a STORE on both sides, and snow ploughed out of the
            // system is a loss — matching legacy massbal.c:685-692. A pack
            // present at the start is water the model was given; a pack still
            // standing at the end is water it still holds.
            double total_in = runoff_rainfall + runoff_runon + runoff_init_store +
                              runoff_init_snow;
            double total_out = runoff_evap + runoff_infil + runoff_runoff +
                               runoff_lid_drain +
                               runoff_final_store + runoff_snowremov +
                               runoff_final_snow;
            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_coupling_out +
                               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 NegativeSourceStats {
        long   clamp_events     = 0;   
        double shortfall_mass   = 0.0; 
        long   age_clamp_events = 0;   
        int    first_node       = -1;  
        bool   first_clamp_recorded = false;
        bool   api_warned       = false; 
        void reset() { *this = NegativeSourceStats{}; }
    } negsrc;

    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;  
        double infil_out             = 0.0;
        double infil_to_aquifer      = 0.0;
        double aquifer_in            = 0.0;
        bool   active                = false;
 
        // Cumulative marcher statistics (published by SurfaceRouter2D at
        // finalize from ISurfaceSolver::run_stats). Printed as the "2D Solver
        // Statistics" report block. -1 = not populated.
        long   solver_nsteps         = -1;   
        long   solver_nrhs           = 0;    
        double solver_avg_h          = 0.0;  
        double solver_last_h         = 0.0;  
 
        // Marcher telemetry (guarded by n_tiers > 0 / non-negative fractions).
        double solver_active_min     = -1.0; 
        double solver_active_mean    = -1.0; 
        double solver_active_max     = -1.0; 
        long   solver_tier_cells[8]  = {0};  
        int    solver_n_tiers        = 0;    

        double error() const {
            double total_in  = rainfall_in + coupling_1d_to_2d_in + outfall_in
                               + boundary_in + aquifer_in + init_storage;
            double total_out = coupling_2d_to_1d_out + outfall_out + boundary_out
                               + evap_out + infil_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;
        }
 
        // ---------------------------------------------------------------
        // FV 1D solver statistics (published by SWMMEngine::end() from
        // INetworkSolver::run_stats, printed as the "FV Solver Statistics"
        // report block). The 1D counterpart of mass_balance_2d.solver_*.
        //
        // -1 = not populated: the run was not FLOW_ROUTING FV, or a backend
        // that carries no counters was selected. The report block is skipped
        // on the sentinel rather than printing a row of zeros, which would
        // read as "the solver did nothing" instead of "nobody counted".
        // ---------------------------------------------------------------
        // Slot-storage share (FV slot program R0). Peak instantaneous
        // system share slot/stored and the time that share exceeded 1 %;
        // the run-level integrated share is Σ links.stat_slot_vol_dt /
        // Σ links.stat_vol_dt, summed at report time. 0 under DW.
        double slot_peak_share   = 0.0;
        double slot_time_above_s = 0.0;
 
        long   fv_nsteps        = -1;   
        long   fv_nflux         = 0;    
        double fv_avg_h         = 0.0;  
        double fv_last_h        = 0.0;  
        double fv_min_h         = 0.0;  
        double fv_active_min    = -1.0; 
        double fv_active_mean   = -1.0; 
        double fv_active_max    = -1.0; 
        long   fv_tier_cells[8] = {0};  
        int    fv_n_tiers       = 0;    
 
        // dt-argmin attribution (slot program R0): who owned the binding
        // CFL element, counted per census / re-tier.
        long   fv_dt_argmin_pressurized = 0;
        long   fv_dt_argmin_band        = 0;
        long   fv_dt_argmin_free        = 0;
        long   fv_dt_argmin_node        = 0;
    } routing_stats;
 
    // =========================================================================
    // Virtual-junction diagnostics (refactored engine only)
    // =========================================================================

    struct VJDiag {
        std::vector<int>       node_idx;   
        std::vector<int>       up_link;    
        std::vector<int>       dn_link;    
        std::vector<double>    resid_max;  
        std::vector<double>    resid_sum;  
        std::vector<long long> resid_n;    
        void clear() {
            node_idx.clear(); up_link.clear(); dn_link.clear();
            resid_max.clear(); resid_sum.clear(); resid_n.clear();
        }
    } vj_diag;
 
    // =========================================================================
    // Street inlet performance diagnostics
    // =========================================================================

    struct InletDiag {
        std::vector<int>     host_node;         
        std::vector<int>     up_link;           
        std::vector<uint8_t> is_sag;            
        std::vector<int>     num_inlets;
        std::vector<int>     flow_periods;      
        std::vector<int>     capture_periods;   
        std::vector<int>     backflow_periods;  
        std::vector<double>  peak_flow;         
        std::vector<double>  peak_flow_capture; 
        std::vector<double>  avg_flow_capture;  
        std::vector<double>  bypass_freq;       
 
        int count() const { return static_cast<int>(host_node.size()); }
 
        void resize(int n) {
            auto un = static_cast<std::size_t>(n);
            host_node.assign(un, -1);
            up_link.assign(un, -1);
            is_sag.assign(un, static_cast<uint8_t>(0));
            num_inlets.assign(un, 1);
            flow_periods.assign(un, 0);
            capture_periods.assign(un, 0);
            backflow_periods.assign(un, 0);
            peak_flow.assign(un, 0.0);
            peak_flow_capture.assign(un, 0.0);
            avg_flow_capture.assign(un, 0.0);
            bypass_freq.assign(un, 0.0);
        }
 
        void clear() {
            host_node.clear(); up_link.clear(); is_sag.clear(); num_inlets.clear();
            flow_periods.clear(); capture_periods.clear(); backflow_periods.clear();
            peak_flow.clear(); peak_flow_capture.clear(); avg_flow_capture.clear();
            bypass_freq.clear();
        }
    } inlet_diag;
 
    // =========================================================================
    // 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();
        deferred_section_rows.clear();
        pending_gw_nodes.clear();
        pending_link_nodes.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();
 
        // Clear inflow-related stores that aren't reset by their owning solvers
        rdii_decay = RDIIDecayData{};
 
        // E3: transport.ard component config — its apply hook resets it, but
        // a reopen WITHOUT the component would otherwise inherit the previous
        // model's dispersion.
        ard_config = ArdConfigData{};
 
        // A1a: same stale-on-reopen hygiene for water age.
        water_age_config = WaterAgeConfigData{};
        water_age_state.clear();
 
        // H1: and for heat.
        heat_config = HeatConfigData{};
        heat_state.clear();
        lid_layer_state.clear();
        bed_state.clear();      // H6b
 
        // Virtual-junction diagnostics
        vj_diag.clear();
        inlet_diag.clear();
 
        // 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 static_cast<int>(tables.count()); }
 
    // =========================================================================
    // Kind-aware table lookups (case-insensitive, legacy hash.c parity).
    // Legacy keeps TSERIES and CURVE in separate hash tables, so the same
    // name may denote both a curve and a timeseries; every consumer knows
    // which kind it wants.
    // =========================================================================

    int find_timeseries(std::string_view name) const noexcept {
        return tables.find_by_kind(name, /*want_timeseries=*/true);
    }

    int find_curve(std::string_view name) const noexcept {
        return tables.find_by_kind(name, /*want_timeseries=*/false);
    }

    int find_table_any(std::string_view name) const noexcept {
        const int ts = find_timeseries(name);
        if (ts >= 0) return ts;
        return find_curve(name);
    }
};
 
} /* namespace openswmm */
 
#endif /* OPENSWMM_ENGINE_SIMULATION_CONTEXT_HPP */