[Numerics] Create intermediate class for sparse Jacobians

include/cantera/numerics/AdaptivePreconditioner.h

@@ -10,10 +10,7 @@
 #ifndef ADAPTIVEPRECONDITIONER_H
 #define ADAPTIVEPRECONDITIONER_H
 
-#include "cantera/numerics/SystemJacobian.h"
-#include "cantera/numerics/eigen_sparse.h"
-#include "cantera/base/global.h"
-#include <iostream>
+#include "cantera/numerics/EigenSparseJacobian.h"
 
 namespace Cantera
 {
@@ -23,7 +20,7 @@ namespace Cantera
 //! the preconditioner by a threshold value. It also neglects pressure
 //! dependence and third body contributions in its formation and has a
 //! finite difference approximation for temperature.
-class AdaptivePreconditioner : public SystemJacobian
+class AdaptivePreconditioner : public EigenSparseJacobian
 {
 public:
     AdaptivePreconditioner();
@@ -31,7 +28,7 @@ public:
     void initialize(size_t networkSize) override;
 
     void reset() override {
-        m_precon_matrix.setZero();
+        m_matrix.setZero();
         m_jac_trips.clear();
     };
 
@@ -39,10 +36,7 @@ public:
     void setup() override; // deprecated
     void factorize() override;
     void solve(const size_t stateSize, double* rhs_vector, double* output) override;
-    void setValue(size_t row, size_t col, double value) override;
     void stateAdjustment(vector<double>& state) override;
-    void updatePreconditioner() override;
-    void updateTransient(double rdt, int* mask) override;
 
     int info() const override {
         return static_cast<int>(m_solver.info());
@@ -51,20 +45,6 @@ public:
     //! Prune preconditioner elements
     void prunePreconditioner();
 
-    //! Return semi-analytical Jacobian of an AdaptivePreconditioner object.
-    //! @ingroup derivGroup
-    Eigen::SparseMatrix<double> jacobian() {
-        Eigen::SparseMatrix<double> jacobian_mat(m_dim, m_dim);
-        jacobian_mat.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
-        return jacobian_mat;
-    }
-
-    //! Return the internal preconditioner matrix
-    Eigen::SparseMatrix<double> matrix() {
-        updatePreconditioner();
-        return m_precon_matrix;
-    }
-
     //! Get the threshold value for setting elements
     double threshold() { return m_threshold; }
 
@@ -95,12 +75,6 @@ public:
         m_solver.setFillfactor(fillFactor);
     }
 
-    //! Print preconditioner contents
-    void printPreconditioner() override;
-
-    //! Print jacobian contents
-    void printJacobian();
-
 protected:
     //! ILUT fill factor
     double m_fill_factor = 0;
@@ -108,15 +82,6 @@ protected:
     //! ILUT drop tolerance
     double m_drop_tol = 0;
 
-    //! Vector of triples representing the jacobian used in preconditioning
-    vector<Eigen::Triplet<double>> m_jac_trips;
-
-    //! Storage of appropriately sized identity matrix for making the preconditioner
-    Eigen::SparseMatrix<double> m_identity;
-
-    //! Container that is the sparse preconditioner
-    Eigen::SparseMatrix<double> m_precon_matrix;
-
     //! Solver used in solving the linear system
     Eigen::IncompleteLUT<double> m_solver;
 

include/cantera/numerics/EigenSparseJacobian.h

@@ -0,0 +1,54 @@
+//! @file EigenSparseJacobian.h
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#ifndef EIGENSPARSEJACOBIAN_H
+#define EIGENSPARSEJACOBIAN_H
+
+#include "cantera/numerics/SystemJacobian.h"
+#include "cantera/numerics/eigen_sparse.h"
+
+namespace Cantera
+{
+
+//! System Jacobians that use Eigen sparse matrices for storage
+class EigenSparseJacobian : public SystemJacobian
+{
+public:
+    EigenSparseJacobian() = default;
+    void initialize(size_t networkSize) override;
+    void setValue(size_t row, size_t col, double value) override;
+    void updatePreconditioner() override;
+    void updateTransient(double rdt, int* mask) override;
+
+    //! Return underlying Jacobian matrix
+    //! @ingroup derivGroup
+    Eigen::SparseMatrix<double> jacobian();
+
+    //! Return the internal preconditioner matrix
+    Eigen::SparseMatrix<double> matrix() {
+        updatePreconditioner();
+        return m_matrix;
+    }
+
+    //! Print preconditioner contents
+    void printPreconditioner() override;
+
+    //! Print jacobian contents
+    void printJacobian();
+
+protected:
+    //! Vector of triples representing the jacobian used in preconditioning
+    vector<Eigen::Triplet<double>> m_jac_trips;
+
+    //! Storage of appropriately sized identity matrix for making the preconditioner
+    Eigen::SparseMatrix<double> m_identity;
+
+    //! Container that is the sparse preconditioner
+    Eigen::SparseMatrix<double> m_matrix;
+};
+
+}
+
+#endif

interfaces/cython/cantera/_onedim.pyx

@@ -1023,7 +1023,7 @@ cdef class Sim1D:
 
     @linear_solver.setter
     def linear_solver(self, SystemJacobian precon):
-        self.sim.setLinearSolver(precon.pbase)
+        self.sim.setLinearSolver(precon._base)
 
     def set_initial_guess(self, *args, **kwargs):
         """

interfaces/cython/cantera/jacobians.pxd

@@ -14,9 +14,16 @@ cdef extern from "cantera/numerics/SystemJacobian.h" namespace "Cantera":
         string type()
         void setPreconditionerSide(string) except +translate_exception
 
+cdef extern from "cantera/numerics/EigenSparseJacobian.h" namespace "Cantera":
+    cdef cppclass CxxEigenSparseJacobian "Cantera::EigenSparseJacobian" \
+        (CxxSystemJacobian):
+        CxxEigenSparseJacobian() except +translate_exception
+        void printPreconditioner()
+        CxxSparseMatrix matrix() except +translate_exception
+
 cdef extern from "cantera/numerics/AdaptivePreconditioner.h" namespace "Cantera":
     cdef cppclass CxxAdaptivePreconditioner "Cantera::AdaptivePreconditioner" \
-        (CxxSystemJacobian):
+        (CxxEigenSparseJacobian):
         CxxAdaptivePreconditioner() except +translate_exception
         void setThreshold(double threshold)
         double threshold()
@@ -24,8 +31,6 @@ cdef extern from "cantera/numerics/AdaptivePreconditioner.h" namespace "Cantera"
         double ilutFillFactor()
         void setIlutDropTol(double droptol)
         double ilutDropTol()
-        void printPreconditioner()
-        CxxSparseMatrix matrix() except +translate_exception
 
 cdef extern from "cantera/oneD/MultiJac.h" namespace "Cantera":
     cdef cppclass CxxMultiJac "Cantera::MultiJac" (CxxSystemJacobian):
@@ -39,12 +44,16 @@ cdef class SystemJacobian:
     @staticmethod
     cdef wrap(shared_ptr[CxxSystemJacobian])
     cdef set_cxx_object(self)
-    cdef shared_ptr[CxxSystemJacobian] pbase
+    cdef shared_ptr[CxxSystemJacobian] _base
 
-cdef class AdaptivePreconditioner(SystemJacobian):
+cdef class EigenSparseJacobian(SystemJacobian):
     cdef set_cxx_object(self)
-    cdef CxxAdaptivePreconditioner* preconditioner
+    cdef CxxEigenSparseJacobian* sparse_jac
+
+cdef class AdaptivePreconditioner(EigenSparseJacobian):
+    cdef set_cxx_object(self)
+    cdef CxxAdaptivePreconditioner* adaptive
 
 cdef class BandedJacobian(SystemJacobian):
     cdef set_cxx_object(self)
-    cdef CxxMultiJac* jacobian
+    cdef CxxMultiJac* band_jac

interfaces/cython/cantera/jacobians.pyx

@@ -17,7 +17,11 @@ cdef class SystemJacobian:
 
     def __cinit__(self, *args, init=True, **kwargs):
         if init:
-            self.pbase = newSystemJacobian(stringify(self._type))
+            self._cinit(*args, **kwargs)
+
+    def _cinit(self, *args, **kwargs):
+        self._base = newSystemJacobian(stringify(self._type))
+        self.set_cxx_object()
 
     @staticmethod
     cdef wrap(shared_ptr[CxxSystemJacobian] base):
@@ -39,7 +43,7 @@ cdef class SystemJacobian:
         cls = _class_registry.get(cxx_type, SystemJacobian)
         cdef SystemJacobian jac
         jac = cls(init=False)
-        jac.pbase = base
+        jac._base = base
         jac.set_cxx_object()
         return jac
 
@@ -53,20 +57,36 @@ cdef class SystemJacobian:
         by all solver types.
         """
         def __get__(self):
-            return pystr(self.pbase.get().preconditionerSide())
+            return pystr(self._base.get().preconditionerSide())
 
         def __set__(self, side):
-            self.pbase.get().setPreconditionerSide(stringify(side))
+            self._base.get().setPreconditionerSide(stringify(side))
 
-cdef class AdaptivePreconditioner(SystemJacobian):
+
+cdef class EigenSparseJacobian(SystemJacobian):
+    _type = "EigenSparse"
+
+    cdef set_cxx_object(self):
+        self.sparse_jac = <CxxEigenSparseJacobian*>self._base.get()
+
+    def print_contents(self):
+        self.sparse_jac.printPreconditioner()
+
+    property matrix:
+        """
+        Property to retrieve the latest internal preconditioner matrix.
+        """
+        def __get__(self):
+            cdef CxxSparseMatrix smat = self.sparse_jac.matrix()
+            return get_from_sparse(smat, smat.rows(), smat.cols())
+
+
+cdef class AdaptivePreconditioner(EigenSparseJacobian):
     _type = "Adaptive"
     linear_solver_type = "GMRES"
 
-    def __cinit__(self, *args, **kwargs):
-        self.preconditioner = <CxxAdaptivePreconditioner*>(self.pbase.get())
-
     cdef set_cxx_object(self):
-        self.preconditioner = <CxxAdaptivePreconditioner*>self.pbase.get()
+        self.adaptive = <CxxAdaptivePreconditioner*>self._base.get()
 
     property threshold:
         """
@@ -84,10 +104,10 @@ cdef class AdaptivePreconditioner(SystemJacobian):
         Default is 0.0.
         """
         def __get__(self):
-            return self.preconditioner.threshold()
+            return self.adaptive.threshold()
 
         def __set__(self, val):
-            self.preconditioner.setThreshold(val)
+            self.adaptive.setThreshold(val)
 
     property ilut_fill_factor:
         """
@@ -104,10 +124,10 @@ cdef class AdaptivePreconditioner(SystemJacobian):
         Default is the state size divided by 4.
         """
         def __set__(self, val):
-            self.preconditioner.setIlutFillFactor(val)
+            self.adaptive.setIlutFillFactor(val)
 
         def __get__(self):
-            return self.preconditioner.ilutFillFactor()
+            return self.adaptive.ilutFillFactor()
 
     property ilut_drop_tol:
         """
@@ -122,21 +142,10 @@ cdef class AdaptivePreconditioner(SystemJacobian):
         Default is 1e-10.
         """
         def __set__(self, val):
-            self.preconditioner.setIlutDropTol(val)
+            self.adaptive.setIlutDropTol(val)
 
         def __get__(self):
-            return self.preconditioner.ilutDropTol()
-
-    def print_contents(self):
-        self.preconditioner.printPreconditioner()
-
-    property matrix:
-        """
-        Property to retrieve the latest internal preconditioner matrix.
-        """
-        def __get__(self):
-            cdef CxxSparseMatrix smat = self.preconditioner.matrix()
-            return get_from_sparse(smat, smat.rows(), smat.cols())
+            return self.adaptive.ilutDropTol()
 
 
 cdef class BandedJacobian(SystemJacobian):
@@ -144,4 +153,4 @@ cdef class BandedJacobian(SystemJacobian):
     linear_solver_type = "direct"
 
     cdef set_cxx_object(self):
-        self.jacobian = <CxxMultiJac*>self.pbase.get()
+        self.band_jac = <CxxMultiJac*>self._base.get()

interfaces/cython/cantera/reactor.pyx

@@ -1984,7 +1984,7 @@ cdef class ReactorNet:
         """Preconditioner associated with integrator"""
         def __set__(self, SystemJacobian precon):
             # set preconditioner
-            self.net.setPreconditioner(precon.pbase)
+            self.net.setPreconditioner(precon._base)
             # set problem type as default of preconditioner
             self.linear_solver_type = precon.linear_solver_type
 

src/numerics/AdaptivePreconditioner.cpp

@@ -14,11 +14,6 @@ AdaptivePreconditioner::AdaptivePreconditioner()
     setPreconditionerSide("right");
 }
 
-void AdaptivePreconditioner::setValue(size_t row, size_t col, double value)
-{
-    m_jac_trips.emplace_back(static_cast<int>(row), static_cast<int>(col), value);
-}
-
 void AdaptivePreconditioner::stateAdjustment(vector<double>& state) {
     // Only keep positive composition based on given tol
     for (size_t i = 0; i < state.size(); i++) {
@@ -28,20 +23,9 @@ void AdaptivePreconditioner::stateAdjustment(vector<double>& state) {
 
 void AdaptivePreconditioner::initialize(size_t networkSize)
 {
+    EigenSparseJacobian::initialize(networkSize);
     // don't use legacy rate constants
     use_legacy_rate_constants(false);
-    // reset arrays in case of re-initialization
-    m_jac_trips.clear();
-    // set dimensions of preconditioner from network
-    m_dim = networkSize;
-    // reserve some space for vectors making up SparseMatrix
-    m_jac_trips.reserve(3 * networkSize);
-    // reserve space for preconditioner
-    m_precon_matrix.resize(m_dim, m_dim);
-    // creating sparse identity matrix
-    m_identity.resize(m_dim, m_dim);
-    m_identity.setIdentity();
-    m_identity.makeCompressed();
     // setting default ILUT parameters
     if (m_drop_tol == 0) {
         setIlutDropTol(1e-12);
@@ -53,7 +37,6 @@ void AdaptivePreconditioner::initialize(size_t networkSize)
     m_init = true;
 }
 
-
 void AdaptivePreconditioner::setup()
 {
     warn_deprecated("AdaptivePreconditioner::setup",
@@ -62,39 +45,15 @@ void AdaptivePreconditioner::setup()
     factorize();
 }
 
-void AdaptivePreconditioner::updatePreconditioner()
+void AdaptivePreconditioner::factorize()
 {
-    // set precon to jacobian
-    m_precon_matrix.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
-    // make into preconditioner as P = (I - gamma * J_bar)
-    m_precon_matrix = m_identity - m_gamma * m_precon_matrix;
-    // prune by threshold if desired
     if (m_prune_precon) {
         prunePreconditioner();
     }
-    factorize();
-}
-
-void AdaptivePreconditioner::updateTransient(double rdt, int* mask)
-{
-    // set matrix to steady Jacobian
-    m_precon_matrix.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
-    // update transient diagonal terms
-    Eigen::VectorXd diag = Eigen::Map<Eigen::VectorXi>(mask, m_dim).cast<double>();
-    m_precon_matrix -= rdt * diag.matrix().asDiagonal();
-    // prune by threshold if desired
-    // if (m_prune_precon) {
-    //     prunePreconditioner();
-    // }
-    factorize();
-}
-
-void AdaptivePreconditioner::factorize()
-{
     // compress sparse matrix structure
-    m_precon_matrix.makeCompressed();
+    m_matrix.makeCompressed();
     // analyze and factorize
-    m_solver.compute(m_precon_matrix);
+    m_solver.compute(m_matrix);
     // check for errors
     if (m_solver.info() != Eigen::Success) {
         throw CanteraError("AdaptivePreconditioner::factorize",
@@ -104,8 +63,8 @@ void AdaptivePreconditioner::factorize()
 
 void AdaptivePreconditioner::prunePreconditioner()
 {
-    for (int k=0; k<m_precon_matrix.outerSize(); ++k) {
-        for (Eigen::SparseMatrix<double>::InnerIterator it(m_precon_matrix, k); it;
+    for (int k=0; k<m_matrix.outerSize(); ++k) {
+        for (Eigen::SparseMatrix<double>::InnerIterator it(m_matrix, k); it;
             ++it) {
             if (std::abs(it.value()) < m_threshold && it.row() != it.col()) {
                 it.valueRef() = 0;
@@ -128,19 +87,4 @@ void AdaptivePreconditioner::solve(const size_t stateSize, double* rhs_vector, d
     }
 }
 
-void AdaptivePreconditioner::printPreconditioner() {
-    std::stringstream ss;
-    Eigen::IOFormat HeavyFmt(Eigen::FullPrecision, 0, ", ", ";\n", "[", "]", "[", "]");
-    ss << Eigen::MatrixXd(m_precon_matrix).format(HeavyFmt);
-    writelog(ss.str());
-}
-
-void AdaptivePreconditioner::printJacobian() {
-    std::stringstream ss;
-    Eigen::SparseMatrix<double> jacobian(m_dim, m_dim);
-    jacobian.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
-    ss << Eigen::MatrixXd(jacobian);
-    writelog(ss.str());
-}
-
 }

src/numerics/EigenSparseJacobian.cpp

@@ -0,0 +1,77 @@
+//! @file EigenSparseJacobian.cpp
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#include "cantera/numerics/EigenSparseJacobian.h"
+#include "cantera/base/global.h"
+
+namespace Cantera
+{
+
+void EigenSparseJacobian::initialize(size_t networkSize)
+{
+    // reset arrays in case of re-initialization
+    m_jac_trips.clear();
+    // set dimensions of preconditioner from network
+    m_dim = networkSize;
+    // reserve some space for vectors making up SparseMatrix
+    m_jac_trips.reserve(3 * networkSize);
+    // reserve space for preconditioner
+    m_matrix.resize(m_dim, m_dim);
+    // creating sparse identity matrix
+    m_identity.resize(m_dim, m_dim);
+    m_identity.setIdentity();
+    m_identity.makeCompressed();
+    // update initialized status
+    m_init = true;
+}
+
+void EigenSparseJacobian::setValue(size_t row, size_t col, double value)
+{
+    m_jac_trips.emplace_back(static_cast<int>(row), static_cast<int>(col), value);
+}
+
+void EigenSparseJacobian::updatePreconditioner()
+{
+    // set precon to jacobian
+    m_matrix.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
+    // make into preconditioner as P = (I - gamma * J_bar)
+    m_matrix = m_identity - m_gamma * m_matrix;
+    // prune by threshold if desired
+    factorize();
+}
+
+void EigenSparseJacobian::updateTransient(double rdt, int* mask)
+{
+    // set matrix to steady Jacobian
+    m_matrix.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
+    // update transient diagonal terms
+    Eigen::VectorXd diag = Eigen::Map<Eigen::VectorXi>(mask, m_dim).cast<double>();
+    m_matrix -= rdt * diag.matrix().asDiagonal();
+    factorize();
+}
+
+Eigen::SparseMatrix<double> EigenSparseJacobian::jacobian()
+{
+    Eigen::SparseMatrix<double> jacobian_mat(m_dim, m_dim);
+    jacobian_mat.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
+    return jacobian_mat;
+}
+
+void EigenSparseJacobian::printPreconditioner() {
+    std::stringstream ss;
+    Eigen::IOFormat HeavyFmt(Eigen::FullPrecision, 0, ", ", ";\n", "[", "]", "[", "]");
+    ss << Eigen::MatrixXd(m_matrix).format(HeavyFmt);
+    writelog(ss.str());
+}
+
+void EigenSparseJacobian::printJacobian() {
+    std::stringstream ss;
+    Eigen::SparseMatrix<double> jacobian(m_dim, m_dim);
+    jacobian.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
+    ss << Eigen::MatrixXd(jacobian);
+    writelog(ss.str());
+}
+
+}