ArdEngine.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_TRANSPORT_ARD_ENGINE_HPP
#define OPENSWMM_ENGINE_TRANSPORT_ARD_ENGINE_HPP
 
#include <fstream>
#include <string>
#include <vector>
 
#include "../../../hydraulics/fv/FvKernels.hpp"
#include "../../../hydraulics/fv/NetworkMeshData.hpp"
#include "../../../data/ChannelThermalData.hpp"
 
namespace openswmm {
struct SimulationContext;
}
 
namespace openswmm::transport {
 
class ArdEngine {
public:
    bool init(SimulationContext& ctx);

    void step(SimulationContext& ctx, double dt);

    bool initialized() const noexcept { return initialized_; }
    std::string checkpointState() const;
    // Decode only into a fresh staging engine; caller installs after every
    // physical/output record has passed validation. Detailed writers excluded.
    void restoreCheckpointState(const std::string&, SimulationContext&);
    const std::vector<std::string>& warnings() const noexcept { return warnings_; }

    void absorbTreatedNodeConc(SimulationContext& ctx);

    std::vector<double> totalMass(const SimulationContext& ctx) const;
    std::vector<ChannelThermalRegion> thermalRegions() const;
    std::string applyThermalIncrements(SimulationContext&,
                                       const std::vector<ChannelThermalIncrement>&);
    std::vector<ChannelSpeciesRegion> speciesRegions(const SimulationContext&) const;
    std::string applySpeciesIncrements(SimulationContext&,
                                       const std::vector<ChannelSpeciesIncrement>&);
 
private:
    // Continuity-consistent face fluxes from the link/node fields of the
    // current routing step (plan §3.2): within each conduit the face flux
    // ramps linearly so every cell receives an equal share of the conduit's
    // volume change, and the two end faces carry the flux the node exchange
    // actually saw. Splice faces (virtual junctions) take the mean of the
    // two conduits' end values.
    void projectHydraulics(SimulationContext& ctx, double dt);
 
    // One advection substep of length dt_sub: kernel reconstruction + FCT on
    // interior faces, donor-upwinded node boundary fluxes, cell mass/area
    // update, node CSTR update (external loads prorated by dt_sub / dt_step),
    // then the implicit dispersion solve over the updated cells (E3,
    // sequential Lie split: one full advection step, then one full
    // dispersion step — NOT Strang, which would halve the advection).
    void substep(SimulationContext& ctx, double dt_sub, double load_frac);
 
    // E3: refresh cell_disp_ from the routing step's hydraulics. Per-conduit
    // overrides always win; otherwise VALUE mode broadcasts the global
    // coefficient and FISCHER mode evaluates D = 0.011 v²B²/(Y·U*),
    // U* = √(g·Y·S) from the link fields (all internal ft units).
    void updateDispersion(SimulationContext& ctx);
 
    void publish(SimulationContext& ctx);
 
    fv::NetworkMeshData mesh_;
    fv::NetworkStateData state_;
 
    // Projected per-face volumetric flux for the current routing step and the
    // per-cell start-of-step areas it is consistent with.
    std::vector<double> face_q_;
 
    // Node stores: species mass [node * ns + s] and water volume the mass
    // sits in (tracked from the same fluxes, so node concentration is
    // self-consistent between substeps).
    std::vector<double> node_mass_;
    std::vector<double> node_vol_;
 
    // Lateral loads at VIRTUAL junctions (plans/VJ_LATERAL_INFLOW_PLAN_
    // 2026-09-04.md §E4). A virtual junction owns no faces — its conduits
    // were spliced into one interior face — so a node store there could
    // never discharge and the per-step resync zeroed it: the load was
    // destroyed. Its water and mass go straight into the two cells adjoining
    // the splice, half each, mirroring the hydraulic solver's cell_qlat_
    // split. Rebuilt every substep; has_vj_src_ gates the whole path so a
    // deck without a fed virtual junction executes no new arithmetic.
    std::vector<double> cell_src_vol_;  
    std::vector<double> cell_src_mass_; 
    bool has_vj_src_ = false;
 
    // Kernel scratch (sized once in init; see SpeciesKernelView).
    std::vector<double> f_mass_, f_sstar_, f_phi_l_, f_phi_r_, f_phi_flux_, cell_slope_, lo_flux_,
        anti_flux_, td_, anew_, rplus_, rminus_, cell_u_;
    std::vector<fv::kernels::FaceState> f_state_l_, f_state_r_;
    std::vector<fv::kernels::FaceFlux> f_flux_;
    std::vector<char> cell_active_;
    std::vector<int> active_faces_;
 
    // E5a: transport boundaries/sources resolved onto the mesh at init.
    // Boundaries: the node's external inflow water carries the species at
    // bc_now_ (evaluated once per routing step: VALUE or timeseries).
    // Sources: internal-rate mass distributed over the conduit's wet cells.
    void updateTransportRows(SimulationContext& ctx);
    std::vector<int> bc_node_, bc_srow_, bc_ts_;
    std::vector<double> bc_value_, bc_now_;
    std::vector<int> src_crow_, src_srow_, src_ts_;
    std::vector<double> src_value_, src_len_, src_now_;
 
    // E3 dispersion state. disp_active_ gates the whole path: when false the
    // substep passes no per-cell array and dispersionSolve early-outs — the
    // pre-E3 behavior, bit-identical for existing ARD decks.
    bool disp_active_ = false;
    int disp_mode_ = 0;                   
    double disp_global_ft_ = 0.0;         
    std::vector<double> conduit_disp_ft_; 
    std::vector<double> cell_disp_;       
 
    // E5b: per-cell CSV sidecar ([TRANSPORT_OPTIONS] DETAILED_OUTPUT).
    // Written every routing step (documented decision — a detail feature
    // for short diagnostic runs; cadence control can come later if sizes
    // demand it). Columns: time_s, element, kind (L link cell / N node
    // store), cell index, species name, concentration.
    void writeDetailRows(SimulationContext& ctx);
    std::ofstream detail_out_;
    bool detail_active_ = false;
    double detail_time_s_ = 0.0;

    int age_row_ = -1;

    int nm_ = 0;

    void applyHeatFluxes(SimulationContext& ctx, double dt);

    void applyBedSoluteExchange(SimulationContext& ctx, double dt);

    int cellLink(std::size_t cell) const noexcept;

    int temp_row_ = -1;
 
    std::vector<std::string> warnings_;
    bool initialized_ = false;
    bool warned_cfl_clamp_ = false; 
};
 
} // namespace openswmm::transport
 
#endif // OPENSWMM_ENGINE_TRANSPORT_ARD_ENGINE_HPP