GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/common/doubleexpintegrator.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 34 37 91.9%
Functions: 9 13 69.2%
Branches: 18 18 100.0%
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-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 // SPDX-FileCopyrightInfo: Copyright John D. Cook
7 // SPDX-License-Identifier: BSD-2-Clause
8 // Modified after the original version at
9 // https://www.codeproject.com/Articles/31550/Fast-Numerical-Integration licensed under
10 // BSD-2-Clause. All Changes are licensed under GPL-3.0-or-later.
11 //
12 /*****************************************************************************
13 * This version is modified after the original version by John D. Cook *
14 * see https://www.codeproject.com/ *
15 * Articles/31550/Fast-Numerical-Integration *
16 * which is licensed under BSD-2-clause, which reads as follows: *
17 * Copyright John D. Cook *
18 * *
19 * Redistribution and use in source and binary forms, with or without *
20 * modification, are permitted provided that the following *
21 * conditions are met: *
22 * 1. Redistributions of source code must retain the above *
23 * copyright notice, this list of conditions *
24 and the following disclaimer. *
25 * 2. Redistributions in binary form must reproduce the above *
26 * copyright notice, this list of conditions *
27 * and the following disclaimer *
28 * in the documentation and/or other materials *
29 * provided with the distribution. *
30 * *
31 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS *
32 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT *
33 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS *
34 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE *
35 * COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, *
36 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL *
37 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE *
38 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS *
39 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, *
40 * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE *
41 * OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, *
42 * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
43 * *
44 *****************************************************************************/
45 /*!
46 * \file
47 * \ingroup Core
48 * \brief A double exponential integrator
49 */
50 #ifndef DUMUX_COMMON_DOUBLEEXP_INTEGRATOR_HH
51 #define DUMUX_COMMON_DOUBLEEXP_INTEGRATOR_HH
52
53 #include <cmath>
54 #include <limits>
55 #include <type_traits>
56
57 #include <dumux/common/doubleexpintegrationconstants.hh>
58
59 namespace Dumux {
60
61 /*!
62 * \brief Numerical integration in one dimension using the double exponential method of M. Mori.
63 */
64 template<class Scalar>
65 class DoubleExponentialIntegrator
66 {
67 public:
68 /*!
69 * \brief Integrate an analytic function over a finite interval
70 * \param f the integrand (invocable with a single scalar)
71 * \param a lower limit of integration
72 * \param b upper limit of integration
73 * \param targetAbsoluteError desired bound on error
74 * \param numFunctionEvaluations number of function evaluations used
75 * \param errorEstimate estimate for error in integration
76 * \return The value of the integral
77 */
78 template<class Function,
79 typename std::enable_if_t<std::is_invocable_r_v<Scalar, Function, Scalar>>...>
80 ✗ static Scalar integrate(const Function& f,
81 const Scalar a,
82 const Scalar b,
83 const Scalar targetAbsoluteError,
84 int& numFunctionEvaluations,
85 Scalar& errorEstimate)
86 {
87 // Apply the linear change of variables x = ct + d
88 // $$\int_a^b f(x) dx = c \int_{-1}^1 f( ct + d ) dt$$
89 // c = (b-a)/2, d = (a+b)/2
90
91 2093485 const Scalar c = 0.5*(b - a);
92 2093485 const Scalar d = 0.5*(a + b);
93 2093485 return integrateCore_(f, c, d, targetAbsoluteError, numFunctionEvaluations, errorEstimate,
94 ✗ doubleExponentialIntegrationAbcissas, doubleExponentialIntegrationWeights);
95 }
96
97 /*!
98 * \brief Integrate an analytic function over a finite interval.
99 * \note This version overloaded to not require arguments passed in for
100 * function evaluation counts or error estimates.
101 * \param f the integrand (invocable with a single scalar)
102 * \param a lower integral bound
103 * \param b upper integral bound
104 * \param targetAbsoluteError desired absolute error in the result
105 * \return The value of the integral.
106 */
107 template<class Function,
108 typename std::enable_if_t<std::is_invocable_r_v<Scalar, Function, Scalar>>...>
109 2123726 static Scalar integrate(const Function& f,
110 const Scalar a,
111 const Scalar b,
112 const Scalar targetAbsoluteError)
113 {
114 int numFunctionEvaluations;
115 Scalar errorEstimate;
116 2123726 return integrate(f, a, b, targetAbsoluteError, numFunctionEvaluations, errorEstimate);
117 }
118
119 private:
120 // Integrate f(cx + d) with the given integration constants
121 template<class Function>
122 2093485 static Scalar integrateCore_(const Function& f,
123 const Scalar c, // slope of change of variables
124 const Scalar d, // intercept of change of variables
125 Scalar targetAbsoluteError,
126 int& numFunctionEvaluations,
127 Scalar& errorEstimate,
128 const double* abcissas,
129 const double* weights)
130 {
131 2093485 targetAbsoluteError /= c;
132
133 // Offsets to where each level's integration constants start.
134 // The last element is not a beginning but an end.
135 static const int offsets[] = {1, 4, 7, 13, 25, 49, 97, 193};
136 static const int numLevels = sizeof(offsets)/sizeof(int) - 1;
137
138 2093485 Scalar newContribution = 0.0;
139 2093485 Scalar integral = 0.0;
140 2093485 Scalar h = 1.0;
141 2093485 errorEstimate = std::numeric_limits<Scalar>::max();
142 2093485 Scalar previousDelta, currentDelta = std::numeric_limits<Scalar>::max();
143
144 2093485 integral = f(c*abcissas[0] + d) * weights[0];
145 int i;
146
2/2
✓ Branch 0 taken 6189732 times.
✓ Branch 1 taken 2063244 times.
 8373940 for (i = offsets[0]; i != offsets[1]; ++i)
147 6280455 integral += weights[i]*(f(c*abcissas[i] + d) + f(-c*abcissas[i] + d));
148
149
2/2
✓ Branch 0 taken 12199089 times.
✓ Branch 1 taken 2015783 times.
 14414638 for (int level = 1; level != numLevels; ++level)
150 {
151 12372895 h *= 0.5;
152 12372895 newContribution = 0.0;
153
2/2
✓ Branch 0 taken 381689136 times.
✓ Branch 1 taken 12199089 times.
 399272284 for (i = offsets[level]; i != offsets[level+1]; ++i)
154 386899389 newContribution += weights[i]*(f(c*abcissas[i] + d) + f(-c*abcissas[i] + d));
155 12372895 newContribution *= h;
156
157 // difference in consecutive integral estimates
158 12372895 previousDelta = currentDelta;
159 using std::abs;
160
2/2
✓ Branch 0 taken 2063244 times.
✓ Branch 1 taken 10135845 times.
 12372895 currentDelta = abs(0.5*integral - newContribution);
161 12372895 integral = 0.5*integral + newContribution;
162
163 // Once convergence kicks in, error is approximately squared at each step.
164 // Determine whether we're in the convergent region by looking at the trend in the error.
165
2/2
✓ Branch 0 taken 2063244 times.
✓ Branch 1 taken 10135845 times.
 12372895 if (level == 1)
166 ✗ continue; // previousDelta meaningless, so cannot check convergence.
167
168 // Exact comparison with zero is harmless here. Could possibly be replaced with
169 // a small positive upper limit on the size of currentDelta, but determining
170 // that upper limit would be difficult. At worse, the loop is executed more
171 // times than necessary. But no infinite loop can result since there is
172 // an upper bound on the loop variable.
173
2/2
✓ Branch 0 taken 10133526 times.
✓ Branch 1 taken 2319 times.
 10279410 if (currentDelta == 0.0)
174 break;
175
176 using std::log;
177 10277090 const Scalar rate = log( currentDelta )/log( previousDelta ); // previousDelta != 0 or would have been kicked out previously
178
179
4/4
✓ Branch 0 taken 1158 times.
✓ Branch 1 taken 10132368 times.
✓ Branch 2 taken 1130 times.
✓ Branch 3 taken 28 times.
 10277090 if (rate > 1.9 && rate < 2.1)
180 {
181 // If convergence theory applied perfectly, r would be 2 in the convergence region.
182 // r close to 2 is good enough. We expect the difference between this integral estimate
183 // and the next one to be roughly delta^2.
184 2701 errorEstimate = currentDelta*currentDelta;
185 }
186 else
187 {
188 // Not in the convergence region. Assume only that error is decreasing.
189 10274389 errorEstimate = currentDelta;
190 }
191
192
2/2
✓ Branch 0 taken 10088384 times.
✓ Branch 1 taken 45142 times.
 10277090 if (errorEstimate < 0.1*targetAbsoluteError)
193 break;
194 }
195
196 2093485 numFunctionEvaluations = 2*i - 1;
197 2093485 errorEstimate *= c;
198 2093485 return c*integral;
199 }
200 };
201
202 } // end namespace Dumux
203
204 #endif
205