From bbeda078620013b75ed8a62f2f0a8a6154d5d8df Mon Sep 17 00:00:00 2001 From: Kutalmis Bercin Date: Mon, 14 Sep 2020 17:25:55 +0100 Subject: [PATCH] BUG: polynomialEqns: fix discriminant limit 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. --- src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C | 2 +- .../primitives/polynomialEqns/quadraticEqn/quadraticEqn.C | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C index 17732e91dc..0b2c1f0825 100644 --- a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C +++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C @@ -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; diff --git a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C index 910f9375ce..1ef4e1cc6f 100644 --- a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C +++ b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C @@ -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;