From 5e7a766159e0a2f8a1267b8968cf2802befe6476 Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Tue, 24 Jul 2018 15:54:32 +0100
Subject: [PATCH] BUG: cubicEqn, quadraticEqn: Correction to repeated roots
Also extended the cubic equation test routine and modified the error
methods so that they more accurately generate the round of error of
evaluation.
This resolves bug report https://bugs.openfoam.org/view.php?id=3015
---
applications/test/cubicEqn/Test-cubicEqn.C | 27 +++++++++++++++----
.../polynomialEqns/cubicEqn/cubicEqn.C | 4 +--
.../polynomialEqns/cubicEqn/cubicEqnI.H | 5 +++-
.../polynomialEqns/linearEqn/linearEqnI.H | 2 +-
.../quadraticEqn/quadraticEqn.C | 2 +-
.../quadraticEqn/quadraticEqnI.H | 4 ++-
6 files changed, 33 insertions(+), 11 deletions(-)
diff --git a/applications/test/cubicEqn/Test-cubicEqn.C b/applications/test/cubicEqn/Test-cubicEqn.C
index d0ab6066662..aa0579a4194 100644
--- a/applications/test/cubicEqn/Test-cubicEqn.C
+++ b/applications/test/cubicEqn/Test-cubicEqn.C
@@ -51,7 +51,7 @@ void test(const Type& polynomialEqn, const scalar tol)
case roots::real:
v[i] = polynomialEqn.value(r[i]);
e[i] = polynomialEqn.error(r[i]);
- ok = ok && mag(v[i]) < tol*mag(e[i]);
+ ok = ok && mag(v[i]) <= tol*mag(e[i]);
break;
case roots::posInf:
v[i] = + inf;
@@ -79,8 +79,10 @@ void test(const Type& polynomialEqn, const scalar tol)
int main()
{
- const int nTests = 1000000;
- for (int t = 0; t < nTests; ++ t)
+ const scalar tol = 5;
+
+ const label nTests = 1000000;
+ for (label t = 0; t < nTests; ++ t)
{
test
(
@@ -91,11 +93,26 @@ int main()
randomScalar(1e-50, 1e+50),
randomScalar(1e-50, 1e+50)
),
- 100
+ tol
);
}
+ Info << nTests << " random cubics tested" << endl;
- Info << nTests << " cubics tested" << endl;
+ const label coeffMin = -9, coeffMax = 10, nCoeff = coeffMax - coeffMin;
+ for (label a = coeffMin; a < coeffMax; ++ a)
+ {
+ for (label b = coeffMin; b < coeffMax; ++ b)
+ {
+ for (label c = coeffMin; c < coeffMax; ++ c)
+ {
+ for (label d = coeffMin; d < coeffMax; ++ d)
+ {
+ test(cubicEqn(a, b, c, d), tol);
+ }
+ }
+ }
+ }
+ Info<< nCoeff*nCoeff*nCoeff*nCoeff << " integer cubics tested" << endl;
return 0;
}
diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C
index c46b794c38e..3aad6740018 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C
+++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C
@@ -81,7 +81,7 @@ Foam::Roots<3> Foam::cubicEqn::roots() const
// This is assumed not to over- or under-flow. If it does, all bets are off.
const scalar p = c*a - b*b/3;
- const scalar q = b*b*b*(2.0/27.0) - b*c*a/3 + d*a*a;
+ const scalar q = b*b*b*scalar(2)/27 - b*c*a/3 + d*a*a;
const scalar disc = p*p*p/27 + q*q/4;
// How many roots of what types are available?
@@ -96,7 +96,7 @@ Foam::Roots<3> Foam::cubicEqn::roots() const
if (oneReal)
{
- const Roots<1> r = linearEqn(- a, b/3).roots();
+ const Roots<1> r = linearEqn(a, b/3).roots();
return Roots<3>(r.type(0), r[0]);
}
else if (twoReal)
diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H
index da5fc3cba14..cb528b6be4a 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H
@@ -114,7 +114,10 @@ inline Foam::scalar Foam::cubicEqn::derivative(const scalar x) const
inline Foam::scalar Foam::cubicEqn::error(const scalar x) const
{
- return mag(SMALL*x*(x*(x*3*a() + 2*b()) + c()));
+ return
+ SMALL*magSqr(x)*(mag(x*a()) + mag(b()))
+ + SMALL*mag(x)*(mag(x*(x*a() + b())) + mag(c()))
+ + SMALL*(mag(x*(x*(x*a() + b()) + c())) + mag(d()));
}
diff --git a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H
index 17bfdf452da..f5c4f2b7c94 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H
@@ -82,7 +82,7 @@ inline Foam::scalar Foam::linearEqn::derivative(const scalar x) const
inline Foam::scalar Foam::linearEqn::error(const scalar x) const
{
- return mag(SMALL*x*a());
+ return SMALL*(mag(x*a()) + mag(b()));
}
diff --git a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C
index d84c2928c25..d2ac0ea8a12 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C
+++ b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C
@@ -71,7 +71,7 @@ Foam::Roots<2> Foam::quadraticEqn::roots() const
if (oneReal)
{
- const Roots<1> r = linearEqn(- a, b/2).roots();
+ const Roots<1> r = linearEqn(a, b/2).roots();
return Roots<2>(r, r);
}
else if (twoReal)
diff --git a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H
index 7672a74d724..15ffe4e4e22 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H
+++ b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H
@@ -100,7 +100,9 @@ inline Foam::scalar Foam::quadraticEqn::derivative(const scalar x) const
inline Foam::scalar Foam::quadraticEqn::error(const scalar x) const
{
- return mag(SMALL*x*(x*2*a() + b()));
+ return
+ SMALL*mag(x)*(mag(x*a()) + mag(b()))
+ + SMALL*(mag(x*(x*a() + b())) + mag(c()));
}
--
GitLab