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

ENH: improve stability in polynomialEqns

Browse Source 
- replaces floating-point equal comparisons in
    `linearEqn`, `quadraticEqn`, and `cubicEqn`,
  - ensures `quadraticEqn` and `cubicEqn` can return `complex` roots,
  - reorders if-branches in `quadraticEqn` and `cubicEqn` to avoid
    zero-equal comparison,
  - adds Kahan's cancellation-avoiding algorithm into `quadraticEqn` and
    `cubicEqn` for the numerically-sensitive discriminant computation,

  - adds/improves `polynomialEqns` tests:
    * adds Test-linearEqn.C
    * adds Test-quadraticEqn.C
    * improves Test-cubicEqn.C
This commit is contained in:
Kutalmis Bercin
2020-01-14 12:32:36 +00:00
committed by Andrew Heather
parent 97bdd5bc04
commit 8ca724fffa
16 changed files with 1184 additions and 287 deletions
Download Patch File Download Diff File Expand all files Collapse all files

118
applications/test/cubicEqn/Test-cubicEqn.C
View File

@ -1,118 +0,0 @@
#include <ctime>
#include <random>
#include "cubicEqn.H"
#include "IOstreams.H"
#include "stringList.H"
using namespace Foam;
scalar randomScalar(const scalar min, const scalar max)
{
static_assert
(
sizeof(long) == sizeof(scalar),
"Scalar and long are not the same size"
);
static std::default_random_engine generator(std::time(0));
static std::uniform_int_distribution<long>
distribution
(
std::numeric_limits<long>::min(),
std::numeric_limits<long>::max()
);
scalar x;
do
{
long i = distribution(generator);
x = reinterpret_cast<scalar&>(i);
}
while (min > mag(x) || mag(x) > max || !std::isfinite(x));
return x;
};
template <class Type>
void test(const Type& polynomialEqn, const scalar tol)
{
Roots<Type::nComponents - 1> r = polynomialEqn.roots();
const scalar nan = std::numeric_limits<scalar>::quiet_NaN();
const scalar inf = std::numeric_limits<scalar>::infinity();
FixedList<label, Type::nComponents - 1> t;
FixedList<scalar, Type::nComponents - 1> v(nan);
FixedList<scalar, Type::nComponents - 1> e(nan);
bool ok = true;
forAll(r, i)
{
t[i] = r.type(i);
switch (t[i])
{
case roots::real:
v[i] = polynomialEqn.value(r[i]);
e[i] = polynomialEqn.error(r[i]);
ok = ok && mag(v[i]) <= tol*mag(e[i]);
break;
case roots::posInf:
v[i] = + inf;
e[i] = nan;
break;
case roots::negInf:
v[i] = - inf;
e[i] = nan;
break;
default:
v[i] = e[i] = nan;
break;
}
}
if (!ok)
{
Info<< "Coeffs: " << polynomialEqn << endl
<< " Types: " << t << endl
<< " Roots: " << r << endl
<< "Values: " << v << endl
<< "Errors: " << e << endl << endl;
}
}
int main()
{
const scalar tol = 5;
const label nTests = 1000000;
for (label t = 0; t < nTests; ++ t)
{
test
(
cubicEqn
(
randomScalar(1e-50, 1e+50),
randomScalar(1e-50, 1e+50),
randomScalar(1e-50, 1e+50),
randomScalar(1e-50, 1e+50)
),
tol
);
}
Info << nTests << " random 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;
}

0
applications/test/cubicEqn/Make/files → applications/test/polynomialEqns/cubicEqn/Make/files
View File

0
applications/test/cubicEqn/Make/options → applications/test/polynomialEqns/cubicEqn/Make/options
View File

339
applications/test/polynomialEqns/cubicEqn/Test-cubicEqn.C Normal file
View File

