Commit b2cdc3e7 authored by Alessandro Bolis's avatar Alessandro Bolis
Browse files

Fourier basis moved to a different sign convenction.

First immaginary plane solution set manually to zero.
Efficiency of the Skew-symmetric advection term. 


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4068 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent 401c4617
......@@ -50,10 +50,10 @@ int main(int argc, char *argv[])
//----------------------------------------------
// Define Expansion
lz = vSession->GetParameter("Lz");
int nplanes = vSession->GetParameter("HomModesZ");
lz = vSession->GetParameter("LZ");
bool useFFT = false;
bool deal = false;
int nplanes = 8;
const LibUtilities::PointsKey Pkey(nplanes,LibUtilities::eFourierEvenlySpaced);
const LibUtilities::BasisKey Bkey(LibUtilities::eFourier,nplanes,Pkey);
Exp = MemoryManager<MultiRegions::DisContField3DHomogeneous1D>::
......
......@@ -39,8 +39,8 @@ int main(int argc, char *argv[])
//----------------------------------------------
// Define Expansion
int bc_val = 0;
int nplanes = 8;
NekDouble lz = vSession->GetParameter("Lz");
int nplanes = vSession->GetParameter("HomModesZ");
NekDouble lz = vSession->GetParameter("LZ");
bool useFFT = false;
bool deal = false;
const LibUtilities::PointsKey Pkey(nplanes,LibUtilities::eFourierEvenlySpaced);
......
......@@ -62,8 +62,8 @@ namespace Nektar
for(int i=2;i<m_N;i++)
{
m_FFTW_w[i]=2.0*m_FFTW_w[0];
m_FFTW_w_inv[i]=0.5*m_FFTW_w_inv[0];
m_FFTW_w[i]=m_FFTW_w[0];
m_FFTW_w_inv[i]=m_FFTW_w_inv[0];
}
}
......@@ -115,17 +115,17 @@ namespace Nektar
Vmath::Vmul(m_N, m_wsp, 1, m_FFTW_w, 1, coef, 1);
if(m_N>2)
{
int sign = -1;
for(int i=2;i<(m_N-1);i=i+2)
{
coef[i]=coef[i]*sign;
sign=sign*(-1);
coef[i+1]=coef[i+1]*sign;
}
}
//if(m_N>2)
//{
//int sign = 1;
//for(int i=2;i<(m_N-1);i=i+2)
//{
//
// coef[i]=coef[i]*sign;
// sign=sign*(-1);
// coef[i+1]=coef[i+1]*sign;
//}
//}
return;
}
......@@ -137,17 +137,17 @@ namespace Nektar
Vmath::Vmul(m_N, coef, 1, m_FFTW_w_inv, 1, coef, 1);
if(m_N>2)
{
int sign = -1;
for(int i=2;i<(m_N-1);i=i+2)
{
coef[i]=coef[i]*sign;
sign=sign*(-1);
coef[i+1]=coef[i+1]*sign;
}
}
//if(m_N>2)
//{
// int sign = -1;
// for(int i=2;i<(m_N-1);i=i+2)
// {
//
// coef[i]=coef[i]*sign;
// sign=sign*(-1);
// coef[i+1]=coef[i+1]*sign;
// }
//}
m_wsp[halfN]=coef[1];
......
......@@ -570,10 +570,10 @@ namespace Nektar
for(i = 0; i < numPoints; ++i)
{
m_bdata[ 2*p *numPoints+i] = cos(p*M_PI*z[i]);
m_bdata[(2*p+1)*numPoints+i] = sin(p*M_PI*z[i]);
m_bdata[(2*p+1)*numPoints+i] = -sin(p*M_PI*z[i]);
m_dbdata[ 2*p *numPoints+i] = -p*M_PI*sin(p*M_PI*z[i]);
m_dbdata[(2*p+1)*numPoints+i] = p*M_PI*cos(p*M_PI*z[i]);
m_dbdata[(2*p+1)*numPoints+i] = -p*M_PI*cos(p*M_PI*z[i]);
}
}
......@@ -586,10 +586,10 @@ namespace Nektar
for(i = 0; i < numPoints; ++i)
{
m_bdata[i] = cos(M_PI*z[i]);
m_bdata[numPoints+i] = sin(M_PI*z[i]);
m_bdata[numPoints+i] = -sin(M_PI*z[i]);
m_dbdata[i] = -M_PI*sin(M_PI*z[i]);
m_dbdata[numPoints+i] = M_PI*cos(M_PI*z[i]);
m_dbdata[numPoints+i] = -M_PI*cos(M_PI*z[i]);
}
for (p=1; p < numModes/2; ++p)
......@@ -615,8 +615,8 @@ namespace Nektar
//Fourier Imaginary Half Mode
case eFourierHalfModeIm:
m_bdata[0] = sin(M_PI*z[0]);
m_dbdata[0] = M_PI*cos(M_PI*z[0]);
m_bdata[0] = -sin(M_PI*z[0]);
m_dbdata[0] = -M_PI*cos(M_PI*z[0]);
break;
......
......@@ -76,13 +76,13 @@ namespace Nektar
unsigned int npts = m_pointsKey.GetNumPoints();
if(npts==1)
{
m_weights[0] = 2.0; //midpoint rule
m_weights[0] = 1.0; //midpoint rule
}
else
{
for(unsigned int i=0; i<npts; ++i)
{
m_weights[i] = 2.0/(double)npts;
m_weights[i] = 1.0/(double)npts;
}
}
}
......
......@@ -708,6 +708,7 @@ namespace Nektar
int i,j;
int bndcnt=0;
Array<OneD, NekDouble> gamma(contNcoeffs, 0.0);
for(i = 0; i < m_bndCondExpansions.num_elements(); ++i)
{
if(m_bndConditions[i]->GetBoundaryConditionType() != SpatialDomains::eDirichlet)
......@@ -724,6 +725,7 @@ namespace Nektar
bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
}
}
m_locToGloMap->UniversalAssemble(gamma);
// Add weak boundary conditions to forcing
......
......@@ -163,16 +163,27 @@ namespace Nektar
{
HomogeneousFwdTrans(inarray,fce,(flags.isSet(eUseGlobal))?eGlobal:eLocal);
}
bool smode = false;
if (m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe ||
m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm )
{
smode = true;
}
for(n = 0; n < m_planes.num_elements(); ++n)
{
beta = 2*M_PI*(m_transposition->GetK(n))/m_lhom;
new_factors = factors;
new_factors[StdRegions::eFactorLambda] += beta*beta;
if(n != 1 || m_transposition->GetK(n) != 0 || smode)
{
beta = 2*M_PI*(m_transposition->GetK(n))/m_lhom;
new_factors = factors;
new_factors[StdRegions::eFactorLambda] += beta*beta;
m_planes[n]->HelmSolve(fce + cnt,
m_planes[n]->HelmSolve(fce + cnt,
e_out = outarray + cnt1,
flags, new_factors, varcoeff, dirForcing);
}
cnt += m_planes[n]->GetTotPoints();
cnt1 += m_planes[n]->GetNcoeffs();
......
......@@ -199,7 +199,7 @@ namespace Nektar
for(n = 0; n < m_planes.num_elements(); ++n)
{
m_planes[n]->EvaluateBoundaryConditions(time,0.5*m_lhom*(1.0+local_z[n]));
m_planes[n]->EvaluateBoundaryConditions(time,0.5*m_lhom*(1.0+local_z[n]));
}
// Fourier transform coefficient space boundary values
......@@ -238,12 +238,15 @@ namespace Nektar
for(n = 0; n < m_planes.num_elements(); ++n)
{
beta = 2*M_PI*(m_transposition->GetK(n))/m_lhom;
new_factors = factors;
new_factors[StdRegions::eFactorLambda] += beta*beta;
m_planes[n]->HelmSolve(fce + cnt,
e_out = outarray + cnt1,
flags, new_factors, varcoeff, dirForcing);
if(n != 1 || m_transposition->GetK(n) != 0)
{
beta = 2*M_PI*(m_transposition->GetK(n))/m_lhom;
new_factors = factors;
new_factors[StdRegions::eFactorLambda] += beta*beta;
m_planes[n]->HelmSolve(fce + cnt,
e_out = outarray + cnt1,
flags, new_factors, varcoeff, dirForcing);
}
cnt += m_planes[n]->GetTotPoints();
cnt1 += m_planes[n]->GetNcoeffs();
......
......@@ -857,26 +857,31 @@ namespace Nektar
NekDouble beta;
//Half Mode
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe ||
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)
{
......@@ -953,25 +958,30 @@ namespace Nektar
NekDouble beta;
//HalfMode
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe ||
m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
{
beta = 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 = 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 = -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)
{
out_d = outarray;
......
......@@ -773,7 +773,7 @@ namespace Nektar
//along y
for(int i = 0; i < m_ny; i++)
{
beta = sign*2*M_PI*(i/2)/m_lhom_y;
beta = -sign*2*M_PI*(i/2)/m_lhom_y;
for(int j = 0; j < m_nz; j++)
{
......@@ -787,7 +787,7 @@ namespace Nektar
sign = -1.0;
for(int i = 0; i < m_nz; i++)
{
beta = sign*2*M_PI*(i/2)/m_lhom_z;
beta = -sign*2*M_PI*(i/2)/m_lhom_z;
Vmath::Smul(m_ny*n_points_line,beta,tmp1 = temparray + i*m_ny*n_points_line,1,tmp2 = temparray2 + (i-int(sign))*m_ny*n_points_line,1);
sign = -1.0*sign;
}
......@@ -872,7 +872,7 @@ namespace Nektar
//along y
for(int i = 0; i < m_ny; i++)
{
beta = sign*2*M_PI*(i/2)/m_lhom_y;
beta = -sign*2*M_PI*(i/2)/m_lhom_y;
for(int j = 0; j < m_nz; j++)
{
......@@ -894,7 +894,7 @@ namespace Nektar
//along z
for(int i = 0; i < m_nz; i++)
{
beta = sign*2*M_PI*(i/2)/m_lhom_z;
beta = -sign*2*M_PI*(i/2)/m_lhom_z;
Vmath::Smul(m_ny*n_points_line,beta,tmp1 = temparray + i*m_ny*n_points_line,1,tmp2 = temparray2 + (i-int(sign))*m_ny*n_points_line,1);
sign = -1.0*sign;
}
......
......@@ -240,8 +240,9 @@
<CONDITIONS>
<PARAMETERS>
<P> Lambda = 1 </P>
<P> Lz = 1 </P>
<P> Lambda = 1 </P>
<P> HomModesZ = 8 </P>
<P> LZ = 1 </P>
</PARAMETERS>
<VARIABLES>
......@@ -254,17 +255,17 @@
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="u" VALUE="sin(PI*x)*sin(PI*y)*sin(2*PI*z/Lz)" />
<D VAR="u" VALUE="sin(PI*x)*sin(PI*y)*sin(2*PI*z/LZ)" />
</REGION>
</BOUNDARYCONDITIONS>
<FUNCTION NAME="Forcing">
<E VAR="u"
VALUE="-(Lambda + 2*PI*PI + 4*PI*PI/(Lz*Lz))*sin(PI*x)*sin(PI*y)*sin(2*PI*z/Lz)" />
VALUE="-(Lambda + 2*PI*PI + 4*PI*PI/(LZ*LZ))*sin(PI*x)*sin(PI*y)*sin(2*PI*z/LZ)" />
</FUNCTION>
<FUNCTION NAME="ExactSolution">
<E VAR="u" VALUE="sin(PI*x)*sin(PI*y)*sin(2*PI*z/Lz)" />
<E VAR="u" VALUE="sin(PI*x)*sin(PI*y)*sin(2*PI*z/LZ)" />
</FUNCTION>
</CONDITIONS>
......
......@@ -2,5 +2,5 @@ HDGHelmholtz3DHomo1D
==========================================
helmholtz3D_homo1D_7modes_8nz.xml
L infinity error: 0.0250442
L 2 error : 0.00550986
L 2 error : 0.00389606
----------------------------------------
......@@ -2,5 +2,5 @@ Helmholtz3DHomo1D
==========================================
helmholtz3D_homo1D_7modes_8nz.xml
L infinity error: 0.00918763
L 2 error: 0.00598335
L 2 error: 0.00423087
----------------------------------------
......@@ -33,19 +33,19 @@ L 2 error (variable u) : 1.62207e-08
L inf error (variable u) : 5.66566e-08
----------------------------------------
Test_Helmholtz_3DHomo1D_MVM.xml
L 2 error (variable u) : 3.19231e-05
L 2 error (variable u) : 2.2573e-05
L inf error (variable u) : 2.50889e-05
----------------------------------------
Test_Helmholtz_3DHomo2D_MVM.xml
L 2 error (variable u) : 2.00743e-06
L 2 error (variable u) : 1.00372e-06
L inf error (variable u) : 6.32758e-06
----------------------------------------
Test_Helmholtz_3DHomo1D_FFT.xml
L 2 error (variable u) : 3.19231e-05
L 2 error (variable u) : 2.2573e-05
L inf error (variable u) : 2.50889e-05
----------------------------------------
Test_Helmholtz_3DHomo2D_FFT.xml
L 2 error (variable u) : 2.00743e-06
L 2 error (variable u) : 1.00372e-06
L inf error (variable u) : 6.32758e-06
----------------------------------------
Test_SteadyAdvDiffReact2D_modal.xml
......@@ -133,7 +133,7 @@ L 2 error (variable u) : 1.5458e-08
L inf error (variable u) : 3.64252e-08
----------------------------------------
Test_UnsteadyAdvectionDiffusion_3DHomo1D_MVM.xml
L 2 error (variable u) : 6.99031e-08
L 2 error (variable u) : 4.9429e-08
L inf error (variable u) : 1.53887e-07
----------------------------------------
Test_UnsteadyAdvectionDiffusion_3DHomo2D_MVM.xml
......@@ -141,7 +141,7 @@ L 2 error (variable u) : 5.48207e-09
L inf error (variable u) : 1.39886e-08
----------------------------------------
Test_UnsteadyAdvectionDiffusion_3DHomo1D_FFT.xml
L 2 error (variable u) : 5.99725e-08
L 2 error (variable u) : 4.24069e-08
L inf error (variable u) : 1.02378e-07
----------------------------------------
Test_UnsteadyAdvectionDiffusion_3DHomo2D_FFT.xml
......
......@@ -151,7 +151,7 @@ L 2 error (variable u) : 1.27846e-15
L inf error (variable u) : 6.75498e-15
L 2 error (variable v) : 1.73027e-15
L inf error (variable v) : 7.82707e-15
L 2 error (variable w) : 1.15463e-05
L 2 error (variable w) : 8.16448e-06
L inf error (variable w) : 1.11177e-05
----------------------------------------
Test_2DFlow_lineforcing_bcfromfile.xml
......@@ -207,7 +207,7 @@ L 2 error (variable v) : 1.08091e-05
L inf error (variable v) : 3.95269e-05
L 2 error (variable w) : 1.90154e-05
L inf error (variable w) : 3.94908e-05
L 2 error (variable p) : 8.55847e-05
L 2 error (variable p) : 6.05175e-05
L inf error (variable p) : 0.000513632
----------------------------------------
Test_ChanFlow_m3_SKS.xml
......@@ -225,7 +225,7 @@ L 2 error (variable v) : 1.10749e-05
L inf error (variable v) : 3.96307e-05
L 2 error (variable w) : 1.92288e-05
L inf error (variable w) : 4.03102e-05
L 2 error (variable p) : 8.91694e-05
L 2 error (variable p) : 6.30523e-05
L inf error (variable p) : 0.000565207
----------------------------------------
Test_Hex_channel_m3.xml
......@@ -265,7 +265,7 @@ L 2 error (variable v) : 1.03201e-10
L inf error (variable v) : 1.55026e-09
L 2 error (variable w) : 7.74576e-09
L inf error (variable w) : 2.59194e-08
L 2 error (variable p) : 8.18866e-08
L 2 error (variable p) : 9.34515e-08
L inf error (variable p) : 3.42898e-07
----------------------------------------
Test_Prism_channel_m6.xml
......
......@@ -133,7 +133,6 @@ namespace Nektar
int VelDim = vel_loc.num_elements();
int nqtot = pFields[0]->GetTotPoints();
Array<OneD, Array<OneD, NekDouble> > velocity(VelDim);
Array<OneD, NekDouble > Deriv;
ASSERTL1(nConvectiveFields == pInarray.num_elements(),"Number of convective fields and Inarray are not compatible");
......
......@@ -81,41 +81,69 @@ namespace Nektar
// ToDo: here we should add a check that V has right dimension
int nPointsTot = pFields[0]->GetNpoints();
Array<OneD, NekDouble> gradV0,gradV1,gradV2, gradVV0, gradVV1, gradVV2, Up;
Array<OneD, NekDouble> gradV0,gradV1,gradV2, tmp, Up;
gradV0 = Array<OneD, NekDouble> (nPointsTot);
gradVV0 = Array<OneD, NekDouble> (nPointsTot);
tmp = Array<OneD, NekDouble> (nPointsTot);
/*//////////////////////////////////////////////////////
//FILE *pFile0;
if(pVelocityComponent==0)
{
pFile0= fopen("u_U_V_W.txt", "w");
}
else if(pVelocityComponent==1)
{
pFile0= fopen("v_U_V_W.txt", "w");
}
else
{
pFile0= fopen("w_U_V_W.txt", "w");
}
for(int k=0; k < nPointsTot ; ++k)
{
fprintf(pFile0, "%i %10.20lf\t %10.20lf\t %10.20lf\t %10.20lf\t \n ", k, pU[k], pV[0][k], pV[1][k], pV[2][k]);
}
fclose(pFile0);
Array<OneD, NekDouble> dudx,dudy,dudz,duudx,duvdy,duwdz;
dudx = Array<OneD, NekDouble> (nPointsTot);
dudy = Array<OneD, NekDouble> (nPointsTot);
dudz = Array<OneD, NekDouble> (nPointsTot);
duudx = Array<OneD, NekDouble> (nPointsTot);
duvdy = Array<OneD, NekDouble> (nPointsTot);
duwdz = Array<OneD, NekDouble> (nPointsTot);
//////////////////////////////////////////////////////*/
// Evaluate V\cdot Grad(u)
switch(ndim)
{
case 1:
pFields[0]->PhysDeriv(pU,gradV0);
Vmath::Vmul(nPointsTot,gradV0,1,pV[0],1,pOutarray,1);
pFields[0]->PhysDeriv(pU,gradV0);
Vmath::Vmul(nPointsTot,gradV0,1,pV[0],1,pOutarray,1);
Vmath::Vmul(nPointsTot,pU,1,pV[0],1,gradV0,1);
pFields[0]->PhysDeriv(gradV0,gradVV0);
Vmath::Vadd(nPointsTot,gradVV0,1,pOutarray,1,pOutarray,1);
pFields[0]->PhysDeriv(gradV0,tmp);
Vmath::Vadd(nPointsTot,tmp,1,pOutarray,1,pOutarray,1);
Vmath::Smul(nPointsTot,0.5,pOutarray,1,pOutarray,1);
break;
case 2:
gradV1 = Array<OneD, NekDouble> (nPointsTot);
gradVV1 = Array<OneD, NekDouble> (nPointsTot);
pFields[0]->PhysDeriv(pU,gradV0,gradV1);
Vmath::Vmul (nPointsTot,gradV0,1,pV[0],1,pOutarray,1);
Vmath::Vvtvp(nPointsTot,gradV1,1,pV[1],1,pOutarray,1,pOutarray,1);
gradV1 = Array<OneD, NekDouble> (nPointsTot);
pFields[0]->PhysDeriv(pU,gradV0,gradV1);
Vmath::Vmul (nPointsTot,gradV0,1,pV[0],1,pOutarray,1);
Vmath::Vvtvp(nPointsTot,gradV1,1,pV[1],1,pOutarray,1,pOutarray,1);
Vmath::Vmul(nPointsTot,pU,1,pV[0],1,gradV0,1);
Vmath::Vmul(nPointsTot,pU,1,pV[1],1,gradV1,1);
pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],gradV0,gradVV0);
pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],gradV1,gradVV1);
Vmath::Vadd(nPointsTot,gradVV0,1,pOutarray,1,pOutarray,1);
Vmath::Vadd(nPointsTot,gradVV1,1,pOutarray,1,pOutarray,1);
pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],gradV0,tmp);
Vmath::Vadd(nPointsTot,tmp,1,pOutarray,1,pOutarray,1);
pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],gradV1,tmp);
Vmath::Vadd(nPointsTot,tmp,1,pOutarray,1,pOutarray,1);