Commit e5b0b504 authored by Gabriele Rocco's avatar Gabriele Rocco
Browse files

Cleaning of the stability solver; fixed HOPBCs with MultipleModes; Extension...

Cleaning of the stability solver; fixed HOPBCs with MultipleModes; Extension to Adjoint and TG; post-processing tools


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@3856 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent 416de0bc
......@@ -93,8 +93,7 @@ int main(int argc, char *argv[])
// 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);
E = new StdRegions::StdPointExp(Bkey);
//-----------------------------------------------
......@@ -152,40 +151,20 @@ int main(int argc, char *argv[])
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;
cout << "z[0]= " << z[0] << " coefficient= " << E->GetCoeffs()[0] << 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;
cout << "z[0]= " << z[0] << " physical value= " << E->GetPhys()[0] << 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;
}
//-------------------------------------------
......
......@@ -820,7 +820,7 @@ namespace Nektar
NekDouble sign = -1.0;
NekDouble beta;
//Half Mode
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe ||
m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
{
......@@ -918,6 +918,7 @@ namespace Nektar
NekDouble sign = -1.0;
NekDouble beta;
//HalfMode
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe ||
m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
{
......
......@@ -72,7 +72,9 @@ namespace Nektar
m_period = m_session->GetParameter("TimeStep")* m_session->GetParameter("NumSteps");
m_nfields = m_equ[0]->UpdateFields().num_elements() - 1;
if(m_session->DefinesSolverInfo("SingleMode") && m_session->GetSolverInfo("SingleMode")=="ModifiedBasis")
//if(m_session->DefinesSolverInfo("SingleMode") && m_session->GetSolverInfo("SingleMode")=="ModifiedBasis")
if(m_session->DefinesSolverInfo("ModeType") &&
(m_session->GetSolverInfo("ModeType")=="SingleMode"|| m_session->GetSolverInfo("ModeType")=="HalfMode") )
{
for(int i = 0; i < m_nfields; ++i)
{
......
......@@ -109,10 +109,18 @@ namespace Nektar
m_spacedim = m_graph->GetSpaceDimension();
// Setting parameteres for homogenous problems
m_HomoDirec = 0;
m_useFFT = false;
m_dealiasing = false;
m_SingleMode =false;
m_HomoDirec = 0;
m_useFFT = false;
m_dealiasing = false;
m_SingleMode = false;
m_HalfMode = false;
m_MultipleModes = false;
m_session->MatchSolverInfo("ModeType","SingleMode",m_SingleMode,false);
m_session->MatchSolverInfo("ModeType","HalfMode",m_HalfMode,false);
m_session->MatchSolverInfo("ModeType","MultipleModes",m_MultipleModes,false);
m_HomogeneousType = eNotHomogeneous;
......@@ -125,49 +133,33 @@ namespace Nektar
(HomoStr == "1D")||(HomoStr == "Homo1D"))
{
m_HomogeneousType = eHomogeneous1D;
//m_session->LoadParameter("HomModesZ",m_npointsZ);
m_session->LoadParameter("LZ",m_LhomZ);
m_HomoDirec = 1;
if(m_session->DefinesSolverInfo("SingleMode") &&
m_session->GetSolverInfo("solvertype") != "CoupledLinearisedNS")
{
if(m_session->GetSolverInfo("SingleMode")=="SpecifiedMode")
{
m_SingleMode=true;
if(m_session->DefinesParameter("NumMode"))
{
//read mode from session file
m_NumMode=m_session->GetParameter("NumMode");
}
else
{
//first mode by default
m_NumMode=1;
}
//number of plane to create in case of single modes analysis.
m_npointsZ=2+2*m_NumMode;
}
else if(m_session->GetSolverInfo("SingleMode")=="ModifiedBasis")
{
m_npointsZ=2;
}
else if(m_session->GetSolverInfo("SingleMode")=="HalfMode")
//Stability Analysis flags
if(m_session->DefinesSolverInfo("ModeType"))
{
if(m_SingleMode)
{
m_npointsZ=2;
}
else if(m_HalfMode)
{
m_npointsZ=1;
}
else
{
ASSERTL0(false, "SolverInfo Single Mode not valid");
}
}
else if(m_MultipleModes)
{
m_npointsZ = m_session->GetParameter("HomModesZ");
}
else
{
ASSERTL0(false, "SolverInfo ModeType not valid");
}
}
else
{
m_npointsZ = m_session->GetParameter("HomModesZ");
......@@ -288,8 +280,9 @@ namespace Nektar
{
if(m_HomogeneousType == eHomogeneous1D)
{
//Modified basis for stability analysis
if(m_session->DefinesSolverInfo("SingleMode")&& m_session->GetSolverInfo("SingleMode")=="ModifiedBasis")
//FourierSingleMode basis for stability analysis
if(m_SingleMode)
{
const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierSingleModeSpaced);
......@@ -306,7 +299,7 @@ namespace Nektar
}
}
//Half mode stability analysis
else if(m_session->DefinesSolverInfo("SingleMode")&& m_session->GetSolverInfo("SingleMode")=="HalfMode")
else if(m_HalfMode)
{
const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierSingleModeSpaced);
......@@ -344,7 +337,7 @@ namespace Nektar
::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_dealiasing,m_graph,m_session->GetVariable(i),m_checkIfSystemSingular[i]);
if(m_session->DefinesSolverInfo("SingleMode")&& m_session->GetSolverInfo("SingleMode")=="SpecifiedMode")
if(m_MultipleModes)
{
m_fields[i]->SetWaveSpace(false);
}
......@@ -595,8 +588,7 @@ namespace Nektar
}
EvaluateFunction(fieldStr, m_forces, "BodyForce");
if(m_session->DefinesSolverInfo("SingleMode")&&
(m_session->GetSolverInfo("SingleMode")=="ModifiedBasis" ||m_session->GetSolverInfo("SingleMode")=="HalfMode"))
if(m_SingleMode || m_HalfMode)
{
for(int i=0; i< v_GetForceDimension(); ++i)
{
......@@ -1944,7 +1936,8 @@ namespace Nektar
out << "\tUsing MVM " << endl;
}
if(m_SingleMode==true)
//if(m_SingleMode==true)
if(m_MultipleModes==true)
{
out << "\tSelected Mode : " << m_NumMode << endl;
......
......@@ -405,7 +405,12 @@ namespace Nektar
int m_checksteps; ///< Number of steps between checkpoints
int m_spacedim; ///< Spatial dimension (> expansion dim)
int m_expdim; ///< Dimension of the expansion
bool m_SingleModeBasis; ///< Flag for the SingleMode Basis
bool m_SingleModeBasis; ///< Flag for the SingleMode Basis
bool m_SingleMode; ///< flag to determine if use single mode or not
bool m_HalfMode; ///< flag to determine if use half mode or not
bool m_MultipleModes; ///< flag to determine if use multiple mode or not
/// Type of projection, i.e. Galerkin or DG.
enum MultiRegions::ProjectionType m_projectionType;
......@@ -437,7 +442,6 @@ namespace Nektar
bool m_dealiasing; ///< flag to determine if use dealising or not
bool m_SingleMode; ///< flag to determine if use single mode or not
enum HomogeneousType m_HomogeneousType;
......
......@@ -2251,7 +2251,9 @@ namespace Nektar
(strcmp(LibUtilities::BasisTypeMap[basis[j]], "Modified_C") == 0) ||
(strcmp(LibUtilities::BasisTypeMap[basis[j]], "GLL_Lagrange") == 0) ||
(strcmp(LibUtilities::BasisTypeMap[basis[j]], "Fourier") == 0) ||
(strcmp(LibUtilities::BasisTypeMap[basis[j]], "FourierSingleMode") == 0))
(strcmp(LibUtilities::BasisTypeMap[basis[j]], "FourierSingleMode") == 0)||
(strcmp(LibUtilities::BasisTypeMap[basis[j]], "FourierHalfModeRe") == 0) ||
(strcmp(LibUtilities::BasisTypeMap[basis[j]], "FourierHalfModeIm") == 0))
{
check++;
}
......@@ -3062,6 +3064,80 @@ namespace Nektar
}
}
break;
case eFourierHalfModeRe:
{
switch (shape)
{
case eSegment:
{
const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eFourierSingleModeSpaced);
LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe,nummodes,pkey);
returnval.push_back(bkey);
}
break;
case eQuadrilateral:
{
const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eFourierSingleModeSpaced);
LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe,nummodes,pkey);
returnval.push_back(bkey);
returnval.push_back(bkey);
}
break;
case eHexahedron:
{
const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eFourierSingleModeSpaced);
LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeRe,nummodes,pkey);
returnval.push_back(bkey);
returnval.push_back(bkey);
returnval.push_back(bkey);
}
break;
default:
{
ASSERTL0(false,"Expansion not defined in switch for this shape");
}
break;
}
}
break;
case eFourierHalfModeIm:
{
switch (shape)
{
case eSegment:
{
const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eFourierSingleModeSpaced);
LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm,nummodes,pkey);
returnval.push_back(bkey);
}
break;
case eQuadrilateral:
{
const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eFourierSingleModeSpaced);
LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm,nummodes,pkey);
returnval.push_back(bkey);
returnval.push_back(bkey);
}
break;
case eHexahedron:
{
const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eFourierSingleModeSpaced);
LibUtilities::BasisKey bkey(LibUtilities::eFourierHalfModeIm,nummodes,pkey);
returnval.push_back(bkey);
returnval.push_back(bkey);
returnval.push_back(bkey);
}
break;
default:
{
ASSERTL0(false,"Expansion not defined in switch for this shape");
}
break;
}
}
break;
case eChebyshev:
{
......@@ -3233,6 +3309,23 @@ namespace Nektar
returnval.push_back(bkey1);
}
break;
case eFourierHalfModeRe:
{
const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierSingleModeSpaced);
LibUtilities::BasisKey bkey1(LibUtilities::eFourierHalfModeRe,nummodes_x,pkey1);
returnval.push_back(bkey1);
}
break;
case eFourierHalfModeIm:
{
const LibUtilities::PointsKey pkey1(nummodes_x,LibUtilities::eFourierSingleModeSpaced);
LibUtilities::BasisKey bkey1(LibUtilities::eFourierHalfModeIm,nummodes_x,pkey1);
returnval.push_back(bkey1);
}
break;
case eChebyshev:
{
......@@ -3241,6 +3334,8 @@ namespace Nektar
returnval.push_back(bkey1);
}
break;
default:
{
......@@ -3269,6 +3364,22 @@ namespace Nektar
}
break;
case eFourierHalfModeRe:
{
const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierSingleModeSpaced);
LibUtilities::BasisKey bkey2(LibUtilities::eFourierHalfModeRe,nummodes_y,pkey2);
returnval.push_back(bkey2);
}
break;
case eFourierHalfModeIm:
{
const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eFourierSingleModeSpaced);
LibUtilities::BasisKey bkey2(LibUtilities::eFourierHalfModeIm,nummodes_y,pkey2);
returnval.push_back(bkey2);
}
break;
case eChebyshev:
{
const LibUtilities::PointsKey pkey2(nummodes_y,LibUtilities::eGaussGaussChebyshev);
......@@ -3301,6 +3412,22 @@ namespace Nektar
returnval.push_back(bkey3);
}
break;
case eFourierHalfModeRe:
{
const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierSingleModeSpaced);
LibUtilities::BasisKey bkey3(LibUtilities::eFourierHalfModeRe,nummodes_z,pkey3);
returnval.push_back(bkey3);
}
break;
case eFourierHalfModeIm:
{
const LibUtilities::PointsKey pkey3(nummodes_z,LibUtilities::eFourierSingleModeSpaced);
LibUtilities::BasisKey bkey3(LibUtilities::eFourierHalfModeIm,nummodes_z,pkey3);
returnval.push_back(bkey3);
}
break;
case eChebyshev:
{
......
......@@ -69,6 +69,8 @@ namespace Nektar
eGLL_Lagrange_SEM,
eFourier,
eFourierSingleMode,
eFourierHalfModeRe,
eFourierHalfModeIm,
eChebyshev,
eFourierChebyshev,
eChebyshevFourier,
......@@ -87,6 +89,8 @@ namespace Nektar
"GLL_LAGRANGE_SEM",
"FOURIER",
"FOURIERSINGLEMODE",
"FOURIERHALFMODERE",
"FOURIERHALFMODEIM",
"CHEBYSHEV",
"FOURIER-CHEBYSHEV",
"CHEBYSHEV-FOURIER",
......
......@@ -69,6 +69,7 @@ int main(int argc, char* argv[])
{
quiet = true;
}
// 1D Projection Tests
Execute("StdProject1D", "1 6 7","Seg Ortho. Basis, P=6, Q=7");
......
......@@ -83,6 +83,22 @@ namespace Nektar
protected:
//Storage of the base flow
Array<OneD, MultiRegions::ExpListSharedPtr> m_base;
//number of slices
int m_slices;
//period length
NekDouble m_period;
//interpolation vector
Array<OneD, Array<OneD, NekDouble> > m_interp;
//auxiliary variables
LibUtilities::NektarFFTSharedPtr m_FFT;
Array<OneD,NekDouble> m_tmpIN;
Array<OneD,NekDouble> m_tmpOUT;
bool m_useFFTW;
bool m_SingleModeBasis; ///< Flag for the SingleMode Basis
bool m_SingleMode; ///< flag to determine if use single mode or not
bool m_HalfMode; ///< flag to determine if use half mode or not
bool m_MultipleModes; ///< flag to determine if use multiple mode or not
AdjointAdvection(
const LibUtilities::SessionReaderSharedPtr& pSession,
......@@ -128,6 +144,36 @@ namespace Nektar
int pVelocityComponent,
NekDouble m_time,
Array<OneD, NekDouble> &pWk);
///Parameter for homogeneous expansions
enum HomogeneousType
{
eHomogeneous1D,
eHomogeneous2D,
eHomogeneous3D,
eNotHomogeneous
};
bool m_useFFT; ///< flag to determine if use or not the FFT for transformations
enum HomogeneousType m_HomogeneousType;
NekDouble m_LhomX; ///< physical length in X direction (if homogeneous)
NekDouble m_LhomY; ///< physical length in Y direction (if homogeneous)
NekDouble m_LhomZ; ///< physical length in Z direction (if homogeneous)
int m_npointsX; ///< number of points in X direction (if homogeneous)
int m_npointsY; ///< number of points in Y direction (if homogeneous)
int m_npointsZ; ///< number of points in Z direction (if homogeneous)
int m_HomoDirec; ///< number of homogenous directions
int m_NumMode; ///< Mode to use in case of single mode analysis
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions;
};
......
......@@ -78,47 +78,38 @@ namespace Nektar
(HomoStr == "1D")||(HomoStr == "Homo1D"))
{
m_HomogeneousType = eHomogeneous1D;
// m_npointsZ = m_session->GetParameter("HomModesZ");
m_LhomZ = m_session->GetParameter("LZ");
m_HomoDirec = 1;
m_SingleMode =false;
if(m_session->DefinesSolverInfo("SingleMode"))
m_HalfMode =false;
m_MultipleModes =false;
m_session->MatchSolverInfo("ModeType","SingleMode",m_SingleMode,false);
m_session->MatchSolverInfo("ModeType","HalfMode",m_HalfMode,false);
m_session->MatchSolverInfo("ModeType","MultipleModes",m_MultipleModes,false);
if(m_session->DefinesSolverInfo("ModeType"))
{
if(m_session->GetSolverInfo("SingleMode")=="SpecifiedMode")
{
m_SingleMode=true;
if(m_session->DefinesParameter("NumMode"))
{
//read mode from session file
m_NumMode=m_session->GetParameter("NumMode");
}
else
{
//first mode by default
m_NumMode=1;
}
//number of plane to create in case of single modes analysis.
m_npointsZ=2+2*m_NumMode;
}
else if(m_session->GetSolverInfo("SingleMode")=="ModifiedBasis")
if(m_SingleMode)
{
m_npointsZ=2;
}
else if(m_session->GetSolverInfo("SingleMode")=="HalfMode")
else if(m_HalfMode)
{
m_npointsZ=1;
}
else
else if(m_MultipleModes)
{
ASSERTL0(false, "SolverInfo Single Mode not valid");
m_npointsZ = m_session->GetParameter("HomModesZ");
}
else
{
ASSERTL0(false, "SolverInfo ModeType not valid");
}
}
else
{
......@@ -345,7 +336,6 @@ namespace Nektar
int nvariables = m_session->GetVariables().size();
int i;
m_base = Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
m_base_aux=Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
if (m_projectionType == MultiRegions::eGalerkin)
{
switch (m_expdim)
......@@ -382,9 +372,7 @@ namespace Nektar
{
if(m_HomogeneousType == eHomogeneous1D)
{
if(m_session->DefinesSolverInfo("SingleMode")&&
(m_session->GetSolverInfo("SingleMode")=="ModifiedBasis"))
if(m_SingleMode)
{
const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierSingleModeSpaced);
const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,m_npointsZ,PkeyZ);
......@@ -396,8 +384,7 @@ namespace Nektar
}
}
else if(m_session->DefinesSolverInfo("SingleMode")&&
(m_session->GetSolverInfo("SingleMode")=="HalfMode"))
else if(m_HalfMode)
{
//1 plane field (half mode expansion)