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
@@ -47,13 +48,13 @@ reactor-integration
 
 ## Adding New Features to Cantera
 
-- [](CONTRIBUTING)
-- [](style-guidelines)
-- [](vscode-tips)
-- [](writing-tests)
-- [](running-tests)
-- [](writing-examples)
-- [](doc-formatting)
+-   [](CONTRIBUTING)
+-   [](style-guidelines)
+-   [](vscode-tips)
+-   [](writing-tests)
+-   [](running-tests)
+-   [](writing-examples)
+-   [](doc-formatting)
 
 ```{toctree}
 :caption: Adding New Features to Cantera

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

@@ -9,6 +9,7 @@
 #include "cantera/base/ct_defs.h"
 #include "cantera/base/global.h"
 #include "cantera/base/ctexceptions.h"
+#include "cantera/numerics/eigen_sparse.h"
 
 namespace Cantera
 {
@@ -132,9 +133,41 @@ class FlowDevice
         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;
 

include/cantera/zeroD/IdealGasConstPressureMoleReactor.h

@@ -43,10 +43,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; };
+
 protected:
     void setThermo(ThermoPhase& thermo) override;
 

include/cantera/zeroD/IdealGasMoleReactor.h

@@ -38,15 +38,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:
     void setThermo(ThermoPhase& thermo) override;
 

include/cantera/zeroD/MoleReactor.h

@@ -7,6 +7,7 @@
 #define CT_MOLEREACTOR_H
 
 #include "Reactor.h"
+#include "cantera/zeroD/Wall.h"
 
 namespace Cantera
 {
@@ -38,6 +39,8 @@ class MoleReactor : public Reactor
 
     string componentName(size_t k) 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
@@ -60,6 +63,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

@@ -70,12 +70,6 @@ class Reactor : public ReactorBase
         m_chem = cflag;
     }
 
-    //! Returns `true` if changes in the reactor composition due to chemical reactions
-    //! are enabled.
-    bool chemistryEnabled() const {
-        return m_chem;
-    }
-
     void setEnergy(int eflag=1) override {
         if (eflag > 0) {
             m_energy = true;
@@ -84,11 +78,6 @@ class Reactor : public ReactorBase
         }
     }
 
-    //! Returns `true` if solution of the energy equation is enabled.
-    bool energyEnabled() const {
-        return m_energy;
-    }
-
     //! Number of equations (state variables) for this reactor
     size_t neq() {
         if (!m_nv) {
@@ -187,18 +176,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
@@ -280,8 +302,6 @@ class Reactor : public ReactorBase
 
     vector<double> m_wdot; //!< Species net molar production rates
     vector<double> m_uk; //!< Species molar internal energies
-    bool m_chem = false;
-    bool m_energy = true;
     size_t m_nv = 0;
     size_t m_nv_surf; //!!< Number of variables associated with reactor surfaces
 
@@ -292,6 +312,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

@@ -266,6 +266,34 @@ class ReactorBase
     //! Set the ReactorNet that this reactor belongs to.
     void setNetwork(ReactorNet* net);
 
+    //! 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; };
+
+    //! Returns `true` if solution of the energy equation is enabled.
+    virtual bool energyEnabled() const {
+        return m_energy;
+    }
+
+    //! Returns `true` if changes in the reactor composition due to chemical reactions
+    //! are enabled.
+    bool chemistryEnabled() const {
+        return m_chem;
+    }
+
 protected:
     //! Specify the mixture contained in the reactor. Note that a pointer to
     //! this substance is stored, and as the integration proceeds, the state of
@@ -282,6 +310,12 @@ class ReactorBase
     //! 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 = 1.0; //!< Current volume of the reactor [m^3]
     double m_enthalpy = 0.0; //!< Current specific enthalpy of the reactor [J/kg]
@@ -304,6 +338,12 @@ class ReactorBase
 
     //! Composite thermo/kinetics/transport handler
     shared_ptr<Solution> m_solution;
+
+    //! A bool that enables the energy equation
+    bool m_energy = true;
+
+    //! A bool that enables the chemical kinetics equations
+    bool m_chem = false;
 };
 }