zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

invIncGammaRatio_P: Additional divide by zero protection

Browse Source
This commit is contained in:
Will Bainbridge
2023-04-04 16:42:34 +01:00
parent 9ea964a525
commit 0080177d88
2 changed files with 13 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
applications/test/incGamma/Test-incGamma.C
View File

@ -51,6 +51,8 @@ scalar grade(const scalar f)
int main(int argc, char *argv[])
{
argList args(argc, argv);
const scalarField as({0.25, 0.5, 1, 2, 4});
static const scalar nXs = 10000;
@ -92,7 +94,15 @@ int main(int argc, char *argv[])
{
const scalar P = incGammaRatio_P(as[ai], xs[xi]);
const scalar PStar =
incGammaRatio_P(as[ai], invIncGammaRatio_P(as[ai], P));
incGammaRatio_P
(
as[ai],
invIncGammaRatio_P
(
as[ai],
min(max(P, small), 1 - small)
)
);
scalarPs[ai][xi] = P;
scalarPErrors[ai][xi] = mag(P - PStar);
}

4
src/OpenFOAM/primitives/Scalar/scalar/invIncGamma.C
View File

@ -349,8 +349,8 @@ Foam::scalar Foam::invIncGammaRatio_P(const scalar a, const scalar P)
for (label iter = 0; iter < nIter; ++ iter)
{
const scalar dP = incGammaRatio_P(a, x) - P;
const scalar t = dP/R(a, x);
const scalar dP = incGammaRatio_P(a, max(x, vSmall)) - P;
const scalar t = dP/max(R(a, x), vSmall);
const scalar w = (a - 1 - x)/2;
bool halley = mag(t) < 0.1 && mag(w*t) < 0.1;