Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Implement adaptive preconditioning for reactor networks
Adds PreconditionerBase and derives AdaptivePreconditioner from this base class. PreconditionerBase is intended to be an abstract class so other types of preconditioners can be extended. AdaptivePreconditioner is a specific type of preconditioning based on a tolerance. The preconditioned integrator works with all Sundials versions supported by Cantera. The ILUT algorithm is used to limit fill-in when factorizing the preconditioner. Adds access to CVODES integrator statistics in C++ and Python Adds reactor classes "IdealGasConstPressureMoleReactor" and "IdealGasMoleReactor" that use a mole-based state vector to increase the sparsity of the preconditioner.
- Loading branch information
1 parent
d60e59e
commit c361a5e
Showing
34 changed files
with
2,052 additions
and
37 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,138 @@ | ||
| /** | ||
| * @file AdaptivePreconditioner.h Declarations for the class | ||
| * AdaptivePreconditioner which is a child class of PreconditionerBase | ||
| * for preconditioners used by sundials | ||
| */ | ||
|
|
||
| // 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 ADAPTIVEPRECONDITIONER_H | ||
| #define ADAPTIVEPRECONDITIONER_H | ||
|
|
||
| #include "cantera/numerics/PreconditionerBase.h" | ||
| #include "cantera/numerics/eigen_sparse.h" | ||
| #include "cantera/base/global.h" | ||
| #include <iostream> | ||
|
|
||
| namespace Cantera | ||
| { | ||
|
|
||
| //! AdaptivePreconditioner a preconditioner designed for use with large | ||
| //! mechanisms that leverages sparse solvers. It does this by pruning | ||
| //! the preconditioner by a threshold value. It also neglects pressure | ||
| //! dependence and thirdbody contributions in its formation and has a | ||
| //! finite difference approximation for temperature. | ||
| class AdaptivePreconditioner : public PreconditionerBase | ||
| { | ||
| public: | ||
| AdaptivePreconditioner() {} | ||
|
|
||
| void initialize(size_t networkSize); | ||
|
|
||
| void reset() { | ||
| m_precon_matrix.setZero(); | ||
| m_jac_trips.clear(); | ||
| }; | ||
|
|
||
| void setup(); | ||
|
|
||
| void solve(const size_t stateSize, double *rhs_vector, double* output); | ||
|
|
||
| PreconditionerType preconditionerType() { return PreconditionerType::LEFT_PRECONDITION; } | ||
|
|
||
| void setValue(size_t row, size_t col, double value); | ||
|
|
||
| virtual void stateAdjustment(vector_fp& state); | ||
|
|
||
| //! Transform Jacobian vector and write into | ||
| //! preconditioner | ||
| void transformJacobianToPreconditioner(); | ||
|
|
||
| //! Prune preconditioner elements | ||
| void prunePreconditioner(); | ||
|
|
||
| //! Function used to return semi-analytical jacobian matrix | ||
| Eigen::SparseMatrix<double> getJacobian() { | ||
| Eigen::SparseMatrix<double> jacobian(m_dim, m_dim); | ||
| jacobian.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end()); | ||
| return jacobian; | ||
| } | ||
|
|
||
| //! Get the threshold value for setting elements | ||
| double threshold() { return m_threshold; } | ||
|
|
||
| //! Get ilut fill factor | ||
| double ilutFillFactor() { return m_fill_factor; } | ||
|
|
||
| //! Get ilut drop tolerance | ||
| double ilutDropTol() { return m_drop_tol; } | ||
|
|
||
| //! Set the threshold value to compare elements against | ||
| //! @param threshold double value used in setting by threshold | ||
| void setThreshold(double threshold) { | ||
| m_threshold = threshold; | ||
| m_prune_precon = (threshold <= 0) ? false : true; | ||
| } | ||
|
|
||
| //! Set drop tolerance for ILUT | ||
| //! @param droptol double value used in setting solver drop tolerance | ||
| void setIlutDropTol(double droptol) { | ||
| m_drop_tol = droptol; | ||
| m_solver.setDroptol(droptol); | ||
| } | ||
|
|
||
| //! Set the fill factor for ILUT solver | ||
| //! @param fillFactor fill in factor for ILUT solver | ||
| void setIlutFillFactor(int fillFactor) { | ||
| m_fill_factor = fillFactor; | ||
| m_solver.setFillfactor(fillFactor); | ||
| } | ||
|
|
||
| //! Print preconditioner contents | ||
| void printPreconditioner() { | ||
| std::stringstream ss; | ||
| Eigen::IOFormat HeavyFmt(Eigen::FullPrecision, 0, ", ", ";\n", "[", "]", "[", "]"); | ||
| ss << Eigen::MatrixXd(m_precon_matrix).format(HeavyFmt); | ||
| writelog(ss.str()); | ||
| } | ||
|
|
||
| //! Print jacobian contents | ||
| void 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()); | ||
| } | ||
|
|
||
| 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 | ||
| std::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; | ||
|
|
||
| //! Minimum value a non-diagonal element must be to be included in | ||
| //! the preconditioner | ||
| double m_threshold = 1e-8; | ||
|
|
||
| //! Bool set whether to prune the matrix or not | ||
| double m_prune_precon = true; | ||
| }; | ||
|
|
||
| } | ||
|
|
||
| #endif |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.