diff --git a/CHANGELOG.md b/CHANGELOG.md
index aa20ca35ebffd1347aae13ae514914e375c05c95..e87f05c3803e6f943aeb81f6e836c6df98197c1b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -79,6 +79,7 @@ v4.4.1
- Fixed typo in eIMEXGear part (!854)
- Added regression tests for IMEXOrder1, IMEXOrder2, IMEXOrder3, MCNAB,
IMEXGear, CNAB, 2nd order IMEX-DIRK, 3rd order IMEX-DIRK (!854)
+- Fix bug due to subtractive cancellation in polylib routines (!778)
**FieldConvert:**
- Fix issue with field ordering in the interppointdatatofld module (!754)
diff --git a/library/LibUtilities/Polylib/Polylib.cpp b/library/LibUtilities/Polylib/Polylib.cpp
index a08bcd6db42199738524b64b3860f358f51d1d40..cf725c7a9e526c2a26479f4d674a13a8a8ec619e 100644
--- a/library/LibUtilities/Polylib/Polylib.cpp
+++ b/library/LibUtilities/Polylib/Polylib.cpp
@@ -30,7 +30,42 @@
namespace Polylib {
+ /// The following function is used to circumvent/reduce "Subtractive Cancellation"
+ /// The expression 1/dz is replaced by optinvsub(.,.)
+ /// Added on 26 April 2017
+
+ double optdiff(double xl, double xr)
+ {
+ double m_xln, m_xrn;
+ int m_expn;
+ int m_digits = static_cast<int>(fabs(floor(log10(DBL_EPSILON)))-1);
+
+ if (fabs(xl-xr)<1.e-4){
+
+ m_expn = static_cast<int>(floor(log10(fabs(xl-xr))));
+ m_xln = xl*powl(10.0L,-m_expn)-floor(xl*powl(10.0L,-m_expn)); // substract the digits overlap part
+ m_xrn = xr*powl(10.0L,-m_expn)-floor(xl*powl(10.0L,-m_expn)); // substract the common digits overlap part
+ m_xln = round(m_xln*powl(10.0L,m_digits+m_expn)); // git rid of rubbish
+ m_xrn = round(m_xrn*powl(10.0L,m_digits+m_expn));
+
+ return powl(10.0L,-m_digits)*(m_xln-m_xrn);
+ }else{
+ return (xl-xr);
+ }
+ }
+ double laginterp(double z, int j, const double *zj, int np)
+ {
+ double temp = 1.0;
+ for (int i=0; i<np; i++)
+ {
+ if (j != i)
+ {
+ temp *=optdiff(z,zj[i])/(zj[j]-zj[i]);
+ }
+ }
+ return temp;
+ }
/// Define whether to use polynomial deflation (1) or tridiagonal solver (0).
@@ -1354,29 +1389,15 @@ namespace Polylib {
{
-
-
- double zi, dz, p, pd, h;
-
-
+ double zi, dz;
zi = *(zgj+i);
- dz = z - zi;
-
- if (fabs(dz) < EPS) return 1.0;
-
-
-
- jacobd (1, &zi, &pd , np, alpha, beta);
-
- jacobfd(1, &z , &p, NULL , np, alpha, beta);
-
- h = p/(pd*dz);
+ dz = z-zi;
+ if (fabs(dz) < EPS) return 1.0;
-
- return h;
+ return laginterp(z, i, zgj, np);
}
@@ -1431,34 +1452,15 @@ namespace Polylib {
{
-
- double zi, dz, p, pd, h;
-
-
+ double zi, dz;
zi = *(zgrj+i);
- dz = z - zi;
+ dz = z-zi;
if (fabs(dz) < EPS) return 1.0;
-
-
- jacobfd (1, &zi, &p , NULL, np-1, alpha, beta + 1);
-
- // need to use this routine in caes zi = -1 or 1
-
- jacobd (1, &zi, &pd, np-1, alpha, beta + 1);
-
- h = (1.0 + zi)*pd + p;
-
- jacobfd (1, &z, &p, NULL, np-1, alpha, beta + 1);
-
- h = (1.0 + z )*p/(h*dz);
-
-
-
- return h;
+ return laginterp(z, i, zgrj, np);
}
@@ -1515,34 +1517,15 @@ namespace Polylib {
{
-
- double zi, dz, p, pd, h;
-
-
+ double zi, dz;
zi = *(zgrj+i);
- dz = z - zi;
+ dz = z-zi;
if (fabs(dz) < EPS) return 1.0;
-
-
- jacobfd (1, &zi, &p , NULL, np-1, alpha+1, beta );
-
- // need to use this routine in caes z = -1 or 1
-
- jacobd (1, &zi, &pd, np-1, alpha+1, beta );
-
- h = (1.0 - zi)*pd - p;
-
- jacobfd (1, &z, &p, NULL, np-1, alpha+1, beta);
-
- h = (1.0 - z )*p/(h*dz);
-
-
-
- return h;
+ return laginterp(z, i, zgrj, np);
}
@@ -1598,35 +1581,15 @@ namespace Polylib {
{
- double one = 1., two = 2.;
-
- double zi, dz, p, pd, h;
-
-
+ double zi, dz;
zi = *(zglj+i);
- dz = z - zi;
+ dz = z-zi;
if (fabs(dz) < EPS) return 1.0;
-
-
-
- jacobfd(1, &zi, &p , NULL, np-2, alpha + one, beta + one);
-
- // need to use this routine in caes z = -1 or 1
-
- jacobd (1, &zi, &pd, np-2, alpha + one, beta + one);
-
- h = (one - zi*zi)*pd - two*zi*p;
-
- jacobfd(1, &z, &p, NULL, np-2, alpha + one, beta + one);
-
- h = (one - z*z)*p/(h*dz);
-
-
-
- return h;
+
+ return laginterp(z, i, zglj, np);
}
diff --git a/solvers/CompressibleFlowSolver/Tests/CylinderSubsonic_NS_WeakDG_LDG_SEM_VariableMu.tst b/solvers/CompressibleFlowSolver/Tests/CylinderSubsonic_NS_WeakDG_LDG_SEM_VariableMu.tst
index b57cbb9083c3e03bcfac215270b924e26c1c8265..beca267237b7ae7165a4f36be5893d30546d395e 100644
--- a/solvers/CompressibleFlowSolver/Tests/CylinderSubsonic_NS_WeakDG_LDG_SEM_VariableMu.tst
+++ b/solvers/CompressibleFlowSolver/Tests/CylinderSubsonic_NS_WeakDG_LDG_SEM_VariableMu.tst
@@ -16,7 +16,7 @@
<metric type="Linf" id="2">
<value variable="rho" tolerance="1e-6">0.530796</value>
<value variable="rhou" tolerance="1e-6">92.5827</value>
- <value variable="rhov" tolerance="4e-4">58.1854</value>
+ <value variable="rhov" tolerance="2e-4">58.1851</value>
<value variable="E" tolerance="1e-6">346912</value>
</metric>
</metrics>
diff --git a/solvers/IncNavierStokesSolver/Tests/Womersley_PipeFlow.tst b/solvers/IncNavierStokesSolver/Tests/Womersley_PipeFlow.tst
index 61cef7e0f7a7d146a2d44b0a4e149686793acad8..ce366a9a1d220f2b9f671bc09cbda06f157ba53b 100644
--- a/solvers/IncNavierStokesSolver/Tests/Womersley_PipeFlow.tst
+++ b/solvers/IncNavierStokesSolver/Tests/Womersley_PipeFlow.tst
@@ -12,7 +12,7 @@
<value variable="u" tolerance="1e-8">0.00192614</value>
<value variable="v" tolerance="1e-8">0.00171205</value>
<value variable="w" tolerance="1e-8">0.033844</value>
- <value variable="p" tolerance="2e-7">0.0721992</value>
+ <value variable="p" tolerance="1e-8">0.0721993</value>
</metric>
<metric type="Linf" id="2">
<value variable="u" tolerance="1e-8">0.0102797</value>