doubleexpintegrator.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
// SPDX-FileCopyrightInfo: Copyright John D. Cook
// SPDX-License-Identifier: BSD-2-Clause
// Modified after the original version at
// https://www.codeproject.com/Articles/31550/Fast-Numerical-Integration licensed under
// BSD-2-Clause. All Changes are licensed under GPL-3.0-or-later.
//
/*****************************************************************************
 *   This version is modified after the original version by John D. Cook     *
 *   see https://www.codeproject.com/                                        *
 *       Articles/31550/Fast-Numerical-Integration                           *
 *   which is licensed under BSD-2-clause, which reads as follows:           *
 *   Copyright John D. Cook                                                  *
 *                                                                           *
 *   Redistribution and use in source and binary forms, with or without      *
 *   modification, are permitted provided that the following                 *
 *   conditions are met:                                                     *
 *   1. Redistributions of source code must retain the above                 *
 *      copyright notice, this list of conditions                            *
        and the following disclaimer.                                        *
 *   2. Redistributions in binary form must reproduce the above              *
 *      copyright notice, this list of conditions                            *
 *      and the following disclaimer                                         *
 *   in the documentation and/or other materials                             *
 *   provided with the distribution.                                         *
 *                                                                           *
 *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS     *
 *   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT       *
 *   LIMITED TO, THE IMPLIED  WARRANTIES OF MERCHANTABILITY AND FITNESS      *
 *   FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE          *
 *   COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,              *
 *   INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL              *
 *   DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE       *
 *   GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS           *
 *   INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,            *
 *   WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE    *
 *   OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,       *
 *   EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                      *
 *                                                                           *
 *****************************************************************************/
#ifndef DUMUX_COMMON_DOUBLEEXP_INTEGRATOR_HH
#define DUMUX_COMMON_DOUBLEEXP_INTEGRATOR_HH
 
#include <cmath>
#include <limits>
#include <type_traits>
 
#include <dumux/common/doubleexpintegrationconstants.hh>
 
namespace Dumux {
 
template<class Scalar>
class DoubleExponentialIntegrator
{
public:
    template<class Function,
             typename std::enable_if_t<std::is_invocable_r_v<Scalar, Function, Scalar>>...>
    static Scalar integrate(const Function& f,
                            const Scalar a,
                            const Scalar b,
                            const Scalar targetAbsoluteError,
                            int& numFunctionEvaluations,
                            Scalar& errorEstimate)
    {
        // Apply the linear change of variables x = ct + d
        // $$\int_a^b f(x) dx = c \int_{-1}^1 f( ct + d ) dt$$
        // c = (b-a)/2, d = (a+b)/2
 
        const Scalar c = 0.5*(b - a);
        const Scalar d = 0.5*(a + b);
        return integrateCore_(f, c, d, targetAbsoluteError, numFunctionEvaluations, errorEstimate,
                              doubleExponentialIntegrationAbcissas, doubleExponentialIntegrationWeights);
    }
 
    template<class Function,
             typename std::enable_if_t<std::is_invocable_r_v<Scalar, Function, Scalar>>...>
    static Scalar integrate(const Function& f,
                            const Scalar a,
                            const Scalar b,
                            const Scalar targetAbsoluteError)
    {
        int numFunctionEvaluations;
        Scalar errorEstimate;
        return integrate(f, a, b, targetAbsoluteError, numFunctionEvaluations, errorEstimate);
    }
 
private:
    // Integrate f(cx + d) with the given integration constants
    template<class Function>
    static Scalar integrateCore_(const Function& f,
                                 const Scalar c,   // slope of change of variables
                                 const Scalar d,   // intercept of change of variables
                                 Scalar targetAbsoluteError,
                                 int& numFunctionEvaluations,
                                 Scalar& errorEstimate,
                                 const double* abcissas,
                                 const double* weights)
    {
        targetAbsoluteError /= c;
 
        // Offsets to where each level's integration constants start.
        // The last element is not a beginning but an end.
        static const int offsets[] = {1, 4, 7, 13, 25, 49, 97, 193};
        static const int numLevels = sizeof(offsets)/sizeof(int) - 1;
 
        Scalar newContribution = 0.0;
        Scalar integral = 0.0;
        Scalar h = 1.0;
        errorEstimate = std::numeric_limits<Scalar>::max();
        Scalar previousDelta, currentDelta = std::numeric_limits<Scalar>::max();
 
        integral = f(c*abcissas[0] + d) * weights[0];
        int i;
        for (i = offsets[0]; i != offsets[1]; ++i)
            integral += weights[i]*(f(c*abcissas[i] + d) + f(-c*abcissas[i] + d));
 
        for (int level = 1; level != numLevels; ++level)
        {
            h *= 0.5;
            newContribution = 0.0;
            for (i = offsets[level]; i != offsets[level+1]; ++i)
                newContribution += weights[i]*(f(c*abcissas[i] + d) + f(-c*abcissas[i] + d));
            newContribution *= h;
 
            // difference in consecutive integral estimates
            previousDelta = currentDelta;
            using std::abs;
            currentDelta = abs(0.5*integral - newContribution);
            integral = 0.5*integral + newContribution;
 
            // Once convergence kicks in, error is approximately squared at each step.
            // Determine whether we're in the convergent region by looking at the trend in the error.
            if (level == 1)
                continue; // previousDelta meaningless, so cannot check convergence.
 
            // Exact comparison with zero is harmless here.  Could possibly be replaced with
            // a small positive upper limit on the size of currentDelta, but determining
            // that upper limit would be difficult.  At worse, the loop is executed more
            // times than necessary.  But no infinite loop can result since there is
            // an upper bound on the loop variable.
            if (currentDelta == 0.0)
                break;
 
            using std::log;
            const Scalar rate = log( currentDelta )/log( previousDelta );  // previousDelta != 0 or would have been kicked out previously
 
            if (rate > 1.9 && rate < 2.1)
            {
                // If convergence theory applied perfectly, r would be 2 in the convergence region.
                // r close to 2 is good enough. We expect the difference between this integral estimate
                // and the next one to be roughly delta^2.
                errorEstimate = currentDelta*currentDelta;
            }
            else
            {
                // Not in the convergence region.  Assume only that error is decreasing.
                errorEstimate = currentDelta;
            }
 
            if (errorEstimate < 0.1*targetAbsoluteError)
                break;
        }
 
        numFunctionEvaluations = 2*i - 1;
        errorEstimate *= c;
        return c*integral;
    }
};
 
} // end namespace Dumux
 
#endif