Commit 7cfb55c5 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'feature/turbulence' of localhost:nektar

parents 297a8ca0 5c1e28fa
......@@ -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] = m_FFTW_w[0];
m_FFTW_w_inv[0] = m_N;
m_FFTW_w_inv[1] = m_FFTW_w_inv[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;
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;
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;
}
}
}
......@@ -257,6 +257,20 @@ namespace Nektar
}
}
/**
*
*/
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);
}
/**
* Computes the matrix vector product
......
......@@ -229,6 +229,10 @@ namespace Nektar
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate);
/// 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(
const Array<OneD, const NekDouble> &inarray,
......
......@@ -80,34 +80,57 @@ namespace Nektar
int i,n,nel;
bool False = false;
ContField2DSharedPtr plane_zero;
ContField2DSharedPtr plane_two;
SpatialDomains::BoundaryConditions bcs(m_session, graph2D);
// note that nzplanes can be larger than nzmodes
m_planes[0] = plane_zero = MemoryManager<ContField2D>::AllocateSharedPtr(pSession,graph2D,variable,False,CheckIfSingularSystem);
// Plane zero (k=0 - cos) - singularaty check required for Poisson
// problems
plane_zero = MemoryManager<ContField2D>::AllocateSharedPtr(
pSession, graph2D, variable, false,
CheckIfSingularSystem);
m_exp = MemoryManager<StdRegions::StdExpansionVector>::AllocateSharedPtr();
nel = m_planes[0]->GetExpSize();
plane_two = MemoryManager<ContField2D>::AllocateSharedPtr(
pSession, graph2D, variable, false,
false);
for(i = 0; i < nel; ++i)
{
(*m_exp).push_back(m_planes[0]->GetExp(i));
}
m_exp = MemoryManager<StdRegions::StdExpansionVector>
::AllocateSharedPtr();
for(n = 1; n < m_planes.num_elements(); ++n)
for(n = 0; n < m_planes.num_elements(); ++n)
{
m_planes[n] = MemoryManager<ContField2D>::AllocateSharedPtr(*plane_zero,graph2D,variable,False,CheckIfSingularSystem);
// 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_exp)[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);
}
......@@ -149,6 +172,24 @@ namespace Nektar
};
/**
*
*/
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)
{
m_planes[n]->SmoothField(tmp = field + cnt);
cnt += m_planes[n]->GetTotPoints();
}
}
void ContField3DHomogeneous1D::v_HelmSolve(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
......
......@@ -64,7 +64,10 @@ namespace Nektar
/// Destructor.
MULTI_REGIONS_EXPORT virtual ~ContField3DHomogeneous1D();
MULTI_REGIONS_EXPORT virtual void v_SmoothField(
Array<OneD,NekDouble> &field);
protected:
private:
......
......@@ -660,6 +660,34 @@ namespace Nektar
}
}
/**
* This function smooth a field after some calculaitons which have
* 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.");
}
/**
* This function assembles the block diagonal matrix
* \f$\underline{\boldsymbol{M}}^e\f$, which is the
......
......@@ -245,6 +245,9 @@ namespace Nektar
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate = eLocal);
/// Smooth a field across elements
inline void SmoothField(Array<OneD,NekDouble> &field);
/// Solve helmholtz problem
inline void HelmSolve(
const Array<OneD, const NekDouble> &inarray,
......@@ -1037,9 +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_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);
......@@ -1365,8 +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,13 +71,23 @@ namespace Nektar
{
m_FFT = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW", m_homogeneousBasis->GetNumPoints());
}
if(m_dealiasing)
{
ASSERTL0(m_comm->GetColumnComm()->GetSize() == 1,"Remove dealiasing if you want to run in parallel");
SetPaddingBase();
}
}
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");
}
}
}
/**
......@@ -89,14 +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_dealiasing(In.m_dealiasing),
m_padsize(In.m_padsize),
MatBwdPAD(In.MatBwdPAD),
MatFwdPAD(In.MatFwdPAD),
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());
}
......@@ -136,88 +145,76 @@ namespace Nektar
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate)
{
// inarray1 = first term of the product
// inarray2 = second term of the product
// 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 npoints = outarray.num_elements(); // number of total physical points
int nplanes = m_planes.num_elements(); // number of planes == number of Fourier modes = number of Fourier coeff
int npencils = npoints/nplanes; // number of pencils = numebr of physical points per plane
Array<OneD, NekDouble> V1(npoints);
Array<OneD, NekDouble> V2(npoints);
Array<OneD, NekDouble> V1V2(npoints);
Array<OneD, NekDouble> ShufV1(npoints);
Array<OneD, NekDouble> ShufV2(npoints);
Array<OneD, NekDouble> ShufV1V2(npoints);
if(m_WaveSpace)
{
V1 = inarray1;
V2 = inarray2;
}
else
{
HomogeneousFwdTrans(inarray1,V1,coeffstate);
HomogeneousFwdTrans(inarray2,V2,coeffstate);
}
m_transposition->Transpose(V1,ShufV1,false,LibUtilities::eXYtoZ);
m_transposition->Transpose(V2,ShufV2,false,LibUtilities::eXYtoZ);
/////////////////////////////////////////////////////////////////////////////
// Creating padded vectors for each pencil
Array<OneD, NekDouble> PadV1_pencil_coeff(m_padsize,0.0);
Array<OneD, NekDouble> PadV2_pencil_coeff(m_padsize,0.0);
Array<OneD, NekDouble> PadRe_pencil_coeff(m_padsize,0.0);
Array<OneD, NekDouble> PadV1_pencil_phys(m_padsize,0.0);
Array<OneD, NekDouble> PadV2_pencil_phys(m_padsize,0.0);
Array<OneD, NekDouble> PadRe_pencil_phys(m_padsize,0.0);
NekVector<NekDouble> PadIN_V1(m_padsize,PadV1_pencil_coeff,eWrapper);
NekVector<NekDouble> PadOUT_V1(m_padsize,PadV1_pencil_phys,eWrapper);
NekVector<NekDouble> PadIN_V2(m_padsize,PadV2_pencil_coeff,eWrapper);
NekVector<NekDouble> PadOUT_V2(m_padsize,PadV2_pencil_phys,eWrapper);
NekVector<NekDouble> PadIN_Re(m_padsize,PadRe_pencil_phys,eWrapper);
NekVector<NekDouble> PadOUT_Re(m_padsize,PadRe_pencil_coeff,eWrapper);
//Looping on the pencils
for(int i = 0 ; i< npencils ; i++)
{
//Copying the i-th pencil pf lenght N into a bigger
//pencil of lenght 2N We are in Fourier space
Vmath::Vcopy(nplanes,&(ShufV1[i*nplanes]),1,&(PadV1_pencil_coeff[0]),1);
Vmath::Vcopy(nplanes,&(ShufV2[i*nplanes]),1,&(PadV2_pencil_coeff[0]),1);
//Moving to physical space using the padded system
PadOUT_V1 = (*MatBwdPAD)*PadIN_V1;
PadOUT_V2 = (*MatBwdPAD)*PadIN_V2;
//Perfroming the vectors multiplication in physical space on the padded system
Vmath::Vmul(m_padsize,PadV1_pencil_phys,1,PadV2_pencil_phys,1,PadRe_pencil_phys,1);
//Moving back the result (V1*V2)_phys in Fourier space, padded system
PadOUT_Re = (*MatFwdPAD)*PadIN_Re;
//Copying the first half of the padded pencil in the full vector (Fourier space)
Vmath::Vcopy(nplanes,&(PadRe_pencil_coeff[0]),1,&(ShufV1V2[i*nplanes]),1);
}
if(m_WaveSpace)
{
m_transposition->Transpose(ShufV1V2,outarray,false,LibUtilities::eZtoXY);
}
else
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_proc = m_comm->GetColumnComm()->GetSize();
int num_dfts_per_proc = num_points_per_plane / num_proc
+ (num_points_per_plane % num_proc > 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);