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 <vector>
 
#include "../../../data/ReactionData.hpp"
 
namespace openswmm::transport {

struct RxStepReport {
    bool ok            = true;
    int  substeps      = 0;   
    int  newton_iters  = 0;   
    std::string error;        
};

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<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:
    static RxStepReport step(const ReactionData& rx, bool tank, double dt,
                             double* species, double* hydvar, RxWorkspace& ws,
                             const double* pollutants = nullptr);
};
 
}  // namespace openswmm::transport
 
#endif  // OPENSWMM_ENGINE_TRANSPORT_REACTION_INTEGRATOR_HPP