Build jacobian structure defined for walls and flow devices

Author: anthony-walker

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
 {
@@ -114,7 +115,47 @@ 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.0.
+    */
+    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.0.
+    */
+    virtual void buildNetworkJacobian(vector<Eigen::Triplet<double>>& jacVector) {
+        throw NotImplementedError(type() + "::buildNetworkJacobian");
+    }
+
+    /*! Specify the jacobian terms have been calculated and should not be recalculated.
+     * @warning This function is an experimental part of the %Cantera API and may be changed
+     * or removed without notice.
+     * @since New in %Cantera 3.0.
+    */
+    void jacobianCalculated() { m_jac_calculated = true; };
+
+    /*! Specify that jacobian terms have not been calculated and should be recalculated.
+     * @warning This function is an experimental part of the %Cantera API and may be changed
+     * or removed without notice.
+     * @since New in %Cantera 3.0.
+    */
+    void jacobianNotCalculated() { m_jac_calculated = false; };
+
 protected:
+    //! 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/IdealGasMoleReactor.h

@@ -45,6 +45,12 @@ class IdealGasMoleReactor : public MoleReactor
 
     bool preconditionerSupported() const override {return true;};
 
+    double moleDerivative(size_t index) override;
+
+    double moleRadiationDerivative(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
 {
@@ -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

@@ -232,7 +232,27 @@ class Reactor : public ReactorBase
         throw NotImplementedError(type() + "::buildJacobian");
     }
 
-    // virtual void jacobian()
+    //! 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 jac_vector 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.
     //!
@@ -325,6 +345,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

@@ -248,10 +248,46 @@ 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 heat
+     * transfer equation based on the reactor specific equation of state.
+     * will be adjusted to the networks indexing system outside of the reactor.
+     * @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.0.
+    */
+    virtual double moleDerivative(size_t index) {
+        throw NotImplementedError("Reactor::moleDerivative");
+    }
+
+    /*! Calculate the derivative of T with respect to the ith species in the heat
+     * transfer radiation equation based on the reactor specific equation of state.
+     * will be adjusted to the networks indexing system outside of the reactor.
+     * @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.0.
+    */
+    virtual double moleRadiationDerivative(size_t index) {
+        throw NotImplementedError("Reactor::moleRadiationDerivative");
+    }
+
+    //! 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:
     //! 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]

include/cantera/zeroD/ReactorNet.h

@@ -236,6 +236,20 @@ class ReactorNet : public FuncEval
     //! reactor network.
     size_t globalComponentIndex(const string& component, size_t reactor=0);
 
+    //! Return the index corresponding to the start of the reactor specific state
+    //! vector in the reactor with index *reactor* in the global state vector for the
+    //! reactor network.
+    size_t globalStartIndex(Reactor* curr_reactor) {
+        for (size_t i = 0; i < m_reactors.size(); i++) {
+            if (curr_reactor == m_reactors[i]) {
+                return m_start[i];
+            }
+        }
+        throw CanteraError("ReactorNet::globalStartIndex: ",
+                curr_reactor->name(), " not found in network.");
+        return npos;
+    }
+
     //! Return the name of the i-th component of the global state vector. The
     //! name returned includes both the name of the reactor and the specific
     //! component, for example `'reactor1: CH4'`.
@@ -396,6 +410,10 @@ class ReactorNet : public FuncEval
     //! "left hand side" of each governing equation
     vector<double> m_LHS;
     vector<double> m_RHS;
+
+    //! derivative settings
+    bool m_jac_skip_walls = false;
+    bool m_jac_skip_flow_devices = false;
 };
 }
 

include/cantera/zeroD/Wall.h

@@ -8,6 +8,7 @@
 
 #include "cantera/base/ctexceptions.h"
 #include "cantera/zeroD/ReactorBase.h"
