[Numerics] Add EigenSparseDirectJacobian implementation

include/cantera/numerics/AdaptivePreconditioner.h

@@ -26,12 +26,6 @@ public:
     AdaptivePreconditioner();
 
     void initialize(size_t networkSize) override;
-
-    void reset() override {
-        m_matrix.setZero();
-        m_jac_trips.clear();
-    };
-
     const string type() const override { return "Adaptive"; }
     void setup() override; // deprecated
     void factorize() override;

include/cantera/numerics/EigenSparseDirectJacobian.h

@@ -0,0 +1,28 @@
+//! @file EigenSparseDirectJacobian.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 EIGENSPARSEDIRECTJACOBIAN_H
+#define EIGENSPARSEDIRECTJACOBIAN_H
+
+#include "cantera/numerics/EigenSparseJacobian.h"
+
+namespace Cantera
+{
+
+class EigenSparseDirectJacobian : public EigenSparseJacobian
+{
+public:
+    EigenSparseDirectJacobian() = default;
+    const string type() const override { return "eigen-sparse-direct"; }
+    void factorize() override;
+    void solve(const size_t stateSize, double* rhs_vector, double* output) override;
+
+protected:
+    Eigen::SparseLU<Eigen::SparseMatrix<double>> m_solver;
+};
+
+}
+
+#endif

include/cantera/numerics/EigenSparseJacobian.h

@@ -18,6 +18,7 @@ class EigenSparseJacobian : public SystemJacobian
 public:
     EigenSparseJacobian() = default;
     void initialize(size_t networkSize) override;
+    void reset() override;
     void setValue(size_t row, size_t col, double value) override;
     void updatePreconditioner() override;
     void updateTransient(double rdt, int* mask) override;

interfaces/cython/cantera/jacobians.pxd

@@ -50,6 +50,9 @@ cdef class EigenSparseJacobian(SystemJacobian):
     cdef set_cxx_object(self)
     cdef CxxEigenSparseJacobian* sparse_jac
 
+cdef class EigenSparseDirectJacobian(EigenSparseJacobian):
+    pass
+
 cdef class AdaptivePreconditioner(EigenSparseJacobian):
     cdef set_cxx_object(self)
     cdef CxxAdaptivePreconditioner* adaptive

interfaces/cython/cantera/jacobians.pyx

@@ -64,7 +64,7 @@ cdef class SystemJacobian:
 
 
 cdef class EigenSparseJacobian(SystemJacobian):
-    _type = "EigenSparse"
+    _type = "eigen-sparse"
 
     cdef set_cxx_object(self):
         self.sparse_jac = <CxxEigenSparseJacobian*>self._base.get()
@@ -81,6 +81,15 @@ cdef class EigenSparseJacobian(SystemJacobian):
             return get_from_sparse(smat, smat.rows(), smat.cols())
 
 
+cdef class EigenSparseDirectJacobian(EigenSparseJacobian):
+    """
+    A system matrix solver that uses Eigen's sparse direct (LU) algorithm. Wraps C++
+    class :ct:`EigenSparseDirectJacobian`.
+    """
+
+    _type = "eigen-sparse-direct"
+
+
 cdef class AdaptivePreconditioner(EigenSparseJacobian):
     _type = "Adaptive"
     linear_solver_type = "GMRES"

src/numerics/EigenSparseDirectJacobian.cpp

@@ -0,0 +1,34 @@
+//! @file EigenSparseDirectJacobian.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/EigenSparseDirectJacobian.h"
+#include "cantera/numerics/eigen_dense.h"
+
+namespace Cantera
+{
+
+void EigenSparseDirectJacobian::factorize()
+{
+    m_matrix.makeCompressed();
+    // analyze and factorize
+    m_solver.compute(m_matrix);
+    // check for errors
+    if (m_solver.info() != Eigen::Success) {
+        throw CanteraError("EigenSparseDirectJacobian::factorize",
+                           "error code: {}", static_cast<int>(m_solver.info()));
+    }
+}
+
+void EigenSparseDirectJacobian::solve(const size_t stateSize, double* b, double* x)
+{
+    MappedVector(x, m_dim) = m_solver.solve(MappedVector(b, m_dim));
+    // check for errors
+    if (m_solver.info() != Eigen::Success) {
+        throw CanteraError("EigenSparseDirectJacobian::solve",
+                           "error code: {}", static_cast<int>(m_solver.info()));
+    }
+}
+
+}

src/numerics/EigenSparseJacobian.cpp

@@ -27,6 +27,12 @@ void EigenSparseJacobian::initialize(size_t networkSize)
     m_init = true;
 }
 
+void EigenSparseJacobian::reset()
+{
+    m_matrix.setZero();
+    m_jac_trips.clear();
+}
+
 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);

src/numerics/SystemJacobianFactory.cpp

@@ -5,6 +5,7 @@
 
 #include "cantera/numerics/SystemJacobianFactory.h"
 #include "cantera/numerics/AdaptivePreconditioner.h"
+#include "cantera/numerics/EigenSparseDirectJacobian.h"
 #include "cantera/oneD/MultiJac.h"
 
 namespace Cantera
@@ -31,6 +32,7 @@ SystemJacobianFactory::SystemJacobianFactory()
 {
     reg("Adaptive", []() { return new AdaptivePreconditioner(); });
     reg("banded-direct", []() { return new MultiJac(); });
+    reg("eigen-sparse-direct", []() { return new EigenSparseDirectJacobian(); });
 }
 
 shared_ptr<SystemJacobian> newSystemJacobian(const string& type)