diff --git a/opm/autodiff/AutoDiffMatrix.hpp b/opm/autodiff/AutoDiffMatrix.hpp index 04b925838..c06a535ca 100644 --- a/opm/autodiff/AutoDiffMatrix.hpp +++ b/opm/autodiff/AutoDiffMatrix.hpp @@ -491,7 +491,9 @@ namespace Opm M mat(n, n); mat.reserve(Eigen::ArrayXi::Ones(n, 1)); for (M::Index i = 0; i < n; ++i) { - mat.insert(i, i) = d[i]; + if (d[i] != 0.0) { + mat.insert(i, i) = d[i]; + } } return mat; diff --git a/tests/test_block.cpp b/tests/test_block.cpp index 51da5f0dd..9d4345ac0 100644 --- a/tests/test_block.cpp +++ b/tests/test_block.cpp @@ -71,6 +71,52 @@ namespace { // Note: Investigate implementing this operator as // return A.cwiseNotEqual(B).count() == 0; } + + + bool operator==(const AutoDiffMatrix& lhs, + const AutoDiffMatrix& rhs) + { + Eigen::SparseMatrix lhs_s, rhs_s; + + lhs.toSparse(lhs_s); + rhs.toSparse(rhs_s); + + return lhs_s == rhs_s; + } + + void checkClose(const AutoDiffBlock& lhs, const AutoDiffBlock& rhs, double tolerance) { + + BOOST_CHECK(lhs.value().isApprox(rhs.value(), tolerance)); + + auto lhs_d = lhs.derivative(); + auto rhs_d = rhs.derivative(); + + Eigen::SparseMatrix lhs_s, rhs_s; + + //If lhs has no derivatives, make sure all rhs derivatives are zero + if (lhs_d.size() == 0) { + for (size_t i=0; i& J = x.derivative(); BOOST_REQUIRE(J[0].nonZeros() == v.size()); - const Eigen::Diagonal& d = J[0].diagonal(); - BOOST_REQUIRE((d.array() == 1.0).all()); + const ADB::M& d = J[0]; + for (int i=0; i::const_iterator b = J.begin() + 1, e = J.end(); b != e; ++b) { @@ -131,8 +186,9 @@ BOOST_AUTO_TEST_CASE(FunctionInitialisation) std::vector jacs(num_blocks); for (std::vector::size_type j = 0; j < num_blocks; ++j) { - jacs[j] = ADB::M(blocksizes[FirstVar], blocksizes[j]); - jacs[j].insert(0,0) = -1.0; + Eigen::SparseMatrix sm(blocksizes[FirstVar], blocksizes[j]); + sm.insert(0,0) = -1.0; + jacs[j] = ADB::M(sm); } ADB::V v_copy(v); @@ -212,24 +268,18 @@ BOOST_AUTO_TEST_CASE(AssignAddSubtractOperators) z += y; ADB sum = x + y; const double tolerance = 1e-14; - BOOST_CHECK(z.value().isApprox(sum.value(), tolerance)); - BOOST_CHECK(z.derivative()[0].isApprox(sum.derivative()[0], tolerance)); - BOOST_CHECK(z.derivative()[1].isApprox(sum.derivative()[1], tolerance)); + checkClose(z, sum, tolerance); + z -= y; - BOOST_CHECK(z.value().isApprox(x.value(), tolerance)); - BOOST_CHECK(z.derivative()[0].isApprox(x.derivative()[0], tolerance)); - BOOST_CHECK(z.derivative()[1].isApprox(x.derivative()[1], tolerance)); + checkClose(z, x, tolerance); // Testing the case when the left hand side has empty() jacobian. ADB yconst = ADB::constant(vy); z = yconst; z -= x; ADB diff = yconst - x; - BOOST_CHECK(z.value().isApprox(diff.value(), tolerance)); - BOOST_CHECK(z.derivative()[0].isApprox(diff.derivative()[0], tolerance)); - BOOST_CHECK(z.derivative()[1].isApprox(diff.derivative()[1], tolerance)); + checkClose(z, diff, tolerance); + z += x; - BOOST_CHECK(z.value().isApprox(yconst.value(), tolerance)); - BOOST_CHECK(z.derivative()[0].isApprox(Eigen::Matrix::Zero())); - BOOST_CHECK(z.derivative()[1].isApprox(Eigen::Matrix::Zero())); + checkClose(z, yconst, tolerance); }