Extend Domain1D (#12)

- Add solver matrices to Domain1D class
chengcli · May 13, 2024 · ac614f7 · ac614f7
1 parent 18e42c6
commit ac614f7
Show file tree

Hide file tree

Showing 4 changed files with 24 additions and 6 deletions.
diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h
@@ -1430,9 +1430,10 @@ class Kinetics
         }
 
         if ((il  + 1) * nWavelengths() > m_actinicFlux.lock()->size()) {
+            auto max_levels = m_actinicFlux.lock()->size() / nWavelengths();
             throw CanteraError("Kinetics::setActinicFluxLevel",
-                               "Requested start index {} is out of bounds for actinic flux vector of length {}",
-                               il, nWavelengths());
+                               "Requested end level {} is out of bounds for actinic flux vector of length {}",
+                               (il + 1), max_levels);
         }
 
         if (il != m_actinicFluxLevel) {

diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h
@@ -7,6 +7,8 @@
 #define CT_DOMAIN1D_H
 
 #include "cantera/base/ctexceptions.h"
+#include "cantera/numerics/eigen_dense.h"
+#include "cantera/numerics/eigen_sparse.h"
 
 namespace Cantera
 {
@@ -365,7 +367,7 @@ class Domain1D
         return m_solution;
     }
 
-    size_t size() const {
+    virtual size_t size() const {
         return m_nv*m_points;
     }
 
@@ -521,6 +523,17 @@ class Domain1D
         m_state = data;
     }
 
+    virtual void update(double const *y) {}
+
+    //----------------------------------------
+    //        Tri-diagonal matrix solver
+    //----------------------------------------
+    Domain1D* forwardSweep(double const *residual);
+    Domain1D* backwardSweep(double *result);
+    virtual void modifyLeft(Eigen::VectorXd &delta) {}
+    virtual void modifyRight(Eigen::SparseMatrix<double> &a,
+                             Eigen::VectorXd &delta) {}
+
 protected:
     //! Retrieve meta data
     virtual AnyMap getMeta() const;
@@ -563,6 +576,10 @@ class Domain1D
 
     //! Composite thermo/kinetics/transport handler
     shared_ptr<Solution> m_solution;
+
+    //! tri-diagonal matrices 
+    vector<Eigen::SparseMatrix<double>> m_A, m_B, m_C;
+    vector<Eigen::VectorXd> m_D;
 };
 }
 

diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h
@@ -51,7 +51,7 @@ class OneDim
      * @param x1         Final solution satisfying F(x1) = 0.
      * @param loglevel   Controls amount of diagnostic output.
      */
-    int solve(double* x0, double* x1, int loglevel);
+    virtual int solve(double* x0, double* x1, int loglevel);
 
     //! Number of domains.
     size_t nDomains() const {
@@ -63,7 +63,7 @@ class OneDim
         return *m_dom[i];
     }
 
-    size_t domainIndex(const string& name);
+    size_t domainIndex(const string& name) const;
 
     //! Check that the specified domain index is in range.
     //! Throws an exception if n is greater than nDomains()-1

diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp
@@ -37,7 +37,7 @@ OneDim::~OneDim()
 {
 }
 
-size_t OneDim::domainIndex(const string& name)
+size_t OneDim::domainIndex(const string& name) const
 {
     for (size_t n = 0; n < m_dom.size(); n++) {
         if (domain(n).id() == name) {