diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index d0aa5270864..7fa964cd045 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/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) - diff --git a/test/zeroD/test_zeroD.cpp b/test/zeroD/test_zeroD.cpp index 25ac0aaf808..bcd8a0823e0 100644 --- a/test/zeroD/test_zeroD.cpp +++ b/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 wallJacMat; + wallJacMat.resize(network.neq(), network.neq()); + // manually get wall jacobian elements + vector> 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::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::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::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");