@ -0,0 +1,339 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2018 OpenFOAM Foundation
Copyright (C) 2020 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 <http://www.gnu.org/licenses/>.
Application
Test-cubicEqn
Description
Tests for \c cubicEqn constructors, and member functions.
\*---------------------------------------------------------------------------*/
#include <ctime>
#include <random>
#include "cubicEqn.H"
#include "vector.H"
#include "scalarList.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Total number of unit tests
unsigned nTest_ = 0;
// Total number of failed unit tests
unsigned nFail_ = 0;
// Compare two floating point types, and print output.
// Do ++nFail_ if values of two objects are not equal within a given tolerance.
// The function is converted from PEP-485.
void cmp
(
const word& msg,
const scalar& x,
const scalar& y,
const scalar relTol = 1e-8, //<! are values the same within 8 decimals
const scalar absTol = 0 //<! useful for cmps near zero
)
{
Info<< msg << x << endl;
unsigned nFail = 0;
if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
{
++nFail;
}
if (nFail)
{
Info<< nl
<< " #### Fail in " << nFail << " comps ####" << nl << endl;
++nFail_;
}
++nTest_;
}
// Create each constructor of cubicEqn, and print output
void test_constructors()
{
#if 0
{
Info<< "# Construct null:" << nl;
const cubicEqn T();
Info<< T << endl;
}
#endif
{
Info<< "# Construct initialized to zero:" << nl;
const cubicEqn T(Zero);
Info<< T << endl;
}
{
Info<< "# Construct from components:" << nl;
const cubicEqn T(1, 2, 3, 4);
Info<< T << endl;
}
}
// Execute each member function of cubicEqn, and print output
void test_member_funcs()
{
const cubicEqn T(1, -6, -2, 12);
const scalar x = 9.7;
Info<< "# Operands: " << nl
<< " cubicEqn = " << T << nl
<< " x = " << x << endl;
{
Info<< "# Access:" << nl;
cmp(" a = ", 1, T.a());
cmp(" b = ", -6, T.b());
cmp(" c = ", -2, T.c());
cmp(" d = ", 12, T.d());
}
{
Info<< "# Evaluate:" << nl;
cmp(" T.value(x) = ", T.value(x), 340.73299999999983);
cmp(" T.derivative(x) = ", T.derivative(x), 163.87);
cmp
(
" T.error(x) = ",
T.error(x),
2.1854789999999994e-13,
1e-8,
1e-8
);
const vector rts(6, 1.41421356, -1.41421356);
forAll(rts, i)
{
cmp(" T.roots() = ", T.roots()[i], rts[i]);
}
}
{
Info<< "# Special cases for the root evaluation: " << nl;
{
Info<< " Three distinct real roots:" << nl;
const cubicEqn Tc(1.0, -4.0, 1.0, 6.0);
const vector rts(3, 2, -1);
forAll(rts, i)
{
cmp(" root = ", Tc.roots()[i], rts[i]);
}
}
{
Info<< " One real root and one complex conjugate-pair root:" << nl;
const cubicEqn Tc(1.0, 2.0, 3.0, 4.0);
const vector rts(-1.65062919, -0.1746854, 1.54686889);
forAll(rts, i)
{
cmp(" root = ", Tc.roots()[i], rts[i], 1e-8, 1e-8);
}
}
{
Info<< " Two identical real roots and a distinct real root:" << nl;
const cubicEqn Tc(1.0, -5.0, 8.0, -4.0);
const vector rts(2, 2, 1);
forAll(rts, i)
{
cmp(" root = ", Tc.roots()[i], rts[i]);
}
}
{
Info<< " Three identical real roots:" << nl;
const cubicEqn Tc(1.0, -6.0, 12.0, -8.0);
const vector rts(2, 2, 2);
forAll(rts, i)
{
cmp(" root = ", Tc.roots()[i], rts[i]);
}
}
{
Info<< " Arbitrary roots with non-unity leading term:" << nl;
const cubicEqn Tc(-8.59, 6.77, -0.2, 8.0);
const vector rts(1.31167947, -0.26177687, 0.80093099);
forAll(rts, i)
{
cmp(" root = ", Tc.roots()[i], rts[i]);
}
}
}
}
scalar randomScalar(const scalar min, const scalar max)
{
static_assert
(
sizeof(long) == sizeof(scalar),
"Scalar and long are not the same size"
);
static std::default_random_engine generator(std::time(0));
static std::uniform_int_distribution<long>
distribution
(
std::numeric_limits<long>::min(),
std::numeric_limits<long>::max()
);
scalar x = 0;
do
{
x = scalar(distribution(generator));
}
while (min > mag(x) || mag(x) > max || !std::isfinite(x));
return x;
};
template <class Type>
void test(const Type& polynomialEqn, const scalar tol)
{
Roots<Type::nComponents - 1> r = polynomialEqn.roots();
const scalar nan = std::numeric_limits<scalar>::quiet_NaN();
const scalar inf = std::numeric_limits<scalar>::infinity();
FixedList<label, Type::nComponents - 1> t;
FixedList<scalar, Type::nComponents - 1> v(nan);
FixedList<scalar, Type::nComponents - 1> e(nan);
bool ok = true;
forAll(r, i)
{
t[i] = r.type(i);
switch (t[i])
{
case roots::real:
v[i] = polynomialEqn.value(r[i]);
e[i] = polynomialEqn.error(r[i]);
ok = ok && mag(v[i]) <= tol*mag(e[i]);
break;
case roots::posInf:
v[i] = + inf;
e[i] = nan;
break;
case roots::negInf:
v[i] = - inf;
e[i] = nan;
break;
default:
v[i] = e[i] = nan;
break;
}
}
if (!ok)
{
Info<< "Coeffs: " << polynomialEqn << endl
<< " Types: " << t << endl
<< " Roots: " << r << endl
<< "Values: " << v << endl
<< "Errors: " << e << endl << endl;
++nFail_;
}
}
// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * //
int main()
{
Info<< nl << " ## Test constructors: ##" << nl;
test_constructors();
Info<< nl << " ## Test member functions: ##" << nl;
test_member_funcs();
// Pre-v2006 tests
const scalar tol = 5;
const label nTests = 1000000;
for (label t = 0; t < nTests; ++ t)
{
test
(
cubicEqn
(
randomScalar(1e-50, 1e+50),
randomScalar(1e-50, 1e+50),
randomScalar(1e-50, 1e+50),
randomScalar(1e-50, 1e+50)
),
tol
);
++nTest_;
}
Info << nTests << " scalar cubic equations were 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);
++nTest_;
}
}
}
}
Info<< nCoeff*nCoeff*nCoeff*nCoeff
<< " label cubic equations were tested." << endl;
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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

