Commit 66afe343 authored by Michael Turner's avatar Michael Turner

Merge remote-tracking branch 'upstream/master' into fix/missing-2d-configs

parents f5033c84 03d884a1
......@@ -23,6 +23,7 @@ v5.0.0
- Refactor ParseUtils to be more consistent (!843)
- Added support for using the distance to a specific region (e.g. outlet) in the
function definitions for the Absorption Forcing (!769)
- Improve performance of DisContField2D::v_ExtractTracePhys (!824)
**NekMesh**:
- Add feature to read basic 2D geo files as CAD (!731)
......@@ -49,9 +50,18 @@ v5.0.0
- Enable output to multiple files (!844)
- Allow using xml file without expansion tag in FieldConvert (!849)
**CompressibleFlowSolver**
- Add 3D regression tests (!567)
**Documentation**:
- Added the developer-guide repository as a submodule (!751)
v4.4.2
------
**Library**
- Fix evaluation of points (e.g. HistoryPoints, Interpolation to pts) close to
the interface of two elements (!836)
v4.4.1
------
**Library**
......@@ -69,7 +79,7 @@ v4.4.1
- Fix deadlock with HDF5 input (!786)
- Fix missing entriess in LibUtilities::kPointsTypeStr (!792)
- Fix compiler warnings with CommDataType (!793)
- Fix ability to set default implementation in Collections and added an option
- Fix ability to set default implementation in Collections and added an option
to set eNoCollections in FieldConvert as default (!789)
- Fix performance issue in ProcessIsoContour in relation to memory consumption
(!821)
......@@ -77,6 +87,11 @@ v4.4.1
- Fix available classes being listed multiple times (!817)
- Fix Intel compiler warnings (!837)
- Fix overwriting and backup of chk/fld files on slow file systes (!741)
- Fix DriverAdaptive with second order IMEX (!850)
- 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)
......
......@@ -21,3 +21,11 @@ ADD_NEKTAR_TEST(NodalDemo_Prism_Interp_P7)
ADD_NEKTAR_TEST(NodalDemo_Tet_Deriv_P8)
ADD_NEKTAR_TEST(NodalDemo_Tet_Integral_P6)
ADD_NEKTAR_TEST(NodalDemo_Tet_Interp_P7)
ADD_NEKTAR_TEST(TimeIntegrationDemoIMEXGear)
ADD_NEKTAR_TEST(TimeIntegrationDemoIMEXOrder1)
ADD_NEKTAR_TEST(TimeIntegrationDemoIMEXOrder2)
ADD_NEKTAR_TEST(TimeIntegrationDemoIMEXOrder3)
ADD_NEKTAR_TEST(TimeIntegrationDemoDIRKIMEXOrder2)
ADD_NEKTAR_TEST(TimeIntegrationDemoDIRKIMEXOrder3)
ADD_NEKTAR_TEST(TimeIntegrationDemoCNAB)
ADD_NEKTAR_TEST(TimeIntegrationDemoMCNAB)
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 7</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.10135</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 4</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.0616905</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 5</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.061634</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 6</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.315262</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 1</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.224678</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 2</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.0798986</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 3</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.0610636</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 8</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.134315</value>
</metric>
</metrics>
</test>
......@@ -350,7 +350,16 @@ void Interpolator::Interpolate(
// Obtain Element and LocalCoordinate to interpolate
int elmtid = m_expInField[0]->GetExpIndex(Scoords, Lcoords,
NekConstants::kNekZeroTol);
NekConstants::kGeomFactorsTol);
// we use kGeomFactorsTol as tolerance, while StdPhysEvaluate has
// kNekZeroTol hardcoded, so we need to limit Lcoords to not produce
// a ton of warnings
for(int j = 0; j < nInDim; ++j)
{
Lcoords[j] = std::max(Lcoords[j], -1.0);
Lcoords[j] = std::min(Lcoords[j], 1.0);
}
if (elmtid >= 0)
{
......@@ -431,7 +440,16 @@ void Interpolator::Interpolate(
// Obtain Element and LocalCoordinate to interpolate
int elmtid = m_expInField[0]->GetExpIndex(coords, Lcoords,
NekConstants::kNekZeroTol);
NekConstants::kGeomFactorsTol);
// we use kGeomFactorsTol as tolerance, while StdPhysEvaluate has
// kNekZeroTol hardcoded, so we need to limit Lcoords to not produce
// a ton of warnings
for(int j = 0; j < nInDim; ++j)
{
Lcoords[j] = std::max(Lcoords[j], -1.0);
Lcoords[j] = std::min(Lcoords[j], 1.0);
}
if (elmtid >= 0)
{
......
......@@ -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);
}
......
......@@ -742,7 +742,7 @@ namespace Nektar
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,twothirdth);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,twothirdth);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,twothirdth);
m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
......
......@@ -1488,28 +1488,43 @@ namespace Nektar
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
// Loop over elemente and collect forward expansion
int nexp = GetExpSize();
int n, e, offset, phys_offset;
Array<OneD,NekDouble> e_tmp;
Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
if (basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
{
Vmath::Zero(outarray.num_elements(), outarray, 1);
Array<OneD, NekDouble> tracevals(
m_locTraceToTraceMap->GetNFwdLocTracePts());
m_locTraceToTraceMap->FwdLocTracesFromField(inarray,tracevals);
m_locTraceToTraceMap->
InterpLocEdgesToTrace(0,tracevals,outarray);
m_traceMap->UniversalTraceAssemble(outarray);
}
else
{
// Loop over elemente and collect forward expansion
int nexp = GetExpSize();
int n, e, offset, phys_offset;
Array<OneD,NekDouble> e_tmp;
Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
ASSERTL1(outarray.num_elements() >= m_trace->GetNpoints(),
"input array is of insufficient length");
ASSERTL1(outarray.num_elements() >= m_trace->GetNpoints(),
"input array is of insufficient length");
// use m_trace tmp space in element to fill values
for(n = 0; n < nexp; ++n)
{
phys_offset = GetPhys_Offset(n);
for(e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
// use m_trace tmp space in element to fill values
for(n = 0; n < nexp; ++n)
{
offset = m_trace->GetPhys_Offset(
elmtToTrace[n][e]->GetElmtId());
(*m_exp)[n]->GetEdgePhysVals(e, elmtToTrace[n][e],
inarray + phys_offset,
e_tmp = outarray + offset);
phys_offset = GetPhys_Offset(n);
for(e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
{
offset = m_trace->GetPhys_Offset(
elmtToTrace[n][e]->GetElmtId());
(*m_exp)[n]->GetEdgePhysVals(e, elmtToTrace[n][e],
inarray + phys_offset,
e_tmp = outarray + offset);
}
}
}
}
......
......@@ -84,36 +84,25 @@ namespace Nektar
int nCoeffs = fields[0]->GetNcoeffs();
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
Array<OneD, NekDouble> qcoeffs(nCoeffs);
Array<OneD, Array<OneD, NekDouble> > fluxvector(nDim);
Array<OneD, NekDouble> tmp(nCoeffs);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > flux (nDim);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > qfield(nDim);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > qfieldStd(nDim);
for (j = 0; j < nDim; ++j)
{
qfield[j] =
Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
qfieldStd[j] =
Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
flux[j] =
Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
for (i = 0; i < nConvectiveFields; ++i)
{
qfield[j][i] = Array<OneD, NekDouble>(nPts, 0.0);
qfieldStd[j][i] = Array<OneD, NekDouble>(nPts, 0.0);
flux[j][i] = Array<OneD, NekDouble>(nTracePts, 0.0);
}
}
for (k = 0; k < nDim; ++k)
{
fluxvector[k] = Array<OneD, NekDouble>(nPts, 0.0);
}
// Compute q_{\eta} and q_{\xi}
// Obtain numerical fluxes
......@@ -123,12 +112,12 @@ namespace Nektar
{
for (i = 0; i < nConvectiveFields; ++i)
{
fields[i]->IProductWRTDerivBase(j, inarray[i], qcoeffs);
Vmath::Neg (nCoeffs, qcoeffs, 1);
fields[i]->AddTraceIntegral (flux[j][i], qcoeffs);
fields[i]->IProductWRTDerivBase(j, inarray[i], tmp);
Vmath::Neg (nCoeffs, tmp, 1);
fields[i]->AddTraceIntegral (flux[j][i], tmp);
fields[i]->SetPhysState (false);
fields[i]->MultiplyByElmtInvMass(qcoeffs, qcoeffs);
fields[i]->BwdTrans (qcoeffs, qfield[j][i]);
fields[i]->MultiplyByElmtInvMass(tmp, tmp);
fields[i]->BwdTrans (tmp, qfield[j][i]);
}
}
// Compute u from q_{\eta} and q_{\xi}
......@@ -198,7 +187,6 @@ namespace Nektar
}
}
Array<OneD, NekDouble> tmp = Array<OneD, NekDouble>(nCoeffs, 0.0);
Array<OneD, Array<OneD, NekDouble> > qdbase(nDim);
for (i = 0; i < nConvectiveFields; ++i)
......
......@@ -226,12 +226,18 @@ void FilterHistoryPoints::v_Initialise(
if (m_isHomogeneous1D)
{
idList[i] = pFields[0]->GetPlane(0)->GetExpIndex(gloCoord,locCoords,
NekConstants::kNekZeroTol);
NekConstants::kGeomFactorsTol);
}
else
{
idList[i] = pFields[0]->GetExpIndex(gloCoord,locCoords,
NekConstants::kNekZeroTol);
NekConstants::kGeomFactorsTol);
}
for(int j = 0; j < 3; ++j)
{
locCoords[j] = std::max(locCoords[j], -1.0);
locCoords[j] = std::min(locCoords[j], 1.0);
}
// Save Local coordinates for later
......
......@@ -248,6 +248,7 @@ namespace Nektar
LibUtilities::Timer timer;
bool doCheckTime = false;
int step = m_initialStep;
int stepCounter = 0;
NekDouble intTime = 0.0;
NekDouble lastCheckTime = 0.0;
NekDouble cpuTime = 0.0;
......@@ -283,7 +284,7 @@ namespace Nektar
}
fields = m_intScheme->TimeIntegrate(
step, m_timestep, m_intSoln, m_ode);
stepCounter, m_timestep, m_intSoln, m_ode);
timer.Stop();
m_time += m_timestep;
......@@ -395,6 +396,7 @@ namespace Nektar
// Step advance
++step;
++stepCounter;
}
// Print out summary statistics
......
......@@ -308,7 +308,7 @@ namespace Nektar
}
}
v_GetLocCoords(gloCoord, locCoord);
resid = v_GetLocCoords(gloCoord, locCoord);
if (locCoord[0] >= -(1+tol) && locCoord[0] <= 1+tol
&& locCoord[1] >= -(1+tol) && locCoord[1] <= 1+tol
......@@ -325,12 +325,12 @@ namespace Nektar
{
if(locCoord[i] <-(1+tol))
{
locCoord[i] = -(1+tol);
locCoord[i] = -1;
}
if(locCoord[i] > (1+tol))
{
locCoord[i] = 1+tol;
locCoord[i] = 1;
}
}
......
......@@ -204,12 +204,12 @@ namespace Nektar
{
if(locCoord[i] <-(1+tol))
{
locCoord[i] = -(1+tol);
locCoord[i] = -1;
}
if(locCoord[i] > (1+tol))
{
locCoord[i] = 1+tol;
locCoord[i] = 1;
}
}
......
......@@ -167,12 +167,12 @@ namespace Nektar
{
if(locCoord[i] <-(1+tol))
{
locCoord[i] = -(1+tol);
locCoord[i] = -1;
}
if(locCoord[i] > (1+tol))
{
locCoord[i] = 1+tol;
locCoord[i] = 1;
}
}
......
......@@ -92,8 +92,8 @@ namespace Nektar
NekDouble val;
DNekMatSharedPtr I = m_base[0]->GetI(Lcoord);
ASSERTL2(Lcoord[0] >= -1,"Lcoord[0] < -1");
ASSERTL2(Lcoord[0] <= 1,"Lcoord[0] > 1");
ASSERTL2(Lcoord[0] >= -1 - NekConstants::kNekZeroTol,"Lcoord[0] < -1");
ASSERTL2(Lcoord[0] <= 1 + NekConstants::kNekZeroTol,"Lcoord[0] > 1");
val = Blas::Ddot(nquad, I->GetPtr(), 1, physvals, 1);
......
......@@ -145,12 +145,12 @@ namespace Nektar
Array<OneD, NekDouble> eta = Array<OneD, NekDouble>(3);
Array<OneD, DNekMatSharedPtr> I(3);
WARNINGL2(coords[0] >= -1,"coord[0] < -1");
WARNINGL2(coords[0] <= 1,"coord[0] > 1");
WARNINGL2(coords[1] >= -1,"coord[1] < -1");
WARNINGL2(coords[1] <= 1,"coord[1] > 1");
WARNINGL2(coords[2] >= -1,"coord[2] < -1");
WARNINGL2(coords[2] <= 1,"coord[2] > 1");
WARNINGL2(coords[0] >= -1 - NekConstants::kNekZeroTol,"coord[0] < -1");
WARNINGL2(coords[0] <= 1 + NekConstants::kNekZeroTol,"coord[0] > 1");
WARNINGL2(coords[1] >= -1 - NekConstants::kNekZeroTol,"coord[1] < -1");
WARNINGL2(coords[1] <= 1 + NekConstants::kNekZeroTol,"coord[1] > 1");
WARNINGL2(coords[2] >= -1 - NekConstants::kNekZeroTol,"coord[2] < -1");
WARNINGL2(coords[2] <= 1 + NekConstants::kNekZeroTol,"coord[2] > 1");
// Obtain local collapsed corodinate from
// cartesian coordinate.
......
......@@ -97,6 +97,8 @@ IF( NEKTAR_SOLVER_COMPRESSIBLE_FLOW )
ADD_NEKTAR_TEST(Nozzle_AxiSym_NoSwirl)
ADD_NEKTAR_TEST(Nozzle_AxiSym_Swirl)
ADD_NEKTAR_TEST(Nozzle_Quasi1D_P6)
ADD_NEKTAR_TEST (hump3D_GLL)
ADD_NEKTAR_TEST (hump3D_SEM)
IF (NEKTAR_USE_MPI)
#ADD_NEKTAR_TEST(Perturbation_M05_square_CBC_par LENGTHY)
......