Commit 633fc740 authored by Ale's avatar Ale
Browse files

Turbulence, acceleration term, smoothing, FFT dealiasing and other bugs.

parent edd440c0
......@@ -57,10 +57,10 @@ namespace Nektar
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[1] = 0.0;
m_FFTW_w_inv[0] = m_N;
m_FFTW_w_inv[1] = m_FFTW_w_inv[0];
m_FFTW_w_inv[1] = 0.0;
for(int i=2;i<m_N;i++)
{
......
......@@ -256,6 +256,20 @@ 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);
}
/**
......
......@@ -228,6 +228,8 @@ 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 forwarder for MultiplyByInvMassMatrix().
MULTI_REGIONS_EXPORT virtual void v_MultiplyByInvMassMatrix(
......
......@@ -80,35 +80,45 @@ 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);
plane_two = MemoryManager<ContField2D>::AllocateSharedPtr(pSession,graph2D,variable,False,false);
m_exp = MemoryManager<StdRegions::StdExpansionVector>::AllocateSharedPtr();
nel = m_planes[0]->GetExpSize();
for(i = 0; i < nel; ++i)
{
(*m_exp).push_back(m_planes[0]->GetExp(i));
}
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();
SetupBoundaryConditions(HomoBasis,lhom,bcs,variable);
}
......@@ -147,6 +157,22 @@ namespace Nektar
m_planes[n]->GlobalToLocal();
}
};
/**
*
*/
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(
......
......@@ -64,6 +64,8 @@ namespace Nektar
/// Destructor.
MULTI_REGIONS_EXPORT virtual ~ContField3DHomogeneous1D();
MULTI_REGIONS_EXPORT virtual void v_SmoothField(Array<OneD,NekDouble> &field);
protected:
......
......@@ -659,6 +659,27 @@ namespace Nektar
e_outarray = outarray+m_coeff_offset[i]);
}
}
/**
* 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
......
......@@ -244,6 +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);
/// Solve helmholtz problem
inline void HelmSolve(
......@@ -1038,8 +1041,9 @@ namespace Nektar
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);
......@@ -1369,6 +1373,14 @@ namespace Nektar
/**
*
*/
inline void ExpList::SmoothField(Array<OneD,NekDouble> &field)
{
v_SmoothField(field);
}
/**
*
*/
inline void ExpList::BwdTrans (const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate)
......
......@@ -74,8 +74,16 @@ namespace Nektar
if(m_dealiasing)
{
ASSERTL0(m_comm->GetColumnComm()->GetSize() == 1,"Remove dealiasing if you want to run in parallel");
SetPaddingBase();
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");
}
}
}
......@@ -90,10 +98,9 @@ namespace Nektar
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),
MatBwdPAD(In.MatBwdPAD),
MatFwdPAD(In.MatFwdPAD),
m_tmpIN(In.m_tmpIN),
m_tmpOUT(In.m_tmpOUT),
m_transposition(In.m_transposition)
......@@ -136,88 +143,66 @@ 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 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 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);
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);
/////////////////////////////////////////////////////////////////////////////
// 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++)
//Looping on the pencils
for(int i = 0 ; i < num_dfts_per_proc ; 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);
Vmath::Vcopy(N,&(ShufV1[i*N]),1,&(ShufV1_PAD_coef[0]),1);
Vmath::Vcopy(N,&(ShufV2[i*N]),1,&(ShufV2_PAD_coef[0]),1);
//Moving to physical space using the padded system
PadOUT_V1 = (*MatBwdPAD)*PadIN_V1;
PadOUT_V2 = (*MatBwdPAD)*PadIN_V2;
m_FFT_deal->FFTBwdTrans(ShufV1_PAD_coef,ShufV1_PAD_phys);
m_FFT_deal->FFTBwdTrans(ShufV2_PAD_coef,ShufV2_PAD_phys);
//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);
Vmath::Vmul(m_padsize,ShufV1_PAD_phys,1,ShufV2_PAD_phys,1,ShufV1V2_PAD_phys,1);
//Moving back the result (V1*V2)_phys in Fourier space, padded system
PadOUT_Re = (*MatFwdPAD)*PadIN_Re;
m_FFT_deal->FFTFwdTrans(ShufV1V2_PAD_phys,ShufV1V2_PAD_coef);
//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
{
m_transposition->Transpose(ShufV1V2,V1V2,false,LibUtilities::eZtoXY);
//Moving the results in physical space for the output
HomogeneousBwdTrans(V1V2,outarray,coeffstate);
Vmath::Vcopy(N,&(ShufV1V2_PAD_coef[0]),1,&(ShufV1V2[i*N]),1);
}
m_transposition->Transpose(ShufV1V2,V1V2,false,LibUtilities::eZtoXY);
//Moving the results in physical space for the output
HomogeneousBwdTrans(V1V2,outarray,coeffstate);
}
/**
* Forward transform
*/
......@@ -1037,28 +1022,6 @@ namespace Nektar
v_PhysDeriv(edir,inarray,out_d);
}
/*
* Setting the Padding base for dealisaing
*/
void ExpListHomogeneous1D::SetPaddingBase(void)
{
NekDouble size = 1.5*m_homogeneousBasis->GetNumPoints();
m_padsize = int(size);
const LibUtilities::PointsKey Ppad(m_padsize,LibUtilities::eFourierEvenlySpaced);
const LibUtilities::BasisKey Bpad(LibUtilities::eFourier,m_padsize,Ppad);
m_paddingBasis = LibUtilities::BasisManager()[Bpad];
StdRegions::StdSegExp StdSeg(m_paddingBasis->GetBasisKey());
StdRegions::StdMatrixKey matkey1(StdRegions::eFwdTrans,StdSeg.DetExpansionType(),StdSeg);
StdRegions::StdMatrixKey matkey2(StdRegions::eBwdTrans,StdSeg.DetExpansionType(),StdSeg);
MatFwdPAD = StdSeg.GetStdMatrix(matkey1);
MatBwdPAD = StdSeg.GetStdMatrix(matkey2);
}
Array<OneD, unsigned int> ExpListHomogeneous1D::v_GetZIDs(void)
{
return m_transposition->GetPlanesIDs();
......
......@@ -113,9 +113,7 @@ namespace Nektar
const Array<OneD, NekDouble> &inarray2,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate = eLocal);
MULTI_REGIONS_EXPORT void SetPaddingBase(void);
LibUtilities::BasisSharedPtr GetHomogeneousBasis(void)
{
return m_homogeneousBasis;
......@@ -142,6 +140,9 @@ namespace Nektar
/// FFT variables
bool m_useFFT;
LibUtilities::NektarFFTSharedPtr m_FFT;
LibUtilities::NektarFFTSharedPtr m_FFT_deal;
Array<OneD,NekDouble> m_tmpIN;
Array<OneD,NekDouble> m_tmpOUT;
......@@ -149,7 +150,6 @@ namespace Nektar
/// quadrature points. Sets up the storage for \a m_coeff and \a
/// m_phys.
LibUtilities::BasisSharedPtr m_homogeneousBasis;
LibUtilities::BasisSharedPtr m_paddingBasis;
NekDouble m_lhom; ///< Width of homogeneous direction
Homo1DBlockMatrixMapShPtr m_homogeneous1DBlockMat;
Array<OneD, ExpListSharedPtr> m_planes;
......@@ -250,8 +250,6 @@ namespace Nektar
//Padding operations variables
bool m_dealiasing;
int m_padsize;
DNekMatSharedPtr MatFwdPAD;
DNekMatSharedPtr MatBwdPAD;
};
inline void ExpListHomogeneous1D::HomogeneousFwdTrans(const Array<OneD, const NekDouble> &inarray,
......
......@@ -130,11 +130,27 @@ namespace Nektar
Vmath::Vvtvp(nPointsTot,grad1,1,pV[1],1,pOutarray,1,pOutarray,1);
Vmath::Vvtvp(nPointsTot,grad2,1,pV[2],1,pOutarray,1,pOutarray,1);
}
else if(pFields[0]->GetWaveSpace() == true && m_dealiasing == true)
{
pFields[0]->HomogeneousBwdTrans(grad0,pOutarray);
pFields[0]->DealiasedProd(pV[0],pOutarray,grad0,m_CoeffState);
pFields[0]->HomogeneousBwdTrans(grad1,pOutarray);
pFields[0]->DealiasedProd(pV[1],pOutarray,grad1,m_CoeffState);
pFields[0]->HomogeneousBwdTrans(grad2,pOutarray);
pFields[0]->DealiasedProd(pV[2],pOutarray,grad2,m_CoeffState);
Vmath::Vadd(nPointsTot,grad0,1,grad1,1,grad0,1);
Vmath::Vadd(nPointsTot,grad0,1,grad2,1,grad0,1);
pFields[0]->HomogeneousFwdTrans(grad0,pOutarray);
}
else
{
ASSERTL0(false,"Dealiasing can't be used in combination with Wave-Space time-integration ");
ASSERTL0(false,"Advection term calculation not implented or possible with the current problem set up");
}
break;
break;
default:
ASSERTL0(false,"dimension unknown");
}
......
......@@ -181,7 +181,7 @@ namespace Nektar
}
else
{
ASSERTL0(false,"Dealiasing can't be used in combination with Wave-Space time-integration ");
ASSERTL0(false,"Dealiasing is not allowed in combination with the Skew-Symmetric advection form for efficiency reasons.");
}
break;
default:
......
......@@ -55,7 +55,8 @@ namespace Nektar
//EquationSystem(pSession),
UnsteadySystem(pSession),
m_steadyStateSteps(0),
m_subSteppingScheme(false)
m_subSteppingScheme(false),
m_SmoothAdvection(false)
{
}
......
......@@ -146,6 +146,7 @@ namespace Nektar
LibUtilities::TimeIntegrationSolutionSharedPtr m_integrationSoln;
bool m_subSteppingScheme; // bool to identify if using a substepping scheme
bool m_SmoothAdvection; // bool to identify if advection term smoothing is requested
LibUtilities::TimeIntegrationSchemeSharedPtr m_subStepIntegrationScheme;
LibUtilities::TimeIntegrationSchemeOperators m_subStepIntegrationOps;
......
......@@ -94,6 +94,8 @@ namespace Nektar
m_session->MatchSolverInfo("SubSteppingScheme","True",m_subSteppingScheme,false);
m_session->MatchSolverInfo("ShowTimings","True",m_showTimings,false);
m_session->MatchSolverInfo("SmoothAdvection","True",m_SmoothAdvection,false);
if(m_subSteppingScheme)
{
......@@ -234,6 +236,7 @@ namespace Nektar
// Storage array for high order pressure BCs
m_pressureHBCs = Array<OneD, Array<OneD, NekDouble> > (m_intSteps);
m_acceleration = Array<OneD, Array<OneD, NekDouble> > (m_intSteps+1);
m_HBCnumber = 0;
for(cnt = n = 0; n < PBndConds.num_elements(); ++n)
......@@ -248,11 +251,15 @@ namespace Nektar
if (m_HBCnumber > 0)
{
for(n = 0; n < m_intSteps; ++n)
m_acceleration[0] = Array<OneD, NekDouble>(cnt,0.0);
for(n = 0; n < m_intSteps; ++n)
{
m_pressureHBCs[n] = Array<OneD, NekDouble>(cnt);
m_pressureHBCs[n] = Array<OneD, NekDouble>(cnt,0.0);
m_acceleration[n+1] = Array<OneD, NekDouble>(cnt,0.0);
}
}
m_HBCsize = cnt;
// creating a Map to store the information regarding
// High-Order pressure BCs to improve efficiency
......@@ -410,6 +417,14 @@ namespace Nektar
{
cout << "\t Adv. eval Time : "<< timer.TimePerTest(1) << endl;
}
if(m_SmoothAdvection && m_pressureCalls > 1)
{
for(int i = 0; i < m_nConvectiveFields; ++i)
{
m_pressure->SmoothField(outarray[i]);
}
}
//add the force
if(m_session->DefinesFunction("BodyForce"))
......@@ -551,7 +566,7 @@ namespace Nektar
void VelocityCorrectionScheme::EvaluatePressureBCs(const Array<OneD, const Array<OneD, NekDouble> > &fields, const Array<OneD, const Array<OneD, NekDouble> > &N, const NekDouble Aii_Dt)
{
Array<OneD, NekDouble> tmp;
Array<OneD, NekDouble> tmp,accelerationTerm;
Array<OneD, const SpatialDomains::BoundaryConditionShPtr > PBndConds;
Array<OneD, MultiRegions::ExpListSharedPtr> PBndExp;
int n,cnt;
......@@ -560,14 +575,15 @@ namespace Nektar
PBndConds = m_pressure->GetBndConditions();
PBndExp = m_pressure->GetBndCondExpansions();
int acc_order = 0;
accelerationTerm = Array<OneD, NekDouble>(m_HBCsize,0.0);
// Reshuffle Bc Storage vector
tmp = m_pressureHBCs[nlevels-1];
for(n = nlevels-1; n > 0; --n)
{
m_pressureHBCs[n] = m_pressureHBCs[n-1];
}
m_pressureHBCs[0] = tmp;
// Roll HOPBCs storage
Roll(m_pressureHBCs);
// Roll acceleration term
Roll(m_acceleration);
// Calculate BCs at current level
CalcPressureBCs(fields,N);
......@@ -583,18 +599,39 @@ namespace Nektar
cnt += nq;
}