ReactionData.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_DATA_REACTION_DATA_HPP
#define OPENSWMM_ENGINE_DATA_REACTION_DATA_HPP
 
#include <cstdint>
#include <string>
#include <string_view>
#include <utility>
#include <vector>
 
#include "ReactionTokens.hpp"
#include "MsxSurfaceData.hpp"
 
namespace openswmm {
 
enum class ReactionSolverKind : int { EUL = 0, RK5 = 1, ROS2 = 2, BDF2 = 3 };
enum class ReactionCoupling : int { NONE = 0, FULL = 1 };
enum class ReactionRateUnits : int { SEC = 0, MIN = 1, HR = 2, DAY = 3 };
enum class ReactionAreaUnits : int { FT2 = 0, M2 = 1, CM2 = 2 };
enum class ReactionExprForm : int { NONE = 0, RATE = 1, EQUIL = 2, FORMULA = 3 };
 
struct ReactionData {
    // ---- [REACTION_OPTIONS] ------------------------------------------------
    ReactionSolverKind solver = ReactionSolverKind::RK5;
    ReactionCoupling coupling = ReactionCoupling::NONE;
    ReactionRateUnits rate_units = ReactionRateUnits::HR;
    ReactionAreaUnits area_units = ReactionAreaUnits::FT2;
    bool wall_geometry = false; 
    double timestep = 0.0;      
    double atol = 1.0e-6;
    double rtol = 1.0e-4;
    double default_temp_c = 20.0;
 
    // ---- [REACTION_SPECIES] — index-aligned with the SpeciesRegistry MSX
    //      block (registry index = registry_base + i). ------------------------
    int registry_base = -1;                 
    std::vector<std::string> species_name;  
    std::vector<uint8_t> species_is_wall;   
    std::vector<std::string> species_units; 
    std::vector<double> species_atol;       
    std::vector<double> species_rtol;       
 
    // ---- [REACTION_COEFFICIENTS] -------------------------------------------
    std::vector<std::string> coef_name; 
    std::vector<uint8_t> coef_is_param; 
    std::vector<double> coef_value;     

    struct ParameterOverride {
        bool link;
        int element, coefficient;
        double value;
    };
    std::vector<ParameterOverride> parameter_overrides;
    struct WallArea {
        bool link;
        int element;
        double area;
    }; // declared AREA_UNITS
    std::vector<WallArea> wall_areas;
    struct LidLaw {
        int control = -1, layer = 0, species = -1;
        ReactionExprForm form = ReactionExprForm::NONE;
        std::string expression;
        RxExprSpan compiled;
    };
    std::vector<int> lid_controls;
    std::vector<LidLaw> lid_laws;
    struct LidState {
        std::vector<std::vector<double>> node_mass, node_outlet;
        std::vector<double> node_mobile_old;
        std::vector<int> group_offsets;
        std::vector<double> mass, old_water, surface_conc,
            drain_conc;                                    // conventional units, conc*ft3
        std::vector<double> surface_export, surface_carry; // [subcatchment*species], conc*ft3
        std::vector<uint8_t> subcatchments;
    } lid;

    bool lid_enabled(int control) const {
        for (int c : lid_controls)
            if (c == control)
                return true;
        return false;
    }

    const double* coefficients(bool link, int element, std::vector<double>& scratch) const {
        if (parameter_overrides.empty() || element < 0)
            return coef_value.data();
        scratch = coef_value;
        for (const auto& row : parameter_overrides)
            if (row.link == link && row.element == element)
                scratch[static_cast<std::size_t>(row.coefficient)] = row.value;
        return scratch.data();
    }
 
    struct SourcePattern {
        std::string name;
        std::vector<double> factors;
    };
    struct Source {
        bool mass;
        int node, species, pattern;
        double value;
        int booster = 0;
    }; // 1 FLOWPACED, 2 SETPOINT
    std::vector<double> booster_mass; // actual cumulative injection, concentration*ft3
    double pattern_timestep = 3600.0; 
    std::vector<SourcePattern> source_patterns;
    std::vector<Source> sources;

    std::vector<uint8_t> report_enabled;
    std::vector<int> reported_species;
 
    // ---- [REACTION_TERMS] --------------------------------------------------
    std::vector<std::string> term_name;     
    std::vector<std::string> term_expr_src; 
 
    struct SubsurfaceLaw {
        int scope = 0, cell = -1, zone = 0, layer = -1,
            species = -1; // scope: global/tag/cell; zone: sat/unsat/layer
        std::string tag, expression;
        ReactionExprForm form = ReactionExprForm::NONE;
        RxExprSpan compiled;
    };
    std::vector<SubsurfaceLaw> subsurface_laws;
 
