diff --git a/applications/test/MathFunctions/Make/files b/applications/test/MathFunctions/Make/files
new file mode 100644
index 0000000000..e32f24c285
--- /dev/null
+++ b/applications/test/MathFunctions/Make/files
@@ -0,0 +1,3 @@
+Test-MathFunctions.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-MathFunctions
diff --git a/applications/test/MathFunctions/Make/options b/applications/test/MathFunctions/Make/options
new file mode 100644
index 0000000000..a134fe4eb4
--- /dev/null
+++ b/applications/test/MathFunctions/Make/options
@@ -0,0 +1 @@
+EXE_INC = -I../TestTools
diff --git a/applications/test/MathFunctions/Test-MathFunctions.C b/applications/test/MathFunctions/Test-MathFunctions.C
new file mode 100644
index 0000000000..e1ac201693
--- /dev/null
+++ b/applications/test/MathFunctions/Test-MathFunctions.C
@@ -0,0 +1,179 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+Application
+ Test-MathFunctions
+
+Description
+ Tests for \c Math namespace member functions
+ using \c doubleScalar base type.
+
+\*---------------------------------------------------------------------------*/
+
+#include "MathFunctions.H"
+#include "mathematicalConstants.H"
+#include "IOmanip.H"
+#include "TestTools.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Execute each member functions of Foam::Math::, and print output
+template
+void test_member_funcs(Type)
+{
+ Random rndGen(1234);
+ const label numberOfTests = 10000;
+ const scalar tolerance = 1e-3;
+
+ for (label i = 0; i < numberOfTests; ++i)
+ {
+ Info<< "# Inverse error functions:" << endl;
+ const Type x1 = rndGen.sample01();
+ const Type x2 = rndGen.sample01();
+ Info<< " # Operands:" << nl
+ << " x1 = " << x1 << nl
+ << " x2 = " << x2 << endl;
+
+ cmp
+ (
+ " # erf(erfinv(x1)) = x1 = ",
+ x1,
+ Foam::erf(Foam::Math::erfInv(x1)),
+ tolerance
+ );
+
+ cmp
+ (
+ " # erf(erfinv(-x2)) = -x2 = ",
+ x2,
+ Foam::erf(Foam::Math::erfInv(x2)),
+ tolerance
+ );
+
+ Info<< "# Incomplete gamma functions:" << endl;
+ const Type a = rndGen.position(SMALL, Type(100));
+ const Type x = rndGen.position(Type(0), Type(100));
+ Info<< " # Operands:" << nl
+ << " a = " << a << nl
+ << " x = " << x << endl;
+
+ cmp
+ (
+ " # incGammaRatio_Q(a,x) + incGammaRatio_P(a,x) = 1 = ",
+ Foam::Math::incGammaRatio_Q(a,x) + Foam::Math::incGammaRatio_P(a,x),
+ scalar(1),
+ tolerance
+ );
+
+ cmp
+ (
+ " # incGamma_Q(1, x) = exp(-x) = ",
+ Foam::Math::incGamma_Q(1, x),
+ Foam::exp(-x),
+ tolerance
+ );
+
+ cmp
+ (
+ " # incGamma_Q(0.5, x) = sqrt(pi)*erfc(sqrt(x)) = ",
+ Foam::Math::incGamma_Q(0.5, x),
+ Foam::sqrt(Foam::constant::mathematical::pi)
+ *Foam::erfc(Foam::sqrt(x)),
+ tolerance
+ );
+
+ cmp
+ (
+ " # incGamma_P(1,x) = 1 - exp(x) = ",
+ Foam::Math::incGamma_P(1, x),
+ 1 - Foam::exp(-x),
+ tolerance
+ );
+
+ cmp
+ (
+ " # incGamma_P(0.5, x) = sqrt(pi)*erf(sqrt(x)) = ",
+ Foam::Math::incGamma_P(0.5, x),
+ Foam::sqrt(Foam::constant::mathematical::pi)
+ *Foam::erf(Foam::sqrt(x)),
+ tolerance
+ );
+ }
+}
+
+
+// Do compile-time recursion over the given types
+template
+inline typename std::enable_if::type
+run_tests(const std::tuple& types, const List& typeID){}
+
+
+template
+inline typename std::enable_if::type
+run_tests(const std::tuple& types, const List& typeID)
+{
+ Info<< nl << " ## Test member functions: "<< typeID[I] <<" ##" << nl;
+ test_member_funcs(std::get(types));
+
+ run_tests(types, typeID);
+}
+
+
+// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * //
+
+int main()
+{
+ Info<< setprecision(15);
+
+ const std::tuple types
+ (
+ std::make_tuple(Zero)
+ );
+
+ const List typeID
+ ({
+ "doubleScalar"
+ });
+
+ run_tests(types, typeID);
+
+
+ if (nFail_)
+ {
+ Info<< nl << " #### "
+ << "Failed in " << nFail_ << " tests "
+ << "out of total " << nTest_ << " tests "
+ << "####\n" << endl;
+ return 1;
+ }
+
+ Info<< nl << " #### Passed all " << nTest_ <<" tests ####\n" << endl;
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/TestTools/TestTools.H b/applications/test/TestTools/TestTools.H
index c24cb9fa05..3fed3b7d5d 100644
--- a/applications/test/TestTools/TestTools.H
+++ b/applications/test/TestTools/TestTools.H
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2020 OpenCFD Ltd.
+ Copyright (C) 2020-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -30,6 +30,10 @@ Description
using namespace Foam;
+#include "floatScalar.H"
+#include "doubleScalar.H"
+#include "complex.H"
+#include "Matrix.H"
#include "Random.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 5382d32670..8a66579f1f 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -60,8 +60,6 @@ $(ints)/lists/labelListIOList.C
primitives/Scalar/doubleScalar/doubleScalar.C
primitives/Scalar/floatScalar/floatScalar.C
primitives/Scalar/scalar/scalar.C
-primitives/Scalar/scalar/incGamma.C
-primitives/Scalar/scalar/invIncGamma.C
primitives/Scalar/lists/scalarList.C
primitives/Scalar/lists/scalarIOList.C
primitives/Scalar/lists/scalarListIOList.C
@@ -117,6 +115,12 @@ primitives/functions/Function1/quarterCosineRamp/quarterCosineRamp.C
primitives/functions/Function1/halfCosineRamp/halfCosineRamp.C
primitives/functions/Polynomial/polynomialFunction.C
+/* Math functions */
+math = primitives/functions/Math
+$(math)/erfInv.C
+$(math)/incGamma.C
+$(math)/invIncGamma.C
+
primitives/subModelBase/subModelBase.C
strings = primitives/strings
diff --git a/src/OpenFOAM/primitives/Scalar/scalar/scalar.H b/src/OpenFOAM/primitives/Scalar/scalar/scalar.H
index c947c05df2..58bb002e98 100644
--- a/src/OpenFOAM/primitives/Scalar/scalar/scalar.H
+++ b/src/OpenFOAM/primitives/Scalar/scalar/scalar.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
- Copyright (C) 2017-2020 OpenCFD Ltd.
+ Copyright (C) 2017-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -213,25 +213,9 @@ inline float narrowFloat(const double val)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-// Additional transcendental functions and specialisations
namespace Foam
{
- //- Inverse normalized incomplete gamma function
- scalar invIncGamma(const scalar a, const scalar P);
-
- //- Normalized upper incomplete gamma function
- scalar incGammaRatio_Q(const scalar a, const scalar x);
-
- //- Normalized lower incomplete gamma function
- scalar incGammaRatio_P(const scalar a, const scalar x);
-
- //- Upper incomplete gamma function
- scalar incGamma_Q(const scalar a, const scalar x);
-
- //- Lower incomplete gamma function
- scalar incGamma_P(const scalar a, const scalar x);
-
//- Type to use for extended precision
template<>
class typeOfSolve
diff --git a/src/OpenFOAM/primitives/functions/Math/MathFunctions.H b/src/OpenFOAM/primitives/functions/Math/MathFunctions.H
new file mode 100644
index 0000000000..eb0757cc83
--- /dev/null
+++ b/src/OpenFOAM/primitives/functions/Math/MathFunctions.H
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+Namespace
+ Foam::Math
+
+Description
+ A namespace for various mathematical functions.
+
+ Reference:
+ \verbatim
+ Inverse error function (tag:W):
+ Winitzki, S. (2008).
+ A handy approximation for the error function and its inverse.
+ A lecture note obtained through private communication.
+ URL:https://sites.google.com/site/winitzki/sergei-winitzkis-files
+ (Retrieved on: 16 Feb 2021).
+
+ Incomplete gamma functions (tag:DM):
+ DiDonato, A. R., & Morris Jr, A. H. (1986).
+ Computation of the incomplete gamma
+ function ratios and their inverse.
+ ACM Transactions on Mathematical Software (TOMS), 12(4), 377-393.
+ DOI:10.1145/22721.23109
+ \endverbatim
+
+Note
+ - The algorithm in \c invIncGamma is described in (DM:Sec. 4).
+ - The algorithm in \c incGammaRatio_Q is described in (DM:Sec. 3).
+ - The accuracy parameter \c IND is set to a value of 1.
+
+SourceFiles
+ erfInv.C
+ incGamma.C
+ invIncGamma.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MathFunctions_H
+#define MathFunctions_H
+
+#include "scalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Namespace Math Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace Math
+{
+
+//- Inverse error function of a real-number argument
+// \param y Real-number argument at which to evaluate. Domain: (-1, 1)
+// \return The inverse of error function of y
+scalar erfInv(const scalar y);
+
+
+// Incomplete gamma functions
+
+//- Inverse of regularised lower incomplete gamma function
+// \param a Real-number argument. Domain: (0, infty]
+// \param P Real-number argument. Domain: [0,1]
+scalar invIncGamma(const scalar a, const scalar P);
+
+//- Regularised upper incomplete gamma function
+// \param a Real-number argument. Domain: (0, infty]
+// \param x Real-number argument. Domain: [0, infty]
+scalar incGammaRatio_Q(const scalar a, const scalar x);
+
+//- Regularised lower incomplete gamma function
+// \param a Real-number argument. Domain: (0, infty]
+// \param x Real-number argument. Domain: [0, infty]
+scalar incGammaRatio_P(const scalar a, const scalar x);
+
+//- Upper incomplete gamma function
+// \param a Real-number argument. Domain: (0, infty]
+// \param x Real-number argument. Domain: [0, infty]
+scalar incGamma_Q(const scalar a, const scalar x);
+
+//- Lower incomplete gamma function
+// \param a Real-number argument. Domain: (0, infty]
+// \param x Real-number argument. Domain: [0, infty]
+scalar incGamma_P(const scalar a, const scalar x);
+
+
+} // End namespace Math
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/functions/Math/erfInv.C b/src/OpenFOAM/primitives/functions/Math/erfInv.C
new file mode 100644
index 0000000000..ccce80dbec
--- /dev/null
+++ b/src/OpenFOAM/primitives/functions/Math/erfInv.C
@@ -0,0 +1,68 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+\*---------------------------------------------------------------------------*/
+
+#include "MathFunctions.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * //
+
+Foam::scalar Foam::Math::erfInv(const scalar y)
+{
+ #ifdef FULLDEBUG
+ if (mag(y) >= scalar(1))
+ {
+ WarningInFunction
+ << "The domain of inverse error function argument "
+ << "(i.e. y) should be limited to (-1, 1):" << nl
+ << " y = " << y
+ << endl;
+
+ return std::numeric_limits::infinity();
+ }
+ #endif
+
+ // (W:p. 2) to reduce the max relative error to O(1e-4)
+ constexpr scalar a = 0.147;
+
+ const scalar k =
+ scalar(2)/(a*constant::mathematical::pi) + 0.5*log(scalar(1) - sqr(y));
+
+ const scalar h = log(scalar(1) - sqr(y))/a;
+
+ // (W:Eq. 7)
+ const scalar x = sqrt(-k + sqrt(sqr(k) - h));
+
+ if (y < 0)
+ {
+ return -x;
+ }
+
+ return x;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Scalar/scalar/incGamma.C b/src/OpenFOAM/primitives/functions/Math/incGamma.C
similarity index 64%
rename from src/OpenFOAM/primitives/Scalar/scalar/incGamma.C
rename to src/OpenFOAM/primitives/functions/Math/incGamma.C
index f6fcafaf28..fcee586dda 100644
--- a/src/OpenFOAM/primitives/Scalar/scalar/incGamma.C
+++ b/src/OpenFOAM/primitives/functions/Math/incGamma.C
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019 OpenFOAM Foundation
+ Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -24,38 +25,23 @@ License
along with OpenFOAM. If not, see .
Global
- Foam::incGamma
+ Foam::Math::incGamma
Description
- Calculates the upper and lower incomplete gamma functions as well as their
- normalized versions.
-
- The algorithm is described in detail in DiDonato et al. (1986).
-
- \verbatim
- DiDonato, A. R., & Morris Jr, A. H. (1986).
- Computation of the incomplete gamma function ratios and their inverse.
- ACM Transactions on Mathematical Software (TOMS), 12(4), 377-393.
- \endverbatim
-
- All equation numbers in the following code refer to the above paper.
- The algorithm in function 'incGammaRatio_Q' is described in section 3.
- The accuracy parameter IND is set to a value of 1.
+ Implementation of the incomplete gamma functions.
\*---------------------------------------------------------------------------*/
+#include "MathFunctions.H"
#include "mathematicalConstants.H"
+#include
-using namespace Foam::constant::mathematical;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
namespace Foam
{
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// Eqn. (13)
+// (DM:Eq. 13)
static scalar calcQE11(const scalar a, const scalar x, const int e = 30)
{
scalar a_2n = 0;
@@ -88,7 +74,7 @@ static scalar calcQE11(const scalar a, const scalar x, const int e = 30)
}
-// Eqn. (15)
+// (DM:Eq. 15)
static scalar calcPE15(const scalar a, const scalar x, const int nmax = 20)
{
scalar prod = 1;
@@ -106,7 +92,7 @@ static scalar calcPE15(const scalar a, const scalar x, const int nmax = 20)
}
-// Eq. (16)
+// (DM:Eq. 16)
static scalar calcQE16(const scalar a, const scalar x, const int N = 20)
{
scalar an = 1;
@@ -124,7 +110,7 @@ static scalar calcQE16(const scalar a, const scalar x, const int N = 20)
}
-// Eq. (18)
+// (DM:Eq. 18)
static scalar calcTE18
(
const scalar a,
@@ -135,55 +121,46 @@ static scalar calcTE18
const scalar phi
)
{
- static const scalar D0[] =
- {
- -0.333333333333333E-00,
- 0.833333333333333E-01,
- -0.148148148148148E-01,
- 0.115740740740741E-02,
- 0.352733686067019E-03,
- -0.178755144032922E-03,
- 0.391926317852244E-04,
- -0.218544851067999E-05,
- -0.185406221071516E-05,
- 0.829671134095309E-06,
- -0.176659527368261E-06,
- 0.670785354340150E-08,
- 0.102618097842403E-07,
- -0.438203601845335E-08
- };
+ constexpr scalar D0_0 = -0.333333333333333E-00;
+ constexpr scalar D0_1 = 0.833333333333333E-01;
+ constexpr scalar D0_2 = -0.148148148148148E-01;
+ constexpr scalar D0_3 = 0.115740740740741E-02;
+ constexpr scalar D0_4 = 0.352733686067019E-03;
+ constexpr scalar D0_5 = -0.178755144032922E-03;
+ constexpr scalar D0_6 = 0.391926317852244E-04;
+ // unused: constexpr scalar D0_7 = -0.218544851067999E-05;
+ // unused: constexpr scalar D0_8 = -0.185406221071516E-05;
+ // unused: constexpr scalar D0_9 = 0.829671134095309E-06;
+ // unused: constexpr scalar D0_10 = -0.176659527368261E-06;
+ // unused: constexpr scalar D0_11 = 0.670785354340150E-08;
+ // unused: constexpr scalar D0_12 = 0.102618097842403E-07;
+ // unused: constexpr scalar D0_13 = -0.438203601845335E-08;
- static const scalar D1[] =
- {
- -0.185185185185185E-02,
- -0.347222222222222E-02,
- 0.264550264550265E-02,
- -0.990226337448560E-03,
- 0.205761316872428E-03,
- -0.401877572016461E-06,
- -0.180985503344900E-04,
- 0.764916091608111E-05,
- -0.161209008945634E-05,
- 0.464712780280743E-08,
- 0.137863344691572E-06,
- -0.575254560351770E-07,
- 0.119516285997781E-07
- };
+ constexpr scalar D1_0 = -0.185185185185185E-02;
+ constexpr scalar D1_1 = -0.347222222222222E-02;
+ constexpr scalar D1_2 = 0.264550264550265E-02;
+ constexpr scalar D1_3 = -0.990226337448560E-03;
+ constexpr scalar D1_4 = 0.205761316872428E-03;
+ // unused: constexpr scalar D1_5 = -0.401877572016461E-06;
+ // unused: constexpr scalar D1_6 = -0.180985503344900E-04;
+ // unused: constexpr scalar D1_7 = 0.764916091608111E-05;
+ // unused: constexpr scalar D1_8 = -0.161209008945634E-05;
+ // unused: constexpr scalar D1_9 = 0.464712780280743E-08;
+ // unused: constexpr scalar D1_10 = 0.137863344691572E-06;
+ // unused: constexpr scalar D1_11 = -0.575254560351770E-07;
+ // unused: constexpr scalar D1_12 = 0.119516285997781E-07;
- static const scalar D2[] =
- {
- 0.413359788359788E-02,
- -0.268132716049383E-02,
- 0.771604938271605E-03,
- 0.200938786008230E-05,
- -0.107366532263652E-03,
- 0.529234488291201E-04,
- -0.127606351886187E-04,
- 0.342357873409614E-07,
- 0.137219573090629E-05,
- -0.629899213838006E-06,
- 0.142806142060642E-06
- };
+ constexpr scalar D2_0 = 0.413359788359788E-02;
+ constexpr scalar D2_1 = -0.268132716049383E-02;
+ // unused: constexpr scalar D2_2 = 0.771604938271605E-03;
+ // unused: constexpr scalar D2_3 = 0.200938786008230E-05;
+ // unused: constexpr scalar D2_4 = -0.107366532263652E-03;
+ // unused: constexpr scalar D2_5 = 0.529234488291201E-04;
+ // unused: constexpr scalar D2_6 = -0.127606351886187E-04;
+ // unused: constexpr scalar D2_7 = 0.342357873409614E-07;
+ // unused: constexpr scalar D2_8 = 0.137219573090629E-05;
+ // unused: constexpr scalar D2_9 = -0.629899213838006E-06;
+ // unused: constexpr scalar D2_10 = 0.142806142060642E-06;
const scalar u = 1/a;
scalar z = sqrt(2*phi);
@@ -196,44 +173,66 @@ static scalar calcTE18
if (sigma > (e0/sqrt(a)))
{
const scalar C0 =
- D0[6]*pow6(z) + D0[5]*pow5(z) + D0[4]*pow4(z)
- + D0[3]*pow3(z) + D0[2]*sqr(z) + D0[1]*z + D0[0];
+ D0_6*pow6(z) + D0_5*pow5(z) + D0_4*pow4(z)
+ + D0_3*pow3(z) + D0_2*sqr(z) + D0_1*z + D0_0;
const scalar C1 =
- D1[4]*pow4(z) + D1[3]*pow3(z) + D1[2]*sqr(z) + D1[1]*z + D1[0];
+ D1_4*pow4(z) + D1_3*pow3(z) + D1_2*sqr(z) + D1_1*z + D1_0;
- const scalar C2 = D2[1]*z + D2[0];
+ const scalar C2 = D2_1*z + D2_0;
return C2*sqr(u) + C1*u + C0;
}
else
{
- const scalar C0 = D0[2]*sqr(z) + D0[1]*z + D0[0];
- const scalar C1 = D1[1]*z + D1[0];
- const scalar C2 = D2[1]*z + D2[0];
+ const scalar C0 = D0_2*sqr(z) + D0_1*z + D0_0;
+ const scalar C1 = D1_1*z + D1_0;
+ const scalar C2 = D2_1*z + D2_0;
return C2*sqr(u) + C1*u + C0;
}
}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
} // End namespace Foam
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-Foam::scalar Foam::incGammaRatio_Q(const scalar a, const scalar x)
+Foam::scalar Foam::Math::incGammaRatio_Q(const scalar a, const scalar x)
{
- const scalar BIG = 14;
- const scalar x0 = 17;
- const scalar e0 = 0.025;
+ using namespace Foam::constant::mathematical;
+
+ #ifdef FULLDEBUG
+ if (a <= 0)
+ {
+ WarningInFunction
+ << "The parameter (i.e. a) cannot be negative or zero"
+ << " a = " << a
+ << endl;
+
+ return std::numeric_limits::infinity();
+ }
+
+ if (x < 0)
+ {
+ WarningInFunction
+ << "The parameter (i.e. x) cannot be negative"
+ << " x = " << x
+ << endl;
+
+ return std::numeric_limits::infinity();
+ }
+ #endif
+
+ constexpr scalar BIG = 14;
+ constexpr scalar x0 = 17;
+ constexpr scalar e0 = 0.025;
if (a < 1)
{
if (a == 0.5)
{
- // Eqn. (8)
+ // (DM:Eq. 8)
if (x < 0.25)
{
return 1 - erf(sqrt(x));
@@ -245,7 +244,7 @@ Foam::scalar Foam::incGammaRatio_Q(const scalar a, const scalar x)
}
else if ( x < 1.1)
{
- // Eqn. (12)
+ // (DM:Eq. 12)
scalar alpha = x/2.59;
if (x < 0.5)
@@ -264,12 +263,12 @@ Foam::scalar Foam::incGammaRatio_Q(const scalar a, const scalar x)
if (a > alpha || a == alpha)
{
- // Eqn. (9)
+ // (DM:Eq. 9)
return 1 - (pow(x, a)*(1 - J))/tgamma(a + 1);
}
else
{
- // Eqn. (10)
+ // (DM:Eq. 10)
const scalar L = exp(a*log(x)) - 1;
const scalar H = 1/(tgamma(a + 1)) - 1;
@@ -278,7 +277,7 @@ Foam::scalar Foam::incGammaRatio_Q(const scalar a, const scalar x)
}
else
{
- // Eqn. (11)
+ // (DM:Eq. 11)
const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
return R*calcQE11(a, x);
@@ -290,7 +289,7 @@ Foam::scalar Foam::incGammaRatio_Q(const scalar a, const scalar x)
if (sigma <= e0/sqrt(a))
{
- // Eqn. (19)
+ // (DM:Eq. 19)
const scalar lambda = x/a;
const scalar phi = lambda - 1 - log(lambda);
const scalar y = a*phi;
@@ -319,7 +318,7 @@ Foam::scalar Foam::incGammaRatio_Q(const scalar a, const scalar x)
{
if (sigma <= 0.4)
{
- // Eqn. (17)
+ // (DM:Eq. 17)
const scalar lambda = x/a;
const scalar phi = lambda - 1 - log(lambda);
const scalar y = a*phi;
@@ -344,19 +343,19 @@ Foam::scalar Foam::incGammaRatio_Q(const scalar a, const scalar x)
{
if (x <= max(a, log(10.0)))
{
- // Eqn. (15)
+ // (DM:Eq. 15)
return 1 - calcPE15(a, x);
}
else if (x < x0)
{
- // Eqn. (11)
+ // (DM:Eq. 11)
const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
return R*calcQE11(a, x);
}
else
{
- // Eqn. (16)
+ // (DM:Eq. 16)
return calcQE16(a, x);
}
}
@@ -368,19 +367,19 @@ Foam::scalar Foam::incGammaRatio_Q(const scalar a, const scalar x)
{
if (x <= max(a, log(10.0)))
{
- // Eqn. (15)
+ // (DM:Eq. 15)
return 1 - calcPE15(a, x);
}
else if ( x < x0)
{
- // Eqn. (11)
+ // (DM:Eq. 11)
const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
return R*calcQE11(a, x);
}
else
{
- // Eqn. (16)
+ // (DM:Eq. 16)
return calcQE16(a, x);
}
}
@@ -388,7 +387,7 @@ Foam::scalar Foam::incGammaRatio_Q(const scalar a, const scalar x)
{
if (floor(2*a) == 2*a)
{
- // Eqn. (14)
+ // (DM:Eq. 14)
if (floor(a) == a)
{
scalar sum = 0;
@@ -417,19 +416,19 @@ Foam::scalar Foam::incGammaRatio_Q(const scalar a, const scalar x)
}
else if (x <= max(a, log(10.0)))
{
- // Eqn. (15)
+ // (DM:Eq. 15)
return 1 - calcPE15(a, x);
}
- else if ( x < x0)
+ else if (x < x0)
{
- // Eqn. (11)
+ // (DM:Eq. 11)
const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
return R*calcQE11(a, x);
}
else
{
- // Eqn. (16)
+ // (DM:Eq. 16)
return calcQE16(a, x);
}
}
@@ -437,19 +436,19 @@ Foam::scalar Foam::incGammaRatio_Q(const scalar a, const scalar x)
}
-Foam::scalar Foam::incGammaRatio_P(const scalar a, const scalar x)
+Foam::scalar Foam::Math::incGammaRatio_P(const scalar a, const scalar x)
{
return 1 - incGammaRatio_Q(a, x);
}
-Foam::scalar Foam::incGamma_Q(const scalar a, const scalar x)
+Foam::scalar Foam::Math::incGamma_Q(const scalar a, const scalar x)
{
return incGammaRatio_Q(a, x)*tgamma(a);
}
-Foam::scalar Foam::incGamma_P(const scalar a, const scalar x)
+Foam::scalar Foam::Math::incGamma_P(const scalar a, const scalar x)
{
return incGammaRatio_P(a, x)*tgamma(a);
}
diff --git a/src/OpenFOAM/primitives/Scalar/scalar/invIncGamma.C b/src/OpenFOAM/primitives/functions/Math/invIncGamma.C
similarity index 82%
rename from src/OpenFOAM/primitives/Scalar/scalar/invIncGamma.C
rename to src/OpenFOAM/primitives/functions/Math/invIncGamma.C
index 5ca1412cb1..3a38e3d7c3 100644
--- a/src/OpenFOAM/primitives/Scalar/scalar/invIncGamma.C
+++ b/src/OpenFOAM/primitives/functions/Math/invIncGamma.C
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016 OpenFOAM Foundation
+ Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -24,60 +25,43 @@ License
along with OpenFOAM. If not, see .
Global
- Foam::invIncGamma
+ Foam::Math::invIncGamma
Description
- Function to calculate the inverse of the
- normalized incomplete gamma function.
-
- The algorithm is described in detain in reference:
- \verbatim
- DiDonato, A. R., & Morris Jr, A. H. (1986).
- Computation of the incomplete gamma function ratios and their inverse.
- ACM Transactions on Mathematical Software (TOMS), 12(4), 377-393.
- \endverbatim
-
- All equation numbers in the following code refer to the above text and
- the algorithm in the function 'invIncGamma' is described in section 4.
+ Implementation of the inverse incomplete gamma function.
\*---------------------------------------------------------------------------*/
+#include "MathFunctions.H"
#include "mathematicalConstants.H"
+
using namespace Foam::constant::mathematical;
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
namespace Foam
{
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
static scalar minimaxs(const scalar P)
{
- // Eqn. 32
+ // (DM:Eq. 32)
- static const scalar a[] =
- {
- 3.31125922108741,
- 11.6616720288968,
- 4.28342155967104,
- 0.213623493715853
- };
+ constexpr scalar a_0 = 3.31125922108741;
+ constexpr scalar a_1 = 11.6616720288968;
+ constexpr scalar a_2 = 4.28342155967104;
+ constexpr scalar a_3 = 0.213623493715853;
- static const scalar b[] =
- {
- 6.61053765625462,
- 6.40691597760039,
- 1.27364489782223,
- 0.03611708101884203
- };
+ constexpr scalar b_0 = 6.61053765625462;
+ constexpr scalar b_1 = 6.40691597760039;
+ constexpr scalar b_2 = 1.27364489782223;
+ constexpr scalar b_3 = 0.03611708101884203;
const scalar t = P < 0.5 ? sqrt(-2*log(P)) : sqrt(-2*log(1 - P));
const scalar s =
t
- - (a[0] + t*(a[1] + t*(a[2] + t*a[3])))
- /(1 + t*(b[0] + t*(b[1] + t*(b[2] + t*b[3]))));
+ - (a_0 + t*(a_1 + t*(a_2 + t*a_3)))
+ /(1 + t*(b_0 + t*(b_1 + t*(b_2 + t*b_3))));
return P < 0.5 ? -s : s;
}
@@ -85,12 +69,12 @@ static scalar minimaxs(const scalar P)
static scalar Sn(const scalar a, const scalar x)
{
- //- Eqn. 34
+ // (DM:Eq. 34)
scalar Sn = 1;
scalar Si = 1;
- for (int i=1; i<100; i++)
+ for (int i=1; i<100; ++i)
{
Si *= x/(a + i);
Sn += Si;
@@ -101,15 +85,35 @@ static scalar Sn(const scalar a, const scalar x)
return Sn;
}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
} // End namespace Foam
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-//- Inverse normalized incomplete gamma function
-Foam::scalar Foam::invIncGamma(const scalar a, const scalar P)
+Foam::scalar Foam::Math::invIncGamma(const scalar a, const scalar P)
{
+ #ifdef FULLDEBUG
+ if (a <= 0)
+ {
+ WarningInFunction
+ << "The parameter (i.e. a) cannot be negative or zero"
+ << " a = " << a
+ << endl;
+
+ return std::numeric_limits::infinity();
+ }
+
+ if (P < 0 || P > 1)
+ {
+ WarningInFunction
+ << "The domain of the parameter (i.e. P) should be limited to [0,1]"
+ << " P = " << P
+ << endl;
+
+ return std::numeric_limits::infinity();
+ }
+ #endif
+
const scalar Q = 1 - P;
if (a == 1)
@@ -123,7 +127,7 @@ Foam::scalar Foam::invIncGamma(const scalar a, const scalar P)
if (B > 0.6 || (B >= 0.45 && a >= 0.3))
{
- // Eqn. 21
+ // (DM:Eq. 21)
const scalar u =
(B*Q > 1e-8) ? pow(P*Ga*a, 1/a) : exp((-Q/a) - Eu);
@@ -131,7 +135,7 @@ Foam::scalar Foam::invIncGamma(const scalar a, const scalar P)
}
else if (a < 0.3 && B >= 0.35)
{
- // Eqn. 22
+ // (DM:Eq. 22)
const scalar t = exp(-Eu - B);
const scalar u = t*exp(t);
@@ -139,7 +143,7 @@ Foam::scalar Foam::invIncGamma(const scalar a, const scalar P)
}
else if (B > 0.15 || a >= 0.3)
{
- // Eqn. 23
+ // (DM:Eq. 23)
const scalar y = -log(B);
const scalar u = y - (1 - a)*log(y);
@@ -147,7 +151,7 @@ Foam::scalar Foam::invIncGamma(const scalar a, const scalar P)
}
else if (B > 0.1)
{
- // Eqn. 24
+ // (DM:Eq. 24)
const scalar y = -log(B);
const scalar u = y - (1 - a)*log(y);
@@ -161,7 +165,7 @@ Foam::scalar Foam::invIncGamma(const scalar a, const scalar P)
}
else
{
- // Eqn. 25:
+ // (DM:Eq. 25)
const scalar y = -log(B);
const scalar c1 = (a - 1)*log(y);
const scalar c12 = c1*c1;
@@ -194,7 +198,7 @@ Foam::scalar Foam::invIncGamma(const scalar a, const scalar P)
}
else
{
- // Eqn. 31:
+ // (DM:Eq. 31)
scalar s = minimaxs(P);
const scalar s2 = sqr(s);
@@ -227,7 +231,7 @@ Foam::scalar Foam::invIncGamma(const scalar a, const scalar P)
if (lnB < -2.3*D)
{
- // Eqn. 25:
+ // (DM:Eq. 25)
const scalar y = -lnB;
const scalar c1 = (a - 1)*log(y);
const scalar c12 = c1*c1;
@@ -270,7 +274,7 @@ Foam::scalar Foam::invIncGamma(const scalar a, const scalar P)
}
else
{
- // Eqn. 33
+ // (DM:Eq. 33)
const scalar u =
-lnB + (a - 1)*log(w) - log(1 + (1 - a)/(1 + w));
@@ -285,7 +289,7 @@ Foam::scalar Foam::invIncGamma(const scalar a, const scalar P)
if (w < 0.15*ap1)
{
- // Eqn. 35
+ // (DM:Eq. 35)
const scalar ap2 = a + 2;
const scalar v = log(P) + lgamma(ap1);
z = exp((v + w)/a);
@@ -303,7 +307,7 @@ Foam::scalar Foam::invIncGamma(const scalar a, const scalar P)
}
else
{
- // Eqn. 36
+ // (DM:Eq. 36)
const scalar lnSn = log(Sn(a, z));
const scalar v = log(P) + lgamma(ap1);
z = exp((v + z - lnSn)/a);
diff --git a/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C b/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C
index e021f2f4f6..130e4a9759 100644
--- a/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C
+++ b/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016 OpenFOAM Foundation
+ Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -26,6 +27,7 @@ License
\*---------------------------------------------------------------------------*/
#include "massRosinRammler.H"
+#include "MathFunctions.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -87,7 +89,7 @@ Foam::scalar Foam::distributionModels::massRosinRammler::sample() const
{
const scalar a = 3/n_ + 1;
const scalar P = rndGen_.sample01();
- const scalar x = invIncGamma(a, P);
+ const scalar x = Math::invIncGamma(a, P);
d = d_*pow(x, 1/n_);
} while (d < minValue_ || d > maxValue_);
diff --git a/src/lagrangian/distributionModels/normal/normal.C b/src/lagrangian/distributionModels/normal/normal.C
index dfc63497c7..803c555b2e 100644
--- a/src/lagrangian/distributionModels/normal/normal.C
+++ b/src/lagrangian/distributionModels/normal/normal.C
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
+ Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -27,7 +28,7 @@ License
#include "normal.H"
#include "addToRunTimeSelectionTable.H"
-#include "mathematicalConstants.H"
+#include "MathFunctions.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -99,7 +100,7 @@ Foam::scalar Foam::distributionModels::normal::sample() const
scalar b = erf((maxValue_ - expectation_)/variance_);
scalar y = rndGen_.sample01();
- scalar x = erfInv(y*(b - a) + a)*variance_ + expectation_;
+ scalar x = Math::erfInv(y*(b - a) + a)*variance_ + expectation_;
// Note: numerical approximation of the inverse function yields slight
// inaccuracies
@@ -128,17 +129,4 @@ Foam::scalar Foam::distributionModels::normal::meanValue() const
}
-Foam::scalar Foam::distributionModels::normal::erfInv(const scalar y) const
-{
- scalar k = 2.0/(constant::mathematical::pi*a_) + 0.5*log(1.0 - y*y);
- scalar h = log(1.0 - y*y)/a_;
- scalar x = sqrt(-k + sqrt(k*k - h));
- if (y < 0.0)
- {
- x *= -1.0;
- }
- return x;
-}
-
-
// ************************************************************************* //
diff --git a/src/lagrangian/distributionModels/normal/normal.H b/src/lagrangian/distributionModels/normal/normal.H
index cb689490e9..f98206b25b 100644
--- a/src/lagrangian/distributionModels/normal/normal.H
+++ b/src/lagrangian/distributionModels/normal/normal.H
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2013 OpenFOAM Foundation
+ Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -60,8 +61,7 @@ class normal
:
public distributionModel
{
- // Private data
-
+ // Private Data
//- Distribution minimum
scalar minValue_;
@@ -116,8 +116,6 @@ public:
//- Return the mean value
virtual scalar meanValue() const;
-
- virtual scalar erfInv(const scalar y) const;
};
diff --git a/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C b/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
index 32754939ff..2f65c0cef9 100644
--- a/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
+++ b/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
+ Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -35,25 +36,6 @@ using namespace Foam::constant;
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
-template
-Foam::scalar Foam::BrownianMotionForce::erfInv(const scalar y) const
-{
- const scalar a = 0.147;
- scalar k = 2.0/(mathematical::pi*a) + 0.5*log(1.0 - y*y);
- scalar h = log(1.0 - y*y)/a;
- scalar x = sqrt(-k + sqrt(k*k - h));
-
- if (y < 0.0)
- {
- return -x;
- }
- else
- {
- return x;
- }
-}
-
-
template
Foam::tmp
Foam::BrownianMotionForce::kModel() const
@@ -200,7 +182,7 @@ Foam::forceSuSp Foam::BrownianMotionForce::calcCoupled
// for (direction dir = 0; dir < vector::nComponents; dir++)
// {
// const scalar x = rndGen_.sample01();
- // const scalar eta = sqrt2*erfInv(2*x - 1.0);
+ // const scalar eta = sqrt2*Math::erfInv(2*x - 1.0);
// value.Su()[dir] = f*eta;
// }
diff --git a/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H b/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H
index 751a83d820..9405a4473d 100644
--- a/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H
+++ b/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
+ Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -84,9 +85,6 @@ class BrownianMotionForce
// Private Member Functions
- //- Inverse erf for Gaussian distribution
- scalar erfInv(const scalar y) const;
-
//- Return the k field from the turbulence model
tmp kModel() const;
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/binaryBreakupModels/LuoSvendsen/LuoSvendsen.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/binaryBreakupModels/LuoSvendsen/LuoSvendsen.C
index 66e190345b..479b083bab 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/binaryBreakupModels/LuoSvendsen/LuoSvendsen.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/binaryBreakupModels/LuoSvendsen/LuoSvendsen.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019 OpenFOAM Foundation
- Copyright (C) 2020 OpenCFD Ltd.
+ Copyright (C) 2020-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -30,6 +30,7 @@ License
#include "addToRunTimeSelectionTable.H"
#include "phaseCompressibleTurbulenceModel.H"
#include "tableBounds.H"
+#include "MathFunctions.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -99,19 +100,19 @@ Foam::diameterModels::binaryBreakupModels::LuoSvendsen::LuoSvendsen
Tuple2 gamma2by11
(
z,
- incGammaRatio_Q(2.0/11.0, z)
+ Math::incGammaRatio_Q(2.0/11.0, z)
);
Tuple2 gamma5by11
(
z,
- incGammaRatio_Q(5.0/11.0, z)
+ Math::incGammaRatio_Q(5.0/11.0, z)
);
Tuple2 gamma8by11
(
z,
- incGammaRatio_Q(8.0/11.0, z)
+ Math::incGammaRatio_Q(8.0/11.0, z)
);
gammaUpperReg2by11Table.append(gamma2by11);