+#include "cantera/numerics/eigen_sparse.h"
 
 namespace Cantera
 {
@@ -115,6 +116,42 @@ class WallBase
         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.0.
+    */
+    virtual void buildReactorJacobian(ReactorBase* r, vector<Eigen::Triplet<double>>& jacVector) {
+        throw NotImplementedError("WallBase::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.0.
+    */
+    virtual void buildNetworkJacobian(vector<Eigen::Triplet<double>>& jacVector) {
+        throw NotImplementedError("WallBase::buildNetworkJacobian");
+    }
+
+    /*! Specify the jacobian terms have been calculated and should not be recalculated.
+     * @warning This function is an experimental part of the %Cantera API and may be changed
+     * or removed without notice.
+     * @since New in %Cantera 3.0.
+    */
+    void jacobianCalculated() { m_jac_calculated = true; };
+
+    /*! Specify that jacobian terms have not been calculated and should be recalculated.
+     * @warning This function is an experimental part of the %Cantera API and may be changed
+     * or removed without notice.
+     * @since New in %Cantera 3.0.
+    */
+    void jacobianNotCalculated() { m_jac_calculated = false; };
+
 protected:
     ReactorBase* m_left = nullptr;
     ReactorBase* m_right = nullptr;
@@ -123,6 +160,10 @@ class WallBase
     double m_time = 0.0;
 
     double m_area = 1.0;
+
+    //! a variable to switch on and off so calculations are not doubled by the calling
+    //! reactor or network
+    bool m_jac_calculated = false;
 };
 
 //! Represents a wall between between two ReactorBase objects.
@@ -255,6 +296,10 @@ class Wall : public WallBase
         return m_k;
     }
 
+    virtual void buildReactorJacobian(ReactorBase* r, vector<Eigen::Triplet<double>>& jacVector) override;
+
+    virtual void buildNetworkJacobian(vector<Eigen::Triplet<double>>& jacVector) override;
+
 protected:
 
     //! expansion rate coefficient

interfaces/cython/cantera/reactor.pxd

@@ -202,6 +202,8 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
         void setPreconditioner(shared_ptr[CxxPreconditionerBase] preconditioner)
         void setDerivativeSettings(CxxAnyMap&)
         CxxAnyMap solverStats() except +translate_exception
+        CxxSparseMatrix jacobian() except +translate_exception
+        CxxSparseMatrix finiteDifferenceJacobian() except +translate_exception
 
 cdef extern from "cantera/zeroD/ReactorDelegator.h" namespace "Cantera":
     cdef cppclass CxxReactorAccessor "Cantera::ReactorAccessor":

interfaces/cython/cantera/reactor.pyx

@@ -1797,3 +1797,29 @@ cdef class ReactorNet:
         """
         def __set__(self, settings):
             self.net.setDerivativeSettings(py_to_anymap(settings))
+
+    property jacobian:
+        """
+        Get the system Jacobian or an approximation thereof.
+
+        **Warning**: Depending on the particular implementation, this may return an
+        approximate Jacobian intended only for use in forming a preconditioner for
+        iterative solvers, excluding terms that would generate a fully-dense Jacobian.
+
+        **Warning**: This method is an experimental part of the Cantera API and may be
+        changed or removed without notice.
+        """
+        def __get__(self):
+            return get_from_sparse(self.net.jacobian(), self.n_vars, self.n_vars)
+
+    property finite_difference_jacobian:
+        """
+        Get the system Jacobian, calculated using a finite difference method.
+
+        **Warning:** this property is an experimental part of the Cantera API and
+        may be changed or removed without notice.
+        """
+        def __get__(self):
+            return get_from_sparse(self.net.finiteDifferenceJacobian(),
+                                   self.n_vars, self.n_vars)
+

src/zeroD/IdealGasMoleReactor.cpp

@@ -244,6 +244,31 @@ void IdealGasMoleReactor::buildJacobian(vector<Eigen::Triplet<double>>& jacVecto
                 (specificHeat[j] * qdot - NCv * uk_dnkdnj_sums[j]) * denom);
         }
     }
+
+    // build wall jacobian
+    buildWallJacobian(jacVector);
+
+    // build flow jacobian
+    buildFlowJacobian(jacVector);
+}
+
+double IdealGasMoleReactor::moleDerivative(size_t index)
+{
+    // derivative of temperature transformed by ideal gas law
+    vector<double> moles(m_nsp);
+    getMoles(moles.data());
+    double dTdni = pressure() * m_vol / GasConstant / std::accumulate(moles.begin(), moles.end(), 0);
+    return dTdni;
+}
+
+double IdealGasMoleReactor::moleRadiationDerivative(size_t index)
+{
+    // derivative of temperature transformed by ideal gas law
+    vector<double> moles(m_nsp);
+    getMoles(moles.data());
+    double dT4dni = std::pow(pressure() * m_vol / GasConstant, 4);
+    dT4dni *= std::pow(1 / std::accumulate(moles.begin(), moles.end(), 0), 5);
+    return dT4dni;
 }
 
 }