    // ---- [REACTION_PIPES] / [REACTION_TANKS] — dense per species -----------
    std::vector<ReactionExprForm> pipe_form; 
    std::vector<std::string> pipe_expr_src;  
    std::vector<ReactionExprForm> tank_form;
    std::vector<std::string> tank_expr_src;
 
    struct SurfaceLaw {
        int subcatchment = -1, landuse = -1, species = -1;
        ReactionExprForm form = ReactionExprForm::NONE;
        std::string expression;
        RxExprSpan compiled;
    };
    std::vector<SurfaceLaw> wet_laws, buildup_laws;
    bool wet_enabled(int sc) const noexcept {
        for (auto f : subcatch_form)
            if (f != ReactionExprForm::NONE)
                return true;
        for (const auto& l : wet_laws)
            if (l.subcatchment == sc)
                return true;
        return false;
    }
    // Wet, non-LID subcatchment pools have their own reactor context.
    std::vector<ReactionExprForm> subcatch_form;
    std::vector<std::string> subcatch_expr_src;
    std::vector<RxExprSpan> subcatch_expr;
    bool has_subcatch_reactions() const noexcept {
        for (const auto form : subcatch_form)
            if (form != ReactionExprForm::NONE)
                return true;
        return !wet_laws.empty();
    }
 
    // ---- [REACTION_QUALITY] GLOBAL initial values (R1) ---------------------
    std::vector<double> init_global; 
    std::vector<double> rain_conc;
 
    // ---- [REACTION_QUALITY] NODE/LINK per-element overrides (E-B) ----------
    // Sparse row table over the GLOBAL seed; covered by clear() via the
    // *this = ReactionData{} reset.
    std::vector<uint8_t> init_elem_is_link; 
    std::vector<int> init_elem_idx;         
    std::vector<int> init_elem_species;     
    std::vector<double> init_elem_value;
 
    // ---- Compiled bytecode (R2 — D-L3 flat pool; spans index token_pool) ----
    std::vector<RxToken> token_pool;   
    std::vector<RxExprSpan> term_expr; 
    std::vector<RxExprSpan> pipe_expr; 
    std::vector<RxExprSpan> tank_expr; 
    bool compiled = false;             
    bool warned_react_failure = false;
    bool warned_negative_clamp = false; 
 
    // ---- R4: MSX species element state under QUALITY_SOLVER LEGACY --------
    // [element * n_species + s]; sized lazily by the legacy binding. R4b
    // (2026-09-01): species react per element AND advect between elements
    // (routeLegacyMsx, the CSTR mirror family) — the not-transported
    // warning and its once-per-run flag are gone with the limitation.
    std::vector<double> msx_node_conc;
    std::vector<double> msx_link_conc;
 
    // ---- U2 (2026-09-07): external MSX loads at the node seam ------------
    // [node * n_species + m], a mass RATE in the species' internal unit
    // convention (conc-units × ft3/s — the qual_mass_in shape), assembled
    // every routing step by InflowSolver::computeAll from [INFLOWS] rows
    // that name a species (CONCEN × node external flow, or MASS ÷ LperFT3).
    // Empty until a species inflow row exists. Consumed by every quality
    // engine's node mixing: routeLegacyMsx (LEGACY), ArdEngine stage 1b /
    // 1a' (EULERIAN_ARD) and the LARD node stage.
    std::vector<double> msx_ext_mass_in;

    std::vector<double> msx_coupling_queue;

    struct MsxRoutingLedger {
        bool armed = false;
        std::vector<double> init, inflow, outflow, flooding, seepage, reacted, clamped, final_dry;
        void arm(std::size_t ns) {
            armed = true;
            for (auto* v :
                 {&init, &inflow, &outflow, &flooding, &seepage, &reacted, &clamped, &final_dry})
                v->assign(ns, 0.0);
        }
    };
    MsxRoutingLedger msx_ledger;

    MsxSurfaceData surface;
 
    bool configured = false; 

    int n_species() const noexcept { return static_cast<int>(species_name.size()); }
    int find_species(std::string_view n) const {
        for (std::size_t i = 0; i < species_name.size(); ++i)
            if (species_name[i] == n)
                return static_cast<int>(i);
        return -1;
    }

    int find_coef(std::string_view n) const {
        for (std::size_t i = 0; i < coef_name.size(); ++i)
            if (coef_name[i] == n)
                return static_cast<int>(i);
        return -1;
    }

    int find_term(std::string_view n) const {
        for (std::size_t i = 0; i < term_name.size(); ++i)
            if (term_name[i] == n)
                return static_cast<int>(i);
        return -1;
    }

    void clear() {
        MsxSurfaceData keep = std::move(surface);
        *this = ReactionData{};
        surface = std::move(keep);
    }
};
 
} // namespace openswmm
 
#endif // OPENSWMM_ENGINE_DATA_REACTION_DATA_HPP