3
applications/test/polynomialEqns/linearEqn/Make/files Normal file
View File

@ -0,0 +1,3 @@
Test-linearEqn.C
EXE = $(FOAM_USER_APPBIN)/Test-linearEqn

3
applications/test/polynomialEqns/linearEqn/Make/options Normal file
View File

@ -0,0 +1,3 @@
EXE_INC =
EXE_LIBS =

270
applications/test/polynomialEqns/linearEqn/Test-linearEqn.C Normal file
View File

@ -0,0 +1,270 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 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 <http://www.gnu.org/licenses/>.
Application
Test-linearEqn
Description
Tests for \c linearEqn constructors, and member functions.
\*---------------------------------------------------------------------------*/
#include <ctime>
#include <random>
#include "linearEqn.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Total number of unit tests
unsigned nTest_ = 0;
// Total number of failed unit tests
unsigned nFail_ = 0;
// Compare two floating point types, and print output.
// Do ++nFail_ if values of two objects are not equal within a given tolerance.
// The function is converted from PEP-485.
void cmp
(
const word& msg,
const scalar& x,
const scalar& y,
const scalar relTol = 1e-8, //<! are values the same within 8 decimals
const scalar absTol = 0 //<! useful for cmps near zero
)
{
Info<< msg << x << endl;
unsigned nFail = 0;
if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
{
++nFail;
}
if (nFail)
{
Info<< nl
<< " #### Fail in " << nFail << " comps ####" << nl << endl;
++nFail_;
}
++nTest_;
}
// Create each constructor of linearEqn, and print output
void test_constructors()
{
#if 0
{
Info<< "# Construct null:" << nl;
const linearEqn T();
Info<< T << endl;
}
#endif
{
Info<< "# Construct initialized to zero:" << nl;
const linearEqn T(Zero);
Info<< T << endl;
}
{
Info<< "# Construct from components:" << nl;
const linearEqn T(1, 2);
Info<< T << endl;
}
}
// Execute each member function of linearEqn, and print output
void test_member_funcs()
{
const linearEqn T(1, -6);
const scalar x = 9.7;
Info<< "# Operands: " << nl
<< " linearEqn = " << T << nl
<< " x = " << x << endl;
{
Info<< "# Access:" << nl;
cmp(" a = ", 1, T.a());
cmp(" b = ", -6, T.b());
}
{
Info<< "# Evaluate:" << nl;
cmp(" T.value(x) = ", T.value(x), 3.6999999999999993);
cmp(" T.derivative(x) = ", T.derivative(x), 1);
cmp
(
" T.error(x) = ",
T.error(x),
1.5699999999999999e-15,
1e-8,
1e-8
);
cmp(" T.roots() = ", T.roots()[0], 6);
}
}
scalar randomScalar(const scalar min, const scalar max)
{
static_assert
(
sizeof(long) == sizeof(scalar),
"Scalar and long are not the same size"
);
static std::default_random_engine generator(std::time(0));
static std::uniform_int_distribution<long>
distribution
(
std::numeric_limits<long>::min(),
std::numeric_limits<long>::max()
);
scalar x = 0;
do
{
x = scalar(distribution(generator));
}
while (min > mag(x) || mag(x) > max || !std::isfinite(x));
return x;
};
template <class Type>
void test(const Type& polynomialEqn, const scalar tol)
{
Roots<Type::nComponents - 1> r = polynomialEqn.roots();
const scalar nan = std::numeric_limits<scalar>::quiet_NaN();
const scalar inf = std::numeric_limits<scalar>::infinity();
FixedList<label, Type::nComponents - 1> t;
FixedList<scalar, Type::nComponents - 1> v(nan);
FixedList<scalar, Type::nComponents - 1> e(nan);
bool ok = true;
forAll(r, i)
{
t[i] = r.type(i);
switch (t[i])
{
case roots::real:
v[i] = polynomialEqn.value(r[i]);
e[i] = polynomialEqn.error(r[i]);
ok = ok && mag(v[i]) <= tol*mag(e[i]);
break;
case roots::posInf:
v[i] = + inf;
e[i] = nan;
break;
case roots::negInf:
v[i] = - inf;
e[i] = nan;
break;
default:
v[i] = e[i] = nan;
break;
}
}
if (!ok)
{
Info<< "Coeffs: " << polynomialEqn << endl
<< " Types: " << t << endl
<< " Roots: " << r << endl
<< "Values: " << v << endl
<< "Errors: " << e << endl << endl;
++nFail_;
}
}
// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * //
int main()
{
Info<< nl << " ## Test constructors: ##" << nl;
test_constructors();
Info<< nl << " ## Test member functions: ##" << nl;
test_member_funcs();
// Sanity checks
const scalar tol = 5;
const label nTests = 1000000;
for (label t = 0; t < nTests; ++t)
{
test
(
linearEqn
(
randomScalar(1e-50, 1e+50),
randomScalar(1e-50, 1e+50)
),
tol
);
++nTest_;
}
Info << nTests << " scalar linear equations were 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)
{
test(linearEqn(a, b), tol);
++nTest_;
}
}
Info<< nCoeff*nCoeff << " label linear equations were tested." << endl;
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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

