zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

BUG: polynomialEqns: fix discriminant limit

Browse Source 
It was observed in a MPPICDyMFoam simulation involving a single particle
in a moving mesh that the barocentric trajectory of the particle follows
an unexpected path at some arbitrary instant in time.

The issue was tracked to "hitEqn()" where cubicEqn/quadraticEqn computes
one of the roots wrongly due to the discriminant limit we set, e.g. for:

    0x^3 + 1.4334549e-33 x^2 - 9.0869006e-10 x + 0.0027666538

Although the discriminant limit was carefully selected to avoid various
problems at the time, the new change is required more due to its exposition
to a wider spectrum of applications.
This commit is contained in:
Kutalmis Bercin
2020-09-14 17:25:55 +01:00
parent b42db6cee5
commit bbeda07862
2 changed files with 2 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C
View File

@ -60,7 +60,7 @@ Foam::Roots<3> Foam::cubicEqn::roots() const
const scalar p = -(fma(-a, c, w) + fma(b, b/3.0, -w));
const scalar q = b*b*b*scalar(2)/27 - b*c*a/3 + d*a*a;
const scalar numDiscr = p*p*p/27 + q*q/4;
const scalar discr = (mag(numDiscr) > sqr(SMALL)) ? numDiscr : 0;
const scalar discr = (mag(numDiscr) > VSMALL) ? numDiscr : 0;
// Determine the number and types of the roots
const bool threeReal = discr < 0;

2
src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C
View File

@ -46,7 +46,7 @@ Foam::Roots<2> Foam::quadraticEqn::roots() const
// (JLM:p. 2246) [discriminant = b*b/4 - a*c]
const scalar w = a*c;
const scalar numDiscr = fma(-a, c, w) + fma(b, b/4, -w);
const scalar discr = (mag(numDiscr) > SMALL) ? numDiscr : 0;
const scalar discr = (mag(numDiscr) > VSMALL) ? numDiscr : 0;
// Find how many roots of what types are available
const bool twoReal = discr > 0;