ReactionIntegrator.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_REACTION_INTEGRATOR_HPP
#define OPENSWMM_ENGINE_TRANSPORT_REACTION_INTEGRATOR_HPP
 
#include <string>
#include <span>
#include <vector>
 
#include "../../../data/ReactionData.hpp"
 
namespace openswmm::transport {

struct RxStepReport {
    bool ok = true;
    bool interval_limit = false; 
    int substeps = 0;            
    int newton_iters = 0;        
    std::string error;           
};
// Several aqueous phase contexts may contribute RATE terms to the same
// physical inventory. Each scale is phase aqueous volume / shared capacity.
// All contexts see the same species state at every integrator evaluation.
// Algebraic phase projections retain their separate existing contract.
struct RxRateContext {
    const std::vector<ReactionExprForm>* forms = nullptr;
    const std::vector<RxExprSpan>* expressions = nullptr;
    double* hydvar = nullptr;
    const double* pollutants = nullptr;
    const double* kinetic_scale = nullptr;
};

class RxWorkspace {
public:
    void init(const ReactionData& rx);
    bool initialized() const noexcept { return n_ > 0 || terms_.capacity() > 0; }
 
private:
    friend class ReactionIntegrator;
    int n_ = 0;                              
    std::vector<unsigned char> term_needed_; 
    std::vector<double> terms_;              
    std::vector<double> y_, y0_, ytmp_, rates_, err_;
    std::vector<double> k1_, k2_, k3_, k4_, k5_, k6_; 
    std::vector<double> jac_, lu_;                    
    std::vector<int> piv_;
    std::vector<int> rate_idx_, equil_idx_, formula_idx_;
    std::vector<double> res_, res2_, dy_; 
    std::vector<double> grp_out_;         
 
    // FD-Jacobian / LU reuse across substeps. The Jacobian costs `gn` extra
    // RHS evaluations, each of which re-evaluates every term and rate
    // expression — the dominant per-substep cost once a system has more than
    // a species or two. J depends on y (not on h), so it survives a step-size
    // change; the FACTORED matrix I - c*h*J does not, hence the separate
    // scale key.
    bool jac_valid_ = false; 
    int jac_age_ = 0;        
    bool lu_valid_ = false;  
    double lu_scale_ = 0.0;  
};
 
class ReactionIntegrator {
public:
    // Explicit reactor context. Expressions and forms refer to the same
    // compiled token pool; groundwater never borrows pipe/tank laws.
    static RxStepReport step(const ReactionData& rx, const std::vector<ReactionExprForm>& forms,
                             const std::vector<RxExprSpan>& expressions, double dt, double* species,
                             double* hydvar, RxWorkspace& ws, const double* pollutants = nullptr,
                             const double* kinetic_scale = nullptr,
                             const double* coefficients = nullptr,
                             std::span<const RxRateContext> contexts = {});
    // Opt-in transactional retry for an owned interval. Only a work/floor limit
    // permits bisection; invalid laws and Newton failures retain their error.
    // Tolerances are unchanged, and failure restores every supplied species.
    static RxStepReport
    stepSubcycled(const ReactionData& rx, const std::vector<ReactionExprForm>& forms,
                  const std::vector<RxExprSpan>& expressions, double dt, double* species,
                  double* hydvar, RxWorkspace& ws, const double* pollutants = nullptr,
                  const double* kinetic_scale = nullptr, const double* coefficients = nullptr,
                  std::span<const RxRateContext> contexts = {});
    static RxStepReport step(const ReactionData& rx, bool tank, double dt, double* species,
                             double* hydvar, RxWorkspace& ws, const double* pollutants = nullptr,
                             const double* coefficients = nullptr);
};
 
} // namespace openswmm::transport
 
#endif // OPENSWMM_ENGINE_TRANSPORT_REACTION_INTEGRATOR_HPP