From f9f2d3d1f327df9954faac3c4b9eb444a41a3577 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 13 Dec 2024 22:11:35 -0500 Subject: [PATCH 1/9] [Numerics] Avoid having MultiJac inherit from BandMatrix First step toward allowing Jacobians implemented using different matrix storage and factorization methods. --- include/cantera/oneD/MultiJac.h | 22 +++++++++++++++++++++- src/oneD/MultiJac.cpp | 11 ++++++----- 2 files changed, 27 insertions(+), 6 deletions(-) diff --git a/include/cantera/oneD/MultiJac.h b/include/cantera/oneD/MultiJac.h index e9cd68f1ff..5404b4ad6c 100644 --- a/include/cantera/oneD/MultiJac.h +++ b/include/cantera/oneD/MultiJac.h @@ -20,7 +20,7 @@ namespace Cantera * @ingroup onedUtilsGroup * @ingroup derivGroup */ -class MultiJac : public BandMatrix +class MultiJac { public: //! Constructor. @@ -72,6 +72,22 @@ class MultiJac : public BandMatrix m_age = age; } + double& value(size_t i, size_t j) { + return m_mat.value(i, j); + } + + double value(size_t i, size_t j) const { + return m_mat.value(i, j); + } + + void solve(const double* const b, double* const x) { + m_mat.solve(b, x); + } + + int info() const { + return m_mat.info(); + } + //! Return the transient mask. vector& transientMask() { return m_mask; @@ -85,6 +101,9 @@ class MultiJac : public BandMatrix */ OneDim* m_resid; + BandMatrix m_mat; //!< Underlying matrix storage + size_t m_n; //!< Dimension of the (square) matrix + vector m_r1; //!< Perturbed residual vector double m_rtol = 1e-5; //!< Relative tolerance for perturbing solution components @@ -100,6 +119,7 @@ class MultiJac : public BandMatrix int m_nevals = 0; //!< Number of Jacobian evaluations. int m_age = 100000; //!< Age of the Jacobian (times incrementAge() has been called) }; + } #endif diff --git a/src/oneD/MultiJac.cpp b/src/oneD/MultiJac.cpp index 0c4965146b..370646902e 100644 --- a/src/oneD/MultiJac.cpp +++ b/src/oneD/MultiJac.cpp @@ -10,7 +10,8 @@ namespace Cantera { MultiJac::MultiJac(OneDim& r) - : BandMatrix(r.size(),r.bandwidth(),r.bandwidth()) + : m_mat(r.size(),r.bandwidth(),r.bandwidth()) + , m_n(r.size()) { m_resid = &r; m_r1.resize(m_n); @@ -21,7 +22,7 @@ MultiJac::MultiJac(OneDim& r) void MultiJac::updateTransient(double rdt, integer* mask) { for (size_t n = 0; n < m_n; n++) { - value(n,n) = m_ssdiag[n] - mask[n]*rdt; + m_mat.value(n,n) = m_ssdiag[n] - mask[n]*rdt; } } @@ -29,7 +30,7 @@ void MultiJac::eval(double* x0, double* resid0, double rdt) { m_nevals++; clock_t t0 = clock(); - bfill(0.0); + m_mat.bfill(0.0); size_t ipt=0; for (size_t j = 0; j < m_resid->points(); j++) { @@ -56,7 +57,7 @@ void MultiJac::eval(double* x0, double* resid0, double rdt) size_t mv = m_resid->nVars(i); size_t iloc = m_resid->loc(i); for (size_t m = 0; m < mv; m++) { - value(m+iloc,ipt) = (m_r1[m+iloc] - resid0[m+iloc])*rdx; + m_mat.value(m+iloc,ipt) = (m_r1[m+iloc] - resid0[m+iloc])*rdx; } } } @@ -66,7 +67,7 @@ void MultiJac::eval(double* x0, double* resid0, double rdt) } for (size_t n = 0; n < m_n; n++) { - m_ssdiag[n] = value(n,n); + m_ssdiag[n] = m_mat.value(n,n); } m_elapsed += double(clock() - t0)/CLOCKS_PER_SEC; From f282572feafb2fff53908418e024b59d87488aa9 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 30 Dec 2024 13:37:41 -0500 Subject: [PATCH 2/9] [Numerics] Remove MultiJac argument from MultiNewton methods One step toward allowing alternative Jacobian implementations --- include/cantera/oneD/MultiNewton.h | 43 ++++++++++++++++++++++++++---- src/oneD/MultiNewton.cpp | 15 ++++++----- src/oneD/OneDim.cpp | 2 +- 3 files changed, 48 insertions(+), 12 deletions(-) diff --git a/include/cantera/oneD/MultiNewton.h b/include/cantera/oneD/MultiNewton.h index b90f622c9a..ef573eaf45 100644 --- a/include/cantera/oneD/MultiNewton.h +++ b/include/cantera/oneD/MultiNewton.h @@ -6,7 +6,7 @@ #ifndef CT_MULTINEWTON_H #define CT_MULTINEWTON_H -#include "MultiJac.h" +#include "OneDim.h" namespace Cantera { @@ -37,7 +37,17 @@ class MultiNewton //! Compute the undamped Newton step. The residual function is evaluated //! at `x`, but the Jacobian is not recomputed. - void step(double* x, double* step, OneDim& r, MultiJac& jac, int loglevel); + //! @since Starting in %Cantera 3.2, the Jacobian is accessed via the OneDim object. + void step(double* x, double* step, OneDim& r, int loglevel); + + //! @deprecated Use version without MultiJac argument. To be removed after + //! %Cantera 3.2. + void step(double* x, double* stp, OneDim& r, MultiJac& jac, int loglevel) { + warn_deprecated("MultiNewton::step(..., MultiJac&, ...)", + "Use version without MultiJac argument. " + "To be removed after Cantera 3.2"); + step(x, stp, r, loglevel); + } /** * Return the factor by which the undamped Newton step 'step0' @@ -99,7 +109,6 @@ class MultiNewton * @param[out] s1 norm of the subsequent Newton step after taking the damped Newton * step * @param[in] r domain object, used for evaluating residuals over all domains - * @param[in] jac Jacobian evaluator * @param[in] loglevel controls amount of printed diagnostics * @param[in] writetitle controls if logging title is printed * @@ -109,9 +118,22 @@ class MultiNewton * - `-2` no suitable damping coefficient was found within the maximum iterations. * - `-3` the current solution `x0` is too close to the solution bounds and the * step would exceed the bounds on one or more components. + * + * @since Starting in %Cantera 3.2, the Jacobian is accessed via the OneDim object. */ int dampStep(const double* x0, const double* step0, double* x1, double* step1, - double& s1, OneDim& r, MultiJac& jac, int loglevel, bool writetitle); + double& s1, OneDim& r, int loglevel, bool writetitle); + + //! @deprecated Use version without MultiJac argument. To be removed after + //! %Cantera 3.2. + int dampStep(const double* x0, const double* step0, double* x1, double* step1, + double& s1, OneDim& r, MultiJac& jac, int loglevel, bool writetitle) + { + warn_deprecated("MultiNewton::dampStep(..., MultiJac&, ...)", + "Use version without MultiJac argument. " + "To be removed after Cantera 3.2"); + return dampStep(x0, step0, x1, step1, s1, r, loglevel, writetitle); + } //! Compute the weighted 2-norm of `step`. double norm2(const double* x, const double* step, OneDim& r) const; @@ -129,8 +151,19 @@ class MultiNewton * - `-2` no suitable damping coefficient was found within the maximum iterations. * - `-3` the current solution `x0` is too close to the solution bounds and the * step would exceed the bounds on one or more components. + * + * @since Starting in %Cantera 3.2, the Jacobian is accessed via the OneDim object. */ - int solve(double* x0, double* x1, OneDim& r, MultiJac& jac, int loglevel); + int solve(double* x0, double* x1, OneDim& r, int loglevel); + + //! @deprecated Use version without MultiJac argument. To be removed after + //! %Cantera 3.2. + int solve(double* x0, double* x1, OneDim& r, MultiJac& jac, int loglevel) { + warn_deprecated("MultiNewton::solve(..., MultiJac&, ...)", + "Use version without MultiJac argument. " + "To be removed after Cantera 3.2"); + return solve(x0, x1, r, loglevel); + } //! Set options. //! @param maxJacAge Maximum number of steps that can be taken before requiring diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp index 27ae94bb9c..324809bcb5 100644 --- a/src/oneD/MultiNewton.cpp +++ b/src/oneD/MultiNewton.cpp @@ -181,13 +181,14 @@ double MultiNewton::norm2(const double* x, const double* step, OneDim& r) const return sqrt(sum); } -void MultiNewton::step(double* x, double* step, OneDim& r, MultiJac& jac, int loglevel) +void MultiNewton::step(double* x, double* step, OneDim& r, int loglevel) { r.eval(npos, x, step); for (size_t n = 0; n < r.size(); n++) { step[n] = -step[n]; } + auto& jac = r.jacobian(); try { jac.solve(step, step); } catch (CanteraError&) { @@ -228,7 +229,7 @@ double MultiNewton::boundStep(const double* x0, const double* step0, const OneDi int MultiNewton::dampStep(const double* x0, const double* step0, double* x1, double* step1, double& s1, - OneDim& r, MultiJac& jac, int loglevel, bool writetitle) + OneDim& r, int loglevel, bool writetitle) { // write header if (loglevel > 0 && writetitle) { @@ -259,6 +260,7 @@ int MultiNewton::dampStep(const double* x0, const double* step0, // fbound factor to ensure that the solution remains within bounds. double alpha = fbound*1.0; size_t m; + auto& jac = r.jacobian(); for (m = 0; m < m_maxDampIter; m++) { // step the solution by the damped step size // x_{k+1} = x_k + alpha_k*J(x_k)^-1 F(x_k) @@ -268,7 +270,7 @@ int MultiNewton::dampStep(const double* x0, const double* step0, // compute the next undamped step that would result if x1 is accepted // J(x_k)^-1 F(x_k+1) - step(x1, step1, r, jac, loglevel-1); + step(x1, step1, r, loglevel-1); // compute the weighted norm of step1 s1 = norm2(x1, step1, r); @@ -308,7 +310,7 @@ int MultiNewton::dampStep(const double* x0, const double* step0, } } -int MultiNewton::solve(double* x0, double* x1, OneDim& r, MultiJac& jac, int loglevel) +int MultiNewton::solve(double* x0, double* x1, OneDim& r, int loglevel) { clock_t t0 = clock(); int status = 0; @@ -320,6 +322,7 @@ int MultiNewton::solve(double* x0, double* x1, OneDim& r, MultiJac& jac, int log double rdt = r.rdt(); int nJacReeval = 0; + auto& jac = r.jacobian(); while (true) { // Check whether the Jacobian should be re-evaluated. if (jac.age() > m_maxAge) { @@ -337,13 +340,13 @@ int MultiNewton::solve(double* x0, double* x1, OneDim& r, MultiJac& jac, int log } // compute the undamped Newton step - step(&m_x[0], &m_stp[0], r, jac, loglevel-1); + step(&m_x[0], &m_stp[0], r, loglevel-1); // increment the Jacobian age jac.incrementAge(); // damp the Newton step - status = dampStep(&m_x[0], &m_stp[0], x1, &m_stp1[0], s1, r, jac, loglevel-1, write_header); + status = dampStep(&m_x[0], &m_stp[0], x1, &m_stp1[0], s1, r, loglevel-1, write_header); write_header = false; // Successful step, but not converged yet. Take the damped step, and try diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index 6e8290e393..4a558e4351 100644 --- a/src/oneD/OneDim.cpp +++ b/src/oneD/OneDim.cpp @@ -213,7 +213,7 @@ int OneDim::solve(double* x, double* xnew, int loglevel) m_jac->updateTransient(m_rdt, m_mask.data()); m_jac_ok = true; } - return m_newt->solve(x, xnew, *this, *m_jac, loglevel); + return m_newt->solve(x, xnew, *this, loglevel); } void OneDim::evalSSJacobian(double* x, double* rsd) From 1b390e7f32add3831bb9034ea18666cfb49ec8bf Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 30 Dec 2024 15:09:17 -0500 Subject: [PATCH 3/9] [Numerics] Move 1D Jacobian evaluation into OneDim --- include/cantera/oneD/MultiJac.h | 14 +++++++- include/cantera/oneD/OneDim.h | 14 ++++++++ src/oneD/MultiJac.cpp | 60 +++++++++------------------------ src/oneD/MultiNewton.cpp | 3 +- src/oneD/OneDim.cpp | 56 +++++++++++++++++++++++++++--- 5 files changed, 96 insertions(+), 51 deletions(-) diff --git a/include/cantera/oneD/MultiJac.h b/include/cantera/oneD/MultiJac.h index 5404b4ad6c..264ad80c46 100644 --- a/include/cantera/oneD/MultiJac.h +++ b/include/cantera/oneD/MultiJac.h @@ -7,6 +7,7 @@ #define CT_MULTIJAC_H #include "cantera/numerics/BandMatrix.h" +#include "cantera/numerics/PreconditionerBase.h" #include "OneDim.h" namespace Cantera @@ -20,7 +21,7 @@ namespace Cantera * @ingroup onedUtilsGroup * @ingroup derivGroup */ -class MultiJac +class MultiJac : public PreconditionerBase { public: //! Constructor. @@ -41,18 +42,29 @@ class MultiJac * @param x0 Solution vector at which to evaluate the Jacobian * @param resid0 Residual vector at x0 * @param rdt Reciprocal of the time step + * @deprecated To be removed after %Cantera 3.2. Jacobian evaluation moved to + * OneDim::evalJacobian(). */ void eval(double* x0, double* resid0, double rdt); + void reset() override; + void setValue(size_t row, size_t col, double value) override; + //! Elapsed CPU time spent computing the Jacobian. double elapsedTime() const { return m_elapsed; } + void updateElapsed(double evalTime) { + m_elapsed += evalTime; + } //! Number of Jacobian evaluations. int nEvals() const { return m_nevals; } + void incrementEvals() { + m_nevals++; + } //! Number of times 'incrementAge' has been called since the last evaluation int age() const { diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h index 1a73ab74aa..97162b20f9 100644 --- a/include/cantera/oneD/OneDim.h +++ b/include/cantera/oneD/OneDim.h @@ -195,6 +195,17 @@ class OneDim */ void eval(size_t j, double* x, double* r, double rdt=-1.0, int count = 1); + /** + * Evaluates the Jacobian at x0 using finite differences. + * + * The Jacobian is computed by perturbing each component of `x0`, evaluating the + * residual function, and then estimating the partial derivatives numerically using + * finite differences to determine the corresponding column of the Jacobian. + * + * @param x0 State vector at which to evaluate the Jacobian + */ + void evalJacobian(double* x0); + //! Return a pointer to the domain global point *i* belongs to. /*! * The domains are scanned right-to-left, and the first one with starting @@ -371,6 +382,9 @@ class OneDim size_t m_bw = 0; //!< Jacobian bandwidth size_t m_size = 0; //!< solution vector size + //! Work arrays used during Jacobian evaluation + vector m_work1, m_work2; + //! All domains comprising the system vector> m_dom; diff --git a/src/oneD/MultiJac.cpp b/src/oneD/MultiJac.cpp index 370646902e..5083270fd4 100644 --- a/src/oneD/MultiJac.cpp +++ b/src/oneD/MultiJac.cpp @@ -19,6 +19,19 @@ MultiJac::MultiJac(OneDim& r) m_mask.resize(m_n); } +void MultiJac::reset() { + m_mat.bfill(0.0); + m_age = 10000; +} + +void MultiJac::setValue(size_t row, size_t col, double value) +{ + m_mat(row, col) = value; + if (row == col) { + m_ssdiag[row] = value; + } +} + void MultiJac::updateTransient(double rdt, integer* mask) { for (size_t n = 0; n < m_n; n++) { @@ -28,50 +41,9 @@ void MultiJac::updateTransient(double rdt, integer* mask) void MultiJac::eval(double* x0, double* resid0, double rdt) { - m_nevals++; - clock_t t0 = clock(); - m_mat.bfill(0.0); - size_t ipt=0; - - for (size_t j = 0; j < m_resid->points(); j++) { - size_t nv = m_resid->nVars(j); - for (size_t n = 0; n < nv; n++) { - // perturb x(n); preserve sign(x(n)) - double xsave = x0[ipt]; - double dx; - if (xsave >= 0) { - dx = xsave*m_rtol + m_atol; - } else { - dx = xsave*m_rtol - m_atol; - } - x0[ipt] = xsave + dx; - dx = x0[ipt] - xsave; - double rdx = 1.0/dx; - - // calculate perturbed residual - m_resid->eval(j, x0, m_r1.data(), rdt, 0); - - // compute nth column of Jacobian - for (size_t i = j - 1; i != j+2; i++) { - if (i != npos && i < m_resid->points()) { - size_t mv = m_resid->nVars(i); - size_t iloc = m_resid->loc(i); - for (size_t m = 0; m < mv; m++) { - m_mat.value(m+iloc,ipt) = (m_r1[m+iloc] - resid0[m+iloc])*rdx; - } - } - } - x0[ipt] = xsave; - ipt++; - } - } - - for (size_t n = 0; n < m_n; n++) { - m_ssdiag[n] = m_mat.value(n,n); - } - - m_elapsed += double(clock() - t0)/CLOCKS_PER_SEC; - m_age = 0; + warn_deprecated("MultiJac::eval", "To be removed after Cantera 3.2. " + "Jacobian evaluation moved to OneDim::evalJacobian()."); + m_resid->evalJacobian(x0); } } // namespace diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp index 324809bcb5..1815657610 100644 --- a/src/oneD/MultiNewton.cpp +++ b/src/oneD/MultiNewton.cpp @@ -333,8 +333,7 @@ int MultiNewton::solve(double* x0, double* x1, OneDim& r, int loglevel) } if (forceNewJac) { - r.eval(npos, &m_x[0], &m_stp[0], 0.0, 0); - jac.eval(&m_x[0], &m_stp[0], 0.0); + r.evalJacobian(&m_x[0]); jac.updateTransient(rdt, r.transientMask().data()); forceNewJac = false; } diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index 4a558e4351..b1feaa6a73 100644 --- a/src/oneD/OneDim.cpp +++ b/src/oneD/OneDim.cpp @@ -208,8 +208,7 @@ void OneDim::resize() int OneDim::solve(double* x, double* xnew, int loglevel) { if (!m_jac_ok) { - eval(npos, x, xnew, 0.0, 0); - m_jac->eval(x, xnew, 0.0); + evalJacobian(x); m_jac->updateTransient(m_rdt, m_mask.data()); m_jac_ok = true; } @@ -221,8 +220,7 @@ void OneDim::evalSSJacobian(double* x, double* rsd) double rdt_save = m_rdt; m_jac_ok = false; setSteadyMode(); - eval(npos, x, rsd, 0.0, 0); - m_jac->eval(x, rsd, 0.0); + evalJacobian(x); m_rdt = rdt_save; } @@ -270,6 +268,56 @@ void OneDim::eval(size_t j, double* x, double* r, double rdt, int count) } } +void OneDim::evalJacobian(double* x0) +{ + m_jac->reset(); + clock_t t0 = clock(); + double rtol = 1e-5; + double atol = sqrt(std::numeric_limits::epsilon()); + + m_work1.resize(size()); + m_work2.resize(size()); + eval(npos, x0, m_work1.data(), 0.0, 0); + size_t ipt = 0; + for (size_t j = 0; j < points(); j++) { + size_t nv = nVars(j); + for (size_t n = 0; n < nv; n++) { + // perturb x(n); preserve sign(x(n)) + double xsave = x0[ipt]; + double dx; + if (xsave >= 0) { + dx = xsave*rtol + atol; + } else { + dx = xsave*rtol - atol; + } + x0[ipt] = xsave + dx; + dx = x0[ipt] - xsave; + double rdx = 1.0/dx; + + // calculate perturbed residual + eval(j, x0, m_work2.data(), 0.0, 0); + + // compute nth column of Jacobian + for (size_t i = j - 1; i != j+2; i++) { + if (i != npos && i < points()) { + size_t mv = nVars(i); + size_t iloc = loc(i); + for (size_t m = 0; m < mv; m++) { + m_jac->setValue(m + iloc, ipt, + (m_work2[m+iloc] - m_work1[m+iloc]) * rdx); + } + } + } + x0[ipt] = xsave; + ipt++; + } + } + + m_jac->updateElapsed(double(clock() - t0) / CLOCKS_PER_SEC); + m_jac->incrementEvals(); + m_jac->setAge(0); +} + double OneDim::ssnorm(double* x, double* r) { eval(npos, x, r, 0.0, 0); From 6a8f4f23cdd24265ec23cf06b119296d244790df Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 30 Dec 2024 16:44:57 -0500 Subject: [PATCH 4/9] [Numerics] Make MultiJac an implementation of PreconditionerBase --- .../cantera/numerics/AdaptivePreconditioner.h | 4 ++ include/cantera/numerics/PreconditionerBase.h | 53 +++++++++++++++ include/cantera/oneD/MultiJac.h | 67 +++++-------------- include/cantera/oneD/OneDim.h | 7 +- include/cantera/oneD/Sim1D.h | 1 + src/numerics/PreconditionerFactory.cpp | 2 + src/oneD/Domain1D.cpp | 2 +- src/oneD/MultiJac.cpp | 31 +++++++-- src/oneD/MultiNewton.cpp | 24 +++---- src/oneD/OneDim.cpp | 14 +++- src/oneD/Sim1D.cpp | 7 +- 11 files changed, 139 insertions(+), 73 deletions(-) diff --git a/include/cantera/numerics/AdaptivePreconditioner.h b/include/cantera/numerics/AdaptivePreconditioner.h index 84891c4050..53f1fbb428 100644 --- a/include/cantera/numerics/AdaptivePreconditioner.h +++ b/include/cantera/numerics/AdaptivePreconditioner.h @@ -45,6 +45,10 @@ class AdaptivePreconditioner : public PreconditionerBase void updatePreconditioner() override; + int info() const override { + return static_cast(m_solver.info()); + } + //! Prune preconditioner elements void prunePreconditioner(); diff --git a/include/cantera/numerics/PreconditionerBase.h b/include/cantera/numerics/PreconditionerBase.h index b735f58334..59b9eddcad 100644 --- a/include/cantera/numerics/PreconditionerBase.h +++ b/include/cantera/numerics/PreconditionerBase.h @@ -59,6 +59,13 @@ class PreconditionerBase m_precon_side = preconSide; } + //! Update the diagonal terms in the Jacobian by using the transient mask. + //! @param rdt Reciprocal of the time step [1/s] + //! @param mask Mask for transient terms: 1 if transient, 0 if algebraic. + virtual void updateTransient(double rdt, int* mask) { + throw NotImplementedError("PreconditionerBase::updateTransient"); + } + //! 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 @@ -84,6 +91,10 @@ class PreconditionerBase throw NotImplementedError("PreconditionerBase::initialize"); }; + //! Used to provide system bandwidth for implementations that use banded matrix + //! storage. Ignored if not needed. + virtual void setBandwidth(size_t bw) {} + //! Print preconditioner contents virtual void printPreconditioner() { throw NotImplementedError("PreconditionerBase::printPreconditioner"); @@ -112,6 +123,43 @@ class PreconditionerBase m_atol = atol; } + //! Get latest status of linear solver. Zero indicates success. Meaning of non-zero + //! values is solver dependent. + virtual int info() const { + throw NotImplementedError("PreconditionerBase::info"); + } + + //! Elapsed CPU time spent computing the Jacobian elements. + double elapsedTime() const { + return m_elapsed; + } + void updateElapsed(double evalTime) { + m_elapsed += evalTime; + } + + //! Number of Jacobian evaluations. + int nEvals() const { + return m_nevals; + } + void incrementEvals() { + m_nevals++; + } + + //! Number of times 'incrementAge' has been called since the last evaluation + int age() const { + return m_age; + } + + //! Increment the Jacobian age. + void incrementAge() { + m_age++; + } + + //! Set the Jacobian age. + void setAge(int age) { + m_age = age; + } + protected: //! Dimension of the preconditioner size_t m_dim; @@ -125,6 +173,11 @@ class PreconditionerBase //! Absolute tolerance of the ODE solver double m_atol = 0; + double m_elapsed = 0.0; //!< Elapsed CPU time taken to compute the Jacobian + int m_nevals = 0; //!< Number of Jacobian evaluations. + + int m_age = 100000; //!< Age of the Jacobian (times incrementAge() has been called) + string m_precon_side = "none"; }; diff --git a/include/cantera/oneD/MultiJac.h b/include/cantera/oneD/MultiJac.h index 264ad80c46..907e77b5b2 100644 --- a/include/cantera/oneD/MultiJac.h +++ b/include/cantera/oneD/MultiJac.h @@ -26,8 +26,11 @@ class MultiJac : public PreconditionerBase public: //! Constructor. //! @param r The nonlinear system for which to compute the Jacobian. + //! @deprecated To be removed after %Cantera 3.2. Use default constructor instead. MultiJac(OneDim& r); + MultiJac() = default; + /** * Evaluates the Jacobian at x0 using finite differences. The unperturbed residual * function is resid0, which must be supplied on input. The parameter 'rdt' is the @@ -49,40 +52,9 @@ class MultiJac : public PreconditionerBase void reset() override; void setValue(size_t row, size_t col, double value) override; - - //! Elapsed CPU time spent computing the Jacobian. - double elapsedTime() const { - return m_elapsed; - } - void updateElapsed(double evalTime) { - m_elapsed += evalTime; - } - - //! Number of Jacobian evaluations. - int nEvals() const { - return m_nevals; - } - void incrementEvals() { - m_nevals++; - } - - //! Number of times 'incrementAge' has been called since the last evaluation - int age() const { - return m_age; - } - - //! Increment the Jacobian age. - void incrementAge() { - m_age++; - } - - //! Update the transient terms in the Jacobian by using the transient mask. - void updateTransient(double rdt, integer* mask); - - //! Set the Jacobian age. - void setAge(int age) { - m_age = age; - } + void initialize(size_t nVars) override; + void setBandwidth(size_t bw) override; + void updateTransient(double rdt, integer* mask) override; double& value(size_t i, size_t j) { return m_mat.value(i, j); @@ -96,12 +68,18 @@ class MultiJac : public PreconditionerBase m_mat.solve(b, x); } - int info() const { + void solve(const size_t stateSize, double* b, double* x) override { + m_mat.solve(b, x); + } + + int info() const override { return m_mat.info(); } //! Return the transient mask. + //! @deprecated Unused. To be removed after %Cantera 3.2. vector& transientMask() { + warn_deprecated("MultiJac::transientMask", "To be removed after Cantera 3.2"); return m_mask; } @@ -110,26 +88,17 @@ class MultiJac : public PreconditionerBase /*! * This is a pointer to the residual evaluator. This object isn't owned by * this Jacobian object. + * + * @deprecated Unused. To be removed after %Cantera 3.2. */ - OneDim* m_resid; + OneDim* m_resid = nullptr; BandMatrix m_mat; //!< Underlying matrix storage - size_t m_n; //!< Dimension of the (square) matrix - - vector m_r1; //!< Perturbed residual vector - double m_rtol = 1e-5; //!< Relative tolerance for perturbing solution components - - //! Absolute tolerance for perturbing solution components - double m_atol = sqrt(std::numeric_limits::epsilon()); - - double m_elapsed = 0.0; //!< Elapsed CPU time taken to compute the Jacobian vector m_ssdiag; //!< Diagonal of the steady-state Jacobian - //! Transient mask for transient terms, 1 if transient, 0 if steady-state + //! Transient mask for transient terms, 1 if transient, 0 if steady-state. + //! @deprecated Unused. To be removed after %Cantera 3.2. vector m_mask; - - int m_nevals = 0; //!< Number of Jacobian evaluations. - int m_age = 100000; //!< Age of the Jacobian (times incrementAge() has been called) }; } diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h index 97162b20f9..580b403cd8 100644 --- a/include/cantera/oneD/OneDim.h +++ b/include/cantera/oneD/OneDim.h @@ -10,6 +10,7 @@ #include "Domain1D.h" #include "MultiJac.h" +#include "cantera/numerics/PreconditionerBase.h" namespace Cantera { @@ -43,6 +44,10 @@ class OneDim //! @ingroup derivGroup MultiJac& jacobian(); + shared_ptr getJacobian() { + return m_jac; + } + //! Return a reference to the Newton iterator. MultiNewton& newton(); @@ -374,7 +379,7 @@ class OneDim shared_ptr> m_state; //!< Solution vector - unique_ptr m_jac; //!< Jacobian evaluator + shared_ptr m_jac; //!< Jacobian evaluator unique_ptr m_newt; //!< Newton iterator double m_rdt = 0.0; //!< reciprocal of time step bool m_jac_ok = false; //!< if true, Jacobian is current diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h index b351934561..86ce8d56ec 100644 --- a/include/cantera/oneD/Sim1D.h +++ b/include/cantera/oneD/Sim1D.h @@ -331,6 +331,7 @@ class Sim1D : public OneDim void getInitialSoln(); //! Get the Jacobian element @f$ J_{ij} = \partial f_i / \partial x_j \f$ + //! @deprecated To be removed after %Cantera 3.2. double jacobian(int i, int j); //! Evaluate the Jacobian in steady-state mode. diff --git a/src/numerics/PreconditionerFactory.cpp b/src/numerics/PreconditionerFactory.cpp index 0ffa42dd86..5911eb0c1d 100644 --- a/src/numerics/PreconditionerFactory.cpp +++ b/src/numerics/PreconditionerFactory.cpp @@ -5,6 +5,7 @@ #include "cantera/numerics/PreconditionerFactory.h" #include "cantera/numerics/AdaptivePreconditioner.h" +#include "cantera/oneD/MultiJac.h" namespace Cantera { @@ -30,6 +31,7 @@ std::mutex PreconditionerFactory::precon_mutex; PreconditionerFactory::PreconditionerFactory() { reg("Adaptive", []() { return new AdaptivePreconditioner(); }); + reg("banded-direct", []() { return new MultiJac(); }); } shared_ptr newPreconditioner(const string& precon) diff --git a/src/oneD/Domain1D.cpp b/src/oneD/Domain1D.cpp index c5c2790d4f..7f4f1e43ec 100644 --- a/src/oneD/Domain1D.cpp +++ b/src/oneD/Domain1D.cpp @@ -113,7 +113,7 @@ void Domain1D::setSteadyTolerances(double rtol, double atol, size_t n) void Domain1D::needJacUpdate() { if (m_container) { - m_container->jacobian().setAge(10000); + m_container->getJacobian()->setAge(10000); m_container->saveStats(); } } diff --git a/src/oneD/MultiJac.cpp b/src/oneD/MultiJac.cpp index 5083270fd4..f6116bee79 100644 --- a/src/oneD/MultiJac.cpp +++ b/src/oneD/MultiJac.cpp @@ -10,20 +10,33 @@ namespace Cantera { MultiJac::MultiJac(OneDim& r) - : m_mat(r.size(),r.bandwidth(),r.bandwidth()) - , m_n(r.size()) { + warn_deprecated("MultiJac::MultiJac(OneDim&)", + "Non-default constructor to be removed after Cantera 3.2."); m_resid = &r; - m_r1.resize(m_n); - m_ssdiag.resize(m_n); - m_mask.resize(m_n); + initialize(m_dim); + setBandwidth(r.bandwidth()); } -void MultiJac::reset() { +void MultiJac::reset() +{ m_mat.bfill(0.0); m_age = 10000; } +void MultiJac::initialize(size_t nVars) +{ + m_dim = nVars; + m_mat.resize(m_dim, m_mat.nSubDiagonals(), m_mat.nSuperDiagonals()); + m_ssdiag.resize(m_dim); + m_mask.resize(m_dim); +} + +void MultiJac::setBandwidth(size_t bw) +{ + m_mat.resize(m_dim, bw, bw); +} + void MultiJac::setValue(size_t row, size_t col, double value) { m_mat(row, col) = value; @@ -34,7 +47,7 @@ void MultiJac::setValue(size_t row, size_t col, double value) void MultiJac::updateTransient(double rdt, integer* mask) { - for (size_t n = 0; n < m_n; n++) { + for (size_t n = 0; n < m_dim; n++) { m_mat.value(n,n) = m_ssdiag[n] - mask[n]*rdt; } } @@ -43,6 +56,10 @@ void MultiJac::eval(double* x0, double* resid0, double rdt) { warn_deprecated("MultiJac::eval", "To be removed after Cantera 3.2. " "Jacobian evaluation moved to OneDim::evalJacobian()."); + if (!m_resid) { + throw CanteraError("MultiJac::eval", "Can only be used in combination with " + "(deprecated) constructor that takes a OneDim object."); + } m_resid->evalJacobian(x0); } diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp index 1815657610..7f317cb125 100644 --- a/src/oneD/MultiNewton.cpp +++ b/src/oneD/MultiNewton.cpp @@ -188,13 +188,13 @@ void MultiNewton::step(double* x, double* step, OneDim& r, int loglevel) step[n] = -step[n]; } - auto& jac = r.jacobian(); + auto jac = r.getJacobian(); try { - jac.solve(step, step); + jac->solve(r.size(), step, step); } catch (CanteraError&) { - if (jac.info() > 0) { + if (jac->info() > 0) { // Positive value for "info" indicates the row where factorization failed - size_t row = static_cast(jac.info() - 1); + size_t row = static_cast(jac->info() - 1); // Find the domain, grid point, and solution component corresponding // to this row for (size_t n = 0; n < r.nDomains(); n++) { @@ -260,7 +260,7 @@ int MultiNewton::dampStep(const double* x0, const double* step0, // fbound factor to ensure that the solution remains within bounds. double alpha = fbound*1.0; size_t m; - auto& jac = r.jacobian(); + auto jac = r.getJacobian(); for (m = 0; m < m_maxDampIter; m++) { // step the solution by the damped step size // x_{k+1} = x_k + alpha_k*J(x_k)^-1 F(x_k) @@ -280,7 +280,7 @@ int MultiNewton::dampStep(const double* x0, const double* step0, writelog("\n {:<4d} {:<9.3e} {:<9.3e} {:>6.3f} {:>6.3f} {:>6.3f} {:<5d} {:d}/{:d}", m, alpha, fbound, log10(ss+SmallNumber), log10(s0+SmallNumber), log10(s1+SmallNumber), - jac.nEvals(), jac.age(), m_maxAge); + jac->nEvals(), jac->age(), m_maxAge); } // If the norm of s1 is less than the norm of s0, then accept this @@ -322,10 +322,10 @@ int MultiNewton::solve(double* x0, double* x1, OneDim& r, int loglevel) double rdt = r.rdt(); int nJacReeval = 0; - auto& jac = r.jacobian(); + auto jac = r.getJacobian(); while (true) { // Check whether the Jacobian should be re-evaluated. - if (jac.age() > m_maxAge) { + if (jac->age() > m_maxAge) { if (loglevel > 1) { writelog("\n Maximum Jacobian age reached ({}), updating it.", m_maxAge); } @@ -334,7 +334,7 @@ int MultiNewton::solve(double* x0, double* x1, OneDim& r, int loglevel) if (forceNewJac) { r.evalJacobian(&m_x[0]); - jac.updateTransient(rdt, r.transientMask().data()); + jac->updateTransient(rdt, r.transientMask().data()); forceNewJac = false; } @@ -342,7 +342,7 @@ int MultiNewton::solve(double* x0, double* x1, OneDim& r, int loglevel) step(&m_x[0], &m_stp[0], r, loglevel-1); // increment the Jacobian age - jac.incrementAge(); + jac->incrementAge(); // damp the Newton step status = dampStep(&m_x[0], &m_stp[0], x1, &m_stp1[0], s1, r, loglevel-1, write_header); @@ -354,14 +354,14 @@ int MultiNewton::solve(double* x0, double* x1, OneDim& r, int loglevel) copy(x1, x1 + m_n, m_x.begin()); } else if (status == 1) { // convergence if (rdt == 0) { - jac.setAge(0); // for efficient sensitivity analysis + jac->setAge(0); // for efficient sensitivity analysis } break; } else if (status < 0) { // If dampStep fails, first try a new Jacobian if an old one was // being used. If it was a new Jacobian, then return -1 to signify // failure. - if (jac.age() > 1) { + if (jac->age() > 1) { forceNewJac = true; if (nJacReeval > 3) { break; diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index b1feaa6a73..d37f786517 100644 --- a/src/oneD/OneDim.cpp +++ b/src/oneD/OneDim.cpp @@ -7,6 +7,7 @@ #include "cantera/numerics/Func1.h" #include "cantera/oneD/MultiNewton.h" #include "cantera/base/AnyMap.h" +#include "cantera/numerics/PreconditionerFactory.h" #include #include @@ -86,7 +87,14 @@ void OneDim::addDomain(shared_ptr d) MultiJac& OneDim::jacobian() { - return *m_jac; + warn_deprecated("OneDim::jacobian", + "Replaced by getJacobian. To be removed after Cantera 3.2."); + auto multijac = dynamic_pointer_cast(m_jac); + if (multijac) { + return *multijac; + } else { + throw CanteraError("OneDim::jacobian", "Active Jacobian is not a MultiJac"); + } } MultiNewton& OneDim::newton() { @@ -201,7 +209,9 @@ void OneDim::resize() m_mask.resize(size()); // delete the current Jacobian evaluator and create a new one - m_jac = make_unique(*this); + m_jac = newPreconditioner("banded-direct"); + m_jac->initialize(size()); + m_jac->setBandwidth(bandwidth()); m_jac_ok = false; } diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index 4ef88354ca..0a0ce46444 100644 --- a/src/oneD/Sim1D.cpp +++ b/src/oneD/Sim1D.cpp @@ -892,6 +892,7 @@ size_t Sim1D::maxGridPoints(size_t dom) double Sim1D::jacobian(int i, int j) { + warn_deprecated("Sim1D::jacobian", "To be removed after Cantera 3.2."); return OneDim::jacobian().value(i,j); } @@ -910,6 +911,10 @@ void Sim1D::solveAdjoint(const double* b, double* lambda) D->forceFullUpdate(false); } + auto multijac = dynamic_pointer_cast(m_jac); + if (!multijac) { + throw CanteraError("Sim1D::solveAdjoint", "Banded (MultiJac) required"); + } // Form J^T size_t bw = bandwidth(); BandMatrix Jt(size(), bw, bw); @@ -917,7 +922,7 @@ void Sim1D::solveAdjoint(const double* b, double* lambda) size_t j1 = (i > bw) ? i - bw : 0; size_t j2 = (i + bw >= size()) ? size() - 1: i + bw; for (size_t j = j1; j <= j2; j++) { - Jt(j,i) = m_jac->value(i,j); + Jt(j,i) = multijac->value(i,j); } } From 8182131da164cab058e1158e8fece9203081d041 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 1 Jan 2025 15:35:42 -0500 Subject: [PATCH 5/9] [Numerics] Enable using AdaptivePreconditioner as OneDim linear solver --- .../cantera/numerics/AdaptivePreconditioner.h | 9 ++-- include/cantera/numerics/PreconditionerBase.h | 10 ++++ include/cantera/oneD/MultiJac.h | 4 ++ include/cantera/oneD/OneDim.h | 7 +++ interfaces/cython/cantera/_onedim.pxd | 3 ++ interfaces/cython/cantera/_onedim.pyx | 13 +++++ interfaces/cython/cantera/preconditioners.pxd | 14 ++++++ interfaces/cython/cantera/preconditioners.pyx | 47 ++++++++++++++++++- src/numerics/AdaptivePreconditioner.cpp | 47 ++++++++++++++----- src/oneD/MultiJac.cpp | 1 + src/oneD/OneDim.cpp | 13 ++++- src/zeroD/ReactorNet.cpp | 2 +- test/zeroD/test_zeroD.cpp | 6 +-- 13 files changed, 151 insertions(+), 25 deletions(-) diff --git a/include/cantera/numerics/AdaptivePreconditioner.h b/include/cantera/numerics/AdaptivePreconditioner.h index 53f1fbb428..534daff662 100644 --- a/include/cantera/numerics/AdaptivePreconditioner.h +++ b/include/cantera/numerics/AdaptivePreconditioner.h @@ -35,15 +35,14 @@ class AdaptivePreconditioner : public PreconditionerBase 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& state) override; - void updatePreconditioner() override; + void updateTransient(double rdt, int* mask) override; int info() const override { return static_cast(m_solver.info()); diff --git a/include/cantera/numerics/PreconditionerBase.h b/include/cantera/numerics/PreconditionerBase.h index 59b9eddcad..955368f226 100644 --- a/include/cantera/numerics/PreconditionerBase.h +++ b/include/cantera/numerics/PreconditionerBase.h @@ -34,6 +34,8 @@ class PreconditionerBase virtual ~PreconditionerBase() {} + virtual const string type() const = 0; + //! 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 @@ -75,6 +77,8 @@ class PreconditionerBase }; //! Perform preconditioner specific post-reactor setup operations such as factorize. + //! @deprecated To be removed after %Cantera 3.2. Use updatePreconditioner() and + //! factorize() instead. virtual void setup() { throw NotImplementedError("PreconditionerBase::setup"); }; @@ -161,6 +165,12 @@ class PreconditionerBase } protected: + //! Factorize the system matrix. This method should be called at the end of + //! implementations of updatePreconditioner() and updateTransient(). + virtual void factorize() { + throw NotImplementedError("preconditionerBase::factorize"); + } + //! Dimension of the preconditioner size_t m_dim; diff --git a/include/cantera/oneD/MultiJac.h b/include/cantera/oneD/MultiJac.h index 907e77b5b2..2f02feb996 100644 --- a/include/cantera/oneD/MultiJac.h +++ b/include/cantera/oneD/MultiJac.h @@ -51,10 +51,14 @@ class MultiJac : public PreconditionerBase void eval(double* x0, double* resid0, double rdt); void reset() override; + const string type() const override { return "banded-direct"; } void setValue(size_t row, size_t col, double value) override; void initialize(size_t nVars) override; void setBandwidth(size_t bw) override; void updateTransient(double rdt, integer* mask) override; + void factorize() override { + m_mat.factor(); + } double& value(size_t i, size_t j) { return m_mat.value(i, j); diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h index 580b403cd8..2bdf560c14 100644 --- a/include/cantera/oneD/OneDim.h +++ b/include/cantera/oneD/OneDim.h @@ -51,6 +51,13 @@ class OneDim //! Return a reference to the Newton iterator. MultiNewton& newton(); + //! Set the linear solver used to hold the Jacobian matrix and solve linear systems + //! as part of each Newton iteration. The default is a direct, banded solver. + void setLinearSolver(shared_ptr solver); + + //! Get the type of the linear solver being used. + shared_ptr linearSolver() const { return m_jac; } + /** * Solve F(x) = 0, where F(x) is the multi-domain residual function. * diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index 37c326d186..834af4228f 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -9,6 +9,7 @@ from .solutionbase cimport * from .kinetics cimport * from .func1 cimport * from .thermo cimport * +from .preconditioners cimport * cdef extern from "cantera/oneD/DomainFactory.h" namespace "Cantera": @@ -114,6 +115,8 @@ cdef extern from "cantera/oneD/Sim1D.h": void restoreSteadySolution() except +translate_exception void setMaxTimeStepCount(int) int maxTimeStepCount() + void setLinearSolver(shared_ptr[CxxPreconditionerBase]) except +translate_exception + shared_ptr[CxxPreconditionerBase] linearSolver() void getInitialSoln() except +translate_exception void solve(int, cbool) except +translate_exception void refine(int) except +translate_exception diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index 3ba626eb73..15fc04b784 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -1012,6 +1012,19 @@ cdef class Sim1D: def __set__(self, nmax): self.sim.setMaxTimeStepCount(nmax) + @property + def linear_solver(self): + """ + Get/Set the the linear solver used to hold the Jacobian matrix and solve linear + systems as part of each Newton iteration. The default is a banded, direct + solver. + """ + return PreconditionerBase.wrap(self.sim.linearSolver()) + + @linear_solver.setter + def linear_solver(self, PreconditionerBase precon): + self.sim.setLinearSolver(precon.pbase) + def set_initial_guess(self, *args, **kwargs): """ Store arguments for initial guess and prepare storage for solution. diff --git a/interfaces/cython/cantera/preconditioners.pxd b/interfaces/cython/cantera/preconditioners.pxd index ebe80f6a02..e66a5cd3a1 100644 --- a/interfaces/cython/cantera/preconditioners.pxd +++ b/interfaces/cython/cantera/preconditioners.pxd @@ -11,6 +11,7 @@ cdef extern from "cantera/numerics/PreconditionerBase.h" namespace "Cantera": cdef cppclass CxxPreconditionerBase "Cantera::PreconditionerBase": CxxPreconditionerBase() string preconditionerSide() + string type() void setPreconditionerSide(string) except +translate_exception cdef extern from "cantera/numerics/AdaptivePreconditioner.h" namespace "Cantera": @@ -26,12 +27,25 @@ cdef extern from "cantera/numerics/AdaptivePreconditioner.h" namespace "Cantera" void printPreconditioner() CxxSparseMatrix matrix() except +translate_exception +cdef extern from "cantera/oneD/MultiJac.h" namespace "Cantera": + cdef cppclass CxxMultiJac "Cantera::MultiJac" (CxxPreconditionerBase): + CxxMultiJac() except +translate_exception + cdef extern from "cantera/numerics/PreconditionerFactory.h" namespace "Cantera": cdef shared_ptr[CxxPreconditionerBase] newPreconditioner(string) except\ +translate_exception cdef class PreconditionerBase: + @staticmethod + cdef wrap(shared_ptr[CxxPreconditionerBase]) + cdef set_cxx_object(self) cdef shared_ptr[CxxPreconditionerBase] pbase cdef class AdaptivePreconditioner(PreconditionerBase): + cdef set_cxx_object(self) cdef CxxAdaptivePreconditioner* preconditioner + +cdef class BandedJacobian(PreconditionerBase): + cdef set_cxx_object(self) + cdef CxxMultiJac* preconditioner + diff --git a/interfaces/cython/cantera/preconditioners.pyx b/interfaces/cython/cantera/preconditioners.pyx index ebe5bcc300..d7df649719 100644 --- a/interfaces/cython/cantera/preconditioners.pyx +++ b/interfaces/cython/cantera/preconditioners.pyx @@ -4,6 +4,10 @@ from ._utils cimport stringify, pystr, py_to_anymap, anymap_to_py from .kinetics cimport get_from_sparse +# dictionary to store reaction rate classes +cdef dict _preconditioner_class_registry = {} + + cdef class PreconditionerBase: """ Common base class for preconditioners. @@ -11,8 +15,36 @@ cdef class PreconditionerBase: precon_type = "PreconditionerBase" precon_linear_solver_type = "GMRES" - def __cinit__(self, *args, **kwargs): - self.pbase = newPreconditioner(stringify(self.precon_type)) + def __cinit__(self, *args, init=True, **kwargs): + if init: + self.pbase = newPreconditioner(stringify(self.precon_type)) + + @staticmethod + cdef wrap(shared_ptr[CxxPreconditionerBase] base): + """ + Wrap a C++ PreconditionerBase object with a Python object of the correct + derived type. + """ + # Ensure all derived types are registered + if not _preconditioner_class_registry: + def register_subclasses(cls): + for c in cls.__subclasses__(): + _preconditioner_class_registry[getattr(c, "precon_type")] = c + register_subclasses(c) + + # Update global class registry + register_subclasses(PreconditionerBase) + + cxx_type = pystr(base.get().type()) + cls = _preconditioner_class_registry.get(cxx_type, PreconditionerBase) + cdef PreconditionerBase precon + precon = cls(init=False) + precon.pbase = base + precon.set_cxx_object() + return precon + + cdef set_cxx_object(self): + pass property side: """ @@ -33,6 +65,9 @@ cdef class AdaptivePreconditioner(PreconditionerBase): def __cinit__(self, *args, **kwargs): self.preconditioner = (self.pbase.get()) + cdef set_cxx_object(self): + self.preconditioner = self.pbase.get() + property threshold: """ The threshold of the preconditioner is used to remove or prune any off diagonal @@ -102,3 +137,11 @@ cdef class AdaptivePreconditioner(PreconditionerBase): def __get__(self): cdef CxxSparseMatrix smat = self.preconditioner.matrix() return get_from_sparse(smat, smat.rows(), smat.cols()) + + +cdef class BandedJacobian(PreconditionerBase): + precon_type = "banded-direct" + precon_linear_solver_type = "direct" + + cdef set_cxx_object(self): + self.preconditioner = self.pbase.get() \ No newline at end of file diff --git a/src/numerics/AdaptivePreconditioner.cpp b/src/numerics/AdaptivePreconditioner.cpp index 05145a3ca8..7cb5b2c4b9 100644 --- a/src/numerics/AdaptivePreconditioner.cpp +++ b/src/numerics/AdaptivePreconditioner.cpp @@ -44,10 +44,10 @@ void AdaptivePreconditioner::initialize(size_t networkSize) m_identity.makeCompressed(); // setting default ILUT parameters if (m_drop_tol == 0) { - setIlutDropTol(1e-10); + setIlutDropTol(1e-12); } if (m_drop_tol == 0) { - setIlutFillFactor(static_cast(m_dim) / 4); + setIlutFillFactor(static_cast(m_dim) / 2); } // update initialized status m_init = true; @@ -56,29 +56,50 @@ void AdaptivePreconditioner::initialize(size_t networkSize) void AdaptivePreconditioner::setup() { - // make into preconditioner as P = (I - gamma * J_bar) + warn_deprecated("AdaptivePreconditioner::setup", + "To be removed after Cantera 3.2. Use updatePreconditioner() instead."); updatePreconditioner(); - // compressing sparse matrix structure - m_precon_matrix.makeCompressed(); - // analyze and factorize - m_solver.compute(m_precon_matrix); - // check for errors - if (m_solver.info() != Eigen::Success) { - throw CanteraError("AdaptivePreconditioner::setup", - "error code: {}", static_cast(m_solver.info())); - } + factorize(); } void AdaptivePreconditioner::updatePreconditioner() { // set precon to jacobian m_precon_matrix.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end()); - // convert to preconditioner + // 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(mask, m_dim).cast(); + 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(); + // analyze and factorize + m_solver.compute(m_precon_matrix); + // check for errors + if (m_solver.info() != Eigen::Success) { + throw CanteraError("AdaptivePreconditioner::factorize", + "error code: {}", static_cast(m_solver.info())); + } } void AdaptivePreconditioner::prunePreconditioner() diff --git a/src/oneD/MultiJac.cpp b/src/oneD/MultiJac.cpp index f6116bee79..94f2684c62 100644 --- a/src/oneD/MultiJac.cpp +++ b/src/oneD/MultiJac.cpp @@ -50,6 +50,7 @@ void MultiJac::updateTransient(double rdt, integer* mask) for (size_t n = 0; n < m_dim; n++) { m_mat.value(n,n) = m_ssdiag[n] - mask[n]*rdt; } + factorize(); } void MultiJac::eval(double* x0, double* resid0, double rdt) diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index d37f786517..63127da893 100644 --- a/src/oneD/OneDim.cpp +++ b/src/oneD/OneDim.cpp @@ -96,11 +96,20 @@ MultiJac& OneDim::jacobian() throw CanteraError("OneDim::jacobian", "Active Jacobian is not a MultiJac"); } } + MultiNewton& OneDim::newton() { return *m_newt; } +void OneDim::setLinearSolver(shared_ptr solver) +{ + m_jac = solver; + m_jac->initialize(size()); + m_jac->setBandwidth(bandwidth()); + m_jac_ok = false; +} + void OneDim::setJacAge(int ss_age, int ts_age) { m_ss_jac_age = ss_age; @@ -209,7 +218,9 @@ void OneDim::resize() m_mask.resize(size()); // delete the current Jacobian evaluator and create a new one - m_jac = newPreconditioner("banded-direct"); + if (!m_jac) { + m_jac = newPreconditioner("banded-direct"); + } m_jac->initialize(size()); m_jac->setBandwidth(bandwidth()); m_jac_ok = false; diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 3d1c4f2635..3668706125 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -593,7 +593,7 @@ void ReactorNet::preconditionerSetup(double t, double* y, double gamma) } } // post reactor setup operations - precon->setup(); + precon->updatePreconditioner(); } void ReactorNet::updatePreconditioner(double gamma) diff --git a/test/zeroD/test_zeroD.cpp b/test/zeroD/test_zeroD.cpp index aa34f3363e..d1799f78a7 100644 --- a/test/zeroD/test_zeroD.cpp +++ b/test/zeroD/test_zeroD.cpp @@ -145,7 +145,7 @@ TEST(AdaptivePreconditionerTests, test_adaptive_precon_utils) precon.setGamma(gamma); EXPECT_NEAR(precon.gamma(), gamma, tol); // test setup and getting the matrix - precon.setup(); + precon.updatePreconditioner(); Eigen::SparseMatrix identity(testSize, testSize); identity.setIdentity(); EXPECT_TRUE(precon.matrix().isApprox(identity)); @@ -169,11 +169,11 @@ TEST(AdaptivePreconditionerTests, test_adaptive_precon_utils) testMat.setIdentity(); testMat.fill(thresh * 0.9); EXPECT_TRUE(precon.jacobian().isApprox(testMat)); - precon.setup(); + precon.updatePreconditioner(); EXPECT_TRUE(precon.matrix().isApprox(identity * (thresh * 1.1))); // reset and setup then test again precon.reset(); - precon.setup(); + precon.updatePreconditioner(); EXPECT_TRUE(precon.matrix().isApprox(identity)); } From a84390cec370c24d9a1abd3914789ac66a879b81 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 1 Jan 2025 18:51:16 -0500 Subject: [PATCH 6/9] [1D] Improve Jacobian by using better choice of perturbation size Formulas for perturbation size are based roughly on those used in CVODES. --- src/oneD/OneDim.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index 63127da893..67adb1d5cd 100644 --- a/src/oneD/OneDim.cpp +++ b/src/oneD/OneDim.cpp @@ -293,27 +293,25 @@ void OneDim::evalJacobian(double* x0) { m_jac->reset(); clock_t t0 = clock(); - double rtol = 1e-5; - double atol = sqrt(std::numeric_limits::epsilon()); + const double uround = sqrt(std::numeric_limits::epsilon()); m_work1.resize(size()); m_work2.resize(size()); eval(npos, x0, m_work1.data(), 0.0, 0); size_t ipt = 0; for (size_t j = 0; j < points(); j++) { + Domain1D* dom = pointDomain(j); size_t nv = nVars(j); for (size_t n = 0; n < nv; n++) { // perturb x(n); preserve sign(x(n)) double xsave = x0[ipt]; - double dx; - if (xsave >= 0) { - dx = xsave*rtol + atol; - } else { - dx = xsave*rtol - atol; + double dx = std::max(sqrt(uround) * fabs(xsave), + uround * (dom->rtol(n) * fabs(xsave) + dom->atol(n))); + if (xsave < 0) { + dx = -dx; } x0[ipt] = xsave + dx; - dx = x0[ipt] - xsave; - double rdx = 1.0/dx; + double rdx = 1.0 / (x0[ipt] - xsave); // calculate perturbed residual eval(j, x0, m_work2.data(), 0.0, 0); @@ -324,8 +322,10 @@ void OneDim::evalJacobian(double* x0) size_t mv = nVars(i); size_t iloc = loc(i); for (size_t m = 0; m < mv; m++) { - m_jac->setValue(m + iloc, ipt, - (m_work2[m+iloc] - m_work1[m+iloc]) * rdx); + double delta = m_work2[m+iloc] - m_work1[m+iloc]; + if (std::abs(delta) > 1e-10 || m+iloc == ipt) { + m_jac->setValue(m + iloc, ipt, delta * rdx); + } } } } From 865e1a931d0fb3507d1af61930741b44b8d64ebd Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 1 Jan 2025 21:53:51 -0500 Subject: [PATCH 7/9] Rename PreconditionerBase to SystemJacobian --- .../cantera/numerics/AdaptivePreconditioner.h | 6 +- include/cantera/numerics/Integrator.h | 8 +- include/cantera/numerics/PreconditionerBase.h | 194 +---------------- .../cantera/numerics/PreconditionerFactory.h | 33 +-- include/cantera/numerics/SystemJacobian.h | 202 ++++++++++++++++++ .../cantera/numerics/SystemJacobianFactory.h | 42 ++++ include/cantera/oneD/MultiJac.h | 4 +- include/cantera/oneD/OneDim.h | 10 +- include/cantera/zeroD/ReactorNet.h | 6 +- interfaces/cython/cantera/__init__.pxd | 2 +- interfaces/cython/cantera/_cantera.pyx | 2 +- interfaces/cython/cantera/_onedim.pxd | 6 +- interfaces/cython/cantera/_onedim.pyx | 4 +- .../{preconditioners.pxd => jacobians.pxd} | 27 ++- .../{preconditioners.pyx => jacobians.pyx} | 48 ++--- interfaces/cython/cantera/reactor.pxd | 4 +- interfaces/cython/cantera/reactor.pyx | 4 +- src/numerics/PreconditionerFactory.cpp | 42 ---- src/numerics/SystemJacobianFactory.cpp | 41 ++++ src/oneD/OneDim.cpp | 6 +- src/zeroD/ReactorNet.cpp | 2 +- test/zeroD/test_zeroD.cpp | 4 +- 22 files changed, 369 insertions(+), 328 deletions(-) create mode 100644 include/cantera/numerics/SystemJacobian.h create mode 100644 include/cantera/numerics/SystemJacobianFactory.h rename interfaces/cython/cantera/{preconditioners.pxd => jacobians.pxd} (64%) rename interfaces/cython/cantera/{preconditioners.pyx => jacobians.pyx} (77%) delete mode 100644 src/numerics/PreconditionerFactory.cpp create mode 100644 src/numerics/SystemJacobianFactory.cpp diff --git a/include/cantera/numerics/AdaptivePreconditioner.h b/include/cantera/numerics/AdaptivePreconditioner.h index 534daff662..74d3c1ac14 100644 --- a/include/cantera/numerics/AdaptivePreconditioner.h +++ b/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,7 +10,7 @@ #ifndef ADAPTIVEPRECONDITIONER_H #define ADAPTIVEPRECONDITIONER_H -#include "cantera/numerics/PreconditionerBase.h" +#include "cantera/numerics/SystemJacobian.h" #include "cantera/numerics/eigen_sparse.h" #include "cantera/base/global.h" #include @@ -23,7 +23,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 PreconditionerBase +class AdaptivePreconditioner : public SystemJacobian { public: AdaptivePreconditioner(); diff --git a/include/cantera/numerics/Integrator.h b/include/cantera/numerics/Integrator.h index 43b4df12d1..ed7fd4158b 100644 --- a/include/cantera/numerics/Integrator.h +++ b/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 preconditioner) { + virtual void setPreconditioner(shared_ptr 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 preconditioner() { + virtual shared_ptr 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 m_preconditioner; + shared_ptr m_preconditioner; //! Type of preconditioning used in applyOptions PreconditionerSide m_prec_side = PreconditionerSide::NO_PRECONDITION; // methods for DAE solvers diff --git a/include/cantera/numerics/PreconditionerBase.h b/include/cantera/numerics/PreconditionerBase.h index 955368f226..10ef1e1172 100644 --- a/include/cantera/numerics/PreconditionerBase.h +++ b/include/cantera/numerics/PreconditionerBase.h @@ -1,195 +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() {} - - virtual const string type() const = 0; - - //! 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& 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; - } - - //! Update the diagonal terms in the Jacobian by using the transient mask. - //! @param rdt Reciprocal of the time step [1/s] - //! @param mask Mask for transient terms: 1 if transient, 0 if algebraic. - virtual void updateTransient(double rdt, int* mask) { - throw NotImplementedError("PreconditionerBase::updateTransient"); - } - - //! 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. - //! @deprecated To be removed after %Cantera 3.2. Use updatePreconditioner() and - //! factorize() instead. - 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"); - }; - - //! Used to provide system bandwidth for implementations that use banded matrix - //! storage. Ignored if not needed. - virtual void setBandwidth(size_t bw) {} - - //! 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; - } - - //! Get latest status of linear solver. Zero indicates success. Meaning of non-zero - //! values is solver dependent. - virtual int info() const { - throw NotImplementedError("PreconditionerBase::info"); - } - - //! Elapsed CPU time spent computing the Jacobian elements. - double elapsedTime() const { - return m_elapsed; - } - void updateElapsed(double evalTime) { - m_elapsed += evalTime; - } - - //! Number of Jacobian evaluations. - int nEvals() const { - return m_nevals; - } - void incrementEvals() { - m_nevals++; - } - - //! Number of times 'incrementAge' has been called since the last evaluation - int age() const { - return m_age; - } - - //! Increment the Jacobian age. - void incrementAge() { - m_age++; - } - - //! Set the Jacobian age. - void setAge(int age) { - m_age = age; - } - -protected: - //! Factorize the system matrix. This method should be called at the end of - //! implementations of updatePreconditioner() and updateTransient(). - virtual void factorize() { - throw NotImplementedError("preconditionerBase::factorize"); - } - - //! 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; - - double m_elapsed = 0.0; //!< Elapsed CPU time taken to compute the Jacobian - int m_nevals = 0; //!< Number of Jacobian evaluations. - - int m_age = 100000; //!< Age of the Jacobian (times incrementAge() has been called) - - string m_precon_side = "none"; -}; - -} #endif diff --git a/include/cantera/numerics/PreconditionerFactory.h b/include/cantera/numerics/PreconditionerFactory.h index d8c1b2ac61..24cb067b59 100644 --- a/include/cantera/numerics/PreconditionerFactory.h +++ b/include/cantera/numerics/PreconditionerFactory.h @@ -3,34 +3,13 @@ // 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 PRECONDITIONER_FACTORY_H -#define PRECONDITIONER_FACTORY_H +#ifndef PRECONDITIONERFACTORY_FACTORY_H +#define PRECONDITIONERFACTORY_FACTORY_H -#include "cantera/base/FactoryBase.h" +#pragma message("warning: PreconditionerFactory.h and class PreconditionerFactory " \ + "are deprecated and will be removed after Cantera 3.2. Use " \ + "SystemJacobianFactory.h and class SystemJacobianFactory.") -namespace Cantera -{ - -class PreconditionerBase; - -//! Factory class to create preconditioner objects -class PreconditionerFactory : public Factory -{ -public: - static PreconditionerFactory* factory(); - - //! Delete preconditioner factory - void deleteFactory() override; - -private: - static PreconditionerFactory* s_factory; - static std::mutex precon_mutex; - PreconditionerFactory(); -}; - -//! Create a Preconditioner object of the specified type -shared_ptr newPreconditioner(const string& precon); - -} +#include "SystemJacobianFactory.h" #endif diff --git a/include/cantera/numerics/SystemJacobian.h b/include/cantera/numerics/SystemJacobian.h new file mode 100644 index 0000000000..5e99318efe --- /dev/null +++ b/include/cantera/numerics/SystemJacobian.h @@ -0,0 +1,202 @@ +//! @file SystemJacobian.h Declarations for class SystemJacobian + +// 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 SYSTEMJACOBIAN_H +#define SYSTEMJACOBIAN_H + +#include "cantera/base/ctexceptions.h" +#include "cantera/base/global.h" + +namespace Cantera +{ + +/** + * 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 +}; + +//! Abstract base class representing Jacobian matrices and preconditioners used in +//! nonlinear solvers. +class SystemJacobian +{ +public: + SystemJacobian() {} + + virtual ~SystemJacobian() {} + + virtual const string type() const = 0; + + //! 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("SystemJacobian::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& state) { + throw NotImplementedError("SystemJacobian::stateAdjustment"); + } + + //! Get preconditioner application side for CVODES + string preconditionerSide() const { + return m_precon_side; + } + + virtual void setPreconditionerSide(const string& preconSide) { + m_precon_side = preconSide; + } + + //! Update the diagonal terms in the Jacobian by using the transient mask. + //! @param rdt Reciprocal of the time step [1/s] + //! @param mask Mask for transient terms: 1 if transient, 0 if algebraic. + virtual void updateTransient(double rdt, int* mask) { + throw NotImplementedError("SystemJacobian::updateTransient"); + } + + //! 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("SystemJacobian::solve"); + }; + + //! Perform preconditioner specific post-reactor setup operations such as factorize. + //! @deprecated To be removed after %Cantera 3.2. Use updatePreconditioner() and + //! factorize() instead. + virtual void setup() { + throw NotImplementedError("SystemJacobian::setup"); + }; + + //! Reset parameters as needed + virtual void reset() { + throw NotImplementedError("SystemJacobian::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("SystemJacobian::initialize"); + }; + + //! Used to provide system bandwidth for implementations that use banded matrix + //! storage. Ignored if not needed. + virtual void setBandwidth(size_t bw) {} + + //! Print preconditioner contents + virtual void printPreconditioner() { + throw NotImplementedError("SystemJacobian::printPreconditioner"); + }; + + //! Transform Jacobian vector and write into + //! preconditioner, P = (I - gamma * J) + virtual void updatePreconditioner() { + throw NotImplementedError("SystemJacobian::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; + } + + //! Get latest status of linear solver. Zero indicates success. Meaning of non-zero + //! values is solver dependent. + virtual int info() const { + throw NotImplementedError("SystemJacobian::info"); + } + + //! Elapsed CPU time spent computing the Jacobian elements. + double elapsedTime() const { + return m_elapsed; + } + void updateElapsed(double evalTime) { + m_elapsed += evalTime; + } + + //! Number of Jacobian evaluations. + int nEvals() const { + return m_nevals; + } + void incrementEvals() { + m_nevals++; + } + + //! Number of times 'incrementAge' has been called since the last evaluation + int age() const { + return m_age; + } + + //! Increment the Jacobian age. + void incrementAge() { + m_age++; + } + + //! Set the Jacobian age. + void setAge(int age) { + m_age = age; + } + +protected: + //! Factorize the system matrix. This method should be called at the end of + //! implementations of updatePreconditioner() and updateTransient(). + virtual void factorize() { + throw NotImplementedError("SystemJacobian::factorize"); + } + + //! Dimension of the system + size_t m_dim; + + //! gamma value used in M = I - gamma*J + double m_gamma = 1.0; + + //! bool saying whether or not the system is initialized + bool m_init = false; + + //! Absolute tolerance of the ODE solver + double m_atol = 0; + + double m_elapsed = 0.0; //!< Elapsed CPU time taken to compute the Jacobian + int m_nevals = 0; //!< Number of Jacobian evaluations. + + int m_age = 100000; //!< Age of the Jacobian (times incrementAge() has been called) + + string m_precon_side = "none"; +}; + +//! @deprecated To be removed after Cantera 3.2. Renamed to SystemJacobian. +class PreconditionerBase : public SystemJacobian +{ + PreconditionerBase() { + warn_deprecated("PreconditionerBase::PreconditionerBase", + "To be removed after Cantera 3.2. Renamed to SystemJacobian."); + } +}; + +} +#endif diff --git a/include/cantera/numerics/SystemJacobianFactory.h b/include/cantera/numerics/SystemJacobianFactory.h new file mode 100644 index 0000000000..00b6358130 --- /dev/null +++ b/include/cantera/numerics/SystemJacobianFactory.h @@ -0,0 +1,42 @@ +//! @file SystemJacobianFactory.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 SYSTEMJACOBIANFACTORY_FACTORY_H +#define SYSTEMJACOBIANFACTORY_FACTORY_H + +#include "cantera/base/FactoryBase.h" +#include "cantera/base/global.h" + +namespace Cantera +{ + +class SystemJacobian; + +//! Factory class to create Jacobian objects for use by linear solvers +class SystemJacobianFactory : public Factory +{ +public: + static SystemJacobianFactory* factory(); + void deleteFactory() override; + +private: + static SystemJacobianFactory* s_factory; + static std::mutex jac_mutex; + SystemJacobianFactory(); +}; + +//! Create a SystemJacobian object of the specified type +shared_ptr newSystemJacobian(const string& type); + +//! @deprecated To be removed after %Cantera 3.2. Renamed to newSystemJacobian() +inline shared_ptr newPreconditioner(const string& type) { + warn_deprecated("newPreconditioner", + "To be removed after Cantera 3.2. Renamed to newSystemJacobian"); + return newSystemJacobian(type); +} + +} + +#endif diff --git a/include/cantera/oneD/MultiJac.h b/include/cantera/oneD/MultiJac.h index 2f02feb996..a3e9493c66 100644 --- a/include/cantera/oneD/MultiJac.h +++ b/include/cantera/oneD/MultiJac.h @@ -7,7 +7,7 @@ #define CT_MULTIJAC_H #include "cantera/numerics/BandMatrix.h" -#include "cantera/numerics/PreconditionerBase.h" +#include "cantera/numerics/SystemJacobian.h" #include "OneDim.h" namespace Cantera @@ -21,7 +21,7 @@ namespace Cantera * @ingroup onedUtilsGroup * @ingroup derivGroup */ -class MultiJac : public PreconditionerBase +class MultiJac : public SystemJacobian { public: //! Constructor. diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h index 2bdf560c14..53bf0db77b 100644 --- a/include/cantera/oneD/OneDim.h +++ b/include/cantera/oneD/OneDim.h @@ -10,7 +10,7 @@ #include "Domain1D.h" #include "MultiJac.h" -#include "cantera/numerics/PreconditionerBase.h" +#include "cantera/numerics/SystemJacobian.h" namespace Cantera { @@ -44,7 +44,7 @@ class OneDim //! @ingroup derivGroup MultiJac& jacobian(); - shared_ptr getJacobian() { + shared_ptr getJacobian() { return m_jac; } @@ -53,10 +53,10 @@ class OneDim //! Set the linear solver used to hold the Jacobian matrix and solve linear systems //! as part of each Newton iteration. The default is a direct, banded solver. - void setLinearSolver(shared_ptr solver); + void setLinearSolver(shared_ptr solver); //! Get the type of the linear solver being used. - shared_ptr linearSolver() const { return m_jac; } + shared_ptr linearSolver() const { return m_jac; } /** * Solve F(x) = 0, where F(x) is the multi-domain residual function. @@ -386,7 +386,7 @@ class OneDim shared_ptr> m_state; //!< Solution vector - shared_ptr m_jac; //!< Jacobian evaluator + shared_ptr m_jac; //!< Jacobian evaluator unique_ptr m_newt; //!< Newton iterator double m_rdt = 0.0; //!< reciprocal of time step bool m_jac_ok = false; //!< if true, Jacobian is current diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index f1ffe5b249..fd8e5740e7 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -15,7 +15,7 @@ namespace Cantera class Array2D; class Integrator; -class PreconditionerBase; +class SystemJacobian; //! A class representing a network of connected reactors. /*! @@ -44,7 +44,7 @@ class ReactorNet : public FuncEval //! Set preconditioner used by the linear solver //! @param preconditioner preconditioner object used for the linear solver - void setPreconditioner(shared_ptr preconditioner); + void setPreconditioner(shared_ptr preconditioner); //! Set the initial value of the independent variable (typically time). //! Default = 0.0 s. Restarts integration from this value using the current mixture @@ -341,7 +341,7 @@ class ReactorNet : public FuncEval double m_rtolsens = 1.0e-4; double m_atols = 1.0e-15; double m_atolsens = 1.0e-6; - shared_ptr m_precon; + shared_ptr m_precon; string m_linearSolverType; //! Maximum integrator internal timestep. Default of 0.0 means infinity. diff --git a/interfaces/cython/cantera/__init__.pxd b/interfaces/cython/cantera/__init__.pxd index 1a8aa9a8e6..a513264e0c 100644 --- a/interfaces/cython/cantera/__init__.pxd +++ b/interfaces/cython/cantera/__init__.pxd @@ -10,7 +10,7 @@ from cantera.delegator cimport * from cantera.func1 cimport * from cantera.kinetics cimport * from cantera.mixture cimport * -from cantera.preconditioners cimport * +from cantera.jacobians cimport * from cantera.reaction cimport * from cantera.reactionpath cimport * from cantera.reactor cimport * diff --git a/interfaces/cython/cantera/_cantera.pyx b/interfaces/cython/cantera/_cantera.pyx index 24509c54da..1fb7eaff93 100644 --- a/interfaces/cython/cantera/_cantera.pyx +++ b/interfaces/cython/cantera/_cantera.pyx @@ -43,7 +43,7 @@ from .transport import * from .units import * from .yamlwriter import * from .constants import * -from .preconditioners import * +from .jacobians import * # Custom finder/loader no longer needed, so remove it sys.meta_path.pop() diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index 834af4228f..b97c0fc52d 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -9,7 +9,7 @@ from .solutionbase cimport * from .kinetics cimport * from .func1 cimport * from .thermo cimport * -from .preconditioners cimport * +from .jacobians cimport * cdef extern from "cantera/oneD/DomainFactory.h" namespace "Cantera": @@ -115,8 +115,8 @@ cdef extern from "cantera/oneD/Sim1D.h": void restoreSteadySolution() except +translate_exception void setMaxTimeStepCount(int) int maxTimeStepCount() - void setLinearSolver(shared_ptr[CxxPreconditionerBase]) except +translate_exception - shared_ptr[CxxPreconditionerBase] linearSolver() + void setLinearSolver(shared_ptr[CxxSystemJacobian]) except +translate_exception + shared_ptr[CxxSystemJacobian] linearSolver() void getInitialSoln() except +translate_exception void solve(int, cbool) except +translate_exception void refine(int) except +translate_exception diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index 15fc04b784..97f7e98084 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -1019,10 +1019,10 @@ cdef class Sim1D: systems as part of each Newton iteration. The default is a banded, direct solver. """ - return PreconditionerBase.wrap(self.sim.linearSolver()) + return SystemJacobian.wrap(self.sim.linearSolver()) @linear_solver.setter - def linear_solver(self, PreconditionerBase precon): + def linear_solver(self, SystemJacobian precon): self.sim.setLinearSolver(precon.pbase) def set_initial_guess(self, *args, **kwargs): diff --git a/interfaces/cython/cantera/preconditioners.pxd b/interfaces/cython/cantera/jacobians.pxd similarity index 64% rename from interfaces/cython/cantera/preconditioners.pxd rename to interfaces/cython/cantera/jacobians.pxd index e66a5cd3a1..66fad55989 100644 --- a/interfaces/cython/cantera/preconditioners.pxd +++ b/interfaces/cython/cantera/jacobians.pxd @@ -7,16 +7,16 @@ from .ctcxx cimport * from .kinetics cimport CxxSparseMatrix -cdef extern from "cantera/numerics/PreconditionerBase.h" namespace "Cantera": - cdef cppclass CxxPreconditionerBase "Cantera::PreconditionerBase": - CxxPreconditionerBase() +cdef extern from "cantera/numerics/SystemJacobian.h" namespace "Cantera": + cdef cppclass CxxSystemJacobian "Cantera::SystemJacobian": + CxxSystemJacobian() string preconditionerSide() string type() void setPreconditionerSide(string) except +translate_exception cdef extern from "cantera/numerics/AdaptivePreconditioner.h" namespace "Cantera": cdef cppclass CxxAdaptivePreconditioner "Cantera::AdaptivePreconditioner" \ - (CxxPreconditionerBase): + (CxxSystemJacobian): CxxAdaptivePreconditioner() except +translate_exception void setThreshold(double threshold) double threshold() @@ -28,24 +28,23 @@ cdef extern from "cantera/numerics/AdaptivePreconditioner.h" namespace "Cantera" CxxSparseMatrix matrix() except +translate_exception cdef extern from "cantera/oneD/MultiJac.h" namespace "Cantera": - cdef cppclass CxxMultiJac "Cantera::MultiJac" (CxxPreconditionerBase): + cdef cppclass CxxMultiJac "Cantera::MultiJac" (CxxSystemJacobian): CxxMultiJac() except +translate_exception -cdef extern from "cantera/numerics/PreconditionerFactory.h" namespace "Cantera": - cdef shared_ptr[CxxPreconditionerBase] newPreconditioner(string) except\ +cdef extern from "cantera/numerics/SystemJacobianFactory.h" namespace "Cantera": + cdef shared_ptr[CxxSystemJacobian] newSystemJacobian(string) except\ +translate_exception -cdef class PreconditionerBase: +cdef class SystemJacobian: @staticmethod - cdef wrap(shared_ptr[CxxPreconditionerBase]) + cdef wrap(shared_ptr[CxxSystemJacobian]) cdef set_cxx_object(self) - cdef shared_ptr[CxxPreconditionerBase] pbase + cdef shared_ptr[CxxSystemJacobian] pbase -cdef class AdaptivePreconditioner(PreconditionerBase): +cdef class AdaptivePreconditioner(SystemJacobian): cdef set_cxx_object(self) cdef CxxAdaptivePreconditioner* preconditioner -cdef class BandedJacobian(PreconditionerBase): +cdef class BandedJacobian(SystemJacobian): cdef set_cxx_object(self) - cdef CxxMultiJac* preconditioner - + cdef CxxMultiJac* jacobian diff --git a/interfaces/cython/cantera/preconditioners.pyx b/interfaces/cython/cantera/jacobians.pyx similarity index 77% rename from interfaces/cython/cantera/preconditioners.pyx rename to interfaces/cython/cantera/jacobians.pyx index d7df649719..eae5b6ca5c 100644 --- a/interfaces/cython/cantera/preconditioners.pyx +++ b/interfaces/cython/cantera/jacobians.pyx @@ -5,43 +5,43 @@ from ._utils cimport stringify, pystr, py_to_anymap, anymap_to_py from .kinetics cimport get_from_sparse # dictionary to store reaction rate classes -cdef dict _preconditioner_class_registry = {} +cdef dict _class_registry = {} -cdef class PreconditionerBase: +cdef class SystemJacobian: """ - Common base class for preconditioners. + Common base class for Jacobian matrices used in the solution of nonlinear systems. """ - precon_type = "PreconditionerBase" - precon_linear_solver_type = "GMRES" + _type = "SystemJacobian" + linear_solver_type = "GMRES" def __cinit__(self, *args, init=True, **kwargs): if init: - self.pbase = newPreconditioner(stringify(self.precon_type)) + self.pbase = newSystemJacobian(stringify(self._type)) @staticmethod - cdef wrap(shared_ptr[CxxPreconditionerBase] base): + cdef wrap(shared_ptr[CxxSystemJacobian] base): """ - Wrap a C++ PreconditionerBase object with a Python object of the correct + Wrap a C++ SystemJacobian object with a Python object of the correct derived type. """ # Ensure all derived types are registered - if not _preconditioner_class_registry: + if not _class_registry: def register_subclasses(cls): for c in cls.__subclasses__(): - _preconditioner_class_registry[getattr(c, "precon_type")] = c + _class_registry[getattr(c, "_type")] = c register_subclasses(c) # Update global class registry - register_subclasses(PreconditionerBase) + register_subclasses(SystemJacobian) cxx_type = pystr(base.get().type()) - cls = _preconditioner_class_registry.get(cxx_type, PreconditionerBase) - cdef PreconditionerBase precon - precon = cls(init=False) - precon.pbase = base - precon.set_cxx_object() - return precon + cls = _class_registry.get(cxx_type, SystemJacobian) + cdef SystemJacobian jac + jac = cls(init=False) + jac.pbase = base + jac.set_cxx_object() + return jac cdef set_cxx_object(self): pass @@ -58,9 +58,9 @@ cdef class PreconditionerBase: def __set__(self, side): self.pbase.get().setPreconditionerSide(stringify(side)) -cdef class AdaptivePreconditioner(PreconditionerBase): - precon_type = "Adaptive" - precon_linear_solver_type = "GMRES" +cdef class AdaptivePreconditioner(SystemJacobian): + _type = "Adaptive" + linear_solver_type = "GMRES" def __cinit__(self, *args, **kwargs): self.preconditioner = (self.pbase.get()) @@ -139,9 +139,9 @@ cdef class AdaptivePreconditioner(PreconditionerBase): return get_from_sparse(smat, smat.rows(), smat.cols()) -cdef class BandedJacobian(PreconditionerBase): - precon_type = "banded-direct" - precon_linear_solver_type = "direct" +cdef class BandedJacobian(SystemJacobian): + _type = "banded-direct" + linear_solver_type = "direct" cdef set_cxx_object(self): - self.preconditioner = self.pbase.get() \ No newline at end of file + self.jacobian = self.pbase.get() diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd index c06a9e0b8c..670cd43731 100644 --- a/interfaces/cython/cantera/reactor.pxd +++ b/interfaces/cython/cantera/reactor.pxd @@ -7,7 +7,7 @@ from .ctcxx cimport * from .kinetics cimport * from .func1 cimport * -from .preconditioners cimport * +from .jacobians cimport * cdef extern from "cantera/numerics/Integrator.h" namespace "Cantera": # SUNDIALS integrator @@ -204,7 +204,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": string sensitivityParameterName(size_t) except +translate_exception void setLinearSolverType(string& integratorType) except +translate_exception string linearSolverType() - void setPreconditioner(shared_ptr[CxxPreconditionerBase] preconditioner) + void setPreconditioner(shared_ptr[CxxSystemJacobian] preconditioner) void setDerivativeSettings(CxxAnyMap&) CxxAnyMap solverStats() except +translate_exception diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index 39a05f19d4..e5b9f5c502 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -1982,11 +1982,11 @@ cdef class ReactorNet: property preconditioner: """Preconditioner associated with integrator""" - def __set__(self, PreconditionerBase precon): + def __set__(self, SystemJacobian precon): # set preconditioner self.net.setPreconditioner(precon.pbase) # set problem type as default of preconditioner - self.linear_solver_type = precon.precon_linear_solver_type + self.linear_solver_type = precon.linear_solver_type property linear_solver_type: """ diff --git a/src/numerics/PreconditionerFactory.cpp b/src/numerics/PreconditionerFactory.cpp deleted file mode 100644 index 5911eb0c1d..0000000000 --- a/src/numerics/PreconditionerFactory.cpp +++ /dev/null @@ -1,42 +0,0 @@ -//! @file PreconditionerFactory.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/PreconditionerFactory.h" -#include "cantera/numerics/AdaptivePreconditioner.h" -#include "cantera/oneD/MultiJac.h" - -namespace Cantera -{ - -PreconditionerFactory* PreconditionerFactory::factory() { - std::unique_lock lock(precon_mutex); - if (!s_factory) { - s_factory = new PreconditionerFactory; - } - return s_factory; -}; - -//! Delete preconditioner factory -void PreconditionerFactory::deleteFactory() { - std::unique_lock lock(precon_mutex); - delete s_factory; - s_factory = 0; -}; - -PreconditionerFactory* PreconditionerFactory::s_factory = 0; -std::mutex PreconditionerFactory::precon_mutex; - -PreconditionerFactory::PreconditionerFactory() -{ - reg("Adaptive", []() { return new AdaptivePreconditioner(); }); - reg("banded-direct", []() { return new MultiJac(); }); -} - -shared_ptr newPreconditioner(const string& precon) -{ - return shared_ptr(PreconditionerFactory::factory()->create(precon)); -}; - -} diff --git a/src/numerics/SystemJacobianFactory.cpp b/src/numerics/SystemJacobianFactory.cpp new file mode 100644 index 0000000000..cc5170aeb5 --- /dev/null +++ b/src/numerics/SystemJacobianFactory.cpp @@ -0,0 +1,41 @@ +//! @file SystemJacobianFactory.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/SystemJacobianFactory.h" +#include "cantera/numerics/AdaptivePreconditioner.h" +#include "cantera/oneD/MultiJac.h" + +namespace Cantera +{ + +SystemJacobianFactory* SystemJacobianFactory::factory() { + std::unique_lock lock(jac_mutex); + if (!s_factory) { + s_factory = new SystemJacobianFactory; + } + return s_factory; +}; + +void SystemJacobianFactory::deleteFactory() { + std::unique_lock lock(jac_mutex); + delete s_factory; + s_factory = 0; +}; + +SystemJacobianFactory* SystemJacobianFactory::s_factory = 0; +std::mutex SystemJacobianFactory::jac_mutex; + +SystemJacobianFactory::SystemJacobianFactory() +{ + reg("Adaptive", []() { return new AdaptivePreconditioner(); }); + reg("banded-direct", []() { return new MultiJac(); }); +} + +shared_ptr newSystemJacobian(const string& type) +{ + return shared_ptr(SystemJacobianFactory::factory()->create(type)); +}; + +} diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index 67adb1d5cd..f16edf28c1 100644 --- a/src/oneD/OneDim.cpp +++ b/src/oneD/OneDim.cpp @@ -7,7 +7,7 @@ #include "cantera/numerics/Func1.h" #include "cantera/oneD/MultiNewton.h" #include "cantera/base/AnyMap.h" -#include "cantera/numerics/PreconditionerFactory.h" +#include "cantera/numerics/SystemJacobianFactory.h" #include #include @@ -102,7 +102,7 @@ MultiNewton& OneDim::newton() return *m_newt; } -void OneDim::setLinearSolver(shared_ptr solver) +void OneDim::setLinearSolver(shared_ptr solver) { m_jac = solver; m_jac->initialize(size()); @@ -219,7 +219,7 @@ void OneDim::resize() // delete the current Jacobian evaluator and create a new one if (!m_jac) { - m_jac = newPreconditioner("banded-direct"); + m_jac = newSystemJacobian("banded-direct"); } m_jac->initialize(size()); m_jac->setBandwidth(bandwidth()); diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 3668706125..e879090cd5 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -155,7 +155,7 @@ void ReactorNet::setLinearSolverType(const string& linSolverType) m_integrator_init = false; } -void ReactorNet::setPreconditioner(shared_ptr preconditioner) +void ReactorNet::setPreconditioner(shared_ptr preconditioner) { m_precon = preconditioner; m_integrator_init = false; diff --git a/test/zeroD/test_zeroD.cpp b/test/zeroD/test_zeroD.cpp index d1799f78a7..78e5780ce1 100644 --- a/test/zeroD/test_zeroD.cpp +++ b/test/zeroD/test_zeroD.cpp @@ -4,7 +4,7 @@ #include "cantera/zerodim.h" #include "cantera/base/Interface.h" #include "cantera/numerics/eigen_sparse.h" -#include "cantera/numerics/PreconditionerFactory.h" +#include "cantera/numerics/SystemJacobianFactory.h" #include "cantera/numerics/AdaptivePreconditioner.h" using namespace Cantera; @@ -188,7 +188,7 @@ TEST(AdaptivePreconditionerTests, test_precon_solver_stats) ReactorNet network; network.addReactor(reactor); // setup preconditioner - shared_ptr precon_ptr = newPreconditioner("Adaptive"); + shared_ptr precon_ptr = newSystemJacobian("Adaptive"); network.setPreconditioner(precon_ptr); EXPECT_THROW(network.step(), CanteraError); // take a step From 691c4931c10058fd860efea037636a9f744e0832 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 1 Jan 2025 22:35:26 -0500 Subject: [PATCH 8/9] [Numerics] Create intermediate class for sparse Jacobians --- .../cantera/numerics/AdaptivePreconditioner.h | 41 +--------- .../cantera/numerics/EigenSparseJacobian.h | 54 +++++++++++++ interfaces/cython/cantera/_onedim.pyx | 2 +- interfaces/cython/cantera/jacobians.pxd | 23 ++++-- interfaces/cython/cantera/jacobians.pyx | 63 ++++++++------- interfaces/cython/cantera/reactor.pyx | 2 +- src/numerics/AdaptivePreconditioner.cpp | 68 ++-------------- src/numerics/EigenSparseJacobian.cpp | 77 +++++++++++++++++++ 8 files changed, 194 insertions(+), 136 deletions(-) create mode 100644 include/cantera/numerics/EigenSparseJacobian.h create mode 100644 src/numerics/EigenSparseJacobian.cpp diff --git a/include/cantera/numerics/AdaptivePreconditioner.h b/include/cantera/numerics/AdaptivePreconditioner.h index 74d3c1ac14..9c3948a277 100644 --- a/include/cantera/numerics/AdaptivePreconditioner.h +++ b/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 +#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 @@ class AdaptivePreconditioner : public SystemJacobian void initialize(size_t networkSize) override; void reset() override { - m_precon_matrix.setZero(); + m_matrix.setZero(); m_jac_trips.clear(); }; @@ -39,10 +36,7 @@ class AdaptivePreconditioner : public SystemJacobian 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& state) override; - void updatePreconditioner() override; - void updateTransient(double rdt, int* mask) override; int info() const override { return static_cast(m_solver.info()); @@ -51,20 +45,6 @@ class AdaptivePreconditioner : public SystemJacobian //! Prune preconditioner elements void prunePreconditioner(); - //! Return semi-analytical Jacobian of an AdaptivePreconditioner object. - //! @ingroup derivGroup - Eigen::SparseMatrix jacobian() { - Eigen::SparseMatrix 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 matrix() { - updatePreconditioner(); - return m_precon_matrix; - } - //! Get the threshold value for setting elements double threshold() { return m_threshold; } @@ -95,12 +75,6 @@ class AdaptivePreconditioner : public SystemJacobian 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 @@ class AdaptivePreconditioner : public SystemJacobian //! ILUT drop tolerance double m_drop_tol = 0; - //! Vector of triples representing the jacobian used in preconditioning - vector> m_jac_trips; - - //! Storage of appropriately sized identity matrix for making the preconditioner - Eigen::SparseMatrix m_identity; - - //! Container that is the sparse preconditioner - Eigen::SparseMatrix m_precon_matrix; - //! Solver used in solving the linear system Eigen::IncompleteLUT m_solver; diff --git a/include/cantera/numerics/EigenSparseJacobian.h b/include/cantera/numerics/EigenSparseJacobian.h new file mode 100644 index 0000000000..8760322f26 --- /dev/null +++ b/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 jacobian(); + + //! Return the internal preconditioner matrix + Eigen::SparseMatrix 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> m_jac_trips; + + //! Storage of appropriately sized identity matrix for making the preconditioner + Eigen::SparseMatrix m_identity; + + //! Container that is the sparse preconditioner + Eigen::SparseMatrix m_matrix; +}; + +} + +#endif diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index 97f7e98084..5a9ed1b9e7 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/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): """ diff --git a/interfaces/cython/cantera/jacobians.pxd b/interfaces/cython/cantera/jacobians.pxd index 66fad55989..8656b45fa8 100644 --- a/interfaces/cython/cantera/jacobians.pxd +++ b/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 EigenSparseJacobian(SystemJacobian): + cdef set_cxx_object(self) + cdef CxxEigenSparseJacobian* sparse_jac -cdef class AdaptivePreconditioner(SystemJacobian): +cdef class AdaptivePreconditioner(EigenSparseJacobian): cdef set_cxx_object(self) - cdef CxxAdaptivePreconditioner* preconditioner + cdef CxxAdaptivePreconditioner* adaptive cdef class BandedJacobian(SystemJacobian): cdef set_cxx_object(self) - cdef CxxMultiJac* jacobian + cdef CxxMultiJac* band_jac diff --git a/interfaces/cython/cantera/jacobians.pyx b/interfaces/cython/cantera/jacobians.pyx index eae5b6ca5c..9685da7a31 100644 --- a/interfaces/cython/cantera/jacobians.pyx +++ b/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 EigenSparseJacobian(SystemJacobian): + _type = "EigenSparse" + + cdef set_cxx_object(self): + self.sparse_jac = self._base.get() + + def print_contents(self): + self.sparse_jac.printPreconditioner() -cdef class AdaptivePreconditioner(SystemJacobian): + 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 = (self.pbase.get()) - cdef set_cxx_object(self): - self.preconditioner = self.pbase.get() + self.adaptive = 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) - - def __get__(self): - return self.preconditioner.ilutDropTol() - - def print_contents(self): - self.preconditioner.printPreconditioner() + self.adaptive.setIlutDropTol(val) - 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 = self.pbase.get() + self.band_jac = self._base.get() diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index e5b9f5c502..1da8be2e98 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/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 diff --git a/src/numerics/AdaptivePreconditioner.cpp b/src/numerics/AdaptivePreconditioner.cpp index 7cb5b2c4b9..67aa8db268 100644 --- a/src/numerics/AdaptivePreconditioner.cpp +++ b/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(row), static_cast(col), value); -} - void AdaptivePreconditioner::stateAdjustment(vector& 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& 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(mask, m_dim).cast(); - 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::InnerIterator it(m_precon_matrix, k); it; + for (int k=0; k::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 jacobian(m_dim, m_dim); - jacobian.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end()); - ss << Eigen::MatrixXd(jacobian); - writelog(ss.str()); -} - } diff --git a/src/numerics/EigenSparseJacobian.cpp b/src/numerics/EigenSparseJacobian.cpp new file mode 100644 index 0000000000..6d468adcee --- /dev/null +++ b/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(row), static_cast(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(mask, m_dim).cast(); + m_matrix -= rdt * diag.matrix().asDiagonal(); + factorize(); +} + +Eigen::SparseMatrix EigenSparseJacobian::jacobian() +{ + Eigen::SparseMatrix 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 jacobian(m_dim, m_dim); + jacobian.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end()); + ss << Eigen::MatrixXd(jacobian); + writelog(ss.str()); +} + +} From bb71f0a097f26fc48058fc44a3789e564fe6009d Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Thu, 2 Jan 2025 11:21:57 -0500 Subject: [PATCH 9/9] [Numerics] Add EigenSparseDirectJacobian implementation --- .../cantera/numerics/AdaptivePreconditioner.h | 6 ---- .../numerics/EigenSparseDirectJacobian.h | 28 +++++++++++++++ .../cantera/numerics/EigenSparseJacobian.h | 1 + interfaces/cython/cantera/jacobians.pxd | 3 ++ interfaces/cython/cantera/jacobians.pyx | 6 +++- src/numerics/EigenSparseDirectJacobian.cpp | 34 +++++++++++++++++++ src/numerics/EigenSparseJacobian.cpp | 6 ++++ src/numerics/SystemJacobianFactory.cpp | 2 ++ 8 files changed, 79 insertions(+), 7 deletions(-) create mode 100644 include/cantera/numerics/EigenSparseDirectJacobian.h create mode 100644 src/numerics/EigenSparseDirectJacobian.cpp diff --git a/include/cantera/numerics/AdaptivePreconditioner.h b/include/cantera/numerics/AdaptivePreconditioner.h index 9c3948a277..e5d11703d5 100644 --- a/include/cantera/numerics/AdaptivePreconditioner.h +++ b/include/cantera/numerics/AdaptivePreconditioner.h @@ -26,12 +26,6 @@ class AdaptivePreconditioner : public EigenSparseJacobian 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; diff --git a/include/cantera/numerics/EigenSparseDirectJacobian.h b/include/cantera/numerics/EigenSparseDirectJacobian.h new file mode 100644 index 0000000000..2b2dbedf3b --- /dev/null +++ b/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> m_solver; +}; + +} + +#endif diff --git a/include/cantera/numerics/EigenSparseJacobian.h b/include/cantera/numerics/EigenSparseJacobian.h index 8760322f26..493932ca10 100644 --- a/include/cantera/numerics/EigenSparseJacobian.h +++ b/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; diff --git a/interfaces/cython/cantera/jacobians.pxd b/interfaces/cython/cantera/jacobians.pxd index 8656b45fa8..f27cfd04f7 100644 --- a/interfaces/cython/cantera/jacobians.pxd +++ b/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(SystemJacobian): + pass + cdef class AdaptivePreconditioner(EigenSparseJacobian): cdef set_cxx_object(self) cdef CxxAdaptivePreconditioner* adaptive diff --git a/interfaces/cython/cantera/jacobians.pyx b/interfaces/cython/cantera/jacobians.pyx index 9685da7a31..5966c06a8b 100644 --- a/interfaces/cython/cantera/jacobians.pyx +++ b/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 = self._base.get() @@ -81,6 +81,10 @@ cdef class EigenSparseJacobian(SystemJacobian): return get_from_sparse(smat, smat.rows(), smat.cols()) +cdef class EigenSparseDirectJacobian(SystemJacobian): + _type = "eigen-sparse-direct" + + cdef class AdaptivePreconditioner(EigenSparseJacobian): _type = "Adaptive" linear_solver_type = "GMRES" diff --git a/src/numerics/EigenSparseDirectJacobian.cpp b/src/numerics/EigenSparseDirectJacobian.cpp new file mode 100644 index 0000000000..863ed808eb --- /dev/null +++ b/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(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(m_solver.info())); + } +} + +} diff --git a/src/numerics/EigenSparseJacobian.cpp b/src/numerics/EigenSparseJacobian.cpp index 6d468adcee..6c93a4bdb1 100644 --- a/src/numerics/EigenSparseJacobian.cpp +++ b/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(row), static_cast(col), value); diff --git a/src/numerics/SystemJacobianFactory.cpp b/src/numerics/SystemJacobianFactory.cpp index cc5170aeb5..dfdb102173 100644 --- a/src/numerics/SystemJacobianFactory.cpp +++ b/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 newSystemJacobian(const string& type)