Simplify Reactor evaluation/synchronization logic

doc/sphinx/python/zerodim.rst

@@ -123,6 +123,14 @@ ReactorSurface
 ^^^^^^^^^^^^^^
 .. autoclass:: ReactorSurface(phase, r=None, clone=None, name="(none)", *, A=None)
 
+ExtensibleReactorSurface
+^^^^^^^^^^^^^^^^^^^^^^^^
+.. autoclass:: ExtensibleReactorSurface(phase, r=None, clone=None, name="(none)", *, A=None)
+
+ExtensibleMoleReactorSurface
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+.. autoclass:: ExtensibleMoleReactorSurface(phase, r=None, clone=None, name="(none)", *, A=None)
+
 Flow Controllers
 ----------------
 

include/cantera/base/ct_defs.h

@@ -27,6 +27,7 @@
 #include <algorithm>
 #include <memory>
 #include <functional>
+#include <span>
 
 /**
  * Namespace for the Cantera kernel.
@@ -45,6 +46,7 @@ using std::map;
 using std::set;
 using std::function;
 using std::pair;
+using std::span;
 
 /**
  * @defgroup physConstants Physical Constants

include/cantera/kinetics/Kinetics.h

@@ -274,6 +274,8 @@ class Kinetics
         return m_start[n] + k;
     }
 
+    size_t speciesOffset(const ThermoPhase& phase) const;
+
     //! Return the name of the kth species in the kinetics manager.
     /*!
      * k is an integer from 0 to ktot - 1, where ktot is the number of

include/cantera/numerics/EigenSparseJacobian.h

@@ -23,6 +23,11 @@ class EigenSparseJacobian : public SystemJacobian
     void updatePreconditioner() override;
     void updateTransient(double rdt, int* mask) override;
 
+    //! Set all Jacobian elements. Replaces the existing elements.
+    void setFromTriplets(const vector<Eigen::Triplet<double>>& trips) {
+        m_jac_trips = trips;
+    }
+
     //! Return underlying Jacobian matrix
     //! @ingroup derivGroup
     Eigen::SparseMatrix<double> jacobian();

include/cantera/zeroD/ConstPressureMoleReactor.h

@@ -20,16 +20,15 @@ namespace Cantera
 class ConstPressureMoleReactor : public MoleReactor
 {
 public:
-    using MoleReactor::MoleReactor; // inherit constructors
+    ConstPressureMoleReactor(shared_ptr<Solution> sol, const string& name="(none)");
+    ConstPressureMoleReactor(shared_ptr<Solution> sol, bool clone,
+                             const string& name="(none)");
 
     string type() const override {
         return "ConstPressureMoleReactor";
     };
 
     void getState(double* y) override;
-
-    void initialize(double t0=0.0) override;
-
     void eval(double t, double* LHS, double* RHS) override;
 
     vector<size_t> steadyConstraints() const override {

include/cantera/zeroD/ConstPressureReactor.h

@@ -23,15 +23,15 @@ namespace Cantera
 class ConstPressureReactor : public Reactor
 {
 public:
-    using Reactor::Reactor; // inherit constructors
+    ConstPressureReactor(shared_ptr<Solution> sol, const string& name="(none)");
+    ConstPressureReactor(shared_ptr<Solution> sol, bool clone,
+                         const string& name="(none)");
 
     string type() const override {
         return "ConstPressureReactor";
     }
 
     void getState(double* y) override;
-
-    void initialize(double t0=0.0) override;
     void eval(double t, double* LHS, double* RHS) override;
     vector<size_t> steadyConstraints() const override;
 

include/cantera/zeroD/FlowReactor.h

@@ -16,7 +16,8 @@ namespace Cantera
 class FlowReactor : public IdealGasReactor
 {
 public:
-    using IdealGasReactor::IdealGasReactor; // inherit constructors
+    FlowReactor(shared_ptr<Solution> sol, const string& name="(none)");
+    FlowReactor(shared_ptr<Solution> sol, bool clone, const string& name="(none)");
 
     string type() const override {
         return "FlowReactor";
@@ -36,8 +37,6 @@ class FlowReactor : public IdealGasReactor
     }
 
     void getStateDae(double* y, double* ydot) override;
-    void initialize(double t0=0.0) override;
-    void syncState() override;
     void updateState(double* y) override;
 
     //! Not implemented; FlowReactor implements evalDae() instead.
@@ -74,102 +73,26 @@ class FlowReactor : public IdealGasReactor
     //! Sets the area of the reactor [m²]
     void setArea(double area) override;
 
-    //! The ratio of the reactor's surface area to volume ratio [m^-1]
-    //! @note If the surface area to volume ratio is unspecified by the user,
-    //!       this will be calculated assuming the reactor is a cylinder.
-    double surfaceAreaToVolumeRatio() const;
-
-    //! Set the reactor's surface area to volume ratio [m^-1]
-    void setSurfaceAreaToVolumeRatio(double sa_to_vol) {
-        m_sa_to_vol = sa_to_vol;
-    }
-
-    //! Get the steady state tolerances used to determine the initial state for
-    //! surface coverages
-    double inletSurfaceAtol() const {
-        return m_ss_atol;
-    }
-
-    //! Set the steady state tolerances used to determine the initial state for
-    //! surface coverages
-    void setInletSurfaceAtol(double atol) {
-        m_ss_atol = atol;
-    }
-
-    //! Get the steady state tolerances used to determine the initial state for
-    //! surface coverages
-    double inletSurfaceRtol() const {
-        return m_ss_rtol;
-    }
-
-    //! Set the steady state tolerances used to determine the initial state for
-    //! surface coverages
-    void setInletSurfaceRtol(double rtol) {
-        m_ss_rtol = rtol;
-    }
-
-    //! Get the steady state tolerances used to determine the initial state for
-    //! surface coverages
-    double inletSurfaceMaxSteps() const {
-        return m_max_ss_steps;
-    }
-
-    //! Set the steady state tolerances used to determine the initial state for
-    //! surface coverages
-    void setInletSurfaceMaxSteps(int max_steps) {
-        m_max_ss_steps = max_steps;
-    }
-
-    //! Get the steady state tolerances used to determine the initial state for
-    //! surface coverages
-    double inletSurfaceMaxErrorFailures() const {
-        return m_max_ss_error_fails;
-    }
-
-    //! Set the steady state tolerances used to determine the initial state for
-    //! surface coverages
-    void setInletSurfaceMaxErrorFailures(int max_fails) {
-        m_max_ss_error_fails = max_fails;
-    }
-
     //! Return the index in the solution vector for this reactor of the component named
     //! *nm*. Possible values for *nm* are "density", "speed", "pressure",
-    //! "temperature", the name of a homogeneous phase species, or the name of a surface
-    //! species.
+    //! "temperature" or the name of a homogeneous phase species.
     size_t componentIndex(const string& nm) const override;
 
     string componentName(size_t k) override;
 
-    void updateSurfaceState(double* y) override;
-
 protected:
     //! Density [kg/m^3]. First component of the state vector.
     double m_rho = NAN;
     //! Axial velocity [m/s]. Second component of the state vector.
     double m_u = -1.0;
-    //! Pressure [Pa]. Third component of the state vector.
-    double m_P = NAN;
-    //! Temperature [K]. Fourth component of the state vector.
-    double m_T = NAN;
     //! offset to the species equations
     const size_t m_offset_Y = 4;
     //! reactor area [m^2]
     double m_area = 1.0;
-    //! reactor surface area to volume ratio [m^-1]
-    double m_sa_to_vol = -1.0;
-    //! temporary storage for surface species production rates
-    vector<double> m_sdot_temp;
     //! temporary storage for species partial molar enthalpies
     vector<double> m_hk;
-    //! steady-state relative tolerance, used to determine initial surface coverages
-    double m_ss_rtol = 1e-7;
-    //! steady-state absolute tolerance, used to determine initial surface coverages
-    double m_ss_atol = 1e-14;
-    //! maximum number of steady-state coverage integrator-steps
-    int m_max_ss_steps = 20000;
-    //! maximum number of steady-state integrator error test failures
-    int m_max_ss_error_fails = 10;
 };
+
 }
 
 #endif

include/cantera/zeroD/IdealGasConstPressureMoleReactor.h

@@ -33,13 +33,14 @@ class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor
     void eval(double t, double* LHS, double* RHS) override;
 
     void updateState(double* y) override;
+    void getJacobianScalingFactors(double& f_species, double* f_energy) override;
 
     //! Calculate an approximate Jacobian to accelerate preconditioned solvers
 
     //! Neglects derivatives with respect to mole fractions that would generate a
     //! fully-dense Jacobian. Currently also neglects terms related to interactions
     //! between reactors, for example via inlets and outlets.
-    Eigen::SparseMatrix<double> jacobian() override;
+    void getJacobianElements(vector<Eigen::Triplet<double>>& trips) override;
 
     bool preconditionerSupported() const override { return true; };
 
@@ -50,6 +51,7 @@ class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor
 
 protected:
     vector<double> m_hk; //!< Species molar enthalpies
+    double m_TotalCp; //!< Total heat capacity (@f$ m c_p @f$) [J/K]
 };
 
 }

include/cantera/zeroD/IdealGasMoleReactor.h

@@ -45,12 +45,14 @@ class IdealGasMoleReactor : public MoleReactor
     //! Neglects derivatives with respect to mole fractions that would generate a
     //! fully-dense Jacobian. Currently, also neglects terms related to interactions
     //! between reactors, for example via inlets and outlets.
-    Eigen::SparseMatrix<double> jacobian() override;
+    void getJacobianElements(vector<Eigen::Triplet<double>>& trips) override;
+    void getJacobianScalingFactors(double& f_species, double* f_energy) override;
 
     bool preconditionerSupported() const override {return true;};
 
 protected:
     vector<double> m_uk; //!< Species molar internal energies
+    double m_TotalCv; //!< Total heat capacity (@f$ m c_v @f$) [J/K]
 };
 
 }

include/cantera/zeroD/MoleReactor.h

@@ -20,32 +20,23 @@ namespace Cantera
 class MoleReactor : public Reactor
 {
 public:
-    using Reactor::Reactor; // inherit constructors
+    MoleReactor(shared_ptr<Solution> sol, const string& name="(none)");
+    MoleReactor(shared_ptr<Solution> sol, bool clone, const string& name="(none)");
 
     string type() const override {
         return "MoleReactor";
     }
 
-    void initialize(double t0=0.0) override;
-
     void getState(double* y) override;
-
     void updateState(double* y) override;
-
     void eval(double t, double* LHS, double* RHS) override;
-
     size_t componentIndex(const string& nm) const override;
     string componentName(size_t k) override;
     double upperBound(size_t k) const override;
     double lowerBound(size_t k) const override;
     void resetBadValues(double* y) override;
 
 protected:
-    //! For each surface in the reactor, update vector of triplets with all relevant
-    //! surface jacobian derivatives of species with respect to species
-    //! which are appropriately offset to align with the reactor's state vector.
-    virtual void addSurfaceJacobian(vector<Eigen::Triplet<double>> &triplets);
-
     //! Get moles of the system from mass fractions stored by thermo object
     //! @param y vector for moles to be put into
     void getMoles(double* y);
@@ -54,12 +45,6 @@ class MoleReactor : public Reactor
     //! @param y vector of moles of the system
     void setMassFromMoles(double* y);
 
-    void evalSurfaces(double* LHS, double* RHS, double* sdot) override;
-
-    void updateSurfaceState(double* y) override;
-
-    void getSurfaceInitialConditions(double* y) override;
-
     //! const value for the species start index
     const size_t m_sidx = 2;
 };