| 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 Core | ||
| 10 | * \brief A class for numeric differentiation | ||
| 11 | * | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_NUMERIC_DIFFERENTIATION_HH | ||
| 14 | #define DUMUX_NUMERIC_DIFFERENTIATION_HH | ||
| 15 | |||
| 16 | #include <cmath> | ||
| 17 | #include <cassert> | ||
| 18 | #include <limits> | ||
| 19 | |||
| 20 | namespace Dumux { | ||
| 21 | |||
| 22 | /*! | ||
| 23 | * \ingroup Core | ||
| 24 | * \brief A class for numeric differentiation with respect to a scalar parameter | ||
| 25 | */ | ||
| 26 | class NumericDifferentiation | ||
| 27 | { | ||
| 28 | public: | ||
| 29 | |||
| 30 | /*! | ||
| 31 | * \brief Computes the epsilon used for numeric differentiation | ||
| 32 | * \param value The value of the variable with respect to which we are differentiating | ||
| 33 | * \param baseEps The step width which we are using for differentiation | ||
| 34 | */ | ||
| 35 | template<class Scalar> | ||
| 36 |
1/2✓ Branch 1 taken 502 times.
✗ Branch 2 not taken.
|
357544507 | static Scalar epsilon(const Scalar value, const Scalar baseEps = 1e-10) |
| 37 | { | ||
| 38 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 309586381 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 47887336 times.
|
357473717 | assert(std::numeric_limits<Scalar>::epsilon()*1e4 < baseEps); |
| 39 | // the epsilon value used for the numeric differentiation is | ||
| 40 | // now scaled by the absolute value of the primary variable... | ||
| 41 | using std::abs; | ||
| 42 |
1/2✓ Branch 1 taken 502 times.
✗ Branch 2 not taken.
|
357474220 | return baseEps*(abs(value) + 1.0); |
| 43 | } | ||
| 44 | |||
| 45 | /*! | ||
| 46 | * \brief Computes the derivative of a function with respect to a function parameter | ||
| 47 | * \note Overload using default epsilon computation | ||
| 48 | */ | ||
| 49 | template<class Function, class Scalar, class FunctionEvalType> | ||
| 50 |
0/2✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
70288 | static void partialDerivative(const Function& function, Scalar x0, |
| 51 | FunctionEvalType& derivative, | ||
| 52 | const FunctionEvalType& fx0, | ||
| 53 | const int numericDifferenceMethod = 1) | ||
| 54 |
0/2✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
70288 | { partialDerivative(function, x0, derivative, fx0, epsilon(x0), numericDifferenceMethod); } |
| 55 | |||
| 56 | /*! | ||
| 57 | * \brief Computes the derivative of a function with respect to a function parameter | ||
| 58 | * \param function The function to derive | ||
| 59 | * \param x0 The parameter at which the derivative is ought to be evaluated | ||
| 60 | * \param derivative The partial derivative (output) | ||
| 61 | * \param fx0 The result of the function evaluated at x0 | ||
| 62 | * \param eps The numeric epsilon used in the differentiation | ||
| 63 | * \param numericDifferenceMethod The numeric difference method | ||
| 64 | * (1: forward differences (default), 0: central differences, -1: backward differences, 5: five-point stencil method) | ||
| 65 | */ | ||
| 66 | template<class Function, class Scalar, class FunctionEvalType> | ||
| 67 | 370742724 | static void partialDerivative(const Function& function, Scalar x0, | |
| 68 | FunctionEvalType& derivative, | ||
| 69 | const FunctionEvalType& fx0, | ||
| 70 | const Scalar eps, | ||
| 71 | const int numericDifferenceMethod = 1) | ||
| 72 | { | ||
| 73 | // Five-point stencil numeric difference, | ||
| 74 | // Abramowitz & Stegun, Table 25.2. | ||
| 75 | // The error is proportional to eps^4. | ||
| 76 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 116477799 times.
|
370742724 | if (numericDifferenceMethod == 5) |
| 77 | { | ||
| 78 | 2 | derivative = function(x0 + eps); | |
| 79 | 2 | derivative -= function(x0 - eps); | |
| 80 | 2 | derivative *= 8.0; | |
| 81 | 2 | derivative += function(x0 - 2*eps); | |
| 82 | 2 | derivative -= function(x0 + 2*eps); | |
| 83 | 2 | derivative /= 12*eps; | |
| 84 | 266186487 | return; | |
| 85 | } | ||
| 86 | |||
| 87 | // Forward, central, or backward differences | ||
| 88 | 370742722 | Scalar delta = 0.0; | |
| 89 | |||
| 90 | // we are using forward or central differences, i.e. we need to calculate f(x + \epsilon) | ||
| 91 |
1/2✓ Branch 0 taken 116477799 times.
✗ Branch 1 not taken.
|
370742722 | if (numericDifferenceMethod >= 0) |
| 92 | { | ||
| 93 | 370742720 | delta += eps; | |
| 94 | // calculate the function evaluated with the deflected variable | ||
| 95 | 370742720 | derivative = function(x0 + eps); | |
| 96 | } | ||
| 97 | |||
| 98 | // we are using backward differences, | ||
| 99 | // i.e. we don't need to calculate f(x + \epsilon) | ||
| 100 | // we can recycle the (possibly cached) f(x) | ||
| 101 | 2 | else derivative = fx0; | |
| 102 | |||
| 103 | // we are using backward or central differences, | ||
| 104 | // i.e. we need to calculate f(x - \epsilon) | ||
| 105 |
2/2✓ Branch 0 taken 1538880 times.
✓ Branch 1 taken 114938919 times.
|
370742722 | if (numericDifferenceMethod <= 0) |
| 106 | { | ||
| 107 | 81494920 | delta += eps; | |
| 108 | // subtract the function evaluated with the deflected variable | ||
| 109 | 101269906 | derivative -= function(x0 - eps); | |
| 110 | } | ||
| 111 | |||
| 112 | // we are using forward differences, | ||
| 113 | // i.e. we don't need to calculate f(x - \epsilon) | ||
| 114 | // we can recycle the (possibly cached) f(x) | ||
| 115 | 370741120 | else derivative -= fx0; | |
| 116 | |||
| 117 | // divide difference in residuals by the magnitude of the | ||
| 118 | // deflections between the two function evaluation | ||
| 119 | 370742723 | derivative /= delta; | |
| 120 | } | ||
| 121 | }; | ||
| 122 | |||
| 123 | } // end namespace Dumux | ||
| 124 | |||
| 125 | #endif | ||
| 126 |