Steady-State Solver for ZeroD Reactors

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.
 /*!

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

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.

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":

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')
 

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<long int>(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<long int>(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()) {