Commit e2bd965e authored by Kilian Lackhove's avatar Kilian Lackhove

Merge branch 'feature/ExponentialFilterHex' into 'master'

Implement v_ExponentialFilter for Hex elements

See merge request !862
parents 9ea8a971 05c4f0ef
......@@ -57,6 +57,10 @@ v5.0.0
**CompressibleFlowSolver**
- Add 3D regression tests (!567)
- Introduce forcing for quasi-1D Euler simulations (!771)
- Allow performing axi-symmetric Euler simulations (!771)
- Add ability to use an exponential filtering for stabilization with
seg, quad and hex elements (!771, !862)
**Documentation**:
- Added the developer-guide repository as a submodule (!751)
......
......@@ -2436,5 +2436,69 @@ namespace Nektar
// backward transform to physical space
OrthoExp.BwdTrans(orthocoeffs,array);
}
void StdHexExp::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 qc = m_base[2]->GetNumPoints();
int nmodesA = m_base[0]->GetNumModes();
int nmodesB = m_base[1]->GetNumModes();
int nmodesC = m_base[2]->GetNumModes();
int P = nmodesA - 1;
int Q = nmodesB - 1;
int R = nmodesC - 1;
// Declare orthogonal basis.
LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, nmodesA, pa);
LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, nmodesB, pb);
LibUtilities::BasisKey Bc(LibUtilities::eOrtho_A, nmodesC, pc);
StdHexExp OrthoExp(Ba,Bb,Bc);
// Cutoff
int Pcut = cutoff*P;
int Qcut = cutoff*Q;
int Rcut = cutoff*R;
// Project onto orthogonal space.
Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
OrthoExp.FwdTrans(array,orthocoeffs);
//
NekDouble fac, fac1, fac2, fac3;
int index = 0;
for(int i = 0; i < nmodesA; ++i)
{
for(int j = 0; j < nmodesB; ++j)
{
for(int k = 0; k < nmodesC; ++k, ++index)
{
//to filter out only the "high-modes"
if(i > Pcut || j > Qcut || k > Rcut)
{
fac1 = (NekDouble) (i - Pcut)/( (NekDouble)(P - Pcut) );
fac2 = (NekDouble) (j - Qcut)/( (NekDouble)(Q - Qcut) );
fac3 = (NekDouble) (k - Rcut)/( (NekDouble)(R - Rcut) );
fac = max( max(fac1, fac2), fac3);
fac = pow(fac, exponent);
orthocoeffs[index] *= exp(-alpha*fac);
}
}
}
}
// backward transform to physical space
OrthoExp.BwdTrans(orthocoeffs,array);
}
}
}
......@@ -285,6 +285,11 @@ namespace Nektar
STD_REGIONS_EXPORT virtual void v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,const StdMatrixKey &mkey);
STD_REGIONS_EXPORT virtual void v_ExponentialFilter(
Array<OneD, NekDouble> &array,
const NekDouble alpha,
const NekDouble exponent,
const NekDouble cutoff);
};
typedef std::shared_ptr<StdHexExp> StdHexExpSharedPtr;
......
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 to comment