Commit dd165af7 authored by Spencer Sherwin's avatar Spencer Sherwin Committed by Chris Cantwell
Browse files

A first serial working version of SEM dealiasing with Fourier code (but not both Fourier & SEM)

parent 9f887b1c
......@@ -345,13 +345,13 @@ namespace Nektar
{
int num_dofs;
if(IsForwards)
if(IsForwards) // SJS: I reversed this check on Jan 6th
{
num_dofs = inarray.num_elements();
num_dofs = outarray.num_elements();
}
else
{
num_dofs = outarray.num_elements();
num_dofs = inarray.num_elements();
}
if(m_useFFT)
......@@ -369,8 +369,8 @@ namespace Nektar
}
else
{
Vmath::Vcopy(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),inarray,1,fft_in,1);
//fft_in = inarray;
Vmath::Vcopy(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),
inarray,1,fft_in,1);
}
if(IsForwards)
......@@ -394,8 +394,8 @@ namespace Nektar
}
else
{
Vmath::Vcopy(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),fft_out,1,outarray,1);
//outarray = fft_out;
Vmath::Vcopy(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),
fft_out,1,outarray,1);
}
}
else
......@@ -438,7 +438,6 @@ namespace Nektar
else
{
Vmath::Vcopy(ncols,inarray,1,sortedinarray,1);
//sortedinarray = inarray;
}
// Create NekVectors from the given data arrays
......@@ -455,7 +454,6 @@ namespace Nektar
else
{
Vmath::Vcopy(nrows,sortedinarray,1,outarray,1);
//outarray = sortedinarray;
}
}
......@@ -878,31 +876,31 @@ namespace Nektar
NekDouble beta;
//Half Mode
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe)
{
beta = sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
}
else if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
{
beta = -sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
}
//Fully complex
else
{
for(int i = 0; i < m_planes.num_elements(); i++)
{
beta = -sign*2*M_PI*(m_transposition->GetK(i))/m_lhom;
Vmath::Smul(nP_pts,beta,tmp1 = temparray + i*nP_pts,1,tmp2 = outarray + (i-int(sign))*nP_pts,1);
sign = -1.0*sign;
}
}
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe)
{
beta = sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
}
else if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
{
beta = -sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
}
//Fully complex
else
{
for(int i = 0; i < m_planes.num_elements(); i++)
{
beta = -sign*2*M_PI*(m_transposition->GetK(i))/m_lhom;
Vmath::Smul(nP_pts,beta,tmp1 = temparray + i*nP_pts,1,tmp2 = outarray + (i-int(sign))*nP_pts,1);
sign = -1.0*sign;
}
}
if(m_WaveSpace)
{
......
......@@ -155,12 +155,13 @@ namespace Nektar
case 3:
grad1 = Array<OneD, NekDouble> (nPointsTot);
grad2 = Array<OneD, NekDouble> (nPointsTot);
pFields[0]->PhysDeriv(pU,grad0,grad1,grad2);
if(m_dealiasing == true && pFields[0]->GetWaveSpace() == false)
{
cout << "In first option" << endl;
ASSERTL0(m_specHP_dealiasing == false,"Spectral/hp element dealaising is not set up for this option");
pFields[0]->PhysDeriv(pU,grad0,grad1,grad2);
pFields[0]->DealiasedProd(pV[0],grad0,grad0,m_CoeffState);
pFields[0]->DealiasedProd(pV[1],grad1,grad1,m_CoeffState);
pFields[0]->DealiasedProd(pV[2],grad2,grad2,m_CoeffState);
......@@ -169,53 +170,58 @@ namespace Nektar
}
else if(pFields[0]->GetWaveSpace() == true && m_dealiasing == false)
{
//vector reused to avoid even more memory requirements
//names may be misleading
pFields[0]->HomogeneousBwdTrans(grad0,wkSp);
// Put PU in physical space - In principle could use AdvVel
pFields[0]->HomogeneousBwdTrans(pU, pOutarray);
// take d/dx, d/dy gradients in physical Fourier space
pFields[0]->PhysDeriv(pOutarray,grad0,grad1);
// Take d/dz derivative using wave space field
pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],pU,
pOutarray);
pFields[0]->HomogeneousBwdTrans(pOutarray,grad2);
if(m_specHP_dealiasing) //interpolate spectral/hp gradient field
{
pFields[0]->PhysInterp1DScaled(OneDptscale,wkSp,grad0);
Vmath::Vmul(nPointsTot,grad0,1,AdvVel[0],1,Outarray,1);
pFields[0]->PhysInterp1DScaled(OneDptscale,grad0,wkSp);
Vmath::Vmul(nPointsTot,wkSp,1,AdvVel[0],1,Outarray,1);
}
else
{
Vmath::Vmul(nPointsTot,wkSp,1,AdvVel[0],1,Outarray,1);
Vmath::Vmul(nPointsTot,grad0,1,AdvVel[0],1,Outarray,1);
}
pFields[0]->HomogeneousBwdTrans(grad1,wkSp);
if(m_specHP_dealiasing) //interpolate spectral/hp gradient field
{
pFields[0]->PhysInterp1DScaled(OneDptscale,wkSp,grad1);
Vmath::Vvtvp(nPointsTot,grad1,1,AdvVel[1],1,Outarray,1,
pFields[0]->PhysInterp1DScaled(OneDptscale,grad1,wkSp);
Vmath::Vvtvp(nPointsTot,wkSp,1,AdvVel[1],1,Outarray,1,
Outarray,1);
}
else
{
Vmath::Vvtvp(nPointsTot,wkSp,1,AdvVel[1],1,Outarray,1,
Vmath::Vvtvp(nPointsTot,grad1,1,AdvVel[1],1,Outarray,1,
Outarray,1);
}
pFields[0]->HomogeneousBwdTrans(grad2,wkSp);
if(m_specHP_dealiasing) //interpolate spectral/hp gradient field
{
pFields[0]->PhysInterp1DScaled(OneDptscale,wkSp,grad2);
Vmath::Vvtvp(nPointsTot,grad2,1,AdvVel[2],1,Outarray,1,
grad0,1);
pFields[0]->PhysInterp1DScaled(OneDptscale,grad2,wkSp);
Vmath::Vvtvp(nPointsTot,wkSp,1,AdvVel[2],1,Outarray,1,Outarray,1);
pFields[0]->PhysGalerkinProjection1DScaled(OneDptscale,Outarray,pOutarray);
pFields[0]->HomogeneousFwdTrans(pOutarray,pOutarray);
}
else
{
Vmath::Vvtvp(nPointsTot,wkSp,1,AdvVel[2],1,Outarray,1,grad0,1);
}
if(m_specHP_dealiasing) // Galerkin project solution back to origianl space
{
pFields[0]->PhysGalerkinProjection1DScaled(OneDptscale,Outarray,grad0);
Vmath::Vvtvp(nPointsTot,grad2,1,AdvVel[2],1,Outarray,1,grad0,1);
pFields[0]->HomogeneousFwdTrans(grad0,pOutarray);
}
pFields[0]->HomogeneousFwdTrans(grad0,pOutarray);
}
else if(pFields[0]->GetWaveSpace() == false && m_dealiasing == false)
{
pFields[0]->PhysDeriv(pU,grad0,grad1,grad2);
if(m_specHP_dealiasing) // interpolate gradient field
{
pFields[0]->PhysInterp1DScaled(OneDptscale,grad0,wkSp);
......@@ -238,7 +244,9 @@ namespace Nektar
else if(pFields[0]->GetWaveSpace() == true && m_dealiasing == true)
{
ASSERTL0(m_specHP_dealiasing == false,"Spectral/hp element dealaising is not set up for this option");
cout << "In thrid option" << endl;
pFields[0]->PhysDeriv(pU,grad0,grad1,grad2);
pFields[0]->HomogeneousBwdTrans(grad0, pOutarray);
pFields[0]->DealiasedProd(pV[0], pOutarray, grad0,
m_CoeffState);
......
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