3
applications/test/polynomialEqns/quadraticEqn/Make/files Normal file
View File

@ -0,0 +1,3 @@
Test-quadraticEqn.C
EXE = $(FOAM_USER_APPBIN)/Test-quadraticEqn

3
applications/test/polynomialEqns/quadraticEqn/Make/options Normal file
View File

@ -0,0 +1,3 @@
EXE_INC =
EXE_LIBS =

321
applications/test/polynomialEqns/quadraticEqn/Test-quadraticEqn.C Normal file
View File

@ -0,0 +1,321 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 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 <http://www.gnu.org/licenses/>.
Application
Test-quadraticEqn
Description
Tests for \c quadraticEqn constructors, and member functions.
\*---------------------------------------------------------------------------*/
#include <ctime>
#include <random>
#include "quadraticEqn.H"
#include "vector2D.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Total number of unit tests
unsigned nTest_ = 0;
// Total number of failed unit tests
unsigned nFail_ = 0;
// Compare two floating point types, and print output.
// Do ++nFail_ if values of two objects are not equal within a given tolerance.
// The function is converted from PEP-485.
void cmp
(
const word& msg,
const scalar& x,
const scalar& y,
const scalar relTol = 1e-8, //<! are values the same within 8 decimals
const scalar absTol = 0 //<! useful for cmps near zero
)
{
Info<< msg << x << endl;
unsigned nFail = 0;
if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
{
++nFail;
}
if (nFail)
{
Info<< nl
<< " #### Fail in " << nFail << " comps ####" << nl << endl;
++nFail_;
}
++nTest_;
}
// Create each constructor of quadraticEqn, and print output
void test_constructors()
{
#if 0
{
Info<< "# Construct null:" << nl;
const quadraticEqn T();
Info<< T << endl;
}
#endif
{
Info<< "# Construct initialized to zero:" << nl;
const quadraticEqn T(Zero);
Info<< T << endl;
}
{
Info<< "# Construct from components:" << nl;
const quadraticEqn T(1, 2, 3);
Info<< T << endl;
}
}
// Execute each member function of quadraticEqn, and print output
void test_member_funcs()
{
const quadraticEqn T(1, -6, -2);
const scalar x = 9.7;
Info<< "# Operands: " << nl
<< " quadraticEqn = " << T << nl
<< " x = " << x << endl;
{
Info<< "# Access:" << nl;
cmp(" a = ", 1, T.a());
cmp(" b = ", -6, T.b());
cmp(" c = ", -2, T.c());
}
{
Info<< "# Evaluate:" << nl;
cmp(" T.value(x) = ", T.value(x), 33.88999999999999);
cmp(" T.derivative(x) = ", T.derivative(x), 13.399999999999999);
cmp
(
" T.error(x) = ",
T.error(x),
1.9017999999999998e-14,
1e-8,
1e-8
);
const vector2D rts(6.31662479, -0.31662479);
forAll(rts, i)
{
cmp(" T.roots() = ", T.roots()[i], rts[i]);
}
}
{
Info<< "# Special cases for the root evaluation: " << nl;
{
Info<< " Two distinct real roots:" << nl;
const quadraticEqn Tc(1.0, -11.0, 24.0);
const vector2D rts(8, 3);
forAll(rts, i)
{
cmp(" root = ", Tc.roots()[i], rts[i]);
}
}
{
Info<< " Two identical real roots:" << nl;
const quadraticEqn Tc(1.0, -8.0, 16.0);
const vector2D rts(4, 4);
forAll(rts, i)
{
cmp(" root = ", Tc.roots()[i], rts[i]);
}
}
{
Info<< " One complex conjugate-pair root:" << nl;
const quadraticEqn Tc(2.0, -2.0, 1.0);
const vector2D rts(0.5, -0.5);
forAll(rts, i)
{
cmp(" root = ", Tc.roots()[i], rts[i], 1e-8, 1e-8);
}
}
{
Info<< " Cancellation problem:" << nl;
const quadraticEqn Tc(1.0, 68.5, 0.1);
const vector2D rts(-6.84985401e+01, -1.45988513e-03);
forAll(rts, i)
{
cmp(" root = ", Tc.roots()[i], rts[i], 1e-8, 1e-8);
}
}
}
}
scalar randomScalar(const scalar min, const scalar max)
{
static_assert
(
sizeof(long) == sizeof(scalar),
"Scalar and long are not the same size"
);
static std::default_random_engine generator(std::time(0));
static std::uniform_int_distribution<long>
distribution
(
std::numeric_limits<long>::min(),
std::numeric_limits<long>::max()
);
scalar x = 0;
do
{
x = scalar(distribution(generator));
}
while (min > mag(x) || mag(x) > max || !std::isfinite(x));
return x;
};
template <class Type>
void test(const Type& polynomialEqn, const scalar tol)
{
Roots<Type::nComponents - 1> r = polynomialEqn.roots();
const scalar nan = std::numeric_limits<scalar>::quiet_NaN();
const scalar inf = std::numeric_limits<scalar>::infinity();
FixedList<label, Type::nComponents - 1> t;
FixedList<scalar, Type::nComponents - 1> v(nan);
FixedList<scalar, Type::nComponents - 1> e(nan);
bool ok = true;
forAll(r, i)
{
t[i] = r.type(i);
switch (t[i])
{
case roots::real:
v[i] = polynomialEqn.value(r[i]);
e[i] = polynomialEqn.error(r[i]);
ok = ok && mag(v[i]) <= tol*mag(e[i]);
break;
case roots::posInf:
v[i] = + inf;
e[i] = nan;
break;
case roots::negInf:
v[i] = - inf;
e[i] = nan;
break;
default:
v[i] = e[i] = nan;
break;
}
}
if (!ok)
{
Info<< "Coeffs: " << polynomialEqn << endl
<< " Types: " << t << endl
<< " Roots: " << r << endl
<< "Values: " << v << endl
<< "Errors: " << e << endl << endl;
++nFail_;
}
}
// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * //
int main()
{
Info<< nl << " ## Test constructors: ##" << nl;
test_constructors();
Info<< nl << " ## Test member functions: ##" << nl;
test_member_funcs();
// Sanity checks
const scalar tol = 5;
const label nTests = 1000000;
for (label t = 0; t < nTests; ++t)
{
test
(
quadraticEqn
(
randomScalar(1e-50, 1e+50),
randomScalar(1e-50, 1e+50),
randomScalar(1e-50, 1e+50)
),
tol
);
++nTest_;
}
Info << nTests << " scalar quadratic equations were 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)
{
test(quadraticEqn(a, b, c), tol);
++nTest_;
}
}
}
Info<< nCoeff*nCoeff*nCoeff
<< " label quadratic equations were tested." << endl;
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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -33,132 +34,104 @@ License
Foam::Roots<3> Foam::cubicEqn::roots() const
{
/*
This function solves a cubic equation of the following form:
a*x^3 + b*x^2 + c*x + d = 0
x^3 + B*x^2 + C*x + D = 0
The following two substitutions are used:
x = t - B/3
t = w - P/3/w
This reduces the problem to a quadratic in w^3.
w^6 + Q*w^3 - P = 0
Where Q and P are given in the code below.
The properties of the cubic can be related to the properties of this
quadratic in w^3. If it has a repeated root a zero, the cubic has a tripl
root. If it has a repeated root not at zero, the cubic has two real roots,
one repeated and one not. If it has two complex roots, the cubic has three
real roots. If it has two real roots, then the cubic has one real root and
two complex roots.
This is solved for the most numerically accurate value of w^3. See the
quadratic function for details on how to pick a value. This single value of
w^3 can yield up to three cube roots for w, which relate to the three
solutions for x.
Only a single root, or pair of conjugate roots, is directly evaluated; the
one, or ones with the lowest relative numerical error. Root identities are
then used to recover the remaining roots, possibly utilising a quadratic
and/or linear solution. This seems to be a good way of maintaining the
accuracy of roots at very different magnitudes.
*/
const scalar a = this->a();
const scalar b = this->b();
const scalar c = this->c();
const scalar d = this->d();
if (a == 0)
// Check the leading term in the cubic eqn exists
if (mag(a) < VSMALL)
{
return Roots<3>(quadraticEqn(b, c, d).roots(), roots::nan, 0);
}
// This is assumed not to over- or under-flow. If it does, all bets are off.
const scalar p = c*a - b*b/3;
// (JLM:p. 2246) [p = a*c - b*b/3]
const scalar w = a*c;
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 disc = p*p*p/27 + q*q/4;
const scalar numDiscr = p*p*p/27 + q*q/4;
const scalar discr = (mag(numDiscr) > SMALL) ? numDiscr : 0;
// How many roots of what types are available?
const bool oneReal = disc == 0 && p == 0;
const bool twoReal = disc == 0 && p != 0;
const bool threeReal = disc < 0;
//const bool oneRealTwoComplex = disc > 0;
// Determine the number and types of the roots
const bool threeReal = discr < 0;
const bool oneRealTwoComplex = discr > 0;
const bool twoReal = //p != 0; && discr == 0;
(mag(p) > VSMALL) && !(threeReal || oneRealTwoComplex);
// const bool oneReal = p == 0 && discr == 0;
static const scalar sqrt3 = sqrt(3.0);
scalar x;
scalar x = 0;
if (oneReal)
if (threeReal)
{
const Roots<1> r = linearEqn(a, b/3).roots();
return Roots<3>(r.type(0), r[0]);
}
else if (twoReal)
{
if (q*b > 0)
{
x = - 2*cbrt(q/2) - b/3;
}
else
{
x = cbrt(q/2) - b/3;
const Roots<1> r = linearEqn(- a, x).roots();
return Roots<3>(Roots<2>(r, r), linearEqn(x*x, a*d).roots());
}
}
else if (threeReal)
{
const scalar wCbRe = - q/2, wCbIm = sqrt(- disc);
const scalar wCbRe = -q/2;
const scalar wCbIm = sqrt(-discr);
const scalar wAbs = cbrt(hypot(wCbRe, wCbIm));
const scalar wArg = atan2(wCbIm, wCbRe)/3;
const scalar wRe = wAbs*cos(wArg), wIm = wAbs*sin(wArg);
const scalar wRe = wAbs*cos(wArg);
const scalar wIm = wAbs*sin(wArg);
if (b > 0)
{
x = - wRe - mag(wIm)*sqrt3 - b/3;
x = -wRe - mag(wIm)*sqrt3 - b/3;
}
else
{
x = 2*wRe - b/3;
}
}
else // if (oneRealTwoComplex)
else if (oneRealTwoComplex)
{
const scalar wCb = - q/2 - sign(q)*sqrt(disc);
const scalar wCb = -q/2 - sign(q)*sqrt(discr);
const scalar w = cbrt(wCb);
const scalar t = w - p/(3*w);
if (p + t*b < 0)
{
x = t - b/3;
}
else
{
const scalar xRe = - t/2 - b/3, xIm = sqrt3/2*(w + p/3/w);
x = - a*a*d/(xRe*xRe + xIm*xIm);
const scalar xRe = -t/2 - b/3;
const scalar xIm = sqrt3/2*(w + p/3/w);
// This form of solving for the remaining roots seems more stable
// for this case. This has been determined by trial and error.
return
Roots<3>
(
linearEqn(- a, x).roots(),
quadraticEqn(a*x, x*x + b*x, - a*d).roots()
Roots<1>(roots::real, -a*d/(xRe*xRe + xIm*xIm)),
Roots<2>
(
Roots<1>(roots::complex, xRe),
Roots<1>(roots::complex, xIm)
)
);
}
}
else if (twoReal)
{
if (q*b > 0)
{
x = -2*cbrt(q/2) - b/3;
}
else
{
x = cbrt(q/2) - b/3;
const Roots<1> r(linearEqn(-a, x).roots());
return Roots<3>(Roots<2>(r, r), linearEqn(x*x, a*d).roots());
}
}
else // (oneReal)
{
const Roots<1> r(linearEqn(a, b/3).roots());
return Roots<3>(r.type(0), r[0]);
}
return
Roots<3>
(
linearEqn(- a, x).roots(),
quadraticEqn(- x*x, c*x + a*d, d*x).roots()
linearEqn(-a, x).roots(),
quadraticEqn(-x*x, c*x + a*d, d*x).roots()
);
}

