3.6-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
doubleexpintegrator.hh
Go to the documentation of this file.
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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 * *
19 * This version is modified after the original version by John D. Cook *
20 * see https://www.codeproject.com/ *
21 * Articles/31550/Fast-Numerical-Integration *
22 * which is licensed under BSD-2-clause, which reads as follows: *
23 * Copyright John D. Cook *
24 * *
25 * Redistribution and use in source and binary forms, with or without *
26 * modification, are permitted provided that the following *
27 * conditions are met: *
28 * 1. Redistributions of source code must retain the above *
29 * copyright notice, this list of conditions *
30 and the following disclaimer. *
31 * 2. Redistributions in binary form must reproduce the above *
32 * copyright notice, this list of conditions *
33 * and the following disclaimer *
34 * in the documentation and/or other materials *
35 * provided with the distribution. *
36 * *
37 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS *
38 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT *
39 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS *
40 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE *
41 * COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, *
42 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL *
43 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE *
44 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS *
45 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, *
46 * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE *
47 * OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, *
48 * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
49 * *
50 *****************************************************************************/
56#ifndef DUMUX_COMMON_DOUBLEEXP_INTEGRATOR_HH
57#define DUMUX_COMMON_DOUBLEEXP_INTEGRATOR_HH
58
59#include <cmath>
60#include <limits>
61#include <type_traits>
62
64
65namespace Dumux {
66
70template<class Scalar>
72{
73public:
84 template<class Function,
85 typename std::enable_if_t<std::is_invocable_r_v<Scalar, Function, Scalar>>...>
86 static Scalar integrate(const Function& f,
87 const Scalar a,
88 const Scalar b,
89 const Scalar targetAbsoluteError,
90 int& numFunctionEvaluations,
91 Scalar& errorEstimate)
92 {
93 // Apply the linear change of variables x = ct + d
94 // $$\int_a^b f(x) dx = c \int_{-1}^1 f( ct + d ) dt$$
95 // c = (b-a)/2, d = (a+b)/2
96
97 const Scalar c = 0.5*(b - a);
98 const Scalar d = 0.5*(a + b);
99 return integrateCore_(f, c, d, targetAbsoluteError, numFunctionEvaluations, errorEstimate,
101 }
102
113 template<class Function,
114 typename std::enable_if_t<std::is_invocable_r_v<Scalar, Function, Scalar>>...>
115 static Scalar integrate(const Function& f,
116 const Scalar a,
117 const Scalar b,
118 const Scalar targetAbsoluteError)
119 {
120 int numFunctionEvaluations;
121 Scalar errorEstimate;
122 return integrate(f, a, b, targetAbsoluteError, numFunctionEvaluations, errorEstimate);
123 }
124
125private:
126 // Integrate f(cx + d) with the given integration constants
127 template<class Function>
128 static Scalar integrateCore_(const Function& f,
129 const Scalar c, // slope of change of variables
130 const Scalar d, // intercept of change of variables
131 Scalar targetAbsoluteError,
132 int& numFunctionEvaluations,
133 Scalar& errorEstimate,
134 const double* abcissas,
135 const double* weights)
136 {
137 targetAbsoluteError /= c;
138
139 // Offsets to where each level's integration constants start.
140 // The last element is not a beginning but an end.
141 static const int offsets[] = {1, 4, 7, 13, 25, 49, 97, 193};
142 static const int numLevels = sizeof(offsets)/sizeof(int) - 1;
143
144 Scalar newContribution = 0.0;
145 Scalar integral = 0.0;
146 Scalar h = 1.0;
147 errorEstimate = std::numeric_limits<Scalar>::max();
148 Scalar previousDelta, currentDelta = std::numeric_limits<Scalar>::max();
149
150 integral = f(c*abcissas[0] + d) * weights[0];
151 int i;
152 for (i = offsets[0]; i != offsets[1]; ++i)
153 integral += weights[i]*(f(c*abcissas[i] + d) + f(-c*abcissas[i] + d));
154
155 for (int level = 1; level != numLevels; ++level)
156 {
157 h *= 0.5;
158 newContribution = 0.0;
159 for (i = offsets[level]; i != offsets[level+1]; ++i)
160 newContribution += weights[i]*(f(c*abcissas[i] + d) + f(-c*abcissas[i] + d));
161 newContribution *= h;
162
163 // difference in consecutive integral estimates
164 previousDelta = currentDelta;
165 using std::abs;
166 currentDelta = abs(0.5*integral - newContribution);
167 integral = 0.5*integral + newContribution;
168
169 // Once convergence kicks in, error is approximately squared at each step.
170 // Determine whether we're in the convergent region by looking at the trend in the error.
171 if (level == 1)
172 continue; // previousDelta meaningless, so cannot check convergence.
173
174 // Exact comparison with zero is harmless here. Could possibly be replaced with
175 // a small positive upper limit on the size of currentDelta, but determining
176 // that upper limit would be difficult. At worse, the loop is executed more
177 // times than necessary. But no infinite loop can result since there is
178 // an upper bound on the loop variable.
179 if (currentDelta == 0.0)
180 break;
181
182 using std::log;
183 const Scalar rate = log( currentDelta )/log( previousDelta ); // previousDelta != 0 or would have been kicked out previously
184
185 if (rate > 1.9 && rate < 2.1)
186 {
187 // If convergence theory applied perfectly, r would be 2 in the convergence region.
188 // r close to 2 is good enough. We expect the difference between this integral estimate
189 // and the next one to be roughly delta^2.
190 errorEstimate = currentDelta*currentDelta;
191 }
192 else
193 {
194 // Not in the convergence region. Assume only that error is decreasing.
195 errorEstimate = currentDelta;
196 }
197
198 if (errorEstimate < 0.1*targetAbsoluteError)
199 break;
200 }
201
202 numFunctionEvaluations = 2*i - 1;
203 errorEstimate *= c;
204 return c*integral;
205 }
206};
207
208} // end namespace Dumux
209
210#endif
Exponential integration abcissas.
Adaption of the non-isothermal two-phase two-component flow model to problems with CO2.
Definition: adapt.hh:29
constexpr double doubleExponentialIntegrationWeights[]
Definition: doubleexpintegrationconstants.hh:296
constexpr double doubleExponentialIntegrationAbcissas[]
Definition: doubleexpintegrationconstants.hh:92
Numerical integration in one dimension using the double exponential method of M. Mori.
Definition: doubleexpintegrator.hh:72
static Scalar integrate(const Function &f, const Scalar a, const Scalar b, const Scalar targetAbsoluteError, int &numFunctionEvaluations, Scalar &errorEstimate)
Integrate an analytic function over a finite interval.
Definition: doubleexpintegrator.hh:86
static Scalar integrate(const Function &f, const Scalar a, const Scalar b, const Scalar targetAbsoluteError)
Integrate an analytic function over a finite interval.
Definition: doubleexpintegrator.hh:115