Cpp Wall Jacobian Test

Author: anthony-walker

interfaces/cython/cantera/reactor.pyx

@@ -1822,4 +1822,3 @@ cdef class ReactorNet:
         def __get__(self):
             return get_from_sparse(self.net.finiteDifferenceJacobian(),
                                    self.n_vars, self.n_vars)
-

test/zeroD/test_zeroD.cpp

@@ -190,6 +190,87 @@ TEST(AdaptivePreconditionerTests, test_precon_solver_stats)
     EXPECT_GE(stats["nonlinear_conv_fails"].asInt(), 0);
 }
 
+TEST(JacobianTests, test_wall_jacobian_build)
+{
+    // create first reactor
+    auto sol1 = newSolution("h2o2.yaml");
+    sol1->thermo()->setState_TPY(1000.0, OneAtm, " O2:1.0");
+    IdealGasMoleReactor reactor1;
+    reactor1.insert(sol1);
+    reactor1.setInitialVolume(1.0);
+    // create second reactor
+    auto sol2 = newSolution("h2o2.yaml");
+    sol2->thermo()->setState_TPY(900.0, OneAtm, " O2:1.0");
+    IdealGasConstPressureMoleReactor reactor2;
+    reactor2.insert(sol2);
+    reactor2.setInitialVolume(1.0);
+    // create the wall
+    Wall w;
+    w.install(reactor1, reactor2);
+    w.setArea(2.0);
+    w.setHeatTransferCoeff(3.0);
+    // setup reactor network and integrate
+    ReactorNet network;
+    network.addReactor(reactor1);
+    network.addReactor(reactor2);
+    network.initialize();
+    // create jacobian the size of network
+    Eigen::SparseMatrix<double> wallJacMat;
+    wallJacMat.resize(network.neq(), network.neq());
+    // manually get wall jacobian elements
+    vector<Eigen::Triplet<double>>  wallJac;
+    // build jac for reactor 1 wall only
+    w.buildReactorJacobian(&reactor1, wallJac);
+    wallJacMat.setFromTriplets(wallJac.begin(), wallJac.end());
+    // check appropriate values
+    // double tol = 1e-8;
+    for (int k = 0; k < wallJacMat.outerSize(); k++) {
+        for (Eigen::SparseMatrix<double>::InnerIterator it(wallJacMat, k); it; ++it) {
+            EXPECT_DOUBLE_EQ(it.value(), 6000.0);
+            EXPECT_EQ(it.row(), 0); // check that it is the first row
+            EXPECT_GE(it.col(), reactor1.speciesOffset());
+            EXPECT_LT(it.col(), reactor1.neq());
+        }
+    }
+    // build jac for reactor 2 wall only
+    wallJac.clear();
+    w.buildReactorJacobian(&reactor2, wallJac);
+    wallJacMat.setZero();
+    wallJacMat.setFromTriplets(wallJac.begin(), wallJac.end());
+    // check appropriate values
+    // double tol = 1e-8;
+    for (int k = 0; k < wallJacMat.outerSize(); k++) {
+        for (Eigen::SparseMatrix<double>::InnerIterator it(wallJacMat, k); it; ++it) {
+            EXPECT_DOUBLE_EQ(it.value(), -5400.0);
+            EXPECT_EQ(it.row(), 0); // check that it is the first row
+            EXPECT_GE(it.col(), reactor2.speciesOffset());
+            EXPECT_LT(it.col(), reactor2.neq());
+        }
+    }
+    // build jac for network terms
+    wallJac.clear();
+    w.buildNetworkJacobian(wallJac);
+    wallJacMat.setZero();
+    wallJacMat.setFromTriplets(wallJac.begin(), wallJac.end());
+    // check appropriate values
+    // double tol = 1e-8;
+    for (int k = 0; k < wallJacMat.outerSize(); k++) {
+        for (Eigen::SparseMatrix<double>::InnerIterator it(wallJacMat, k); it; ++it) {
+            if (it.value() < 0) {
+                EXPECT_DOUBLE_EQ(it.value(), -5400.0);
+                EXPECT_EQ(it.row(), 0); // check that it is the first row
+                EXPECT_GE(it.col(), reactor1.neq() + reactor2.speciesOffset());
+                EXPECT_LT(it.col(), reactor1.neq() + reactor2.neq());
+            } else {
+                EXPECT_DOUBLE_EQ(it.value(), 6000.0);
+                EXPECT_EQ(it.row(), reactor1.neq()); // check that it is the first row
+                EXPECT_GE(it.col(), reactor2.speciesOffset());
+                EXPECT_LT(it.col(), reactor1.neq());
+            }
+        }
+    }
+}
+
 int main(int argc, char** argv)
 {
     printf("Running main() from test_zeroD.cpp\n");