98
src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -27,7 +28,67 @@ Class
Foam::cubicEqn
Description
Cubic equation of the form a*x^3 + b*x^2 + c*x + d = 0
Container to encapsulate various operations for
cubic equation of the forms with real coefficients:
\f[
a*x^3 + b*x^2 + c*x + d = 0
x^3 + B*x^2 + C*x + D = 0
\f]
The following two substitutions into the above forms are used:
\f[
x = t - B/3
t = w - P/3/w
\f]
This reduces the problem to a quadratic in \c w^3:
\f[
w^6 + Q*w^3 - P = 0
\f]
where \c Q and \c P are given in the code.
The properties of the cubic can be related to the properties of this
quadratic in \c w^3.
If the quadratic eqn has two identical real roots at zero, three identical
real roots exist in the cubic eqn.
If the quadratic eqn has two identical real roots at non-zero, two identical
and one distinct real roots exist in the cubic eqn.
If the quadratic eqn has two complex roots (a complex conjugate pair),
three distinct real roots exist in the cubic eqn.
If the quadratic eqn has two distinct real roots, one real root and two
complex roots (a complex conjugate pair) exist in the cubic eqn.
The quadratic eqn is solved for the most numerically accurate value
of \c w^3. See the \link quadraticEqn.H \endlink for details on how to
pick a value. This single value of \c w^3 can yield up to three cube
roots for \c w, which relate to the three solutions for \c x.
Only a single root, or pair of conjugate roots, is directly evaluated; the
one, or ones with the lowest relative numerical error. Root identities are
then used to recover the remaining roots, possibly utilising a quadratic
and/or linear solution. This seems to be a good way of maintaining the
accuracy of roots at very different magnitudes.
Reference:
\verbatim
Kahan's algo. to compute 'p' using fused multiply-adds (tag:JML):
Jeannerod, C. P., Louvet, N., & Muller, J. M. (2013).
Further analysis of Kahan's algorithm for the accurate
computation of 2× 2 determinants.
Mathematics of Computation, 82(284), 2245-2264.
DOI:10.1090/S0025-5718-2013-02679-8
\endverbatim
See also
Test-cubicEqn.C
SourceFiles
cubicEqnI.H
@ -79,28 +140,37 @@ public:
// Member Functions
// Access
// Access
inline scalar a() const;
inline scalar b() const;
inline scalar c() const;
inline scalar d() const;
inline scalar a() const;
inline scalar b() const;
inline scalar c() const;
inline scalar d() const;
inline scalar& a();
inline scalar& b();
inline scalar& c();
inline scalar& d();
inline scalar& a();
inline scalar& b();
inline scalar& c();
inline scalar& d();
//- Evaluate at x
// Evaluate
//- Evaluate the cubic equation at x
inline scalar value(const scalar x) const;
//- Evaluate the derivative at x
//- Evaluate the derivative of the cubic equation at x
inline scalar derivative(const scalar x) const;
//- Estimate the error of evaluation at x
//- Estimate the error of evaluation of the cubic equation at x
inline scalar error(const scalar x) const;
//- Get the roots
//- Return the roots of the cubic equation with no particular order
// if discriminant > 0: return three distinct real roots
// if discriminant < 0: return one real root and one complex root
// (one member of the complex conjugate pair)
// if discriminant = 0: return two identical roots,
// and one distinct real root
// if identical zero Hessian: return three identical real roots
// where the discriminant = - 4p^3 - 27q^2.
Roots<3> roots() const;
};

