GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/nonlinear/findscalarroot.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 59 62 95.2%
Functions: 10 10 100.0%
Branches: 59 162 36.4%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup Nonlinear
10 * \brief Root finding algorithms for scalar functions
11 */
12 #ifndef DUMUX_COMMON_SCALAR_ROOT_FINDING_HH
13 #define DUMUX_COMMON_SCALAR_ROOT_FINDING_HH
14
15 #include <cmath>
16 #include <limits>
17 #include <type_traits>
18
19 #include <dumux/common/exceptions.hh>
20 #include <dumux/common/parameters.hh>
21 #include <dumux/common/numericdifferentiation.hh>
22
23 namespace Dumux {
24
25 /*!
26 * \ingroup Nonlinear
27 * \brief Newton's root finding algorithm for scalar functions (secant method)
28 * \param xOld initial guess
29 * \param residual Residual function
30 * \param derivative Derivative of the residual
31 * \param tol Relative shift tolerance
32 * \param maxIter Maximum number of iterations
33 */
34 template<class Scalar, class ResFunc, class DerivFunc,
35 typename std::enable_if_t<std::is_invocable_r_v<Scalar, ResFunc, Scalar> &&
36 std::is_invocable_r_v<Scalar, DerivFunc, Scalar>>...>
37 506 Scalar findScalarRootNewton(Scalar xOld, const ResFunc& residual, const DerivFunc& derivative,
38 const Scalar tol = 1e-13, const int maxIter = 200)
39 {
40 506 Scalar xNew = xOld;
41 506 Scalar r = residual(xNew);
42
43 506 int n = maxIter;
44 506 Scalar relativeShift = std::numeric_limits<Scalar>::max();
45
3/4
✓ Branch 0 taken 1521 times.
✓ Branch 1 taken 504 times.
✓ Branch 2 taken 1521 times.
✗ Branch 3 not taken.
 2041 while (relativeShift > tol && n > 0)
46 {
47
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
 1535 xNew = xOld - r/derivative(xOld);
48 1535 r = residual(xNew);
49
50 using std::abs; using std::max;
51
2/2
✓ Branch 0 taken 732 times.
✓ Branch 1 taken 789 times.
 1535 relativeShift = abs(xOld-xNew)/max(abs(xOld), abs(xNew));
52 1535 xOld = xNew;
53 1535 n--;
54 }
55
56 using std::isfinite;
57
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 504 times.
 506 if (!isfinite(r))
58 ✗ DUNE_THROW(NumericalProblem, "Residual is not finite: " << r << " after " << maxIter - n << " iterations!");
59
60
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 504 times.
 506 if (relativeShift > tol)
61 ✗ DUNE_THROW(NumericalProblem, "Scalar newton solver did not converge after " << maxIter << " iterations!");
62
63 506 return xNew;
64 }
65
66 /*!
67 * \ingroup Nonlinear
68 * \brief Newton's root finding algorithm for scalar functions (secant method)
69 * \note The derivative is numerically computed. If the derivative is know use signature with derivative function.
70 * \param xOld initial guess
71 * \param residual Residual function
72 * \param tol Relative shift tolerance
73 * \param maxIter Maximum number of iterations
74 */
75 template<class Scalar, class ResFunc,
76 typename std::enable_if_t<std::is_invocable_r_v<Scalar, ResFunc, Scalar>>...>
77
1/2
✓ Branch 1 taken 502 times.
✗ Branch 2 not taken.
 503 Scalar findScalarRootNewton(Scalar xOld, const ResFunc& residual,
78 const Scalar tol = 1e-13, const int maxIter = 200)
79 {
80 503 const auto eps = NumericDifferentiation::epsilon(xOld);
81
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
 2016 auto derivative = [&](const auto x){ return (residual(x + eps)-residual(x))/eps; };
82
1/2
✓ Branch 1 taken 502 times.
✗ Branch 2 not taken.
 503 return findScalarRootNewton(xOld, residual, derivative, tol, maxIter);
83 }
84
85 /*!
86 * \ingroup Nonlinear
87 * \brief Brent's root finding algorithm for scalar functions
88 * \note Modified from pseudo-code on wikipedia: https://en.wikipedia.org/wiki/Brent%27s_method
89 * \note See also R.P. Brent "An algorithm with guaranteed convergence for finding a zero of a function", The Computer Journal (1971).
90 * \note This is usually more robust than Newton's method
91 * \param a Lower bound
92 * \param b Upper bound
93 * \param residual Residual function
94 * \param tol Relative shift tolerance
95 * \param maxIter Maximum number of iterations
96 */
97 template<class Scalar, class ResFunc,
98 typename std::enable_if_t<std::is_invocable_r_v<Scalar, ResFunc, Scalar>>...>
99 10829087 Scalar findScalarRootBrent(Scalar a, Scalar b, const ResFunc& residual,
100 const Scalar tol = 1e-13, const int maxIter = 200)
101 {
102 // precompute the residuals (minimize function evaluations)
103
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 10829087 Scalar fa = residual(a);
104 10829087 Scalar fb = residual(b);
105 10829087 Scalar fs = 0.0;
106
107 // check if the root is inside the given interval
108 using std::signbit;
109
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 10829086 times.
 10829087 if (!signbit(fa*fb))
110
11/22
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 1 times.
✗ Branch 29 not taken.
✓ Branch 32 taken 1 times.
✗ Branch 33 not taken.
 5 DUNE_THROW(NumericalProblem, "Brent's algorithm failed: [a,b] does not contain any, or no uniquely defined root!");
111
112 // sort bounds
113 using std::abs; using std::swap;
114
2/2
✓ Branch 0 taken 9799697 times.
✓ Branch 1 taken 1029389 times.
 10829086 if (abs(fa) < abs(fb))
115 {
116 9799697 swap(a, b);
117 9799697 swap(fa, fb);
118 }
119
120 10829086 Scalar c = a;
121 10829086 Scalar fc = fa;
122 10829086 Scalar d = 0.0;
123 10829086 Scalar s = 0.0;
124 10829086 bool flag = true;
125
126
1/2
✓ Branch 0 taken 261812146 times.
✗ Branch 1 not taken.
 261812146 for (int i = 0; i < maxIter; ++i)
127 {
128 // stopping criterion
129 using std::max;
130
4/4
✓ Branch 0 taken 107563990 times.
✓ Branch 1 taken 154248156 times.
✓ Branch 2 taken 10829086 times.
✓ Branch 3 taken 250983060 times.
 369376136 if (abs(b-a) < tol*max(abs(a), abs(b)))
131 10829086 return b;
132
133 // inverse quadratic interpolation
134
4/4
✓ Branch 0 taken 207035948 times.
✓ Branch 1 taken 43947112 times.
✓ Branch 2 taken 34661581 times.
✓ Branch 3 taken 172374367 times.
 250983060 if (fa != fc && fb != fc)
135 {
136 34661581 const auto fab = fa-fb;
137 34661581 const auto fac = fa-fc;
138 34661581 const auto fbc = fb-fc;
139 34661581 s = a*fb*fc/(fab*fac) - b*fa*fc/(fab*fbc) + c*fa*fb/(fac*fbc);
140 34661581 }
141
142 // secant method
143 else
144 {
145 216321479 s = b - fb*(b-a)/(fb-fa);
146 }
147
148 // bisection method
149
2/2
✓ Branch 0 taken 35016364 times.
✓ Branch 1 taken 12148387 times.
 47164751 if ( (s < (3*a + b)*0.25 || s > b)
150
4/4
✓ Branch 0 taken 11426227 times.
✓ Branch 1 taken 23590137 times.
✓ Branch 2 taken 18686252 times.
✓ Branch 3 taken 4903885 times.
 35016364 || (flag && abs(s-b) >= abs(b-c)*0.5)
151
2/2
✓ Branch 0 taken 11238599 times.
✓ Branch 1 taken 187628 times.
 11426227 || (!flag && abs(s-b) >= abs(c-d)*0.5)
152
4/4
✓ Branch 0 taken 5854368 times.
✓ Branch 1 taken 12831884 times.
✓ Branch 2 taken 18519154 times.
✓ Branch 3 taken 167098 times.
 24540620 || (flag && abs(b-c) < tol*max(abs(b), abs(c)))
153
7/8
✓ Branch 0 taken 47164751 times.
✓ Branch 1 taken 203818309 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 18519154 times.
✓ Branch 4 taken 2090330 times.
✓ Branch 5 taken 9148269 times.
✓ Branch 6 taken 75112 times.
✓ Branch 7 taken 11163487 times.
 282831143 || (!flag && abs(c-d) < tol*max(abs(c), abs(d))) )
154 {
155 221300419 s = (a+b)*0.5;
156 221300419 flag = true;
157 }
158 else
159 flag = false;
160
161 // compute residual at new guess
162 250983060 fs = residual(s);
163 250983060 d = c;
164 250983060 c = b;
165 250983060 fc = fb;
166
167 // check on which side of the root s falls
168
2/2
✓ Branch 0 taken 215386353 times.
✓ Branch 1 taken 35596707 times.
 250983060 if (fa*fs < 0.0)
169 {
170 b = s;
171 fb = fs;
172 }
173 else
174 {
175 215386353 a = s;
176 215386353 fa = fs;
177 }
178
179 // sort if necessary
180
2/2
✓ Branch 0 taken 34966868 times.
✓ Branch 1 taken 216016192 times.
 250983060 if (abs(fa) < abs(fb))
181 {
182 34966868 swap(a, b);
183 34966868 swap(fa, fb);
184 }
185 }
186
187 ✗ DUNE_THROW(NumericalProblem, "Brent's algorithm didn't converge after " << maxIter << " iterations!");
188 }
189
190 } // end namespace Dumux
191
192 #endif
193