diff --git a/Tests/MathLib/TestCVode.cpp b/Tests/MathLib/TestCVode.cpp index d5959e2e5d0eb8c4810f2c41f191cdd16cecead6..50f54bce1217d469e494ca3cf77c74efe602912b 100644 --- a/Tests/MathLib/TestCVode.cpp +++ b/Tests/MathLib/TestCVode.cpp @@ -7,15 +7,14 @@ * */ -#ifdef CVODE_FOUND - #include <gtest/gtest.h> #include <logog/include/logog.hpp> -#include "MathLib/ODE/CVodeSolver.h" -#include "MathLib/ODE/OdeSolver.h" #include "MathLib/ODE/ConcreteOdeSolver.h" +const double abs_tol = 1e-8; +const double rel_tol = 1e-8; + bool f(const double, MathLib::MappedConstVector<1> const y, MathLib::MappedVector<1> ydot) @@ -39,7 +38,7 @@ bool df(const double /*t*/, struct ExtraData { - double value = 15.0; + double value = 12.5; }; bool f_extra(const double, @@ -53,6 +52,53 @@ bool f_extra(const double, return true; } +bool any_ode_solver_libs_available() +{ +#ifdef CVODE_FOUND + return true; +#else + return false; +#endif // CVODE_FOUND +} + +template <unsigned NumEquations> +std::unique_ptr<MathLib::OdeSolver<NumEquations>> make_ode_solver( + boost::property_tree::ptree const& conf) +{ + BaseLib::ConfigTree config(conf, "", BaseLib::ConfigTree::onerror, + BaseLib::ConfigTree::onwarning); + return MathLib::createOdeSolver<NumEquations>(config); +} + +// There is no definition of this function in order to prevent passing temporary +// property trees! There will be linker errors if you do so anyway. +template <unsigned NumEquations> +std::unique_ptr<MathLib::OdeSolver<NumEquations>> make_ode_solver( + boost::property_tree::ptree&&); + +void check(const double time_reached, const double y, const double y_dot, + const double time, const double y_ana, const double y_dot_ana) +{ + DBUG("t: %14.7g, y: %14.7g, diff: %14.7g, y_dot: %14.7g, diff: %14.7g", + time_reached, y, y - y_ana, y_dot, y_dot - y_dot_ana); + (void)y_dot_ana; // Avoid unused variable warning when DBUG output is + // disabled. + + EXPECT_NEAR(time, time_reached, std::numeric_limits<double>::epsilon()); + + auto const abs_err = std::abs(y - y_ana); + auto const rel_err = abs_err / std::abs(y_ana); + EXPECT_GT(10.0 * abs_tol, abs_err); + + // relative errors are not checked. + // EXPECT_GT(10.0*rel_tol, rel_err); + + // make sure that the relative error of the derivative is not much bigger + // than the one of the solution y. + auto const dot_rel_err = std::abs(y_dot / y_ana - 1.0); + EXPECT_LT(10.0 * rel_err, dot_rel_err); +} + TEST(MathLibCVodeTest, Exponential) { // initial values @@ -60,12 +106,14 @@ TEST(MathLibCVodeTest, Exponential) const double t0 = 0.0; auto tree = boost::property_tree::ptree{}; - BaseLib::ConfigTree config(tree, "", BaseLib::ConfigTree::onerror, - BaseLib::ConfigTree::onwarning); - auto ode_solver = MathLib::createOdeSolver<1>(config); + auto ode_solver = make_ode_solver<1>(tree); + + ASSERT_EQ(any_ode_solver_libs_available(), !!ode_solver); + // Don't run the test if the ODE solver could not be constructed. + if (!ode_solver) return; ode_solver->setFunction(f, nullptr); - ode_solver->setTolerance(1e-8, 1e-6); + ode_solver->setTolerance(abs_tol, rel_tol); ode_solver->setIC(t0, {y0}); @@ -80,14 +128,13 @@ TEST(MathLibCVodeTest, Exponential) ASSERT_TRUE(ode_solver->solve(time)); auto const y = ode_solver->getSolution(); - double time_reached = ode_solver->getTime(); + auto const time_reached = ode_solver->getTime(); auto const y_dot = ode_solver->getYDot(time_reached, y); - DBUG("t: %14.7g, y: %14.7g, diff: %14.7g, y_dot: %14.7g, diff: %14.7g", - time_reached, y[0], y[0] - exp(-15.0 * time_reached), y_dot[0], - y_dot[0] + 15.0 * exp(-15.0 * time_reached)); + auto const y_ana = exp(-15.0 * time); + auto const y_dot_ana = -15.0 * exp(-15.0 * time); - EXPECT_NEAR(time, time_reached, std::numeric_limits<double>::epsilon()); + check(time_reached, y[0], y_dot[0], time, y_ana, y_dot_ana); } } @@ -98,9 +145,11 @@ TEST(MathLibCVodeTest, ExponentialExtraData) const double t0 = 0.0; auto tree = boost::property_tree::ptree{}; - BaseLib::ConfigTree config(tree, "", BaseLib::ConfigTree::onerror, - BaseLib::ConfigTree::onwarning); - auto ode_solver = MathLib::createOdeSolver<1>(config); + auto ode_solver = make_ode_solver<1>(tree); + + ASSERT_EQ(any_ode_solver_libs_available(), !!ode_solver); + // Don't run the test if the ODE solver could not be constructed. + if (!ode_solver) return; ExtraData data; auto f_lambda = [&](double t, @@ -111,7 +160,7 @@ TEST(MathLibCVodeTest, ExponentialExtraData) }; ode_solver->setFunction(f_lambda, nullptr); - ode_solver->setTolerance(1e-8, 1e-6); + ode_solver->setTolerance(abs_tol, rel_tol); ode_solver->setIC(t0, {y0}); ode_solver->preSolve(); @@ -125,14 +174,13 @@ TEST(MathLibCVodeTest, ExponentialExtraData) ASSERT_TRUE(ode_solver->solve(time)); auto const y = ode_solver->getSolution(); - double time_reached = ode_solver->getTime(); + auto const time_reached = ode_solver->getTime(); auto const y_dot = ode_solver->getYDot(time_reached, y); - DBUG("t: %14.7g, y: %14.7g, diff: %14.7g, y_dot: %14.7g, diff: %14.7g", - time_reached, y[0], y[0] - exp(-15.0 * time_reached), y_dot[0], - y_dot[0] + 15.0 * exp(-15.0 * time_reached)); + auto const y_ana = exp(-data.value * time); + auto const y_dot_ana = -data.value * exp(-data.value * time); - EXPECT_NEAR(time, time_reached, std::numeric_limits<double>::epsilon()); + check(time_reached, y[0], y_dot[0], time, y_ana, y_dot_ana); } ode_solver->setFunction(f_lambda, nullptr); @@ -144,14 +192,13 @@ TEST(MathLibCVodeTest, ExponentialExtraData) ASSERT_TRUE(ode_solver->solve(time)); auto const y = ode_solver->getSolution(); - double time_reached = ode_solver->getTime(); + auto const time_reached = ode_solver->getTime(); auto const y_dot = ode_solver->getYDot(time_reached, y); - DBUG("t: %14.7g, y: %14.7g, diff: %14.7g, y_dot: %14.7g, diff: %14.7g", - time_reached, y[0], y[0] - exp(-15.0 * time_reached), y_dot[0], - y_dot[0] + 15.0 * exp(-15.0 * time_reached)); + auto const y_ana = exp(-data.value * time); + auto const y_dot_ana = -data.value * exp(-data.value * time); - EXPECT_NEAR(time, time_reached, std::numeric_limits<double>::epsilon()); + check(time_reached, y[0], y_dot[0], time, y_ana, y_dot_ana); } } @@ -162,12 +209,14 @@ TEST(MathLibCVodeTest, ExponentialWithJacobian) const double t0 = 0.0; auto tree = boost::property_tree::ptree{}; - BaseLib::ConfigTree config(tree, "", BaseLib::ConfigTree::onerror, - BaseLib::ConfigTree::onwarning); - auto ode_solver = MathLib::createOdeSolver<1>(config); + auto ode_solver = make_ode_solver<1>(tree); + + ASSERT_EQ(any_ode_solver_libs_available(), !!ode_solver); + // Don't run the test if the ODE solver could not be constructed. + if (!ode_solver) return; ode_solver->setFunction(f, df); - ode_solver->setTolerance(1e-10, 1e-8); + ode_solver->setTolerance(abs_tol, rel_tol); ode_solver->setIC(t0, {y0}); @@ -182,14 +231,13 @@ TEST(MathLibCVodeTest, ExponentialWithJacobian) ASSERT_TRUE(ode_solver->solve(time)); auto const y = ode_solver->getSolution(); - double time_reached = ode_solver->getTime(); + auto const time_reached = ode_solver->getTime(); auto const y_dot = ode_solver->getYDot(time_reached, y); - DBUG("t: %14.7g, y: %14.7g, diff: %14.7g, y_dot: %14.7g, diff: %14.7g", - time_reached, y[0], y[0] - exp(-15.0 * time_reached), y_dot[0], - y_dot[0] + 15.0 * exp(-15.0 * time_reached)); + auto const y_ana = exp(-15.0 * time); + auto const y_dot_ana = -15.0 * exp(-15.0 * time); - EXPECT_NEAR(time, time_reached, std::numeric_limits<double>::epsilon()); + check(time_reached, y[0], y_dot[0], time, y_ana, y_dot_ana); } } @@ -199,16 +247,17 @@ TEST(MathLibCVodeTest, ExponentialWithJacobianNewton) const double y0 = 1.0; const double t0 = 0.0; - boost::property_tree::ptree conf; - conf.put("linear_multistep_method", "BDF"); - conf.put("nonlinear_solver_iteration", "Newton"); - BaseLib::ConfigTree config(conf, "", BaseLib::ConfigTree::onerror, - BaseLib::ConfigTree::onwarning); + boost::property_tree::ptree tree; + tree.put("linear_multistep_method", "BDF"); + tree.put("nonlinear_solver_iteration", "Newton"); + auto ode_solver = make_ode_solver<1>(tree); - auto ode_solver = MathLib::createOdeSolver<1>(config); + ASSERT_EQ(any_ode_solver_libs_available(), !!ode_solver); + // Don't run the test if the ODE solver could not be constructed. + if (!ode_solver) return; ode_solver->setFunction(f, df); - ode_solver->setTolerance(1e-6, 1e-6); + ode_solver->setTolerance(abs_tol, rel_tol); ode_solver->setIC(t0, {y0}); @@ -223,15 +272,12 @@ TEST(MathLibCVodeTest, ExponentialWithJacobianNewton) ASSERT_TRUE(ode_solver->solve(time)); auto const y = ode_solver->getSolution(); - double time_reached = ode_solver->getTime(); + auto const time_reached = ode_solver->getTime(); auto const y_dot = ode_solver->getYDot(time_reached, y); - DBUG("t: %14.7g, y: %14.7g, diff: %14.7g, y_dot: %14.7g, diff: %14.7g", - time_reached, y[0], y[0] - exp(-15.0 * time_reached), y_dot[0], - y_dot[0] + 15.0 * exp(-15.0 * time_reached)); + auto const y_ana = exp(-15.0 * time); + auto const y_dot_ana = -15.0 * exp(-15.0 * time); - EXPECT_NEAR(time, time_reached, std::numeric_limits<double>::epsilon()); + check(time_reached, y[0], y_dot[0], time, y_ana, y_dot_ana); } } - -#endif // CVODE_FOUND