Commit 7f1bb27e authored by Gabriele Rocco's avatar Gabriele Rocco
Browse files

Half-Mode method


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@3851 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent 26c572e1
......@@ -49,7 +49,7 @@ int main(int argc, char *argv[])
const LibUtilities::PointsKey PkeyZ(nzpoints,LibUtilities::eFourierSingleModeSpaced);
const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourierSingleMode,nzpoints,PkeyZ);
const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourierHalfModeRe,nzpoints,PkeyZ);
Exp_u = MemoryManager<MultiRegions::ContField3DHomogeneous1D>::AllocateSharedPtr(vSession,BkeyZ,lz,useFFT,deal,graph2D,vSession->GetVariable(0));
Exp_v = MemoryManager<MultiRegions::ContField3DHomogeneous1D>::AllocateSharedPtr(vSession,BkeyZ,lz,useFFT,deal,graph2D,vSession->GetVariable(1));
......@@ -112,8 +112,10 @@ int main(int argc, char *argv[])
//----------------------------------------------
//Taking derivative and printing the error
cout << "Deriv u" << endl;
Exp_u->PhysDeriv(Exp_u->GetPhys(),Exp_u->UpdatePhys(),dump,dump,false);
cout << "Deriv u done" << endl;
cout << "L infinity error: " << Exp_u->Linf(dudx) << endl;
cout << "L 2 error : " << Exp_u->L2 (dudx) << endl;
......
......@@ -3,6 +3,10 @@ SET(LinkLibraries
optimized LibUtilities debug LibUtilities-g
optimized ${TINYXML_LIB} debug ${TINYXML_LIB}
)
SET(StdProject0DSource StdProject0D.cpp)
ADD_NEKTAR_EXECUTABLE(StdProject0D demos StdProject0DSource)
TARGET_LINK_LIBRARIES(StdProject0D ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(StdProject0D)
SET(StdProject1DSource StdProject1D.cpp)
ADD_NEKTAR_EXECUTABLE(StdProject1D demos StdProject1DSource)
......
#include <cstdio>
#include <cstdlib>
#include "StdRegions/StdSegExp.h"
#include "StdRegions/StdPointExp.h"
#include "StdRegions/StdRegions.hpp"
#include "LibUtilities/Foundations/Foundations.hpp"
using namespace Nektar;
// This routine projects a polynomial or trigonmetric functions which
// has energy in all mdoes of the expansions and report an error.
int main(int argc, char *argv[])
{
int i,j;
int order, nq;
LibUtilities::PointsType Qtype;
LibUtilities::BasisType btype;
StdRegions::StdExpansion0D *E;
Array<OneD,NekDouble> sol;
if(argc != 4)
{
fprintf(stderr,"Usage: Project1D Type order nq \n");
fprintf(stderr,"Where type is an integer value which "
"dictates the basis as:\n");
fprintf(stderr,"\t Ortho_A = 1\n");
fprintf(stderr,"\t Modified_A = 4\n");
fprintf(stderr,"\t Fourier = 7\n");
fprintf(stderr,"\t Lagrange = 8\n");
fprintf(stderr,"\t Legendre = 9\n");
fprintf(stderr,"\t Chebyshev = 10\n");
fprintf(stderr,"\t Monomial = 11\n");
fprintf(stderr,"\t FourierSingleMode = 12\n");
fprintf(stderr,"\t FourierHalfModeRe = 13\n");
fprintf(stderr,"\t FourierHalfModeIm = 14\n");
fprintf(stderr,"Note type = 1,2,4,5 are for higher dimensional basis\n");
exit(1);
}
btype = (LibUtilities::BasisType) atoi(argv[1]);
if(btype == LibUtilities::eNoBasisType)
{
NEKERROR(ErrorUtil::efatal,
"No Basis Type requested");
}
// Check to see that only 1D Expansions are used
if((btype == LibUtilities::eOrtho_B)||(btype == LibUtilities::eOrtho_B)||
(btype == LibUtilities::eModified_B)||(btype == LibUtilities::eModified_C))
{
NEKERROR(ErrorUtil::efatal,
"This basis is for 2 or 3D expansions");
}
order = atoi(argv[2]);
nq = atoi(argv[3]);
sol = Array<OneD, NekDouble>(nq);
if(btype== LibUtilities::eFourier)
{
Qtype = LibUtilities::eFourierEvenlySpaced;
}
else if(btype== LibUtilities::eFourierSingleMode)
{
Qtype = LibUtilities::eFourierSingleModeSpaced;
}
else if(btype== LibUtilities::eFourierHalfModeRe || btype== LibUtilities::eFourierHalfModeIm )
{
//Qtype = LibUtilities::eFourierHalfModeSpaced;
Qtype = LibUtilities::eFourierSingleModeSpaced;
}
else
{
Qtype = LibUtilities::eGaussLobattoLegendre;
}
//-----------------------------------------------
// Define a point expansion based on basis definition
const LibUtilities::PointsKey Pkey(nq,Qtype);
const LibUtilities::BasisKey Bkey(btype,order,Pkey);
cout << "StdSegExp" << endl;
E = new StdRegions::StdPointExp(Bkey);
//-----------------------------------------------
//----------------------------------------------
// Define solution to be projected
Array<OneD,NekDouble> z = Array<OneD,NekDouble>(nq);
E->GetCoords(z);
if(btype== LibUtilities::eFourier)
{
for(i = 0; i < nq; ++i)
{
sol[i] = 1.0;
for(j = 0; j < order/2-1; ++j)
{
sol[i] += sin(j*M_PI*z[i]) + cos(j*M_PI*z[i]);
}
}
}
else if(btype== LibUtilities::eFourierSingleMode)
{
for(i = 0; i < nq; ++i)
{
sol[i] = (0.45*sin(M_PI*z[i]) + 0.25*cos(M_PI*z[i]));
}
}
else if(btype== LibUtilities::eFourierHalfModeRe)
{
for(i = 0; i < nq; ++i)
{
sol[i] = (0.45*cos(M_PI*z[i]));
}
}
else if(btype== LibUtilities::eFourierHalfModeIm)
{
for(i = 0; i < nq; ++i)
{
sol[i] = (0.25*sin(M_PI*z[i]));
}
}
else {
for(i = 0; i < nq; ++i)
{
sol[i] = 0.0;
for(j = 0; j < order; ++j)
{
sol[i] += pow(z[i],j);
}
}
}
//---------------------------------------------
// Project onto Expansion
E->FwdTrans(sol,E->UpdateCoeffs());
//---------------------------------------------
cout << "z[0]=" << z[0] << " m_coeffs[0]=" << E->GetCoeffs()[0] << endl;
//cout << "z[1]=" << z[1] << " m_coeffs[1]=" << E->GetCoeffs()[1] << endl;
//-------------------------------------------
// Backward Transform Solution to get projected values
E->BwdTrans(E->GetCoeffs(),E->UpdatePhys());
//-------------------------------------------
cout << "z[0]=" << z[0] << " m_phys[0]=" << E->GetPhys()[0] << endl;
//cout << "z[1]=" << z[1] << " m_phys[1]=" << E->GetPhys()[1] << endl;
//--------------------------------------------
// Calculate L_inf error
cout << "L infinity error: " << E->Linf(sol) << endl;
//--------------------------------------------
//-------------------------------------------
// Evaulate solution at mid point and print error
if(btype != LibUtilities::eFourierSingleMode || (btype != LibUtilities::eFourierHalfModeRe) || (btype != LibUtilities::eFourierHalfModeIm) )
{
if(btype != LibUtilities::eFourier)
{
sol[0] = 1;
}
else
{
sol[0] = order/2;
}
Array<OneD,NekDouble> x = Array<OneD, NekDouble>(1);
x[0] = 0;
NekDouble nsol = E->PhysEvaluate(x);
cout << "error at x = 0: " << nsol - sol[0] << endl;
}
//-------------------------------------------
return 0;
}
......@@ -85,6 +85,12 @@ namespace Nektar
m_K[1] = 1;
}
if(HomoBasis0.GetBasisType() == LibUtilities::eFourierHalfModeRe ||
HomoBasis0.GetBasisType() == LibUtilities::eFourierHalfModeIm)
{
m_K[0]=1;
}
//=================================================================================
}
......
......@@ -605,7 +605,22 @@ namespace Nektar
}
break;
//Fourier Real Half Mode
case eFourierHalfModeRe:
m_bdata[0] = cos(M_PI*z[0]);
m_dbdata[0] = -M_PI*sin(M_PI*z[0]);
break;
//Fourier Imaginary Half Mode
case eFourierHalfModeIm:
m_bdata[0] = sin(M_PI*z[0]);
m_dbdata[0] = M_PI*cos(M_PI*z[0]);
break;
case eChebyshev:
mode = m_bdata.data();
......
......@@ -117,6 +117,8 @@ namespace Nektar
case eDG_HU_Left:
case eDG_HU_Right:
case eFourierSingleMode:
case eFourierHalfModeRe:
case eFourierHalfModeIm:
value = m_nummodes;
break;
......
......@@ -57,6 +57,8 @@ namespace Nektar
eChebyshev, //!< Chebyshev Polynomials \f$ T_p(z_i) = P^{-1/2,-1/2}_p(z_i)\f$
eMonomial, //!< Monomial polynomials \f$ L_p(z_i) = z_i^{p}\f$
eFourierSingleMode, //!< Fourier ModifiedExpansion with just the first mode \f$ \exp(i \pi z_i)\f$
eFourierHalfModeRe, //!< Fourier Modified expansions with just the real part of the first mode \f$ Re[\exp(i \pi z_i)]\f$
eFourierHalfModeIm, //!< Fourier Modified expansions with just the imaginary part of the first mode \f$ Im[\exp(i \pi z_i)]\f$
eG_Lagrange, //!< Lagrange Polynomials using the Gauss points \f$ h_p(z_i) \f$
eDG_DG_Left, //!< Derivative of the left correction function for DG FR \f$ dGL_{p}(z_i) \f$
eDG_DG_Right, //!< Derivative of the Right correction function for DG FR \f$ dGR_{p}(z_i) \f$
......@@ -82,6 +84,8 @@ namespace Nektar
"Chebyshev",
"Monomial",
"FourierSingleMode",
"FourierHalfModeRe",
"FourierHalfModeIm",
"G_Lagrange",
"DG_DG_Left",
"DG_DG_Right",
......
......@@ -54,15 +54,15 @@ namespace Nektar
if(npts==1)
{
m_points[0][0] = 0.0;
m_points[0][0] = 0.25;
}
else
{
ASSERTL0(npts==2, "Fourier points for single mode analysis should be 2");
m_points[0][0] = -1.0;
m_points[0][1] = -0.5;
m_points[0][0] = 0.0;
m_points[0][1] = 0.5;
}
}
......
......@@ -179,12 +179,13 @@ namespace Nektar
const StdRegions::VarCoeffMap &varcoeff,
const Array<OneD, const NekDouble> &dirForcing)
{
int n;
int cnt = 0;
int cnt1 = 0;
NekDouble beta;
StdRegions::ConstFactorMap new_factors;
Array<OneD, NekDouble> e_out;
Array<OneD, NekDouble> fce(inarray.num_elements());
......
......@@ -490,26 +490,55 @@ namespace Nektar
BlkMatrix = MemoryManager<DNekBlkMat>
::AllocateSharedPtr(nrows,ncols,blkmatStorage);
StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
{
StdRegions::StdMatrixKey matkey(StdRegions::eFwdTrans,
StdSeg.DetExpansionType(),
StdSeg);
loc_mat = StdSeg.GetStdMatrix(matkey);
//Half Mode
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
{
StdRegions::StdPointExp StdPoint(m_homogeneousBasis->GetBasisKey());
}
else
{
StdRegions::StdMatrixKey matkey(StdRegions::eBwdTrans,
StdSeg.DetExpansionType(),
StdSeg);
loc_mat = StdSeg.GetStdMatrix(matkey);
if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
{
StdRegions::StdMatrixKey matkey(StdRegions::eFwdTrans,
StdPoint.DetExpansionType(),
StdPoint);
loc_mat = StdPoint.GetStdMatrix(matkey);
}
else
{
StdRegions::StdMatrixKey matkey(StdRegions::eBwdTrans,
StdPoint.DetExpansionType(),
StdPoint);
loc_mat = StdPoint.GetStdMatrix(matkey);
}
}
//other cases
else
{
StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
}
if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
{
StdRegions::StdMatrixKey matkey(StdRegions::eFwdTrans,
StdSeg.DetExpansionType(),
StdSeg);
loc_mat = StdSeg.GetStdMatrix(matkey);
}
else
{
StdRegions::StdMatrixKey matkey(StdRegions::eBwdTrans,
StdSeg.DetExpansionType(),
StdSeg);
loc_mat = StdSeg.GetStdMatrix(matkey);
}
}
// set up array of block matrices.
for(int i = 0; i < num_trans_per_proc; ++i)
......@@ -776,7 +805,8 @@ namespace Nektar
m_planes[i]->PhysDeriv(tmp1 = inarray + i*nP_pts ,tmp2 = out_d0 + i*nP_pts , tmp3 = out_d1 + i*nP_pts );
}
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierSingleMode)
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierSingleMode ||
m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
{
if(m_WaveSpace)
{
......@@ -790,13 +820,26 @@ namespace Nektar
NekDouble sign = -1.0;
NekDouble beta;
for(int i = 0; i < m_planes.num_elements(); i++)
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe ||
m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
{
beta = sign*2*M_PI*(m_transposition->GetK(i))/m_lhom;
beta = sign*2*M_PI*(m_transposition->GetK(0))/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;
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)
......@@ -860,7 +903,8 @@ namespace Nektar
}
else
{
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierSingleMode)
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierSingleMode ||
m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
{
if(m_WaveSpace)
{
......@@ -873,14 +917,25 @@ namespace Nektar
NekDouble sign = -1.0;
NekDouble beta;
for(int i = 0; i < m_planes.num_elements(); i++)
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe ||
m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
{
beta = sign*2*M_PI*(m_transposition->GetK(i))/m_lhom;
beta = 2*M_PI*(m_transposition->GetK(0))/m_lhom;
Vmath::Smul(nP_pts,beta,tmp1 = temparray + i*nP_pts,1,tmp2 = outarray + (i-int(sign))*nP_pts,1);
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;
sign = -1.0*sign;
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)
{
......
......@@ -64,8 +64,10 @@ namespace Nektar
*/
void DriverStandard::v_InitObject(ostream &out)
{
cout << "InitObject" << endl;
Driver::v_InitObject(out);
}
cout << "InitObject done " << endl;
}
void DriverStandard::v_Execute(ostream &out)
......
......@@ -156,6 +156,11 @@ namespace Nektar
m_npointsZ=2;
}
else if(m_session->GetSolverInfo("SingleMode")=="HalfMode")
{
m_npointsZ=1;
}
else
{
ASSERTL0(false, "SolverInfo Single Mode not valid");
......@@ -243,6 +248,7 @@ namespace Nektar
int nvariables = m_session->GetVariables().size();
bool DeclareCoeffPhysArrays = true;
m_fields = Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
m_spacedim = m_graph->GetSpaceDimension()+m_HomoDirec;
m_expdim = m_graph->GetMeshDimension();
......@@ -298,7 +304,33 @@ namespace Nektar
m_fields[i]->SetWaveSpace(false);
}
}
}
//Half mode stability analysis
else if(m_session->DefinesSolverInfo("SingleMode")&& m_session->GetSolverInfo("SingleMode")=="HalfMode")
{
const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierSingleModeSpaced);
const LibUtilities::BasisKey BkeyZR(LibUtilities::eFourierHalfModeRe,m_npointsZ,PkeyZ);
const LibUtilities::BasisKey BkeyZI(LibUtilities::eFourierHalfModeIm,m_npointsZ,PkeyZ);
for(i = 0 ; i < m_fields.num_elements(); i++)
{
if(i==m_fields.num_elements()-2)
{
m_fields[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
::AllocateSharedPtr(m_session,BkeyZI,m_LhomZ,m_useFFT,m_dealiasing,m_graph,m_session->GetVariable(i),m_checkIfSystemSingular[i]);
}
m_fields[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
::AllocateSharedPtr(m_session,BkeyZR,m_LhomZ,m_useFFT,m_dealiasing,m_graph,m_session->GetVariable(i),m_checkIfSystemSingular[i]);
//necessary to perform to HomoBwdTrans for the fields (are complex differently from the base flow)
m_fields[i]->SetWaveSpace(false);
}
}
//normal homogeneous 1D
else
{
......@@ -563,7 +595,8 @@ namespace Nektar
}
EvaluateFunction(fieldStr, m_forces, "BodyForce");
if(m_session->DefinesSolverInfo("SingleMode")&& m_session->GetSolverInfo("SingleMode")=="ModifiedBasis")
if(m_session->DefinesSolverInfo("SingleMode")&&
(m_session->GetSolverInfo("SingleMode")=="ModifiedBasis" ||m_session->GetSolverInfo("SingleMode")=="HalfMode"))
{
for(int i=0; i< v_GetForceDimension(); ++i)
{
......@@ -584,6 +617,7 @@ namespace Nektar
// Zero all physical fields initially.
ZeroPhysFields();
}
/**
......
......@@ -42,19 +42,90 @@ namespace Nektar
namespace StdRegions
{
StdExpansion0D::StdExpansion0D()
{
}
StdExpansion0D::StdExpansion0D(const StdExpansion0D &T):StdExpansion(T)
{
}
StdExpansion0D::~StdExpansion0D()
{
}
StdExpansion0D::StdExpansion0D()
{
}
StdExpansion0D::StdExpansion0D(int numcoeffs, const LibUtilities::BasisKey &Ba):
StdExpansion(numcoeffs,1,Ba)
{
}
StdExpansion0D::StdExpansion0D(const StdExpansion0D &T):StdExpansion(T)
{
}