32
src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -27,7 +28,16 @@ Class
Foam::linearEqn
Description
Linear equation of the form a*x + b = 0
Container to encapsulate various operations for
linear equation of the forms with real coefficients:
\f[
a*x + b = 0
x + B = 0
\f]
See also
Test-linearEqn.C
SourceFiles
linearEqnI.H
@ -72,24 +82,26 @@ public:
// Member Functions
// Access
// Access
inline scalar a() const;
inline scalar b() const;
inline scalar a() const;
inline scalar b() const;
inline scalar& a();
inline scalar& b();
inline scalar& a();
inline scalar& b();
//- Evaluate at x
// Evaluate
//- Evaluate the linear equation at x
inline scalar value(const scalar x) const;
//- Evaluate the derivative at x
//- Evaluate the derivative of the linear equation at x
inline scalar derivative(const scalar x) const;
//- Estimate the error of evaluation at x
//- Estimate the error of evaluation of the linear equation at x
inline scalar error(const scalar x) const;
//- Get the roots
//- Return the real root of the linear equation
inline Roots<1> roots() const;
};

17
src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -90,29 +91,19 @@ inline Foam::scalar Foam::linearEqn::error(const scalar x) const
inline Foam::Roots<1> Foam::linearEqn::roots() const
{
/*
This function solves a linear equation of the following form:
a*x + b = 0
x + B = 0
*/
const scalar a = this->a();
const scalar b = this->b();
if (a == 0)
if (mag(a) < VSMALL)
{
return Roots<1>(roots::nan, 0);
}
if (mag(b/VGREAT) >= mag(a))
else if (mag(b/VGREAT) >= mag(a))
{
return Roots<1>(sign(a) == sign(b) ? roots::negInf : roots::posInf, 0);
}
return Roots<1>(roots::real, - b/a);
return Roots<1>(roots::real, -b/a);
}

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

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -32,59 +33,43 @@ License
Foam::Roots<2> Foam::quadraticEqn::roots() const
{
/*
This function solves a quadraticEqn equation of the following form:
a*x^2 + b*x + c = 0
x^2 + B*x + C = 0
The quadraticEqn formula is as follows:
x = - B/2 +- sqrt(B*B - 4*C)/2
If the sqrt generates a complex number, this provides the result. If not
then the real root with the smallest floating point error is calculated.
x0 = - B/2 - sign(B)*sqrt(B*B - 4*C)/2
The other root is the obtained using an identity.
x1 = C/x0
*/
const scalar a = this->a();
const scalar b = this->b();
const scalar c = this->c();
if (a == 0)
// Check the leading term in the quadratic eqn exists
if (mag(a) < VSMALL)
{
return Roots<2>(linearEqn(b, c).roots(), roots::nan, 0);
}
// This is assumed not to over- or under-flow. If it does, all bets are off.
const scalar disc = b*b/4 - a*c;
// (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;
// How many roots of what types are available?
const bool oneReal = disc == 0;
const bool twoReal = disc > 0;
//const bool twoComplex = disc < 0;
// Find how many roots of what types are available
const bool twoReal = discr > 0;
const bool twoComplex = discr < 0;
//const bool oneReal = discr == 0;
if (oneReal)
if (twoReal)
{
const Roots<1> r = linearEqn(a, b/2).roots();
// (F:Exp. 8.9)
const scalar x = -b/2 - sign(b)*sqrt(discr);
return Roots<2>(linearEqn(-a, x).roots(), linearEqn(-x, c).roots());
}
else if (twoComplex)
{
const Roots<1> xRe(roots::type::complex, -b/2/a);
const Roots<1> xIm(roots::type::complex, sign(b)*sqrt(mag(discr))/a);
return Roots<2>(xRe, xIm);
}
else // (oneReal)
{
const Roots<1> r(linearEqn(a, b/2).roots());
return Roots<2>(r, r);
}
else if (twoReal)
{
const scalar x = - b/2 - sign(b)*sqrt(disc);
return Roots<2>(linearEqn(- a, x).roots(), linearEqn(- x, c).roots());
}
else // if (twoComplex)
{
return Roots<2>(roots::complex, 0);
}
}
// ************************************************************************* //

