From 34b8b965e7d390a9df33ba61315f3120279977fa Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Fri, 13 Dec 2019 17:22:48 +0100
Subject: [PATCH] [test] Add test for scalar root finding algorithms
---
dumux/nonlinear/findscalarroot.hh | 14 +++---
test/CMakeLists.txt | 1 +
test/nonlinear/CMakeLists.txt | 1 +
test/nonlinear/findscalarroot/CMakeLists.txt | 2 +
.../findscalarroot/test_findscalarroot.cc | 50 +++++++++++++++++++
5 files changed, 60 insertions(+), 8 deletions(-)
create mode 100644 test/nonlinear/CMakeLists.txt
create mode 100644 test/nonlinear/findscalarroot/CMakeLists.txt
create mode 100644 test/nonlinear/findscalarroot/test_findscalarroot.cc
diff --git a/dumux/nonlinear/findscalarroot.hh b/dumux/nonlinear/findscalarroot.hh
index 2a6d4b1990..7c230d47fe 100644
--- a/dumux/nonlinear/findscalarroot.hh
+++ b/dumux/nonlinear/findscalarroot.hh
@@ -36,14 +36,13 @@
namespace Dumux {
/*!
- * \file
* \ingroup Common
* \brief Newton's root finding algorithm for scalar functions (secant method)
* \param xOld initial guess
* \param residual Residual function
* \param derivative Derivative of the residual
- * \param maxIter Maximum number of iterations
* \param tol Relative shift tolerance
+ * \param maxIter Maximum number of iterations
* \todo replace std::result_of by std::invoke_result for C++17
*/
template<class Scalar, class ResFunc, class DerivFunc,
@@ -68,7 +67,8 @@ Scalar findScalarRootNewton(Scalar xOld, const ResFunc& residual, const DerivFun
n--;
}
- if (!std::isfinite(r))
+ using std::isfinite;
+ if (!isfinite(r))
DUNE_THROW(NumericalProblem, "Residual is not finite: " << r << " after " << maxIter - n << " iterations!");
if (relativeShift > tol)
@@ -78,14 +78,13 @@ Scalar findScalarRootNewton(Scalar xOld, const ResFunc& residual, const DerivFun
}
/*!
- * \file
* \ingroup Common
* \brief Newton's root finding algorithm for scalar functions (secant method)
* \note The derivative is numerically computed. If the derivative is know use signature with derivative function.
* \param xOld initial guess
* \param residual Residual function
- * \param maxIter Maximum number of iterations
* \param tol Relative shift tolerance
+ * \param maxIter Maximum number of iterations
*/
template<class Scalar, class ResFunc,
typename std::enable_if_t<Dune::Std::is_invocable<ResFunc, Scalar>::value>...>
@@ -98,7 +97,6 @@ Scalar findScalarRootNewton(Scalar xOld, const ResFunc& residual,
}
/*!
- * \file
* \ingroup Common
* \brief Brent's root finding algorithm for scalar functions
* \note Modified from pseudo-code on wikipedia: https://en.wikipedia.org/wiki/Brent%27s_method
@@ -107,8 +105,8 @@ Scalar findScalarRootNewton(Scalar xOld, const ResFunc& residual,
* \param a Lower bound
* \param b Upper bound
* \param residual Residual function
- * \param maxIter Maximum number of iterations
* \param tol Relative shift tolerance
+ * \param maxIter Maximum number of iterations
*/
template<class Scalar, class ResFunc,
typename std::enable_if_t<Dune::Std::is_invocable<ResFunc, Scalar>::value>...>
@@ -123,7 +121,7 @@ Scalar findScalarRootBrent(Scalar a, Scalar b, const ResFunc& residual,
// check if the root is inside the given interval
using std::signbit;
if (!signbit(fa*fb))
- DUNE_THROW(NumericalProblem, "Brent's algorithm failed: [a,b] does not contain the root!");
+ DUNE_THROW(NumericalProblem, "Brent's algorithm failed: [a,b] does not contain any, or no uniquely defined root!");
// sort bounds
using std::abs; using std::swap;
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 416332236d..b5c7428a4b 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -4,5 +4,6 @@ add_subdirectory(freeflow)
add_subdirectory(io)
add_subdirectory(material)
add_subdirectory(multidomain)
+add_subdirectory(nonlinear)
add_subdirectory(porousmediumflow)
add_subdirectory(discretization)
diff --git a/test/nonlinear/CMakeLists.txt b/test/nonlinear/CMakeLists.txt
new file mode 100644
index 0000000000..3640235bdc
--- /dev/null
+++ b/test/nonlinear/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(findscalarroot)
diff --git a/test/nonlinear/findscalarroot/CMakeLists.txt b/test/nonlinear/findscalarroot/CMakeLists.txt
new file mode 100644
index 0000000000..004545b738
--- /dev/null
+++ b/test/nonlinear/findscalarroot/CMakeLists.txt
@@ -0,0 +1,2 @@
+dumux_add_test(SOURCES test_findscalarroot.cc
+ LABELS unit nonlinear)
diff --git a/test/nonlinear/findscalarroot/test_findscalarroot.cc b/test/nonlinear/findscalarroot/test_findscalarroot.cc
new file mode 100644
index 0000000000..88dcbea992
--- /dev/null
+++ b/test/nonlinear/findscalarroot/test_findscalarroot.cc
@@ -0,0 +1,50 @@
+#include <config.h>
+
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/float_cmp.hh>
+#include <dumux/nonlinear/findscalarroot.hh>
+
+int main(int argc, char* argv[]) try
+{
+ using namespace Dumux;
+
+ const auto func = [](const double x){ return x*x - 5.0; };
+ const auto deriv = [](const double x){ return 2*x; };
+ const auto absTol = 1e-15;
+ const auto relTol = 1e-15;
+
+ const auto rootNewton = findScalarRootNewton(5.0, func, deriv, absTol);
+ const auto rootNewtonNumDiff = findScalarRootNewton(5.0, func, absTol);
+ const auto rootBrent = findScalarRootBrent(0.0, 5.0, func, relTol);
+ const auto rootExact = std::sqrt(5.0);
+
+ if (Dune::FloatCmp::ne(rootExact, rootNewton, absTol))
+ DUNE_THROW(Dune::Exception, "Wrong root for Newton algorithm (exact derivative)! Relative error: " << std::abs(rootExact-rootNewton)/rootExact);
+ if (Dune::FloatCmp::ne(rootExact, rootNewtonNumDiff, absTol))
+ DUNE_THROW(Dune::Exception, "Wrong root for Newton algorithm (numeric differentiation)! Relative error: " << std::abs(rootExact-rootNewtonNumDiff)/rootExact);
+ if (Dune::FloatCmp::ne(rootExact, rootBrent, absTol))
+ DUNE_THROW(Dune::Exception, "Wrong root for Brent algorithm! Relative error: " << std::abs(rootExact-rootBrent)/rootExact);
+
+ std::cout << "Positive root of f(x) = (x*x - 5) is\n" << std::setprecision(15)
+ << "-- Exact: " << rootExact << "\n"
+ << "-- Newton (exact derivative): " << rootNewton << "\n"
+ << "-- Newton (numeric derivative): " << rootNewtonNumDiff << "\n"
+ << "-- Brent: " << rootBrent << std::endl;
+
+ return 0;
+
+}
+catch (const Dune::Exception& e)
+{
+ std::cout << e << std::endl;
+ return 1;
+}
+catch (...)
+{
+ std::cout << "Unknown exception thrown!" << std::endl;
+ return 1;
+}
--
GitLab