Complete finite difference jacobian

Author: anthony-walker

include/cantera/zeroD/ReactorNet.h

@@ -299,7 +299,45 @@ class ReactorNet : public FuncEval
     //! @param settings the settings map propagated to all reactors and kinetics objects
     virtual void setDerivativeSettings(AnyMap& settings);
 
+    //! Calculate the system Jacobian using a finite difference method.
+    //!
+    //! This method is used only for informational purposes. Jacobian calculations
+    //! for the full reactor system are handled internally by CVODES.
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    virtual Eigen::SparseMatrix<double> finiteDifferenceJacobian();
+
+    //! A wrapper method to calculate the system jacobian as Eigen::SparseMatrix<double>
+    //! @warning Depending on the particular implementation, this may return an
+    //! approximate Jacobian intended only for use in forming a preconditioner for
+    //! iterative solvers.
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    virtual Eigen::SparseMatrix<double> jacobian() {
+        vector<Eigen::Triplet<double>> jac_trips;
+        // Add before, during, after evals
+        buildJacobian(jac_trips);
+        // construct jacobian from vector
+        Eigen::SparseMatrix<double> jac(m_nv, m_nv);
+        jac.setFromTriplets(jac_trips.begin(), jac_trips.end());
+        return jac;
+    }
+
 protected:
+    //! Calculate the Jacobian of the entire reactor network.
+    //! @param jac_vector vector where jacobian triplets are added
+    //! @param offset offset added to the row and col indices of the elements
+    //! @warning Depending on the particular implementation, this may return an
+    //! approximate Jacobian intended only for use in forming a preconditioner for
+    //! iterative solvers.
+    //! @ingroup derivGroup
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    virtual void buildJacobian(vector<Eigen::Triplet<double>>& jacVector);
+
     //! Check that preconditioning is supported by all reactors in the network
     virtual void checkPreconditionerSupported() const;
 

src/zeroD/ReactorNet.cpp

@@ -545,20 +545,10 @@ void ReactorNet::preconditionerSetup(double t, double* y, double gamma)
     updateState(yCopy.data());
     // Create jacobian triplet vector
     vector<Eigen::Triplet<double>> jac_vector;
-    vector<size_t> jstarts;
-    // Get jacobians and give elements to preconditioners
-    jstarts.push_back(jac_vector.size());
-    for (size_t i = 0; i < m_reactors.size(); i++) {
-        m_reactors[i]->buildJacobian(jac_vector);
-        jstarts.push_back(jac_vector.size());
-    }
+    buildJacobian(jac_vector);
     // Add to preconditioner with offset
-    for (size_t i=0; i < m_reactors.size(); i++) {
-        for (size_t j = jstarts[i]; j < jstarts[i+1]; j++) {
-            auto it = jac_vector[j];
-            precon->setValue(it.row() + m_start[i], it.col() + m_start[i],
-            it.value());
-        }
+    for (auto it : jac_vector) {
+        precon->setValue(it.row(), it.col(), it.value());
     }
     // post reactor setup operations
     precon->setup();
@@ -587,4 +577,72 @@ void ReactorNet::checkPreconditionerSupported() const {
     }
 }
 
+void ReactorNet::buildJacobian(vector<Eigen::Triplet<double>>& jacVector)
+{
+    // network must be initialized for the jacobian
+    if (! m_init) {
+        initialize();
+    }
+    // Create jacobian triplet vector
+    vector<size_t> jstarts;
+    // Get jacobians and give elements to preconditioners
+    jstarts.push_back(jacVector.size());
+    for (size_t i = 0; i < m_reactors.size(); i++) {
+        m_reactors[i]->buildJacobian(jacVector);
+        jstarts.push_back(jacVector.size());
+    }
+    // Add to preconditioner with offset
+    for (size_t i=0; i < m_reactors.size(); i++) {
+        for (size_t j = jstarts[i]; j < jstarts[i+1]; j++) {
+            auto it = jacVector[j];
+            auto newTrip = Eigen::Triplet<double>(it.row() + m_start[i], it.col()
+                + m_start[i], it.value());
+            jacVector[j] = newTrip;
+        }
+    }
+}
+
+Eigen::SparseMatrix<double> ReactorNet::finiteDifferenceJacobian()
+{
+    // network must be initialized for the jacobian
+    if (! m_init) {
+        initialize();
+    }
+
+    // clear former jacobian elements
+    vector<Eigen::Triplet<double>> jac_trips;
+
+    // Get the current state
+    Eigen::ArrayXd yCurrent(m_nv);
+    getState(yCurrent.data());
+
+    Eigen::ArrayXd yPerturbed = yCurrent;
+    Eigen::ArrayXd ydotCurrent(m_nv), ydotPerturbed(m_nv);
+
+    eval(m_time, yCurrent.data(), ydotCurrent.data(), m_sens_params.data());
+    double rel_perturb = std::sqrt(std::numeric_limits<double>::epsilon());
+
+    for (size_t j = 0; j < m_nv; j++) {
+        yPerturbed = yCurrent;
+        double delta_y = std::max(std::abs(yCurrent[j]), 1000 * m_atols) * rel_perturb;
+        yPerturbed[j] += delta_y;
+        ydotPerturbed = 0;
+        eval(m_time, yPerturbed.data(), ydotPerturbed.data(), m_sens_params.data());
+        // d ydot_i/dy_j
+        for (size_t i = 0; i < m_nv; i++) {
+            if (ydotCurrent[i] != ydotPerturbed[i]) {
+                jac_trips.emplace_back(
+                    static_cast<int>(i), static_cast<int>(j),
+                    (ydotPerturbed[i] - ydotCurrent[i]) / delta_y);
+            }
+        }
+    }
+    updateState(yCurrent.data());
+
+    Eigen::SparseMatrix<double> jac(m_nv, m_nv);
+    jac.setFromTriplets(jac_trips.begin(), jac_trips.end());
+    return jac;
+}
+
+
 }

test/python/test_reactor.py

@@ -1200,6 +1200,17 @@ class TestIdealGasMoleReactor(TestMoleReactor):
     reactorClass = ct.IdealGasMoleReactor
     test_preconditioner_unsupported = None
 
+    @pytest.mark.diagnose
+    def test_heat_transfer_network(self):
+        # Result should be the same if (m * cp) / (U * A) is held constant
+        self.make_reactors(T1=300, T2=1000)
+        self.add_wall(U=200, A=1.0)
+        self.net.preconditioner = ct.AdaptivePreconditioner()
+        print(self.net.finite_difference_jacobian)
+        # self.net.advance(1.0)
+        # T1a = self.r1.T
+        # T2a = self.r2.T
+
     def test_adaptive_precon_integration(self):
         # Network one with non-mole reactor
         net1 = ct.ReactorNet()