diff --git a/include/cantera/numerics/DenseMatrix.h b/include/cantera/numerics/DenseMatrix.h index 409041ce33..32afda839a 100644 --- a/include/cantera/numerics/DenseMatrix.h +++ b/include/cantera/numerics/DenseMatrix.h @@ -155,6 +155,8 @@ class DenseMatrix : public Array2D friend int invert(DenseMatrix& A, int nn); }; +int factor(DenseMatrix& A); +int solveFactored(DenseMatrix& A, double* b, size_t nrhs=1, size_t ldb=0); //! Solve Ax = b. Array b is overwritten on exit with x. /*! diff --git a/include/cantera/numerics/Newton.h b/include/cantera/numerics/Newton.h new file mode 100644 index 0000000000..596b04a1d9 --- /dev/null +++ b/include/cantera/numerics/Newton.h @@ -0,0 +1,116 @@ +//! @file Newton.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 CT_NEWTON_H +#define CT_NEWTON_H + +#include "FuncEval.h" +#include "DenseMatrix.h" + +namespace Cantera +{ + +struct Configuration { + vector_fp rtol; + vector_fp atol; + double convtol; + double dt; + size_t jac_maxage; + double jac_rtol; + double jac_atol; +}; + +/** + * A Newton solver. + */ +class Newton +{ +public: + Newton(FuncEval& func); + virtual ~Newton() {}; + Newton(const Newton&) = delete; + Newton& operator=(const Newton&) = delete; + + size_t size() { + return m_nv; + } + + //! Compute the undamped Newton step. The residual function is evaluated + //! at `x`, but the Jacobian is not recomputed. + void step(doublereal* x, doublereal* step); + + //! Compute the weighted 2-norm of `step`. + doublereal weightedNorm(const doublereal* x, const doublereal* step) const; + + int hybridSolve(); + + int solve(double* x); + + //TODO: implement get methods + //nice implementation for steady vs transient below + //! Relative tolerance of the nth component. + // doublereal rtol(size_t n) { + // return (m_rdt == 0.0 ? m_rtol_ss[n] : m_rtol_ts[n]); + // } + + //! Set upper and lower bounds on the nth component + void setBounds(size_t n, doublereal lower, doublereal upper) { + m_lower_bounds[n] = lower; + m_upper_bounds[n] = upper; + } + + void setConstants(vector_int constantComponents) { + m_constantComponents = constantComponents; + } + + void evalJacobian(doublereal* x, doublereal* xdot); + + void getSolution(double* x) { + for (size_t i = 0; i < m_nv; i++) { + x[i] = m_x[i]; + } + } + +protected: + FuncEval* m_residfunc; + + //! number of variables + size_t m_nv; + + DenseMatrix m_jacobian, m_jacFactored; + size_t m_jacAge; + + //! work arrays of size #m_nv used in solve(). + vector_fp m_x, m_x1, m_stp, m_stp1; + + vector_fp m_upper_bounds, m_lower_bounds; + + vector_fp m_xlast, m_xsave; + + //! the indexes of any constant variables + vector_int m_constantComponents; + + //! current timestep reciprocal + doublereal m_rdt = 0; + + Configuration m_directsolve_config; + Configuration m_timestep_config; + Configuration* m_config; +}; + +// //! Returns the weighted Root Mean Square Deviation given a vector of residuals and +// // vectors of the corresponding weights and absolute tolerances for each component. +// double weightedRMS(vector_fp residuals, vector_fp weights, vector_fp atols) { +// size_t n = residuals.size(); +// double square = 0.0; +// for (size_t i = 0; i < n; i++) { +// square += pow(residuals[i]/(weights[i] + atols[i]), 2); +// } +// return sqrt(square/n); +// } + +} + +#endif diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 75e81daae5..cdb0dc0fd0 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -253,6 +253,9 @@ class ReactorNet : public FuncEval //! Retrieve absolute step size limits during advance bool getAdvanceLimits(double* limits); + //! Advance the reactor network to steady state. + double solveSteady(); + protected: //! Estimate a future state based on current derivatives. diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 4c80d496e6..c1e9ca9626 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -987,6 +987,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": double sensitivity(string&, size_t, int) except +translate_exception size_t nparams() string sensitivityParameterName(size_t) except +translate_exception + double solveSteady() except +translate_exception cdef extern from "cantera/zeroD/ReactorDelegator.h" namespace "Cantera": cdef cppclass CxxReactorAccessor "Cantera::ReactorAccessor": diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index a18a107992..00beba56e9 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -1463,6 +1463,13 @@ cdef class ReactorNet: if return_residuals: return residuals[:step + 1] + def solve_steady(self): + """ + Advance the reactor network to steady state via direct solution of + the ODE governing equations. + """ + return self.net.solveSteady() + def __reduce__(self): raise NotImplementedError('ReactorNet object is not picklable') diff --git a/src/numerics/DenseMatrix.cpp b/src/numerics/DenseMatrix.cpp index 2d5b014e01..b453957fe3 100644 --- a/src/numerics/DenseMatrix.cpp +++ b/src/numerics/DenseMatrix.cpp @@ -140,6 +140,94 @@ vector_int& DenseMatrix::ipiv() return m_ipiv; } +int factor(DenseMatrix& A) +{ + int info = 0; + #if CT_USE_LAPACK + ct_dgetrf(A.nRows(), A.nColumns(), A.ptrColumn(0), + A.nRows(), &A.ipiv()[0], info); + if (info > 0) { + if (A.m_printLevel) { + writelogf("solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d U(i,i) is exactly zero. The factorization has" + " been completed, but the factor U is exactly singular, and division by zero will occur if " + "it is used to solve a system of equations.\n", info); + } + if (!A.m_useReturnErrorCode) { + throw CanteraError("solve(DenseMatrix& A, double* b)", + "DGETRF returned INFO = {}. U(i,i) is exactly zero. The factorization has" + " been completed, but the factor U is exactly singular, and division by zero will occur if " + "it is used to solve a system of equations.", info); + } + return info; + } else if (info < 0) { + if (A.m_printLevel) { + writelogf("solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d. The argument i has an illegal value\n", info); + } + + throw CanteraError("solve(DenseMatrix& A, double* b)", + "DGETRF returned INFO = {}. The argument i has an illegal value", info); + } + #else + MappedMatrix Am(&A(0,0), A.nRows(), A.nColumns()); + #ifdef NDEBUG + auto lu = Am.partialPivLu(); + #else + auto lu = Am.fullPivLu(); + if (lu.nonzeroPivots() < static_cast(A.nColumns())) { + throw CanteraError("solve(DenseMatrix& A, double* b)", + "Matrix appears to be rank-deficient: non-zero pivots = {}; columns = {}", + lu.nonzeroPivots(), A.nColumns()); + } + #endif + #endif + return info; +} + +int solveFactored(DenseMatrix& A, double* b, size_t nrhs, size_t ldb) +{ + if (A.nColumns() != A.nRows()) { + if (A.m_printLevel) { + writelogf("solve(DenseMatrix& A, double* b): Can only solve a square matrix\n"); + } + throw CanteraError("solve(DenseMatrix& A, double* b)", "Can only solve a square matrix"); + } + + int info = 0; + if (ldb == 0) { + ldb = A.nColumns(); + } + #if CT_USE_LAPACK + ct_dgetrs(ctlapack::NoTranspose, A.nRows(), nrhs, A.ptrColumn(0), + A.nRows(), &A.ipiv()[0], b, ldb, info); + if (info != 0) { + if (A.m_printLevel) { + writelog("solve(DenseMatrix& A, double* b): DGETRS returned INFO = {}\n", info); + } + if (info < 0 || !A.m_useReturnErrorCode) { + throw CanteraError("solve(DenseMatrix& A, double* b)", + "DGETRS returned INFO = {}", info); + } + } + #else + MappedMatrix Am(&A(0,0), A.nRows(), A.nColumns()); + #ifdef NDEBUG + auto lu = Am.partialPivLu(); + #else + auto lu = Am.fullPivLu(); + if (lu.nonzeroPivots() < static_cast(A.nColumns())) { + throw CanteraError("solve(DenseMatrix& A, double* b)", + "Matrix appears to be rank-deficient: non-zero pivots = {}; columns = {}", + lu.nonzeroPivots(), A.nColumns()); + } + #endif + for (size_t i = 0; i < nrhs; i++) { + MappedVector bm(b + ldb*i, A.nColumns()); + bm = lu.solve(bm); + } + #endif + return info; +} + int solve(DenseMatrix& A, double* b, size_t nrhs, size_t ldb) { if (A.nColumns() != A.nRows()) { diff --git a/src/numerics/Newton.cpp b/src/numerics/Newton.cpp new file mode 100644 index 0000000000..bf0fd95f80 --- /dev/null +++ b/src/numerics/Newton.cpp @@ -0,0 +1,264 @@ +//! @file Newton.cpp: A damped & bounded quasi-Newton solver + +// 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/Newton.h" +#include "cantera/base/utilities.h" + +using namespace std; + +namespace Cantera +{ + +// constants +const double damp_factor = sqrt(2.0); +const double damp_min = 0.1; + +Newton::Newton(FuncEval& func) { + m_residfunc = &func; + m_nv = m_residfunc->neq(); + m_x.resize(m_nv); + m_x1.resize(m_nv); + m_stp.resize(m_nv); + m_stp1.resize(m_nv); + m_upper_bounds.resize(m_nv, 0.0); + m_lower_bounds.resize(m_nv, 0.0); + + m_directsolve_config.rtol.resize(m_nv, 1.0e-4); + m_directsolve_config.atol.resize(m_nv, 1.0e-9); + m_directsolve_config.convtol = 1.0e-14; + m_directsolve_config.dt = 0; + m_directsolve_config.jac_maxage = 5; + m_directsolve_config.jac_rtol = 1.0e-15; + m_directsolve_config.jac_atol = sqrt(std::numeric_limits::epsilon()); + + m_timestep_config.rtol.resize(m_nv, 1.0e-4); + m_timestep_config.atol.resize(m_nv, 1.0e-11); + m_timestep_config.convtol = 1.0e-14; + m_timestep_config.dt = 1.0e-5; + m_timestep_config.jac_maxage = 5; + m_timestep_config.jac_rtol = 1.0e-15; + m_timestep_config.jac_atol = sqrt(std::numeric_limits::epsilon()); + + m_config = &m_directsolve_config; + + m_xlast.resize(m_nv); + m_xsave.resize(m_nv); + + m_jacobian = DenseMatrix(m_nv, m_nv); + m_jacAge = npos; +} + +void Newton::evalJacobian(doublereal* x, doublereal* xdot) { + + // calculate unperturbed residual + m_residfunc->eval(0, x, xdot, 0); + + for (size_t n = 0; n < m_nv; n++) { + // calculate the nth Jacobian column + double xsave = x[n]; + + // calculate the perturbation amount, preserving the sign of x[n] + double dx; + if (xsave >= 0) { + dx = xsave*m_config->jac_rtol + m_config->jac_atol; + } else { + dx = xsave*m_config->jac_rtol - m_config->jac_atol; + } + + // perturb the solution vector + x[n] = xsave + dx; + dx = x[n] - xsave; + + // calculate perturbed residual + vector_fp xdotPerturbed(m_nv); //make this member for speed? + m_residfunc->eval(0, x, xdotPerturbed.data(), 0); + + // compute nth column of Jacobian + for (size_t m = 0; m < m_nv; m++) { + m_jacobian.value(m,n) = (xdotPerturbed[m] - xdot[m])/dx; + } + // restore solution vector + x[n] = xsave; + } + + // timestep components + if (m_rdt > 0) { + for (size_t i = 0; i < m_nv; i++) { + m_jacobian.value(i,i) -= m_rdt; + } + } + + // for constant-valued components: 1 in diagonal position, all 0's in row and column + for (size_t i : m_constantComponents) { + for (size_t j = 0; j < m_nv; j++) { + m_jacobian.value(i,j) = 0; + m_jacobian.value(j,i) = 0; + } + m_jacobian.value(i,i) = 1; + } + + // factor and save jacobian, will be reused for faster step computation + m_jacFactored = m_jacobian; + Cantera::factor(m_jacFactored); +} + +// RMSD (weighted) +doublereal Newton::weightedNorm(const doublereal* x, const doublereal* step) const +{ + double square = 0.0; + for (size_t i = 0; i < m_nv; i++) { + square += pow(step[i]/(x[i] + m_config->atol[i]), 2); + } + return sqrt(square/m_nv); +} + +void Newton::step(doublereal* x, doublereal* step) +{ + m_residfunc->eval(0, x, step, 0); + for (size_t n = 0; n < m_nv; n++) { + step[n] = -step[n] + m_rdt*(x[n]-m_xlast[n]); + } + + DenseMatrix solvejac = m_jacFactored; + try { + Cantera::solveFactored(solvejac, step); + } catch (CanteraError&) { + // int iok = m_jac->info() - 1; + int iok = -1; //TODO: enable error info + if (iok >= 0) { + throw CanteraError("Newton::step", + "Jacobian is singular for component {} (Matrix row {})", + iok, iok); //TODO: add component name + } else { + throw; + } + } +} + +int Newton::hybridSolve() { + int MAX = 17; + int newtonsolves = 0; + int timesteps = 0; + + // initial state + m_residfunc->getState(m_x.data()); + copy(m_x.begin(), m_x.end(), m_xsave.begin()); + + for(int i = 0; i < MAX; i++) { + newtonsolves++; + m_config = &m_directsolve_config; + if(solve(m_x.data()) > 0) { + writelog("\nConverged in {} newton solves, {} timesteps.", newtonsolves, timesteps); + return 1; + } + m_config = &m_timestep_config; + copy(m_xsave.begin(), m_xsave.end(), m_x.begin()); + for(int j = 0; j < MAX; j++) { + if (solve(m_x.data()) < 0) { + writelog("\nTimestep failure after {} newton solves, {} timesteps.", newtonsolves, timesteps); + return -1; + } + timesteps++; + } + copy(m_x.begin(), m_x.end(), m_xsave.begin()); + } + writelog("Failure to converge in {} steps.", MAX); + return 0; +} + +//! input initial state using parameter `x` +// for time integration, input parameter `dt` != 0 +// call without `dt`for direct nonlinear solution +// solution state available in *x* on return +int Newton::solve(double* x) +{ + copy(&x[0], &x[m_nv], m_x.begin()); + m_jacAge = npos; + + if (m_config->dt) { + copy(m_x.begin(), m_x.end(), m_xlast.begin()); + m_rdt = 1/m_config->dt; + } else { + m_rdt = 0; + } + + while (true) { + // Check whether the Jacobian should be re-evaluated. + if (m_jacAge > m_config->jac_maxage) { + evalJacobian(&m_x[0], &m_stp[0]); + m_jacAge = 0; + } + + //compute the undamped Newton step and save its weighted norm + step(&m_x[0], &m_stp[0]); + double step_rms = weightedNorm(&m_x[0], &m_stp[0]); + m_jacAge++; + + // compute the multiplier to keep all components in bounds. + double bound_factor = 1.0; + for (size_t i = 0; i < m_nv; i++) { + double upper_bound = m_upper_bounds[i]; + double lower_bound = m_lower_bounds[i]; + double val = m_x[i]; + if (val > upper_bound + 1.0e-12 || val < lower_bound - 1.0e-12) { + throw CanteraError("Newton::dampStep", "solution out of bounds"); + } + double newval = val + m_stp[i]; + + if (newval > upper_bound) { + bound_factor = max(0.0, min(bound_factor, (upper_bound - val)/(newval - val))); + } else if (newval < lower_bound) { + bound_factor = min(bound_factor, (val - lower_bound)/(val - newval)); + } + } + // if bound factor is very small, then x0 is already close to the boundary and + // step0 points out of the allowed domain. In this case, the Newton + // algorithm fails, so return an error condition. + if (bound_factor < 1.e-10) { + return -1; + throw CanteraError("Newton::dampStep", "solution at limits"); + } + + int m = -1; + // damped step - attempt to find a damping coefficient such that the next + // undamped step would have a RMSD smaller than that of step0 + for (double damp = bound_factor; damp > damp_min; damp /= damp_factor) { + // step the solution by the damped step size + for (size_t j = 0; j < m_nv; j++) { + m_x1[j] = damp*m_stp[j] + m_x[j]; + } + // compute the next undamped step that would result if x1 is accepted, and save its weighted norm + step(&m_x1[0], &m_stp1[0]); + double nextstep_rms = weightedNorm(&m_x1[0], &m_stp1[0]); + + // converged solution criteria + if (nextstep_rms < m_config->convtol) { + copy(m_x1.begin(), m_x1.end(), &x[0]); + return 1; + } + // Also accept it if this step would result in a + // converged solution - successful step, but not converged yet. + // Take the damped step and try again. + if (nextstep_rms < step_rms) { + m = 0; + copy(m_x1.begin(), m_x1.end(), m_x.begin()); + break; + } + } + + if (m < 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 (m_jacAge > 1) { + m_jacAge = npos; + } else { + return -1; + } + } + } +} + +} // end namespace Cantera diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 16c183c26b..6d72852ecb 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -9,6 +9,7 @@ #include "cantera/base/utilities.h" #include "cantera/base/Array.h" #include "cantera/numerics/Integrator.h" +#include "cantera/numerics/Newton.h" using namespace std; @@ -430,4 +431,47 @@ size_t ReactorNet::registerSensitivityParameter( return m_sens_params.size() - 1; } +// ----------- 0DSS ----------- +double ReactorNet::solveSteady() +{ + if (m_reactors.empty()) { + throw CanteraError("ReactorNet::solveSteady", + "no reactors in network!"); + } + + //initialize + m_nv = 0; + m_start.assign(1, 0); + for (size_t n = 0; n < m_reactors.size(); n++) { + Reactor& r = *m_reactors[n]; + if (r.nWalls()) { + throw CanteraError("ReactorNet::solveSteady", + "Wall detected in network - cannot solve for" + "steady state with transient components."); + } + //TODO: add check for constant mass flow rate + r.initialize(m_time); + m_nv += r.neq(); + m_start.push_back(m_nv); + } + + //TODO: better initialization/configuration for this solver + //solve + std::unique_ptr m_newt; + m_newt.reset(new Cantera::Newton(*this)); + + for (int i = 0; i < m_nv; i++) + { + m_newt->setBounds(i, -1.0e-3, 1.01); + } + m_newt->setBounds(0, 1.0e-3, 100); + m_newt->setBounds(1, 0, 5); + m_newt->setBounds(2, -10000000, 10000000); + + m_newt->setConstants({1}); + // m_newt->setConstants({0,1,2,11,12}); + + return m_newt->hybridSolve(); +} + }