Make Reactor Jacobian available in Python

Cantera · Aug 30, 2022 · 3fc279e · 3fc279e
1 parent 9101e60
commit 3fc279e
Show file tree

Hide file tree

Showing 4 changed files with 19 additions and 4 deletions.
diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd
@@ -42,6 +42,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
         string componentName(size_t) except +translate_exception
         size_t neq()
         void getState(double*) except +translate_exception
+        CxxSparseMatrix jacobian() except +translate_exception
         void addSurface(CxxReactorSurface*)
         void setAdvanceLimit(string&, double) except +translate_exception
         void addSensitivityReaction(size_t) except +translate_exception

diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx
@@ -333,6 +333,10 @@ cdef class Reactor(ReactorBase):
         self.reactor.getState(&y[0])
         return y
 
+    property jacobian:
+        def __get__(self):
+            return get_from_sparse(self.reactor.jacobian(), self.n_vars, self.n_vars)
+
     def set_advance_limit(self, name, limit):
         """
         Limit absolute change of component ``name`` during `ReactorNet.advance`.

diff --git a/src/zeroD/IdealGasConstPressureMoleReactor.cpp b/src/zeroD/IdealGasConstPressureMoleReactor.cpp
@@ -134,6 +134,10 @@ void IdealGasConstPressureMoleReactor::eval(double time, double* LHS, double* RH
 
 Eigen::SparseMatrix<double> IdealGasConstPressureMoleReactor::jacobian()
 {
+    if (m_nv == 0) {
+        throw CanteraError("IdealGasConstPressureMoleReactor::jacobian",
+                           "Reactor must be initialized first.");
+    }
     // clear former jacobian elements
     m_jac_trips.clear();
     // Determine Species Derivatives
@@ -174,10 +178,11 @@ Eigen::SparseMatrix<double> IdealGasConstPressureMoleReactor::jacobian()
         yNext[0] += deltaTemp;
         // getting perturbed state
         updateState(yNext.data());
-        eval(m_net->time(), yNext.data(), ydotNext.data());
+        double time = (m_net != nullptr) ? m_net->time() : 0.0;
+        eval(time, yNext.data(), ydotNext.data());
         // reset and get original state
         updateState(yCurrent.data());
-        eval(m_net->time(), yCurrent.data(), ydotCurrent.data());
+        eval(time, yCurrent.data(), ydotCurrent.data());
         // d T_dot/dT
         m_jac_trips.emplace_back(0, 0, (ydotNext[0] - ydotCurrent[0]) / deltaTemp);
         // d omega_dot_j/dT

diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp
@@ -158,6 +158,10 @@ void IdealGasMoleReactor::eval(double time, double* LHS, double* RHS)
 
 Eigen::SparseMatrix<double> IdealGasMoleReactor::jacobian()
 {
+    if (m_nv == 0) {
+        throw CanteraError("IdealGasMoleReactor::jacobian",
+                           "Reactor must be initialized first.");
+    }
     // clear former jacobian elements
     m_jac_trips.clear();
     // Determine Species Derivatives
@@ -186,10 +190,11 @@ Eigen::SparseMatrix<double> IdealGasMoleReactor::jacobian()
         yNext[0] += deltaTemp;
         // getting perturbed state
         updateState(yNext.data());
-        eval(m_net->time(), yNext.data(), ydotNext.data());
+        double time = (m_net != nullptr) ? m_net->time() : 0.0;
+        eval(time, yNext.data(), ydotNext.data());
         // reset and get original state
         updateState(yCurrent.data());
-        eval(m_net->time(), yCurrent.data(), ydotCurrent.data());
+        eval(time, yCurrent.data(), ydotCurrent.data());
         // d T_dot/dT
         m_jac_trips.emplace_back(0, 0, (ydotNext[0] - ydotCurrent[0]) / deltaTemp);
         // d omega_dot_j/dT