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/ChannelThermalData.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 "../data/NodeCoupling2DData.hpp"
#include "../hydraulics/Transect.hpp"
#include "../data/HydrologyData.hpp"
#include "../data/ForcingData.hpp"
#include "../hydrology/Climate.hpp"
 
#include <algorithm>
#include <array>
#include <cstdint>
#include <map>
#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 };

enum class HotstartFormat { AUTO = 0, LEGACY = 1, NATIVE = 2 };

struct HotstartSaveEntry {
    FilePathPair path;
    double datetime = 0.0;
    HotstartFormat format = HotstartFormat::AUTO;
};
 
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()> has_fv_retained_dry_mass;
    std::function<std::string()> get_fv_dry_state;
    std::function<std::string(const std::string&,bool)> set_fv_dry_state;
    std::function<std::vector<double>()> get_fv_retained_pollutant_mass;
    std::function<std::vector<double>()> get_fv_retained_species_mass;
    std::function<std::vector<ChannelThermalRegion>()> get_channel_thermal_regions;
    // Each region appears once. Validate the complete batch before mutation.
    std::function<std::string(const std::vector<ChannelThermalIncrement>&)>
        apply_channel_thermal_increments;
    std::function<std::vector<ChannelSpeciesRegion>()> get_channel_species_regions;
    std::function<std::string(const std::vector<ChannelSpeciesIncrement>&)>
        apply_channel_species_increments;
    std::function<std::string()> get_snow_state;
    // Empty result means valid; false validates without publishing.
    std::function<std::string(const std::string&, bool)> set_snow_state;
    std::function<std::string()> get_surface_state;
    std::function<std::string(const std::string&, bool)> set_surface_state;
    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;
