Commit 5c1e28fa authored by Ale's avatar Ale
Browse files

Gettin last modification before final mergeMerge branch 'feature/turbulence'...

Gettin last modification before final mergeMerge branch 'feature/turbulence' of gitlab.nektar.info:nektar into feature/turbulence
parents 147b0bc3 d71b9447
......@@ -41,106 +41,110 @@ namespace Nektar
{
namespace LibUtilities
{
string NekFFTW::className = GetNektarFFTFactory().RegisterCreatorFunction("NekFFTW", NekFFTW::create);
NekFFTW::NekFFTW(int N):
NektarFFT(N)
{
m_wsp = Array<OneD, NekDouble>(m_N);
phys = Array<OneD,NekDouble>(m_N);
coef = Array<OneD,NekDouble>(m_N);
plan_forward = fftw_plan_r2r_1d(m_N,&phys[0],&coef[0],FFTW_R2HC,FFTW_ESTIMATE);
plan_backward = fftw_plan_r2r_1d(m_N,&coef[0],&phys[0],FFTW_HC2R,FFTW_ESTIMATE);
m_FFTW_w = Array<OneD,NekDouble>(m_N);
string NekFFTW::className
= GetNektarFFTFactory().RegisterCreatorFunction("NekFFTW",
NekFFTW::create);
NekFFTW::NekFFTW(int N)
: NektarFFT(N)
{
m_wsp = Array<OneD, NekDouble>(m_N);
phys = Array<OneD,NekDouble>(m_N);
coef = Array<OneD,NekDouble>(m_N);
plan_forward = fftw_plan_r2r_1d(m_N, &phys[0], &coef[0],
FFTW_R2HC, FFTW_ESTIMATE);
plan_backward = fftw_plan_r2r_1d(m_N, &coef[0], &phys[0],
FFTW_HC2R, FFTW_ESTIMATE);
m_FFTW_w = Array<OneD,NekDouble>(m_N);
m_FFTW_w_inv = Array<OneD,NekDouble>(m_N);
m_FFTW_w[0] = 1.0/(double)m_N;
m_FFTW_w[1] = 0.0;
m_FFTW_w_inv[0] = m_N;
m_FFTW_w[0] = 1.0/(double)m_N;
m_FFTW_w[1] = 0.0;
m_FFTW_w_inv[0] = m_N;
m_FFTW_w_inv[1] = 0.0;
for(int i=2;i<m_N;i++)
{
m_FFTW_w[i]=m_FFTW_w[0];
m_FFTW_w_inv[i]=m_FFTW_w_inv[0];
}
}
// Distructor
NekFFTW::~NekFFTW()
{
}
//=================================================================================
// Forward transformation
void NekFFTW::v_FFTFwdTrans(Array<OneD,NekDouble> &inarray, Array<OneD,NekDouble> &outarray)
{
Vmath::Vcopy(m_N, inarray, 1, phys, 1);
fftw_execute(plan_forward);
Reshuffle_FFTW2Nek(coef);
Vmath::Vcopy(m_N, coef, 1, outarray, 1);
}
//=================================================================================
// Backward transformation
void NekFFTW::v_FFTBwdTrans(Array<OneD,NekDouble> &inarray, Array<OneD,NekDouble> &outarray)
{
Vmath::Vcopy(m_N, inarray, 1, coef, 1);
Reshuffle_Nek2FFTW(coef);
fftw_execute(plan_backward);
Vmath::Vcopy(m_N, phys, 1, outarray, 1);
}
//==================================================================================
// Reshuffle FFTW2Nek
void NekFFTW::Reshuffle_FFTW2Nek(Array<OneD,NekDouble> &coef)
{
int halfN = m_N/2;
m_wsp[1]=coef[halfN];
Vmath::Vcopy(halfN, coef, 1, m_wsp, 2);
for(int i=0;i<(halfN-1);i++)
{
m_wsp[(m_N-1)-2*i]=coef[halfN+1+i];
}
Vmath::Vmul(m_N, m_wsp, 1, m_FFTW_w, 1, coef, 1);
return;
}
//===================================================================================
// Reshuffle Nek2FFTW
void NekFFTW::Reshuffle_Nek2FFTW(Array<OneD,NekDouble> &coef)
{
int halfN = m_N/2;
for(int i=2;i<m_N;i++)
{
m_FFTW_w[i] = m_FFTW_w[0];
m_FFTW_w_inv[i] = m_FFTW_w_inv[0];
}
}
// Distructor
NekFFTW::~NekFFTW()
{
}
// Forward transformation
void NekFFTW::v_FFTFwdTrans(
Array<OneD,NekDouble> &inarray,
Array<OneD,NekDouble> &outarray)
{
Vmath::Vcopy(m_N, inarray, 1, phys, 1);
fftw_execute(plan_forward);
Reshuffle_FFTW2Nek(coef);
Vmath::Vcopy(m_N, coef, 1, outarray, 1);
}
// Backward transformation
void NekFFTW::v_FFTBwdTrans(
Array<OneD,NekDouble> &inarray,
Array<OneD,NekDouble> &outarray)
{
Vmath::Vcopy(m_N, inarray, 1, coef, 1);
Reshuffle_Nek2FFTW(coef);
fftw_execute(plan_backward);
Vmath::Vcopy(m_N, phys, 1, outarray, 1);
}
// Reshuffle FFTW2Nek
void NekFFTW::Reshuffle_FFTW2Nek(Array<OneD,NekDouble> &coef)
{
int halfN = m_N/2;
m_wsp[1] = coef[halfN];
Vmath::Vcopy(halfN, coef, 1, m_wsp, 2);
for(int i = 0; i < (halfN - 1); i++)
{
m_wsp[(m_N-1)-2*i] = coef[halfN+1+i];
}
Vmath::Vmul(m_N, m_wsp, 1, m_FFTW_w, 1, coef, 1);
return;
}
// Reshuffle Nek2FFTW
void NekFFTW::Reshuffle_Nek2FFTW(Array<OneD,NekDouble> &coef)
{
int halfN = m_N/2;
Vmath::Vmul(m_N, coef, 1, m_FFTW_w_inv, 1, coef, 1);
m_wsp[halfN]=coef[1];
Vmath::Vcopy(halfN, coef, 2, m_wsp, 1);
for(int i=0;i<(halfN-1);i++)
{
m_wsp[halfN+1+i]=coef[(m_N-1)-2*i];
}
Vmath::Smul(m_N, 1.0/m_N, m_wsp, 1, coef, 1);
return;
}
//===================================================================================
}//end namespace LibUtilities
}//end of namespace Nektar
m_wsp[halfN] = coef[1];
Vmath::Vcopy(halfN, coef, 2, m_wsp, 1);
for(int i = 0; i < (halfN-1); i++)
{
m_wsp[halfN+1+i] = coef[(m_N-1)-2*i];
}
Vmath::Smul(m_N, 1.0/m_N, m_wsp, 1, coef, 1);
return;
}
}
}
......@@ -256,19 +256,19 @@ namespace Nektar
GlobalToLocal(tmp,outarray);
}
}
/**
/**
*
*/
void ContField2D::v_SmoothField(Array<OneD,NekDouble> &field)
{
int gloNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
Array<OneD,NekDouble> tmp1(gloNcoeffs);
Array<OneD,NekDouble> tmp2(gloNcoeffs);
IProductWRTBase(field,tmp1,eGlobal);
MultiplyByInvMassMatrix(tmp1,tmp2,eGlobal);
BwdTrans(tmp2,field,eGlobal);
int gloNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
Array<OneD,NekDouble> tmp1(gloNcoeffs);
Array<OneD,NekDouble> tmp2(gloNcoeffs);
IProductWRTBase(field,tmp1,eGlobal);
MultiplyByInvMassMatrix(tmp1,tmp2,eGlobal);
BwdTrans(tmp2,field,eGlobal);
}
......
......@@ -228,8 +228,10 @@ namespace Nektar
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate);
MULTI_REGIONS_EXPORT virtual void v_SmoothField(Array<OneD,NekDouble> &field);
/// Template method virtual forwarded for SmoothField().
MULTI_REGIONS_EXPORT virtual void v_SmoothField(
Array<OneD,NekDouble> &field);
/// Template method virtual forwarder for MultiplyByInvMassMatrix().
MULTI_REGIONS_EXPORT virtual void v_MultiplyByInvMassMatrix(
......
......@@ -80,45 +80,58 @@ namespace Nektar
int i,n,nel;
bool False = false;
ContField2DSharedPtr plane_zero;
ContField2DSharedPtr plane_two;
ContField2DSharedPtr plane_two;
SpatialDomains::BoundaryConditions bcs(m_session, graph2D);
// Plane zero (k=0 - cos) - singularaty check required for Poisson problems
plane_zero = MemoryManager<ContField2D>::AllocateSharedPtr(pSession,graph2D,variable,False,CheckIfSingularSystem);
plane_two = MemoryManager<ContField2D>::AllocateSharedPtr(pSession,graph2D,variable,False,false);
m_exp = MemoryManager<StdRegions::StdExpansionVector>::AllocateSharedPtr();
// Plane zero (k=0 - cos) - singularaty check required for Poisson
// problems
plane_zero = MemoryManager<ContField2D>::AllocateSharedPtr(
pSession, graph2D, variable, false,
CheckIfSingularSystem);
plane_two = MemoryManager<ContField2D>::AllocateSharedPtr(
pSession, graph2D, variable, false,
false);
m_exp = MemoryManager<StdRegions::StdExpansionVector>
::AllocateSharedPtr();
for(n = 0; n < m_planes.num_elements(); ++n)
{
//Plane zero and one (k=0 - cos and sin) - singularaty check required for Poisson problems
if(m_transposition->GetK(n) == 0)
{
m_planes[n] = MemoryManager<ContField2D>::AllocateSharedPtr(*plane_zero,graph2D,variable,False,CheckIfSingularSystem);
}
else
{
// For k > 0 singularty check not required anymore - creating another ContField2D to avoid Assembly Map copy
// TODO: We may want to deal with it in a more efficient way in the future.
m_planes[n] = MemoryManager<ContField2D>::AllocateSharedPtr(*plane_two,graph2D,variable,False,false);
}
nel = m_planes[n]->GetExpSize();
// Plane zero and one (k=0 - cos and sin) - singularaty check
// required for Poisson problems
if(m_transposition->GetK(n) == 0)
{
m_planes[n] = MemoryManager<ContField2D>
::AllocateSharedPtr(*plane_zero, graph2D, variable,
false, CheckIfSingularSystem);
}
else
{
// For k > 0 singularty check not required anymore -
// creating another ContField2D to avoid Assembly Map copy
// TODO: We may want to deal with it in a more efficient
// way in the future.
m_planes[n] = MemoryManager<ContField2D>
::AllocateSharedPtr(*plane_two, graph2D, variable,
false, false);
}
nel = m_planes[n]->GetExpSize();
for(i = 0; i < nel; ++i)
{
(*m_exp).push_back(m_planes[n]->GetExp(i));
}
}
nel = GetExpSize();
m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>::AllocateSharedPtr(nel);
SetCoeffPhys();
SetCoeffPhys();
SetupBoundaryConditions(HomoBasis,lhom,bcs,variable);
}
......@@ -157,20 +170,22 @@ namespace Nektar
m_planes[n]->GlobalToLocal();
}
};
/**
/**
*
*/
void ContField3DHomogeneous1D::v_SmoothField(Array<OneD,NekDouble> &field)
void ContField3DHomogeneous1D::v_SmoothField(
Array<OneD,NekDouble> &field)
{
int cnt = 0;
Array<OneD, NekDouble> tmp;
for(int n = 0; n < m_planes.num_elements(); ++n)
int cnt = 0;
Array<OneD, NekDouble> tmp;
for(int n = 0; n < m_planes.num_elements(); ++n)
{
m_planes[n]->SmoothField(tmp = field + cnt);
cnt += m_planes[n]->GetTotPoints();
cnt += m_planes[n]->GetTotPoints();
}
}
......
......@@ -64,9 +64,10 @@ namespace Nektar
/// Destructor.
MULTI_REGIONS_EXPORT virtual ~ContField3DHomogeneous1D();
MULTI_REGIONS_EXPORT virtual void v_SmoothField(Array<OneD,NekDouble> &field);
MULTI_REGIONS_EXPORT virtual void v_SmoothField(
Array<OneD,NekDouble> &field);
protected:
private:
......
......@@ -659,27 +659,34 @@ namespace Nektar
e_outarray = outarray+m_coeff_offset[i]);
}
}
/**
/**
* This function smooth a field after some calculaitons which have
* been done elementally.
* been done elementally.
*
* @param field An array containing the field in physical space
*
*/
void ExpList::v_SmoothField(Array<OneD, NekDouble> &field)
{
// Do nothing unless the method is implemented in the appropriate class, i.e. ContField1D,ContField2D, etc.
// So far it has been implemented just for ContField2D and ContField3DHomogeneous1D
// Block in case users try the smoothing with a modal expansion.
// Maybe a different techique for the smoothing require implementation for modal basis.
ASSERTL0((*m_exp)[0]->GetBasisType(0) == LibUtilities::eGLL_Lagrange,
"Smoothing is currently not allowed unless you are using a nodal base for efficiency reasons."
"The implemented smoothing technique requires the mass matrix inversion which is trivial just for GLL_LAGRANGE_SEM expansions.");
}
void ExpList::v_SmoothField(Array<OneD, NekDouble> &field)
{
// Do nothing unless the method is implemented in the appropriate
// class, i.e. ContField1D,ContField2D, etc.
// So far it has been implemented just for ContField2D and
// ContField3DHomogeneous1D
// Block in case users try the smoothing with a modal expansion.
// Maybe a different techique for the smoothing require
// implementation for modal basis.
ASSERTL0((*m_exp)[0]->GetBasisType(0)
== LibUtilities::eGLL_Lagrange,
"Smoothing is currently not allowed unless you are using "
"a nodal base for efficiency reasons. The implemented "
"smoothing technique requires the mass matrix inversion "
"which is trivial just for GLL_LAGRANGE_SEM expansions.");
}
/**
* This function assembles the block diagonal matrix
......
......@@ -244,9 +244,9 @@ namespace Nektar
const Array<OneD,const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate = eLocal);
/// Smooth a field across elements
inline void SmoothField(Array<OneD,NekDouble> &field);
/// Smooth a field across elements
inline void SmoothField(Array<OneD,NekDouble> &field);
/// Solve helmholtz problem
inline void HelmSolve(
......@@ -1040,10 +1040,12 @@ namespace Nektar
virtual void v_FwdTrans(const Array<OneD,const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate);
virtual void v_FwdTrans_IterPerExp(const Array<OneD,const NekDouble> &inarray, Array<OneD,NekDouble> &outarray);
virtual void v_SmoothField(Array<OneD,NekDouble> &field);
virtual void v_FwdTrans_IterPerExp(
const Array<OneD,const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray);
virtual void v_SmoothField(Array<OneD,NekDouble> &field);
virtual void v_IProductWRTBase(const Array<OneD,const NekDouble> &inarray,Array<OneD, NekDouble> &outarray, CoeffState coeffstate);
......@@ -1369,16 +1371,16 @@ namespace Nektar
{
v_FwdTrans_IterPerExp(inarray,outarray);
}
/**
/**
*
*/
inline void ExpList::SmoothField(Array<OneD,NekDouble> &field)
{
v_SmoothField(field);
}
/**
/**
*
*/
inline void ExpList::BwdTrans (const Array<OneD, const NekDouble> &inarray,
......
......@@ -71,21 +71,23 @@ namespace Nektar
{
m_FFT = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW", m_homogeneousBasis->GetNumPoints());
}
if(m_dealiasing)
{
if(m_useFFT)
{
NekDouble size = 1.5*m_homogeneousBasis->GetNumPoints();
m_padsize = int(size);
m_FFT_deal = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW", m_padsize);
}
else
{
ASSERTL0(false,"Dealiasing available just in combiation with FFTW");
}
}
}
if(m_dealiasing)
{
if(m_useFFT)
{
NekDouble size = 1.5*m_homogeneousBasis->GetNumPoints();
m_padsize = int(size);
m_FFT_deal = LibUtilities::GetNektarFFTFactory()
.CreateInstance("NekFFTW", m_padsize);
}
else
{
ASSERTL0(false, "Dealiasing available just in combination "
"with FFTW");
}
}
}
/**
......@@ -97,13 +99,13 @@ namespace Nektar
m_homogeneous1DBlockMat(In.m_homogeneous1DBlockMat),
m_lhom(In.m_lhom),
m_useFFT(In.m_useFFT),
m_FFT(In.m_FFT),
m_FFT_deal(In.m_FFT_deal),
m_dealiasing(In.m_dealiasing),
m_padsize(In.m_padsize),
m_FFT(In.m_FFT),
m_FFT_deal(In.m_FFT_deal),
m_dealiasing(In.m_dealiasing),
m_padsize(In.m_padsize),
m_tmpIN(In.m_tmpIN),
m_tmpOUT(In.m_tmpOUT),
m_transposition(In.m_transposition)
m_transposition(In.m_transposition)
{
m_planes = Array<OneD, ExpListSharedPtr>(In.m_planes.num_elements());
}
......@@ -146,63 +148,73 @@ namespace Nektar
// inarray1 = first term of the product in full physical space
// inarray2 = second term of the product in full physical space
// dealiased product stored in outarray
int num_dofs = inarray1.num_elements();
int N = m_homogeneousBasis->GetNumPoints();
int num_dofs = inarray1.num_elements();
int N = m_homogeneousBasis->GetNumPoints();
Array<OneD, NekDouble> V1(num_dofs);
Array<OneD, NekDouble> V2(num_dofs);
Array<OneD, NekDouble> V1V2(num_dofs);
HomogeneousFwdTrans(inarray1,V1,coeffstate);
HomogeneousFwdTrans(inarray2,V2,coeffstate);
int num_points_per_plane = num_dofs/m_planes.num_elements();
int num_dfts_per_proc = num_points_per_plane/m_comm->GetColumnComm()->GetSize() + (num_points_per_plane%m_comm->GetColumnComm()->GetSize() > 0);
Array<OneD, NekDouble> ShufV1(num_dfts_per_proc*N,0.0);
Array<OneD, NekDouble> ShufV2(num_dfts_per_proc*N,0.0);
Array<OneD, NekDouble> ShufV1V2(num_dfts_per_proc*N,0.0);
Array<OneD, NekDouble> ShufV1_PAD_coef(m_padsize,0.0);
Array<OneD, NekDouble> ShufV2_PAD_coef(m_padsize,0.0);
Array<OneD, NekDouble> ShufV1_PAD_phys(m_padsize,0.0);
Array<OneD, NekDouble> ShufV2_PAD_phys(m_padsize,0.0);
Array<OneD, NekDouble> ShufV1V2_PAD_coef(m_padsize,0.0);
Array<OneD, NekDouble> ShufV1V2_PAD_phys(m_padsize,0.0);
m_transposition->Transpose(V1,ShufV1,false,LibUtilities::eXYtoZ);
m_transposition->Transpose(V2,ShufV2,false,LibUtilities::eXYtoZ);
//Looping on the pencils
Array<OneD, NekDouble> V1V2(num_dofs);
HomogeneousFwdTrans(inarray1,V1,coeffstate);
HomogeneousFwdTrans(inarray2,V2,coeffstate);