Commit 7f32f027 by Chris Cantwell

Merge branch 'fix/polylib' into 'master'

Fixing the bug about subtractive cancellation

See merge request !778
parents bdde0e5e 9b305890
......@@ -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)
......
......@@ -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);
}
......
......@@ -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>
......
......@@ -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>
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment