Enable use of different Jacobian implementations with 1D solver

include/cantera/numerics/AdaptivePreconditioner.h

@@ -1,6 +1,6 @@
 /**
  *  @file AdaptivePreconditioner.h Declarations for the class
- *   AdaptivePreconditioner which is a child class of PreconditionerBase
+ *   AdaptivePreconditioner which is a child class of SystemJacobian
  *   for preconditioners used by sundials
  */
 
@@ -10,10 +10,7 @@
 #ifndef ADAPTIVEPRECONDITIONER_H
 #define ADAPTIVEPRECONDITIONER_H
 
-#include "cantera/numerics/PreconditionerBase.h"
-#include "cantera/numerics/eigen_sparse.h"
-#include "cantera/base/global.h"
-#include <iostream>
+#include "cantera/numerics/EigenSparseJacobian.h"
 
 namespace Cantera
 {
@@ -23,45 +20,25 @@ 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 PreconditionerBase
+class AdaptivePreconditioner : public EigenSparseJacobian
 {
 public:
     AdaptivePreconditioner();
 
     void initialize(size_t networkSize) override;
-
-    void reset() override {
-        m_precon_matrix.setZero();
-        m_jac_trips.clear();
-    };
-
-    void setup() override;
-
+    const string type() const override { return "Adaptive"; }
+    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;
+    int info() const override {
+        return static_cast<int>(m_solver.info());
+    }
 
     //! 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; }
 
@@ -92,28 +69,13 @@ class AdaptivePreconditioner : public PreconditionerBase
         m_solver.setFillfactor(fillFactor);
     }
 
-    //! Print preconditioner contents
-    void printPreconditioner() override;
-
-    //! Print jacobian contents
-    void printJacobian();
-
 protected:
     //! ILUT fill factor
     double m_fill_factor = 0;
 
     //! 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/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

@@ -0,0 +1,55 @@
+//! @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 reset() 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

include/cantera/numerics/Integrator.h

@@ -8,7 +8,7 @@
 #ifndef CT_INTEGRATOR_H
 #define CT_INTEGRATOR_H
 #include "FuncEval.h"
-#include "PreconditionerBase.h"
+#include "SystemJacobian.h"
 
 #include "cantera/base/global.h"
 #include "cantera/base/AnyMap.h"
@@ -91,7 +91,7 @@ class Integrator
     /*!
      * @param preconditioner preconditioner object used for the linear solver
      */
-    virtual void setPreconditioner(shared_ptr<PreconditionerBase> preconditioner) {
+    virtual void setPreconditioner(shared_ptr<SystemJacobian> preconditioner) {
         m_preconditioner = preconditioner;
         if (preconditioner->preconditionerSide() == "none") {
             m_prec_side = PreconditionerSide::NO_PRECONDITION;
@@ -123,7 +123,7 @@ class Integrator
     }
 
     //! Return preconditioner reference to object
-    virtual shared_ptr<PreconditionerBase> preconditioner() {
+    virtual shared_ptr<SystemJacobian> preconditioner() {
         return m_preconditioner;
     }
 
@@ -296,7 +296,7 @@ class Integrator
     //! Pointer to preconditioner object used in integration which is
     //! set by setPreconditioner and initialized inside of
     //! ReactorNet::initialize()
-    shared_ptr<PreconditionerBase> m_preconditioner;
+    shared_ptr<SystemJacobian> m_preconditioner;
     //! Type of preconditioning used in applyOptions
     PreconditionerSide m_prec_side = PreconditionerSide::NO_PRECONDITION;
     // methods for DAE solvers

include/cantera/numerics/PreconditionerBase.h

@@ -1,132 +1,15 @@
-/**
- *  @file PreconditionerBase.h Declarations for the class
- *   PreconditionerBase which is a virtual base class for
- *   preconditioning systems.
- */
+//! @file PreconditionerBase.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 PRECONDITIONERBASE_H
-#define PRECONDITIONERBASE_H
+#ifndef PRECONDITIONERBASE_FACTORY_H
+#define PRECONDITIONERBASE_FACTORY_H
 
-#include "cantera/base/ctexceptions.h"
+#pragma message("warning: PreconditionerBase.h and class PreconditionerBase are " \
+                "deprecated and will be removed after Cantera 3.2. Use " \
+                "SystemJacobian.h and class SystemJacobian.")
 
-namespace Cantera
-{
+#include "SystemJacobian.h"
 
-/**
- * Specifies the side of the system on which the preconditioner is applied. Not all
- * methods are supported by all integrators.
- */
-enum class PreconditionerSide {
-    NO_PRECONDITION, //! No preconditioning
-    LEFT_PRECONDITION, //! Left side preconditioning
-    RIGHT_PRECONDITION, //! Right side preconditioning
-    BOTH_PRECONDITION //! Left and right side preconditioning
-};
-
-//! PreconditionerBase serves as an abstract type to extend different preconditioners
-class PreconditionerBase
-{
-public:
-    PreconditionerBase() {}
-
-    virtual ~PreconditionerBase() {}
-
-    //! Set a value at the specified row and column of the jacobian triplet vector
-    //! @param row row in the jacobian matrix
-    //! @param col column in the jacobian matrix
-    //! @param value value of the element at row and col
-    virtual void setValue(size_t row, size_t col, double value) {
-        throw NotImplementedError("PreconditionerBase::setValue");
-    }
-
-    //! Adjust the state vector based on the preconditioner, e.g., Adaptive
-    //! preconditioning uses a strictly positive composition when preconditioning which
-    //! is handled by this function
-    //! @param state a vector containing the state to be updated
-    virtual void stateAdjustment(vector<double>& state) {
-        throw NotImplementedError("PreconditionerBase::stateAdjustment");
-    }
-
-    //! Get preconditioner application side for CVODES
-    string preconditionerSide() const {
-        return m_precon_side;
-    }
-
-    virtual void setPreconditionerSide(const string& preconSide) {
-        m_precon_side = preconSide;
-    }
-
-    //! Solve a linear system Ax=b where A is the preconditioner
-    //! @param[in] stateSize length of the rhs and output vectors
-    //! @param[in] rhs_vector right hand side vector used in linear system
-    //! @param[out] output output vector for solution
-    virtual void solve(const size_t stateSize, double* rhs_vector, double* output) {
-        throw NotImplementedError("PreconditionerBase::solve");
-    };
-
-    //! Perform preconditioner specific post-reactor setup operations such as factorize.
-    virtual void setup() {
-        throw NotImplementedError("PreconditionerBase::setup");
-    };
-
-    //! Reset preconditioner parameters as needed
-    virtual void reset() {
-        throw NotImplementedError("PreconditionerBase::reset");
-    };
-
-    //! Called during setup for any processes that need
-    //! to be completed prior to setup functions used in sundials.
-    //! @param networkSize the number of variables in the associated reactor network
-    virtual void initialize(size_t networkSize) {
-        throw NotImplementedError("PreconditionerBase::initialize");
-    };
-
-    //! Print preconditioner contents
-    virtual void printPreconditioner() {
-        throw NotImplementedError("PreconditionerBase::printPreconditioner");
-    };
-
-    //! Transform Jacobian vector and write into
-    //! preconditioner, P = (I - gamma * J)
-    virtual void updatePreconditioner() {
-        throw NotImplementedError("PreconditionerBase::updatePreconditioner");
-    }
-
-    //! Set gamma used in preconditioning
-    //! @param gamma used in M = I - gamma*J
-    virtual void setGamma(double gamma) {
-        m_gamma = gamma;
-    };
-
-    //! Get gamma used in preconditioning
-    virtual double gamma() {
-        return m_gamma;
-    };
-
-    //! Set the absolute tolerance in the solver outside of the network initialization
-    //! @param atol the specified tolerance
-    virtual void setAbsoluteTolerance(double atol) {
-        m_atol = atol;
-    }
-
-protected:
-    //! Dimension of the preconditioner
-    size_t m_dim;
-
-    //! gamma value used in M = I - gamma*J
-    double m_gamma = 1.0;
-
-    //! bool saying whether or not the preconditioner is initialized
-    bool m_init = false;
-
-    //! Absolute tolerance of the ODE solver
-    double m_atol = 0;
-
-    string m_precon_side = "none";
-};
-
-}
 #endif