[T] check results are within tolerances

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