Commit 9c446afc authored by Dave Moxey's avatar Dave Moxey

Merge remote-tracking branch 'upstream/master' into feature/cmake-cleanup

Conflicts:
	library/StdRegions/StdExpansion0D.cpp
	library/StdRegions/StdPointExp.cpp
	solvers/CompressibleFlowSolver/CMakeLists.txt
parents 01fb1258 e38421d7
...@@ -437,6 +437,7 @@ Compressible flow is characterised by abrupt changes in density within the flow ...@@ -437,6 +437,7 @@ Compressible flow is characterised by abrupt changes in density within the flow
\begin{equation}\label{eq:sensor} \begin{equation}\label{eq:sensor}
S_e=\frac{||\rho^p_e-\rho^{p-1}_e||_{L_2}}{||\rho_e^p||_{L_2}} S_e=\frac{||\rho^p_e-\rho^{p-1}_e||_{L_2}}{||\rho_e^p||_{L_2}}
\end{equation} \end{equation}
by default the comparison is made with the $p-1$ solution, but this can be changed by setting the parameter \inltt{SensorOffset}.
An artificial diffusion term is introduced locally to the Euler equations to deal with flow discontinuity and the consequential numerical oscillations. Two models are implemented, a non-smooth and a smooth artificial viscosity model. An artificial diffusion term is introduced locally to the Euler equations to deal with flow discontinuity and the consequential numerical oscillations. Two models are implemented, a non-smooth and a smooth artificial viscosity model.
\subsubsection{Non-smooth artificial viscosity model} \subsubsection{Non-smooth artificial viscosity model}
For the non-smooth artificial viscosity model the added artificial viscosity is constant in each element and discontinuous between the elements. The Euler system is augmented by an added laplacian term on right hand side of equation \ref{eq:euler}. The diffusivity of the system is controlled by a variable viscosity coefficient $\epsilon$. The value of $\epsilon$ is dependent on $\epsilon_0$, which is the maximum viscosity that is dependent on the polynomial order ($p$), the mesh size ($h$) and the maximum wave speed and the local sensor value. Based on pre-defined sensor threshold values, the variable viscosity is set accordingly For the non-smooth artificial viscosity model the added artificial viscosity is constant in each element and discontinuous between the elements. The Euler system is augmented by an added laplacian term on right hand side of equation \ref{eq:euler}. The diffusivity of the system is controlled by a variable viscosity coefficient $\epsilon$. The value of $\epsilon$ is dependent on $\epsilon_0$, which is the maximum viscosity that is dependent on the polynomial order ($p$), the mesh size ($h$) and the maximum wave speed and the local sensor value. Based on pre-defined sensor threshold values, the variable viscosity is set accordingly
...@@ -512,4 +513,53 @@ The polynomial order in each element can be adjusted based on the sensor value t ...@@ -512,4 +513,53 @@ The polynomial order in each element can be adjusted based on the sensor value t
\end{equation} \end{equation}
For now, the threshold values $s_e$, $s_{ds}$, $s_{sm}$ and $s_{fl}$ are determined empirically by looking at the sensor distribution in the domain. Once these values are set, two .txt files are outputted, one that has the composites called VariablePComposites.txt and one with the expansions called VariablePExpansions.txt. These values have to copied into a new .xml file to create the adapted mesh. For now, the threshold values $s_e$, $s_{ds}$, $s_{sm}$ and $s_{fl}$ are determined empirically by looking at the sensor distribution in the domain. Once these values are set, two .txt files are outputted, one that has the composites called VariablePComposites.txt and one with the expansions called VariablePExpansions.txt. These values have to copied into a new .xml file to create the adapted mesh.
\subsection{Quasi-1D nozzle flow}
A quasi-1D inviscid flow (flow with area variation) can be obtained using the
\inltt{Quasi1D} forcing in a 1D solution of the Euler equations:
\begin{lstlisting}[style=XMLStyle]
<FORCING>
<FORCE TYPE="Quasi1D">
<AREAFCN> Area </AREAFCN>
</FORCE>
</FORCING>
\end{lstlisting}
in this case a function named \inltt{Area} must be specified in the \inltt{CONDITIONS} section of the session file.
In this case, it is possible to prescribe the inflow conditions in terms of stagnation properties (density and pressure)
by using the following boundary condition
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="rho" USERDEFINEDTYPE="StagnationInflow" VALUE="rhoStag" />
<D VAR="rhou" USERDEFINEDTYPE="StagnationInflow" VALUE="0" />
<D VAR="E" USERDEFINEDTYPE="StagnationInflow" VALUE="pStag/(Gamma-1)" />
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
\subsection{Axi-symmetric flow}
An axi-symmetric inviscid flow (with symmetry axis on x=0) can be obtained using
the \inltt{AxiSymmetric} forcing in a 2D solution of the Euler equations:
\begin{lstlisting}[style=XMLStyle]
<FORCING>
<FORCE TYPE="AxiSymmetric">
</FORCE>
</FORCING>
\end{lstlisting}
The \inltt{StagnationInflow} boundary condition can also be used in this case.
Also, by defining the geometry with \inltt{<GEOMETRY DIM="2" SPACE="3">} (i.e. a two-dimensional
mesh in three-dimensional space) and adding the \inltt{rhow} variable, we obtain an axi-symmetric
flow with swirl, in which case the \inltt{StagnationInflow} boundary condition allows prescribing \inltt{rhow}:
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="rho" USERDEFINEDTYPE="StagnationInflow" VALUE="rhoStag" />
<D VAR="rhou" USERDEFINEDTYPE="StagnationInflow" VALUE="0" />
<D VAR="rhov" USERDEFINEDTYPE="StagnationInflow" VALUE="0" />
<D VAR="rhow" USERDEFINEDTYPE="StagnationInflow" VALUE="x" />
<D VAR="E" USERDEFINEDTYPE="StagnationInflow" VALUE="pStag/(Gamma-1)" />
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
...@@ -147,7 +147,7 @@ class IProductWRTDerivBase_StdMat : public Operator ...@@ -147,7 +147,7 @@ class IProductWRTDerivBase_StdMat : public Operator
{ {
LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys(); LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
m_dim = PtsKey.size(); m_dim = PtsKey.size();
m_coordim = m_stdExp->GetCoordim(); m_coordim = pCollExp[0]->GetCoordim();
int nqtot = m_stdExp->GetTotPoints(); int nqtot = m_stdExp->GetTotPoints();
int nmodes = m_stdExp->GetNcoeffs(); int nmodes = m_stdExp->GetNcoeffs();
......
...@@ -66,7 +66,7 @@ ...@@ -66,7 +66,7 @@
BOOST_PP_COMMA_IF(BOOST_PP_LESS(n, BOOST_PP_SUB(BOOST_PP_ARRAY_SIZE(data), 1))) BOOST_PP_COMMA_IF(BOOST_PP_LESS(n, BOOST_PP_SUB(BOOST_PP_ARRAY_SIZE(data), 1)))
#define NEKTAR_CREATE_EXPLICIT_INTSTANTIATION(z, n, data) \ #define NEKTAR_CREATE_EXPLICIT_INTSTANTIATION(z, n, data) \
template LIB_UTILITIES_EXPORT \ template LIB_UTILITIES_EXPORT COCKCOCK \
BOOST_PP_ARRAY_ELEM(0, BOOST_PP_ARRAY_ELEM(2, data)) BOOST_PP_ARRAY_ELEM(0, data) BOOST_PP_LPAREN() \ BOOST_PP_ARRAY_ELEM(0, BOOST_PP_ARRAY_ELEM(2, data)) BOOST_PP_ARRAY_ELEM(0, data) BOOST_PP_LPAREN() \
BOOST_PP_REPEAT(BOOST_PP_ARRAY_SIZE(BOOST_PP_ARRAY_ELEM(3, data)), NEKTAR_PRINT_ARRAY, BOOST_PP_ARRAY_ELEM(3, data)) \ BOOST_PP_REPEAT(BOOST_PP_ARRAY_SIZE(BOOST_PP_ARRAY_ELEM(3, data)), NEKTAR_PRINT_ARRAY, BOOST_PP_ARRAY_ELEM(3, data)) \
BOOST_PP_COMMA_IF(BOOST_PP_GREATER(BOOST_PP_ARRAY_SIZE(BOOST_PP_ARRAY_ELEM(3, data)), 0)) \ BOOST_PP_COMMA_IF(BOOST_PP_GREATER(BOOST_PP_ARRAY_SIZE(BOOST_PP_ARRAY_ELEM(3, data)), 0)) \
......
...@@ -1111,6 +1111,12 @@ namespace Nektar ...@@ -1111,6 +1111,12 @@ namespace Nektar
TimeIntegrationSolutionSharedPtr y_out = TimeIntegrationSolutionSharedPtr y_out =
MemoryManager<TimeIntegrationSolution>::AllocateSharedPtr(m_schemeKey,y_0,time,timestep); MemoryManager<TimeIntegrationSolution>::AllocateSharedPtr(m_schemeKey,y_0,time,timestep);
if( GetIntegrationSchemeType() == eExplicit)
{
// ensure initial solution is in correct space
op.DoProjection(y_0,y_out->UpdateSolution(),time);
}
// calculate the initial derivative, if is part of the // calculate the initial derivative, if is part of the
// solution vector of the current scheme // solution vector of the current scheme
if(m_numMultiStepDerivs) if(m_numMultiStepDerivs)
...@@ -1573,8 +1579,12 @@ namespace Nektar ...@@ -1573,8 +1579,12 @@ namespace Nektar
} }
else if( type == eExplicit) else if( type == eExplicit)
{ {
// ensure solution is in correct space // Avoid projecting the same solution twice
op.DoProjection(m_Y,m_Y,m_T); if( ! ((i==0) && m_firstStageEqualsOldSolution) )
{
// ensure solution is in correct space
op.DoProjection(m_Y,m_Y,m_T);
}
op.DoOdeRhs(m_Y, m_F[i], m_T); op.DoOdeRhs(m_Y, m_F[i], m_T);
} }
} }
......
...@@ -41,7 +41,21 @@ namespace Nektar ...@@ -41,7 +41,21 @@ namespace Nektar
{ {
namespace LocalRegions namespace LocalRegions
{ {
void Expansion1D::v_NegateVertexNormal(const int vertex)
{
m_negatedNormals[vertex] = true;
for (int i = 0; i < GetCoordim(); ++i)
{
Vmath::Neg(m_vertexNormals[vertex][i].num_elements(),
m_vertexNormals[vertex][i], 1);
}
}
bool Expansion1D::v_VertexNormalNegated(const int vertex)
{
return m_negatedNormals[vertex];
}
DNekMatSharedPtr Expansion1D::v_GenMatrix(const StdRegions::StdMatrixKey &mkey) DNekMatSharedPtr Expansion1D::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
{ {
DNekMatSharedPtr returnval; DNekMatSharedPtr returnval;
......
...@@ -98,6 +98,8 @@ namespace Nektar ...@@ -98,6 +98,8 @@ namespace Nektar
inline SpatialDomains::Geometry1DSharedPtr GetGeom1D() const; inline SpatialDomains::Geometry1DSharedPtr GetGeom1D() const;
protected: protected:
std::map<int, bool> m_negatedNormals;
virtual DNekMatSharedPtr v_GenMatrix( virtual DNekMatSharedPtr v_GenMatrix(
const StdRegions::StdMatrixKey &mkey); const StdRegions::StdMatrixKey &mkey);
...@@ -110,6 +112,10 @@ namespace Nektar ...@@ -110,6 +112,10 @@ namespace Nektar
const int vert, const int vert,
const Array<OneD, const NekDouble > &primCoeffs, const Array<OneD, const NekDouble > &primCoeffs,
Array<OneD, NekDouble> &coeffs); Array<OneD, NekDouble> &coeffs);
virtual void v_NegateVertexNormal (const int vertex);
virtual bool v_VertexNormalNegated(const int vertex);
private: private:
Expansion2DWeakPtr m_elementLeft; Expansion2DWeakPtr m_elementLeft;
......
...@@ -190,6 +190,31 @@ namespace Nektar ...@@ -190,6 +190,31 @@ namespace Nektar
SetUpPhysNormals(); SetUpPhysNormals();
// Set up information for parallel and periodic problems.
for (int i = 0; i < m_trace->GetExpSize(); ++i)
{
LocalRegions::Expansion0DSharedPtr traceEl =
m_trace->GetExp(i)->as<LocalRegions::Expansion0D>();
int offset = m_trace->GetPhys_Offset(i);
int traceGeomId = traceEl->GetGeom0D()->GetGlobalID();
PeriodicMap::iterator pIt = m_periodicVerts.find(traceGeomId);
if (pIt != m_periodicVerts.end() && !pIt->second[0].isLocal)
{
if (traceGeomId != min(pIt->second[0].id, traceGeomId))
{
traceEl->GetLeftAdjacentElementExp()->NegateVertexNormal(
traceEl->GetLeftAdjacentElementVertex());
}
}
else if (m_traceMap->GetTraceToUniversalMapUnique(offset) < 0)
{
traceEl->GetLeftAdjacentElementExp()->NegateVertexNormal(
traceEl->GetLeftAdjacentElementVertex());
}
}
int cnt, n, e; int cnt, n, e;
// Identify boundary verts // Identify boundary verts
...@@ -807,6 +832,13 @@ namespace Nektar ...@@ -807,6 +832,13 @@ namespace Nektar
{ {
m_negatedFluxNormal[2*i+v] = false; m_negatedFluxNormal[2*i+v] = false;
} }
if(vertExp->GetLeftAdjacentElementExp()->
VertexNormalNegated(vertExp->GetLeftAdjacentElementVertex()))
{
m_negatedFluxNormal[2*i+v] =
(!m_negatedFluxNormal[2*i+v]);
}
} }
} }
......
...@@ -654,6 +654,24 @@ namespace Nektar ...@@ -654,6 +654,24 @@ namespace Nektar
} }
} }
void ExpList::ExponentialFilter(
Array<OneD, NekDouble> &array,
const NekDouble alpha,
const NekDouble exponent,
const NekDouble cutoff)
{
Array<OneD,NekDouble> e_array;
for(int i = 0; i < (*m_exp).size(); ++i)
{
(*m_exp)[i]->ExponentialFilter(
e_array = array+m_phys_offset[i],
alpha,
exponent,
cutoff);
}
}
/** /**
* The coefficients of the function to be acted upon * The coefficients of the function to be acted upon
* should be contained in the \param inarray. The * should be contained in the \param inarray. The
......
...@@ -254,6 +254,12 @@ namespace Nektar ...@@ -254,6 +254,12 @@ namespace Nektar
Array<OneD, NekDouble> &outarray, Array<OneD, NekDouble> &outarray,
CoeffState coeffstate = eLocal); CoeffState coeffstate = eLocal);
MULTI_REGIONS_EXPORT void ExponentialFilter(
Array<OneD, NekDouble> &array,
const NekDouble alpha,
const NekDouble exponent,
const NekDouble cutoff);
/// This function elementally mulplies the coefficient space of /// This function elementally mulplies the coefficient space of
/// Sin my the elemental inverse of the mass matrix. /// Sin my the elemental inverse of the mass matrix.
MULTI_REGIONS_EXPORT void MultiplyByElmtInvMass ( MULTI_REGIONS_EXPORT void MultiplyByElmtInvMass (
......
...@@ -139,8 +139,8 @@ namespace Nektar ...@@ -139,8 +139,8 @@ namespace Nektar
{ {
Array<OneD, NekDouble> muvar(nPts, 0.0); Array<OneD, NekDouble> muvar(nPts, 0.0);
m_ArtificialDiffusionVector(inarray, muvar); m_ArtificialDiffusionVector(inarray, muvar);
int numConvFields = nConvectiveFields; int numConvFields = nConvectiveFields;
if (m_shockCaptureType == "Smooth") if (m_shockCaptureType == "Smooth")
{ {
......
...@@ -1583,6 +1583,15 @@ namespace Nektar ...@@ -1583,6 +1583,15 @@ namespace Nektar
ASSERTL0(false, "This function is not defined in StdExpansion."); ASSERTL0(false, "This function is not defined in StdExpansion.");
} }
void StdExpansion::v_ExponentialFilter(
Array<OneD, NekDouble> &array,
const NekDouble alpha,
const NekDouble exponent,
const NekDouble cutoff)
{
ASSERTL0(false, "This function is not defined in StdExpansion.");
}
void StdExpansion::v_ReduceOrderCoeffs(int numMin, void StdExpansion::v_ReduceOrderCoeffs(int numMin,
const Array<OneD, const NekDouble> &inarray, const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray) Array<OneD, NekDouble> &outarray)
...@@ -1720,6 +1729,17 @@ namespace Nektar ...@@ -1720,6 +1729,17 @@ namespace Nektar
ASSERTL0(false, "Cannot compute vertex normal for this expansion."); ASSERTL0(false, "Cannot compute vertex normal for this expansion.");
} }
void StdExpansion::v_NegateVertexNormal(const int vertex)
{
ASSERTL0(false, "Not implemented.");
}
bool StdExpansion::v_VertexNormalNegated(const int vertex)
{
ASSERTL0(false, "Not implemented.");
return false;
}
const NormalVector & StdExpansion::v_GetFaceNormal(const int face) const const NormalVector & StdExpansion::v_GetFaceNormal(const int face) const
{ {
ASSERTL0(false, "Cannot get face normals for this expansion."); ASSERTL0(false, "Cannot get face normals for this expansion.");
......
...@@ -1000,6 +1000,14 @@ namespace Nektar ...@@ -1000,6 +1000,14 @@ namespace Nektar
v_SVVLaplacianFilter(array,mkey); v_SVVLaplacianFilter(array,mkey);
} }
void ExponentialFilter( Array<OneD, NekDouble> &array,
const NekDouble alpha,
const NekDouble exponent,
const NekDouble cutoff)
{
v_ExponentialFilter(array, alpha, exponent, cutoff);
}
void LaplacianMatrixOp(const int k1, const int k2, void LaplacianMatrixOp(const int k1, const int k2,
const Array<OneD, const NekDouble> &inarray, const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray, Array<OneD,NekDouble> &outarray,
...@@ -1305,6 +1313,16 @@ namespace Nektar ...@@ -1305,6 +1313,16 @@ namespace Nektar
v_ComputeVertexNormal(vertex); v_ComputeVertexNormal(vertex);
} }
void NegateVertexNormal(const int vertex)
{
v_NegateVertexNormal(vertex);
}
bool VertexNormalNegated(const int vertex)
{
return v_VertexNormalNegated(vertex);
}
const NormalVector & GetFaceNormal(const int face) const const NormalVector & GetFaceNormal(const int face) const
{ {
return v_GetFaceNormal(face); return v_GetFaceNormal(face);
...@@ -1793,6 +1811,12 @@ namespace Nektar ...@@ -1793,6 +1811,12 @@ namespace Nektar
STD_REGIONS_EXPORT virtual void v_SVVLaplacianFilter(Array<OneD,NekDouble> &array, STD_REGIONS_EXPORT virtual void v_SVVLaplacianFilter(Array<OneD,NekDouble> &array,
const StdMatrixKey &mkey); const StdMatrixKey &mkey);
STD_REGIONS_EXPORT virtual void v_ExponentialFilter(
Array<OneD, NekDouble> &array,
const NekDouble alpha,
const NekDouble exponent,
const NekDouble cutoff);
STD_REGIONS_EXPORT virtual void v_ReduceOrderCoeffs( STD_REGIONS_EXPORT virtual void v_ReduceOrderCoeffs(
int numMin, int numMin,
const Array<OneD, const NekDouble> &inarray, const Array<OneD, const NekDouble> &inarray,
...@@ -1857,6 +1881,10 @@ namespace Nektar ...@@ -1857,6 +1881,10 @@ namespace Nektar
STD_REGIONS_EXPORT virtual void v_ComputeVertexNormal(const int vertex); STD_REGIONS_EXPORT virtual void v_ComputeVertexNormal(const int vertex);
STD_REGIONS_EXPORT virtual void v_NegateVertexNormal(const int vertex);
STD_REGIONS_EXPORT virtual bool v_VertexNormalNegated(const int vertex);
STD_REGIONS_EXPORT virtual const NormalVector & v_GetFaceNormal(const int face) const; STD_REGIONS_EXPORT virtual const NormalVector & v_GetFaceNormal(const int face) const;
STD_REGIONS_EXPORT virtual const NormalVector & STD_REGIONS_EXPORT virtual const NormalVector &
v_GetSurfaceNormal(const int id) const; v_GetSurfaceNormal(const int id) const;
......
...@@ -69,7 +69,7 @@ namespace Nektar ...@@ -69,7 +69,7 @@ namespace Nektar
{ {
int nquad = GetTotPoints(); int nquad = GetTotPoints();
DNekMatSharedPtr D = m_base[0]->GetD(); DNekMatSharedPtr D = m_base[0]->GetD();
if( inarray.data() == outarray.data()) if( inarray.data() == outarray.data())
{ {
Array<OneD, NekDouble> wsp(nquad); Array<OneD, NekDouble> wsp(nquad);
......
...@@ -1019,6 +1019,8 @@ namespace Nektar ...@@ -1019,6 +1019,8 @@ namespace Nektar
nummodesA = nummodes1; nummodesA = nummodes1;
nummodesB = nummodes2; nummodesB = nummodes2;
break; break;
default:
ASSERTL0(false,"fid must be between 0 and 5");
} }
bool CheckForZeroedModes = false; bool CheckForZeroedModes = false;
......
...@@ -97,7 +97,8 @@ namespace Nektar ...@@ -97,7 +97,8 @@ namespace Nektar
} }
else else
{ {
Blas::Dgemv('N',nquad,m_base[0]->GetNumModes(),1.0, (m_base[0]->GetBdata()).get(), Blas::Dgemv('N',nquad,m_base[0]->GetNumModes(),1.0,
(m_base[0]->GetBdata()).get(),
nquad,&inarray[0],1,0.0,&outarray[0],1); nquad,&inarray[0],1,0.0,&outarray[0],1);
} }
} }
......
...@@ -1074,6 +1074,8 @@ namespace Nektar ...@@ -1074,6 +1074,8 @@ namespace Nektar
nummodesA = m_base[1]->GetNumModes(); nummodesA = m_base[1]->GetNumModes();
nummodesB = m_base[2]->GetNumModes(); nummodesB = m_base[2]->GetNumModes();
break; break;
default:
ASSERTL0(false,"fid must be between 0 and 4");
} }
bool CheckForZeroedModes = false; bool CheckForZeroedModes = false;
......
...@@ -1047,6 +1047,8 @@ namespace Nektar ...@@ -1047,6 +1047,8 @@ namespace Nektar
nummodesA = order1; nummodesA = order1;
nummodesB = order2; nummodesB = order2;
break; break;
default:
ASSERTL0(false,"fid must be between 0 and 4");
} }
bool CheckForZeroedModes = false; bool CheckForZeroedModes = false;
......
...@@ -1177,6 +1177,9 @@ namespace Nektar ...@@ -1177,6 +1177,9 @@ namespace Nektar
case 1: case 1:
case 3: case 3:
numModes = order1; numModes = order1;
break;
default:
ASSERTL0(false,"eid must be between 0 and 3");
} }
bool checkForZeroedModes = false; bool checkForZeroedModes = false;
...@@ -1632,6 +1635,58 @@ namespace Nektar ...@@ -1632,6 +1635,58 @@ namespace Nektar
} }
} }
void StdQuadExp::v_ExponentialFilter(
Array<OneD, NekDouble> &array,
const NekDouble alpha,
const NekDouble exponent,
const NekDouble cutoff)
{
// Generate an orthogonal expansion
int qa = m_base[0]->GetNumPoints();
int qb = m_base[1]->GetNumPoints();
int nmodesA = m_base[0]->GetNumModes();
int nmodesB = m_base[1]->GetNumModes();
int P = nmodesA - 1;
int Q = nmodesB - 1;
// Declare orthogonal basis.
LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, nmodesA, pa);
LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, nmodesB, pb);
StdQuadExp OrthoExp(Ba,Bb);
// Cutoff
int Pcut = cutoff*P;
int Qcut = cutoff*Q;
// Project onto orthogonal space.
Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
OrthoExp.FwdTrans(array,orthocoeffs);
//
NekDouble fac, fac1, fac2;
for(int i = 0; i < nmodesA; ++i)
{
for(int j = 0; j < nmodesB; ++j)
{
//to filter out only the "high-modes"
if(i > Pcut || j > Qcut)
{
fac1 = (NekDouble) (i - Pcut)/( (NekDouble)(P - Pcut) );
fac2 = (NekDouble) (j - Qcut)/( (NekDouble)(Q - Qcut) );
fac = max(fac1, fac2);
fac = pow(fac, exponent);
orthocoeffs[i*nmodesB+j] *= exp(-alpha*fac);
}