Cantera  3.1.0a1
AdaptivePreconditioner.cpp
Go to the documentation of this file.
1 //! @file AdaptivePreconditioner.cpp
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at https://cantera.org/license.txt for license and copyright information.
5 
6 #include "cantera/base/global.h"
8 
9 namespace Cantera
10 {
11 
12 AdaptivePreconditioner::AdaptivePreconditioner()
13 {
14  setPreconditionerSide("right");
15 }
16 
17 void AdaptivePreconditioner::setValue(size_t row, size_t col, double value)
18 {
19  m_jac_trips.emplace_back(static_cast<int>(row), static_cast<int>(col), value);
20 }
21 
22 void AdaptivePreconditioner::stateAdjustment(vector<double>& state) {
23  // Only keep positive composition based on given tol
24  for (size_t i = 0; i < state.size(); i++) {
25  state[i] = std::max(state[i], m_atol);
26  }
27 }
28 
29 void AdaptivePreconditioner::initialize(size_t networkSize)
30 {
31  // don't use legacy rate constants
33  // reset arrays in case of re-initialization
34  m_jac_trips.clear();
35  // set dimensions of preconditioner from network
36  m_dim = networkSize;
37  // reserve some space for vectors making up SparseMatrix
38  m_jac_trips.reserve(3 * networkSize);
39  // reserve space for preconditioner
40  m_precon_matrix.resize(m_dim, m_dim);
41  // creating sparse identity matrix
42  m_identity.resize(m_dim, m_dim);
43  m_identity.setIdentity();
44  m_identity.makeCompressed();
45  // setting default ILUT parameters
46  if (m_drop_tol == 0) {
47  setIlutDropTol(1e-10);
48  }
49  if (m_drop_tol == 0) {
50  setIlutFillFactor(static_cast<int>(m_dim) / 4);
51  }
52  // update initialized status
53  m_init = true;
54 }
55 
56 
58 {
59  // make into preconditioner as P = (I - gamma * J_bar)
61  // compressing sparse matrix structure
62  m_precon_matrix.makeCompressed();
63  // analyze and factorize
64  m_solver.compute(m_precon_matrix);
65  // check for errors
66  if (m_solver.info() != Eigen::Success) {
67  throw CanteraError("AdaptivePreconditioner::setup",
68  "error code: {}", static_cast<int>(m_solver.info()));
69  }
70 }
71 
73 {
74  // set precon to jacobian
75  m_precon_matrix.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
76  // convert to preconditioner
78  // prune by threshold if desired
79  if (m_prune_precon) {
81  }
82 }
83 
85 {
86  for (int k=0; k<m_precon_matrix.outerSize(); ++k) {
87  for (Eigen::SparseMatrix<double>::InnerIterator it(m_precon_matrix, k); it;
88  ++it) {
89  if (std::abs(it.value()) < m_threshold && it.row() != it.col()) {
90  it.valueRef() = 0;
91  }
92  }
93  }
94 }
95 
96 void AdaptivePreconditioner::solve(const size_t stateSize, double* rhs_vector, double*
97  output)
98 {
99  // creating vectors in the form of Ax=b
100  Eigen::Map<Eigen::VectorXd> bVector(rhs_vector, stateSize);
101  Eigen::Map<Eigen::VectorXd> xVector(output, stateSize);
102  // solve for xVector
103  xVector = m_solver.solve(bVector);
104  if (m_solver.info() != Eigen::Success) {
105  throw CanteraError("AdaptivePreconditioner::solve",
106  "error code: {}", static_cast<int>(m_solver.info()));
107  }
108 }
109 
111  std::stringstream ss;
112  Eigen::IOFormat HeavyFmt(Eigen::FullPrecision, 0, ", ", ";\n", "[", "]", "[", "]");
113  ss << Eigen::MatrixXd(m_precon_matrix).format(HeavyFmt);
114  writelog(ss.str());
115 }
116 
118  std::stringstream ss;
119  Eigen::SparseMatrix<double> jacobian(m_dim, m_dim);
120  jacobian.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
121  ss << Eigen::MatrixXd(jacobian);
122  writelog(ss.str());
123 }
124 
125 }
Declarations for the class AdaptivePreconditioner which is a child class of PreconditionerBase for pr...
Eigen::IncompleteLUT< double > m_solver
Solver used in solving the linear system.
void setIlutDropTol(double droptol)
Set drop tolerance for ILUT.
double m_drop_tol
ILUT drop tolerance.
void solve(const size_t stateSize, double *rhs_vector, double *output) override
Solve a linear system Ax=b where A is the preconditioner.
void setValue(size_t row, size_t col, double value) override
Set a value at the specified row and column of the jacobian triplet vector.
double m_prune_precon
Bool set whether to prune the matrix or not.
double m_threshold
Minimum value a non-diagonal element must be to be included in the preconditioner.
void initialize(size_t networkSize) override
Called during setup for any processes that need to be completed prior to setup functions used in sund...
vector< Eigen::Triplet< double > > m_jac_trips
Vector of triples representing the jacobian used in preconditioning.
void setIlutFillFactor(int fillFactor)
Set the fill factor for ILUT solver.
Eigen::SparseMatrix< double > m_identity
Storage of appropriately sized identity matrix for making the preconditioner.
void printPreconditioner() override
Print preconditioner contents.
void setup() override
Perform preconditioner specific post-reactor setup operations such as factorize.
void stateAdjustment(vector< double > &state) override
Adjust the state vector based on the preconditioner, e.g., Adaptive preconditioning uses a strictly p...
void updatePreconditioner() override
Transform Jacobian vector and write into preconditioner, P = (I - gamma * J)
void prunePreconditioner()
Prune preconditioner elements.
void printJacobian()
Print jacobian contents.
Eigen::SparseMatrix< double > m_precon_matrix
Container that is the sparse preconditioner.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
size_t m_dim
Dimension of the preconditioner.
double m_gamma
gamma value used in M = I - gamma*J
bool m_init
bool saying whether or not the preconditioner is initialized
double m_atol
Absolute tolerance of the ODE solver.
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
Eigen::SparseMatrix< double > jacobian()
Return semi-analytical Jacobian of an AdaptivePreconditioner object.
void use_legacy_rate_constants(bool legacy)
Set definition used for rate constant calculation.
Definition: global.cpp:102
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:175
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:564