66
src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -27,7 +28,43 @@ Class
Foam::quadraticEqn
Description
Quadratic equation of the form a*x^2 + b*x + c = 0
Container to encapsulate various operations for
quadratic equation of the forms with real coefficients:
\f[
a*x^2 + b*x + c = 0
x^2 + B*x + C = 0
\f]
The expressions for the roots of \c quadraticEqn are as follows:
\f[
x1 = - (b + sign(b) sqrt(b^2 - 4ac)/(2*a))
x2 = c/(a*x1)
\f]
where \c (b^2 - 4ac) is evaluated by fused multiply-adds to avoid
detrimental cancellation.
Reference:
\verbatim
Cancellation-avoiding quadratic formula (tag:F):
Ford, W. (2014).
Numerical linear algebra with applications: Using MATLAB.
London: Elsevier/Academic Press.
DOI:10.1016/C2011-0-07533-6
Kahan's algo. to compute 'b^2-a*c' using fused multiply-adds (tag:JML):
Jeannerod, C. P., Louvet, N., & Muller, J. M. (2013).
Further analysis of Kahan's algorithm for the accurate
computation of 2× 2 determinants.
Mathematics of Computation, 82(284), 2245-2264.
DOI:10.1090/S0025-5718-2013-02679-8
\endverbatim
See also
Test-quadraticEqn.C
SourceFiles
quadraticEqnI.H
@ -73,26 +110,31 @@ public:
// Member Functions
// Access
// Access
inline scalar a() const;
inline scalar b() const;
inline scalar c() const;
inline scalar a() const;
inline scalar b() const;
inline scalar c() const;
inline scalar& a();
inline scalar& b();
inline scalar& c();
inline scalar& a();
inline scalar& b();
inline scalar& c();
//- Evaluate at x
// Evaluate
//- Evaluate the quadratic equation at x
inline scalar value(const scalar x) const;
//- Evaluate the derivative at x
//- Evaluate the derivative of the quadratic equation at x
inline scalar derivative(const scalar x) const;
//- Estimate the error of evaluation at x
//- Estimate the error of evaluation of the quadratic equation at x
inline scalar error(const scalar x) const;
//- Get the roots
//- Return the roots of the quadratic equation with no particular order
// if discriminant > 0: return two distinct real roots
// if discriminant < 0: return one of the complex conjugate-pair roots
// otherwise : return two identical real roots
Roots<2> roots() const;
};