class SurfaceRouter2D;
struct PendingBoundaryRow;
struct PendingEdgeConveyanceRow;
struct PendingInitialQualityRow;
struct PendingInitialBedQualityRow;
struct PendingBoundaryQualityRow;
struct GwTransportData;
class SurfaceQuality2D;          // S7 (2026-09-19)
struct SubsurfaceConfig;         // G1: the [2D_AQUIFER*] rows
struct SubsurfaceState;          // G1: the running two-zone kernel state
struct SubsurfaceTransportState; // T7.5: and what is dissolved in it
} // 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;

    std::map<std::string, std::vector<std::string>> passthrough_sections;
 
    // =========================================================================
    // 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 runoff_interp_f = 1.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;

    void refreshReportedSpecies() {
        reported_species_names.clear();
        for (int p = 0; p < n_pollutants(); ++p)
            reported_species_names.push_back(pollutant_names.name_of(p));
        if (options.water_age)
            reported_species_names.push_back("__WATER_AGE__");
        if (options.heat_transport)
            reported_species_names.push_back("__TEMPERATURE__");
        reactions.reported_species.clear();
        for (int m = 0; m < reactions.n_species(); ++m) {
            if (!reactions.report_enabled.empty() && !reactions.report_enabled[m])
                continue;
            reactions.reported_species.push_back(m);
            reported_species_names.push_back(reactions.species_name[m]);
        }
    }

    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] = {
        0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0}; 
    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;
    std::vector<LidNodeOutlet> lid_node_outlets;
    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::string>> pending_gw_aquifers;

    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::PendingInitialBedQualityRow>* pending_ibq = nullptr;
        std::vector<twoD::PendingBoundaryQualityRow>* pending_bq = nullptr;
        twoD::Infil2D* infil = nullptr;
        twoD::GwTransportData* gw = nullptr;
        twoD::SurfaceQuality2D* surface_quality = nullptr;
        twoD::SubsurfaceConfig* aquifer = nullptr;
        std::vector<std::string>* aquifer_nodes = nullptr;
        std::vector<std::string>* aquifer_links = nullptr;
        twoD::SubsurfaceState* aquifer_state = nullptr;
        twoD::SurfaceRouter2D* surface_router = nullptr;
        twoD::SubsurfaceTransportState* aquifer_transport = nullptr;
    } twod_io;

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

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

    NodeCoupling2DRows node_coupling_2d;

    std::vector<std::uint8_t> coupled_node;

    std::vector<double> coupled_pond_area;
    std::vector<double> coupled_pond_band;

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

    // One row per node, booked at the same committed seam as routing totals.
    // Read-only external coupling accounting; never a second withdrawal.
    struct OutfallExports {
        std::vector<double> water_m3, reverse_m3, heat_j, pollutant_amounts;
    } outfall_exports;
 
    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;
 
        // RDII generation totals (ft3), distinct from routed RDII and runoff
        // rainfall: sewer service areas need not equal subcatchment areas.
        double rdii_rainfall = 0.0;
        double rdii_produced = 0.0;
 
        // Routing totals
        double routing_dry_weather = 0.0;
        double routing_wet_weather = 0.0;
        double routing_gw_inflow = 0.0;
        double routing_link_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;

        std::array<double, 11> routing_report{};
 
        // 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;
        double last_step_flow_error = 0.0;
        double step_ext_withdrawal = 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_ponded_lost; 
        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;    
        std::vector<double> qual_routing_final_dry;

        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_ponded_lost.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);
            qual_routing_final_dry.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 qpl = std::move(qual_ponded_lost);
            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 qrfd = std::move(qual_routing_final_dry);
            auto qrfqi = std::move(routing_forcing_qual_inflow);
            *this = MassBalance{};
            qual_init_buildup = std::move(qi);
            qual_final_buildup = std::move(qf);
            qual_ponded_lost = std::move(qpl);
            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);
            qual_routing_final_dry = std::move(qrfd);
            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_ponded_lost,
                            &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,
                            &qual_routing_final_dry, &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_link_gw_inflow + // G-X4
                              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;

    // Per-incident-link surface donation limit (ft³/s), indexed by node.
    // Empty for ordinary models; infinity for nodes without a surface outfall.
    std::vector<double> surface_outfall_link_limit;

    struct QualityOnlyInflow {
        std::vector<double> vol, mass, msx, age_vol, temp_vol;
        bool any = false;
        void ensure(int nn, int np, int ns) {
            const auto un = static_cast<std::size_t>(nn > 0 ? nn : 0);
            if (vol.size() != un) {
                vol.assign(un, 0.0);
                age_vol.assign(un, 0.0);
                temp_vol.assign(un, 0.0);
            }
            if (mass.size() != un * static_cast<std::size_t>(np > 0 ? np : 0))
                mass.assign(un * static_cast<std::size_t>(np > 0 ? np : 0), 0.0);
            if (msx.size() != un * static_cast<std::size_t>(ns > 0 ? ns : 0))
                msx.assign(un * static_cast<std::size_t>(ns > 0 ? ns : 0), 0.0);
        }

        void clear() {
            for (auto* v : {&vol, &mass, &msx, &age_vol, &temp_vol})
                std::fill(v->begin(), v->end(), 0.0);
            any = false;
        }
    };
    QualityOnlyInflow surface_outfall_quality; 
    QualityOnlyInflow link_gain_quality;       
 
    bool runtime_coupling_used = false; 
    struct MassBalance2D {
        double init_storage = 0.0;                    
        double final_storage = 0.0;                   
        double external_in = 0.0, external_out = 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 infil_aquifer_pending = 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 + external_in;
            double total_out = coupling_2d_to_1d_out + outfall_out + boundary_out + evap_out +
                               infil_out + final_storage + external_out;
            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;  
        double steady_time = 0.0; 
        long report_steps = 0;    
        double report_time = 0.0; 
        double max_outfall_flow =
            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 record_step(double dt, bool steady, bool first_step) {
            if (steady) {
                steady_time += dt;
                return;
            }
            if (!first_step) {
                min_step = std::min(min_step, dt);
                for (int i = 0; i < N_TIME_BINS; ++i)
                    if (dt >= step_intervals[i + 1]) {
                        step_counts[i]++;
                        break;
                    }
            }
            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;
 
        // LTS macro cycles that ran / were rejected. Distinguishes "tiering
        // never engaged" from "tiering did not help" (the tier histogram
        // alone cannot).
        long fv_macro_cycles = 0;
        long fv_macro_rejected = 0;

        int lstep_mode = 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;       
        int rule_idx;       
        double new_setting; 
        double date;        
    };

    std::vector<ControlLogEntry> control_log;
    std::vector<std::string> control_rule_names;
    std::size_t control_log_dropped = 0;
    static constexpr std::size_t kMaxControlLog = 1'000'000;

    void logControlAction(const ControlLogEntry& e) {
        if (control_log.size() >= kMaxControlLog) {
            ++control_log_dropped;
            return;
        }
        control_log.push_back(e);
    }
 
    // =========================================================================
    // 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() {
        // legacy stats_findMaxStats (stats.c): the slots start at -1.0 (time-
        // step critical, flow turns) or 0.0 (non-convergence), so a value
        // must exceed that magnitude to enter.
        for (int k = 0; k < MAX_STATS; ++k) {
            max_courant_crit[k] = MaxStats{0, -1, -1.0};
            max_flow_turns[k] = MaxStats{1, -1, -1.0};
            max_non_converged[k] = MaxStats{0, -1, 0.0};
        }
 
        // Flow instability over the steps after the report start.
        const long rpt_steps = routing_stats.report_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);
            }
        }
 
        // The rest are over all routed (non-steady) steps.
        const long step_count = routing_stats.n_steps;
        const bool dw = options.routing_model == RoutingModel::DYNWAVE;
        if (dw) {
            for (int j = 0; j < n_nodes(); ++j)
                updateMaxStats(max_non_converged, 0, j,
                               static_cast<double>(
                                   nodes.stat_non_converged_count[static_cast<std::size_t>(j)]) /
                                   static_cast<double>(step_count));
        }
        if (!dw || options.variable_step == 0.0 || step_count == 0)
            return;
        const double inv_steps = 1.0 / static_cast<double>(step_count);
        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);
        }
    }

    void reset() {
        outfall_exports = OutfallExports{};
        runtime_coupling_used = false;
        surface_outfall_link_limit.clear();
        surface_outfall_quality = {};
        link_gain_quality = {};
        state = EngineState::CREATED;
        control_log.clear();
        control_rule_names.clear();
        control_log_dropped = 0;
        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();
        unknown_section_headers.clear();
        pending_gw_nodes.clear();
        pending_gw_aquifers.clear();
        pending_link_nodes.clear();
 
        // Clear SoA stores
        nodes = NodeData{};
        links = LinkData{};
        // Reopening a rain-rule model must not reuse the previous subtype
        // join rows or append its rules to the freshly parsed rule store.
        node_subtypes.clear();
        link_subtypes.clear();
        node_coupling_2d.clear();
        control_rules = ControlRuleStore{};
        subcatches = SubcatchData{};
        gages = GageData{};
        pollutants = PollutantData{};
        tables = TableData{};
        snowpacks = SnowpackStore{};
 
        // Clear name indices
        node_names.clear();
        link_names.clear();
        subcatch_names.clear();
        gage_names.clear();
        pollutant_names.clear();
        snowpack_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), 0.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();
        // (SWMMEngine::step rolls the two halves separately — see
        // save_lat_qual_state / save_hyd_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 save_lat_qual_state() noexcept {
        nodes.save_lat_qual_state();
        links.save_qual_state();
    }

    void save_hyd_state() noexcept {
        nodes.save_hyd_state();
        links.save_hyd_state();
    }

    void reset_state() noexcept {
        runtime_coupling_used = false;
        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.node_has_xy.assign(static_cast<std::size_t>(node_names.size()), 0); // G-X2
        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 */