Jacobian utilities for preconditioner and networks with a test

include/cantera/zeroD/ReactorNet.h

@@ -283,6 +283,24 @@ class ReactorNet : public FuncEval
     //! @param settings the settings map propagated to all reactors and kinetics objects
     virtual void setDerivativeSettings(AnyMap& settings);
 
+    //! Calculate the system Jacobian using a finite difference method.
+    //!
+    //! This method is used only for informational purposes. Jacobian calculations
+    //! for the full reactor system are handled internally by CVODES.
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    Eigen::SparseMatrix<double> finiteDifferenceJacobian();
+
+    //! Method to calculate the system jacobian
+    //! @warning Depending on the particular implementation, this may return an
+    //! approximate Jacobian intended only for use in forming a preconditioner for
+    //! iterative solvers.
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    virtual Eigen::SparseMatrix<double> jacobian();
+
 protected:
     //! Check if surfaces and preconditioning are included, if so throw an error because
     //! they are currently not supported.

interfaces/cython/cantera/preconditioners.pxd

@@ -5,6 +5,7 @@
 #distutils: language=c++
 
 from .ctcxx cimport *
+from .kinetics cimport CxxSparseMatrix
 
 cdef extern from "cantera/numerics/PreconditionerBase.h" namespace "Cantera":
     cdef cppclass CxxPreconditionerBase "Cantera::PreconditionerBase":
@@ -23,6 +24,7 @@ cdef extern from "cantera/numerics/AdaptivePreconditioner.h" namespace "Cantera"
         void setIlutDropTol(double droptol)
         double ilutDropTol()
         void printPreconditioner()
+        CxxSparseMatrix matrix() except +
 
 cdef extern from "cantera/numerics/PreconditionerFactory.h" namespace "Cantera":
     cdef shared_ptr[CxxPreconditionerBase] newPreconditioner(string) except\

interfaces/cython/cantera/preconditioners.pyx

@@ -2,6 +2,7 @@
 # at https://cantera.org/license.txt for license and copyright information.
 
 from ._utils cimport stringify, pystr, dict_to_anymap, anymap_to_dict
+from .kinetics cimport get_from_sparse
 
 cdef class PreconditionerBase:
     """
@@ -93,3 +94,8 @@ cdef class AdaptivePreconditioner(PreconditionerBase):
 
     def print_contents(self):
         self.preconditioner.printPreconditioner()
+
+    property matrix:
+        def __get__(self):
+            cdef CxxSparseMatrix smat = self.preconditioner.matrix()
+            return get_from_sparse(smat, smat.rows(), smat.cols())

interfaces/cython/cantera/reactor.pxd

@@ -165,6 +165,8 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
         void setPreconditioner(shared_ptr[CxxPreconditionerBase] preconditioner)
         void setDerivativeSettings(CxxAnyMap&)
         CxxAnyMap solverStats()
+        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

@@ -1602,3 +1602,28 @@ cdef class ReactorNet:
         """
         def __set__(self, settings):
             self.net.setDerivativeSettings(dict_to_anymap(settings))
+
+    property jacobian:
+        """
+        Get the 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/ReactorNet.cpp

@@ -472,4 +472,54 @@ void ReactorNet::updatePreconditioner(double gamma)
     precon->updatePreconditioner();
 }
 
+Eigen::SparseMatrix<double> ReactorNet::jacobian()
+{
+    // network must be initialized for the jacobian
+    if (! m_init) {
+        initialize();
+    }
+    // create system jacobian from reactor jacobians
+    std::vector<Eigen::Triplet<double>> jac_trips;
+    Eigen::SparseMatrix<double> system_jac(m_nv, m_nv);
+    for (size_t n = 0; n < m_reactors.size(); n++) {
+        Reactor& r = *m_reactors[n];
+        Eigen::SparseMatrix<double> r_jac;
+        r_jac = r.jacobian();
+        size_t st = m_start[n];
+        for (int k=0; k<r_jac.outerSize(); ++k) {
+            for (Eigen::SparseMatrix<double>::InnerIterator it(r_jac, k); it; ++it) {
+                jac_trips.emplace_back(it.row() + st, it.col() + st, it.value());
+            }
+        }
+    }
+    // set jacobian from triples
+    system_jac.setFromTriplets(jac_trips.begin(), jac_trips.end());
+    return system_jac;
+}
+
+Eigen::SparseMatrix<double> ReactorNet::finiteDifferenceJacobian()
+{
+    // network must be initialized for the jacobian
+    if (! m_init) {
+        initialize();
+    }
+    // create system jacobian from reactor jacobians
+    std::vector<Eigen::Triplet<double>> jac_trips;
+    Eigen::SparseMatrix<double> system_fd_jac(m_nv, m_nv);
+    for (size_t n = 0; n < m_reactors.size(); n++) {
+        Reactor& r = *m_reactors[n];
+        Eigen::SparseMatrix<double> r_fd_jac;
+        r_fd_jac = r.finiteDifferenceJacobian();
+        size_t st = m_start[n];
+        for (int k=0; k<r_fd_jac.outerSize(); ++k) {
+            for (Eigen::SparseMatrix<double>::InnerIterator it(r_fd_jac, k); it; ++it) {
+                jac_trips.emplace_back(it.row() + st, it.col() + st, it.value());
+            }
+        }
+    }
+    // set jacobian from triples
+    system_fd_jac.setFromTriplets(jac_trips.begin(), jac_trips.end());
+    return system_fd_jac;
+}
+
 }

test/python/test_reactor.py

@@ -1185,6 +1185,42 @@ def test_adaptive_precon_integration(self):
             self.assertNear(r1.thermo.P, r2.thermo.P, rtol=1e-5)
 
 
+class TestReactorJacobians(utilities.CanteraTest):
+
+    def make_reactor(self, rtype=ct.IdealGasMoleReactor, model="ptcombust.yaml",
+                     gas_phase="gas", interface_phase="Pt_surf", add_surf=True):
+        # create gas and reactor
+        gas = ct.Solution(model, gas_phase)
+        r1 = rtype(gas)
+        # create surface
+        if add_surf:
+            surf = ct.Interface(model, interface_phase, [gas])
+            rsurf = ct.ReactorSurface(surf)
+            rsurf.install(r1)
+        self.net = ct.ReactorNet([r1])
+
+    def test_network_jacobian(self):
+        # create reactors
+        gas = ct.Solution("ptcombust.yaml", "gas")
+        r1 = ct.IdealGasMoleReactor(gas)
+        r2 = ct.IdealGasConstPressureMoleReactor(gas)
+        # create surface
+        surf = ct.Interface("ptcombust.yaml", "Pt_surf", [gas])
+        rsurf1 = ct.ReactorSurface(surf)
+        rsurf2 = ct.ReactorSurface(surf)
+        # install surfaces
+        rsurf1.install(r1)
+        rsurf2.install(r2)
+        # create network
+        net = ct.ReactorNet([r1, r2])
+        # net.advance(0.1)
+        net_jac = net.jacobian
+        r1_jac = r1.jacobian
+        r2_jac = r2.jacobian
+        self.assertArrayNear(r1_jac, net_jac[:r1.n_vars, :r1.n_vars], 1e-6, 1e-12)
+        self.assertArrayNear(r2_jac, net_jac[r1.n_vars:, r1.n_vars:], 1e-6, 1e-12)
+
+
 class TestFlowReactor(utilities.CanteraTest):
     gas_def = """
     phases: