Extensible jacobian and network jacobian interfaces

doc/sphinx/develop/index.md

@@ -1,6 +1,7 @@
 # Develop
 
 (sec-compiling)=
+
 ## Compiling Cantera from Source
 
 If you're interested in contributing new features to Cantera, or you want to try the
@@ -12,9 +13,9 @@ Cantera](compiling/configure-build) on your computer.
 
 The following additional references may also be useful:
 
-- [](compiling/dependencies.md)
-- [](compiling/config-options)
-- [](compiling/special-cases)
+-   [](compiling/dependencies.md)
+-   [](compiling/config-options)
+-   [](compiling/special-cases)
 
 ```{toctree}
 :caption: Compiling Cantera from Source
@@ -35,7 +36,7 @@ compiling/special-cases
 This section is a work in progress.
 ```
 
-- [](reactor-integration)
+-   [](reactor-integration)
 
 ```{toctree}
 :caption: How Cantera Works

include/cantera/base/Delegator.h

@@ -10,6 +10,7 @@
 #include "cantera/base/Units.h"
 #include "cantera/base/ctexceptions.h"
 #include "cantera/base/ExtensionManager.h"
+#include "cantera/numerics/eigen_sparse.h"
 #include <array>
 #include <list>
 
@@ -226,6 +227,20 @@ class Delegator
         *m_funcs_v_d_dp_dp[name] = makeDelegate(func, when, *m_funcs_v_d_dp_dp[name]);
     }
 
+    //! Set delegates for member functions with the signature
+    //! `void(vector<Eigen::Triplet<double>>&)`
+    void setDelegate(const string& name,
+                     const function<void(vector<Eigen::Triplet<double>>&)>& func,
+                     const string& when)
+    {
+        if (!m_funcs_v_vet.count(name)) {
+            throw NotImplementedError("Delegator::setDelegate",
+                "for function '{}' with signature "
+                "'void(vector<Eigen::Triplet<double>>&)'.", name);
+        }
+        *m_funcs_v_vet[name] = makeDelegate(func, when, *m_funcs_v_vet[name]);
+    }
+
     //! Set delegates for member functions with the signature
     //! `void(double*, double*, double*)`
     void setDelegate(
@@ -386,6 +401,16 @@ class Delegator
         m_funcs_v_dp_dp_dp[name] = &target;
     }
 
+    //! Install a function with the signature `void(vector<Eigen::Triplet<double>>&)`
+    //! as being delegatable
+    void install(const string& name,
+                 function<void(vector<Eigen::Triplet<double>>&)>& target,
+                 const function<void(vector<Eigen::Triplet<double>>&)>& base)
+    {
+        target = base;
+        m_funcs_v_vet[name] = &target;
+    }
+
     //! Install a function with the signature `double(void*)` as being delegatable
     void install(const string& name, function<double(void*)>& target,
                  const function<double(void*)>& func)
@@ -514,6 +539,7 @@ class Delegator
     //! - `sz` for `size_t`
     //! - `AM` for `AnyMap`
     //! - `US` for `UnitStack`
+    //! - `VET` for `vector<Eigen::Triplet<double>>`
     //! - prefix `c` for `const` arguments
     //! - suffix `r` for reference arguments
     //! - suffix `p` for pointer arguments
@@ -532,6 +558,7 @@ class Delegator
         function<void(std::array<size_t, 2>, double, double*, double*)>*> m_funcs_v_d_dp_dp;
     map<string,
         function<void(std::array<size_t, 3>, double*, double*, double*)>*> m_funcs_v_dp_dp_dp;
+    map<string, function<void(vector<Eigen::Triplet<double>>&)>*> m_funcs_v_vet;
 
     // Delegates with a return value
     map<string, function<double(void*)>> m_base_d_vp;
@@ -542,6 +569,7 @@ class Delegator
 
     map<string, function<size_t(const string&)>> m_base_sz_csr;
     map<string, function<size_t(const string&)>*> m_funcs_sz_csr;
+
     //! @}
 
     //! Handles to wrappers for the delegated object in external language interfaces.

include/cantera/zeroD/FlowDevice.h

@@ -10,6 +10,7 @@
 #include "cantera/base/global.h"
 #include "cantera/base/ctexceptions.h"
 #include "ConnectorNode.h"
+#include "cantera/numerics/eigen_sparse.h"
 
 namespace Cantera
 {
@@ -140,7 +141,42 @@ class FlowDevice : public ConnectorNode
         m_time = time;
     }
 
+    //! Build the Jacobian terms specific to the flow device for the given connected
+    //! reactor.
+    //! @param r a pointer to the calling reactor
+    //! @param jacVector a vector of triplets to be added to the jacobian for the
+    //! reactor
+    //! @warning This function is an experimental part of the %Cantera API and may be
+    //! changed
+    //! or removed without notice.
+    //! @since New in %Cantera 3.1.
+    //!
+    virtual void buildReactorJacobian(ReactorBase* r,
+        vector<Eigen::Triplet<double>>& jacVector) {
+        throw NotImplementedError(type() + "::buildReactorJacobian");
+    }
+
+    //! Build the Jacobian terms specific to the flow device for the network. These
+    //! terms
+    //! will be adjusted to the networks indexing system outside of the reactor.
+    //! @param jacVector a vector of triplets to be added to the jacobian for the
+    //! reactor
+    //! @warning This function is an experimental part of the %Cantera API and may be
+    //! changed
+    //! or removed without notice.
+    //! @since New in %Cantera 3.1.
+    //!
+    virtual void buildNetworkJacobian(vector<Eigen::Triplet<double>>& jacVector) {
+        throw NotImplementedError(type() + "::buildNetworkJacobian");
+    }
+
 protected:
+    string m_name;  //!< Flow device name.
+    bool m_defaultNameSet = false;  //!< `true` if default name has been previously set.
+    //! a variable to switch on and off so calculations are not doubled by the calling
+    //! reactor or network
+    bool m_jac_calculated = false;
+
     double m_mdot = Undef;
 
     //! Function set by setPressureFunction; used by updateMassFlowRate

include/cantera/zeroD/IdealGasConstPressureMoleReactor.h

@@ -39,10 +39,14 @@ class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor
     //! 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 buildJacobian(vector<Eigen::Triplet<double>>& jacVector) override;
 
     bool preconditionerSupported() const override { return true; };
 
+    double temperature_ddni(size_t index) override;
+
+    size_t speciesOffset() const override { return m_sidx; };
+
     size_t componentIndex(const string& nm) const override;
     string componentName(size_t k) override;
     double upperBound(size_t k) const override;

include/cantera/zeroD/IdealGasMoleReactor.h

@@ -40,15 +40,15 @@ class IdealGasMoleReactor : public MoleReactor
 
     void updateState(double* y) 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;
+    //! Calculate the Jacobian to accelerate preconditioned solvers
+    void buildJacobian(vector<Eigen::Triplet<double>>& jacVector) override;
 
     bool preconditionerSupported() const override {return true;};
 
+    double temperature_ddni(size_t index) override;
+
+    size_t speciesOffset() const override { return m_sidx; };
+
 protected:
     vector<double> m_uk; //!< Species molar internal energies
 };

include/cantera/zeroD/MoleReactor.h

@@ -7,6 +7,7 @@
 #define CT_MOLEREACTOR_H
 
 #include "Reactor.h"
+#include "cantera/zeroD/Wall.h"
 
 namespace Cantera
 {
@@ -40,6 +41,8 @@ class MoleReactor : public Reactor
     double lowerBound(size_t k) const override;
     void resetBadValues(double* y) override;
 
+    size_t energyIndex() const override { return m_eidx; };
+
 protected:
     //! For each surface in the reactor, update vector of triplets with all relevant
     //! surface jacobian derivatives of species with respect to species
@@ -62,6 +65,9 @@ class MoleReactor : public Reactor
 
     //! const value for the species start index
     const size_t m_sidx = 2;
+
+    //! index of state variable associated with energy
+    const size_t m_eidx = 0;
 };
 
 }

include/cantera/zeroD/Reactor.h

@@ -199,18 +199,51 @@ class Reactor : public ReactorBase
     //! @param limit value for step size limit
     void setAdvanceLimit(const string& nm, const double limit);
 
+    //! A wrapper for the Jacobian function to return the Eigen::SparseMatrix<double>
+    //! @warning Depending on the particular implementation, this may return an
+    //! approximate Jacobian intended only for use in forming a preconditioner for
+    //! iterative solvers.
+    //! @ingroup derivGroup
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    virtual Eigen::SparseMatrix<double> jacobian();
+
     //! Calculate the Jacobian of a specific Reactor specialization.
+    //! @param jacVector vector where jacobian triplets are added
     //! @warning Depending on the particular implementation, this may return an
     //! approximate Jacobian intended only for use in forming a preconditioner for
     //! iterative solvers.
     //! @ingroup derivGroup
     //!
     //! @warning  This method is an experimental part of the %Cantera
     //! API and may be changed or removed without notice.
-    virtual Eigen::SparseMatrix<double> jacobian() {
-        throw NotImplementedError("Reactor::jacobian");
+    virtual void buildJacobian(vector<Eigen::Triplet<double>>& jacVector) {
+        throw NotImplementedError(type() + "::buildJacobian");
     }
 
+    //! Calculate the Jacobian of a Reactor specialization for wall contributions.
+    //! @param jacVector vector where jacobian triplets are added
+    //! @warning Depending on the particular implementation, this may return an
+    //! approximate Jacobian intended only for use in forming a preconditioner for
+    //! iterative solvers.
+    //! @ingroup derivGroup
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    virtual void buildWallJacobian(vector<Eigen::Triplet<double>>& jacVector);
+
+    //! Calculate flow contributions to the Jacobian of a Reactor specialization.
+    //! @param jacVector vector where jacobian triplets are added
+    //! @warning Depending on the particular implementation, this may return an
+    //! approximate Jacobian intended only for use in forming a preconditioner for
+    //! iterative solvers.
+    //! @ingroup derivGroup
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    virtual void buildFlowJacobian(vector<Eigen::Triplet<double>>& jacVector);
+
     //! Calculate the reactor-specific Jacobian using a finite difference method.
     //!
     //! This method is used only for informational purposes. Jacobian calculations
@@ -298,6 +331,10 @@ class Reactor : public ReactorBase
 
     //! Vector of triplets representing the jacobian
     vector<Eigen::Triplet<double>> m_jac_trips;
+    //! Boolean to skip walls in jacobian
+    bool m_jac_skip_walls = false;
+    //! Boolean to skip flow devices in jacobian
+    bool m_jac_skip_flow_devices = false;
 };
 }
 

include/cantera/zeroD/ReactorBase.h

@@ -302,12 +302,35 @@ class ReactorBase
         return m_sensParams.size();
     }
 
+    //! Calculate the derivative of T with respect to the ith species in the energy
+    //! conservation equation based on the reactor specific equation of state.
+    //! @param index index of the species the derivative is with respect too
+    //! @warning This function is an experimental part of the %Cantera API and may
+    //! be changed or removed without notice.
+    //! @since New in %Cantera 3.1.
+    //!
+    virtual double temperature_ddni(size_t index) {
+        throw NotImplementedError("Reactor::temperature_ddni");
+    }
+
+    //! Return the index associated with energy of the system
+    virtual size_t energyIndex() const { return m_eidx; };
+
+    //! Return the offset between species and state variables
+    virtual size_t speciesOffset() const { return m_sidx; };
+
 protected:
     explicit ReactorBase(const string& name="(none)");
 
     //! Number of homogeneous species in the mixture
     size_t m_nsp = 0;
 
+    //! species offset in the state vector
+    const size_t m_sidx = 3;
+
+    //! index of state variable associated with energy
+    const size_t m_eidx = 1;
+
     ThermoPhase* m_thermo = nullptr;
     double m_vol = 0.0; //!< Current volume of the reactor [m^3]
     double m_mass = 0.0; //!< Current mass of the reactor [kg]

include/cantera/zeroD/ReactorDelegator.h

@@ -10,6 +10,7 @@
 #include "cantera/base/Delegator.h"
 #include "cantera/zeroD/ReactorSurface.h"
 #include "cantera/thermo/SurfPhase.h"
+#include "cantera/numerics/eigen_sparse.h"
 
 namespace Cantera
 {
@@ -51,6 +52,10 @@ class ReactorAccessor
     //! Set the state of the thermo object for surface *n* to correspond to the
     //! state of that surface
     virtual void restoreSurfaceState(size_t n) = 0;
+
+    //! Public access to the default evaluation function so it can be used in
+    //! replace functions
+    virtual void defaultEval(double t, double* LHS, double* RHS) = 0;
 };
 
 //! Delegate methods of the Reactor class to external functions
@@ -69,7 +74,8 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor
         install("updateState", m_updateState,
             [this](std::array<size_t, 1> sizes, double* y) { R::updateState(y); });
         install("updateSurfaceState", m_updateSurfaceState,
-            [this](std::array<size_t, 1> sizes, double* y) { R::updateSurfaceState(y); });
+            [this](std::array<size_t, 1> sizes, double* y)
+            { R::updateSurfaceState(y); });
         install("getSurfaceInitialConditions", m_getSurfaceInitialConditions,
             [this](std::array<size_t, 1> sizes, double* y) {
                 R::getSurfaceInitialConditions(y);
@@ -84,8 +90,8 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor
         );
         install("evalWalls", m_evalWalls, [this](double t) { R::evalWalls(t); });
         install("evalSurfaces", m_evalSurfaces,
-            [this](std::array<size_t, 3> sizes, double* LHS, double* RHS, double* sdot) {
-                R::evalSurfaces(LHS, RHS, sdot);
+            [this](std::array<size_t, 3> sizes, double* LHS, double* RHS, double* sd) {
+                R::evalSurfaces(LHS, RHS, sd);
             }
         );
         install("componentName", m_componentName,
@@ -94,6 +100,8 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor
             [this](const string& nm) { return R::componentIndex(nm); });
         install("speciesIndex", m_speciesIndex,
             [this](const string& nm) { return R::speciesIndex(nm); });
+        install("buildJacobian", m_build_jacobian,
+            [this](vector<Eigen::Triplet<double>>& jv) { R::buildJacobian(jv); });
     }
 
     // Overrides of Reactor methods
@@ -160,6 +168,14 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor
         return m_speciesIndex(nm);
     }
 
+    void buildJacobian(vector<Eigen::Triplet<double>>& jacVector) override {
+        m_build_jacobian(jacVector);
+    }
+
+    void defaultEval(double t, double* LHS, double* RHS) override {
+        R::eval(t, LHS, RHS);
+    }
+
     // Public access to protected Reactor variables needed by derived classes
 
     void setNEq(size_t n) override {
@@ -204,6 +220,7 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor
     function<string(size_t)> m_componentName;
     function<size_t(const string&)> m_componentIndex;
     function<size_t(const string&)> m_speciesIndex;
+    function<void(vector<Eigen::Triplet<double>>&)> m_build_jacobian;
 };
 
 }