Commit 910fe1b7 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'fix/INSCoupledStabilitySolver' into 'master'

Fix/ins coupled stability solver

This MR fixes the coupled stability solver.

See merge request !508
parents 168ca430 da4e8601
......@@ -320,10 +320,10 @@ Possible values are:
{Equations} & {\texttt{EQTYPE}} &{Dim.}&{Projections} & Alg.\\
\midrule
Steady Stokes (SS)& \texttt{SteadyStokes} & All & CG &VCS \\
Steady Onseen (SO) & \texttt{SteadyOseen} & All & CG& DS \\
Steady Oseen (SO) & \texttt{SteadyOseen} & All & CG& DS \\
Unsteady Stokes (US) & \texttt{UnsteadyStokes} & All & CG &VCS \\
Steady Linearised NS (SLNS) & \texttt{SteadyLinearisedNS} & All & CG & DS \\
Unsteady Linearised NS (ULNS) & \texttt{UnsteadyLinearisedNS} & All & CG & DS \\
Unsteady Linearised NS (ULNS) & \texttt{UnsteadyLinearisedNS} & All & CG & VCS,DS \\
Unsteady NS (UNS) & \texttt{UnsteadyNavierStokes} & All & CG,CG-DG & VCS \\
\bottomrule
\end{tabular}
......@@ -676,30 +676,64 @@ The following parameters can be specified in the \texttt{PARAMETERS} section of
\item \inltt{Re}: sets the Reynolds number
\item \inltt{Kinvis}: sets the kinematic viscosity $\nu$.
\item \inltt{kdim}: sets the dimension of the Krylov subspace $\kappa$. Can be used in: \inltt{ModifiedArnoldi} and \inltt{Arpack}. Default value: 16.
\item \inltt{evtol}: sets the tolerance of the eigenvalues. Can be used in: \inltt{ModifiedArnoldi} and \inltt{Arpack}. Default value: $10^{-6}$.
\item \inltt{evtol}: sets the tolerance of the eigenvalues. Can be used in: \item \inltt{imagShift}: provide an imaginary shift to the direct sovler eigenvalue problem by the specified value.
ltt{ModifiedArnoldi} and \inltt{Arpack}. Default value: $10^{-6}$.
\item \inltt{nits}: sets the maximum number of iterations. Can be used in: \inltt{ModifiedArnoldi} and \inltt{Arpack}. Default value: 500.
\item \inltt{LZ}: sets the length in the spanswise direction $L_z$. Can be used in \inltt{Homogeneous} set to \inltt{1D}. Default value: 1.
\item \inltt{HomModesZ}: sets the number of planes in the homogeneous directions. Can be used in \inltt{Homogeneous} set to \inltt{1D} and \inltt{ModeType} set to \inltt{MultipleModes}.
\item \inltt{N\_slices}: sets the number of temporal slices for Floquet stability analysis.
\item \inltt{period}: sets the periodicity of the base flow.
\item \inltt{realShift}: provide a real shift to the direct sovler eigenvalue problem by the specified value.
\end{itemize}
\subsection{Functions}
\begin{itemize}
\item To be INserted
\end{itemize}
When using the direct solver for stability analysis it is necessary to specify a Forcing function ``StabilityCoupledLNS'' in the form:
\begin{lstlisting}[style=XMLStyle]
<FORCING>
<FORCE TYPE="StabilityCoupledLNS">
</FORCE>
</FORCING>
\end{lstlisting}
This is required since we need to tell the solver to use the existing
field as a forcing function to the direct matrix inverse as part of
the Arnoldi iteration.
\begin{notebox}
Examples of the set up of the direct solver stability analysis (and
other incompressible Navier-Stokes solvers) can be found in the
regression test directory
\inlsh{NEKTAR/solvers/IncNavierStokesSolver/Tests}. See for example the files
\inlsh{PPF\_R15000\_ModifiedArnoldi\_Shift.tst} and \inlsh{PPF\_R15000\_3D.xml}
noting that some parameters are specified in the .tst files.
\end{notebox}
\section{Steady-state solver Session file configuration}
\label{SectionSFD_XML}
In this section, we detail how to use the steady-state solver (that implements the selective frequency damping method, see Sec. \ref{SectionSFD}).
Two cases are detailed here: the execution of the classical SFD method and the \textit{adaptive} SFD method, where the control coefficient $\chi$ and the filter width $\Delta$ of the SFD method are updated all along the solver execution. For the second case, the parameters of the SFD method do not need to be defined by the user (they will be automatically calculated all along the solver execution) but several session files must be defined in a very specific way.
In this section, we detail how to use the steady-state solver (that
implements the selective frequency damping method, see
Sec. \ref{SectionSFD}). Two cases are detailed here: the execution of
the classical SFD method and the \textit{adaptive} SFD method, where
the control coefficient $\chi$ and the filter width $\Delta$ of the
SFD method are updated all along the solver execution. For the second
case, the parameters of the SFD method do not need to be defined by
the user (they will be automatically calculated all along the solver
execution) but several session files must be defined in a very
specific way.
\subsection{Execution of the classical steady-state solver}
\subsubsection{Solver Info}
The definition of \inltt{Eqtype}, \inltt{TimeIntegrationMethod} and \inltt{Projection} is similar as what is explained in \ref{SectionIncNS_SolverInfo_Stab}. The use of the steady-state solver is enforced through the definition of the \inltt{Driver}: is has to be \inltt{SteadyState}. \inltt{EvolutionOperator} does not need to be defined to run the unadapted SFD method (by default, it is \inltt{Nonlinear}).
The definition of \inltt{Eqtype}, \inltt{TimeIntegrationMethod} and
\inltt{Projection} is similar as what is explained in
\ref{SectionIncNS_SolverInfo_Stab}. The use of the steady-state solver
is enforced through the definition of the \inltt{Driver} which has to be
\inltt{SteadyState}. \inltt{EvolutionOperator} does not need to be
defined to run the unadapted SFD method (by default, it is set to
\inltt{Nonlinear}).
\subsubsection{Parameters}
......
......@@ -110,17 +110,17 @@ namespace Nektar
EquationSharedPtr m_expression;
std::string m_fileVariable;
};
typedef std::map<std::pair<std::string,int>, FunctionVariableDefinition>
typedef std::map<std::pair<std::string,int>, FunctionVariableDefinition>
FunctionVariableMap;
typedef std::map<std::string, FunctionVariableMap >
typedef std::map<std::string, FunctionVariableMap >
FunctionMap;
class SessionReader;
typedef boost::shared_ptr<SessionReader> SessionReaderSharedPtr;
/// Reads and parses information from a Nektar++ XML session file.
class SessionReader :
class SessionReader :
public boost::enable_shared_from_this<SessionReader>
{
public:
......@@ -157,9 +157,9 @@ namespace Nektar
* of the object.
*/
LIB_UTILITIES_EXPORT static SessionReaderSharedPtr CreateInstance(
int argc,
char *argv[],
std::vector<std::string> &pFilenames,
int argc,
char *argv[],
std::vector<std::string> &pFilenames,
const CommSharedPtr &pComm = CommSharedPtr())
{
SessionReaderSharedPtr p = MemoryManager<
......@@ -170,9 +170,9 @@ namespace Nektar
}
LIB_UTILITIES_EXPORT SessionReader(
int argc,
char *argv[],
const std::vector<std::string> &pFilenames,
int argc,
char *argv[],
const std::vector<std::string> &pFilenames,
const CommSharedPtr &pComm);
/// Destructor
......@@ -207,29 +207,29 @@ namespace Nektar
const std::string &pName) const;
/// Load an integer parameter
LIB_UTILITIES_EXPORT void LoadParameter(
const std::string &name,
const std::string &name,
int &var) const;
/// Check for and load an integer parameter.
LIB_UTILITIES_EXPORT void LoadParameter(
const std::string &name,
int &var,
const std::string &name,
int &var,
const int &def) const;
/// Load a double precision parameter
LIB_UTILITIES_EXPORT void LoadParameter(
const std::string &name,
const std::string &name,
NekDouble &var) const;
/// Check for and load a double-precision parameter.
LIB_UTILITIES_EXPORT void LoadParameter(
const std::string &name,
NekDouble &var,
const std::string &name,
NekDouble &var,
const NekDouble &def) const;
/// Set an integer parameter
LIB_UTILITIES_EXPORT void SetParameter(
const std::string &name,
const std::string &name,
int &var);
/// Set a double precision parameter
LIB_UTILITIES_EXPORT void SetParameter(
const std::string &name,
const std::string &name,
NekDouble &var);
......@@ -252,35 +252,35 @@ namespace Nektar
const std::string &vValue) const;
/// Check for and load a solver info property.
LIB_UTILITIES_EXPORT void LoadSolverInfo(
const std::string &name,
std::string &var,
const std::string &name,
std::string &var,
const std::string &def = "") const;
/// Check if the value of a solver info property matches.
LIB_UTILITIES_EXPORT void MatchSolverInfo(
const std::string &name,
const std::string &trueval,
bool &var,
const std::string &name,
const std::string &trueval,
bool &var,
const bool &def = false) const;
/// Check if the value of a solver info property matches.
LIB_UTILITIES_EXPORT bool MatchSolverInfo(
const std::string &name,
const std::string &name,
const std::string &trueval) const;
/// Check if the value of a solver info property matches.
template<typename T>
inline bool MatchSolverInfoAsEnum(
const std::string &name,
const std::string &name,
const T &trueval) const;
/// Registers an enumeration value.
LIB_UTILITIES_EXPORT inline static std::string RegisterEnumValue(
std::string pEnum,
std::string pString,
std::string pEnum,
std::string pString,
int pEnumValue);
/// Registers the default string value of a solver info property.
LIB_UTILITIES_EXPORT inline static std::string
LIB_UTILITIES_EXPORT inline static std::string
RegisterDefaultSolverInfo(
const std::string &pName,
const std::string &pName,
const std::string &pValue);
/* ----GlobalSysSolnInfo ----- */
LIB_UTILITIES_EXPORT bool DefinesGlobalSysSolnInfo(
......@@ -298,24 +298,24 @@ namespace Nektar
const std::string &name) const;
/// Checks for and load a geometric info string property.
LIB_UTILITIES_EXPORT void LoadGeometricInfo(
const std::string &name,
std::string &var,
const std::string &name,
std::string &var,
const std::string &def = "") const;
/// Checks for and loads a geometric info boolean property.
LIB_UTILITIES_EXPORT void LoadGeometricInfo(
const std::string &name,
bool &var,
const std::string &name,
bool &var,
const bool &def = false) const;
/// Checks for and loads a geometric info double-precision property.
LIB_UTILITIES_EXPORT void LoadGeometricInfo(
const std::string &name,
NekDouble &var,
const std::string &name,
NekDouble &var,
const NekDouble &def = 0.0) const;
/// Check if the value of a geometric info string property matches.
LIB_UTILITIES_EXPORT void MatchGeometricInfo(
const std::string &name,
const std::string &trueval,
bool &var,
const std::string &name,
const std::string &trueval,
bool &var,
const bool &def = false) const;
/* ------ VARIABLES ------ */
......@@ -335,37 +335,37 @@ namespace Nektar
const std::string &name) const;
/// Checks if a specified function has a given variable defined.
LIB_UTILITIES_EXPORT bool DefinesFunction(
const std::string &name,
const std::string &name,
const std::string &variable,
const int pDomain = 0) const;
/// Returns an EquationSharedPtr to a given function variable.
LIB_UTILITIES_EXPORT EquationSharedPtr GetFunction(
const std::string &name,
const std::string &name,
const std::string &variable,
const int pDomain = 0) const;
/// Returns an EquationSharedPtr to a given function variable index.
LIB_UTILITIES_EXPORT EquationSharedPtr GetFunction(
const std::string &name,
const std::string &name,
const unsigned int &var,
const int pDomain = 0) const;
/// Returns the type of a given function variable.
LIB_UTILITIES_EXPORT enum FunctionType GetFunctionType(
const std::string &name,
const std::string &name,
const std::string &variable,
const int pDomain = 0) const;
/// Returns the type of a given function variable index.
LIB_UTILITIES_EXPORT enum FunctionType GetFunctionType(
const std::string &pName,
const std::string &pName,
const unsigned int &pVar,
const int pDomain = 0) const;
/// Returns the filename to be loaded for a given variable.
LIB_UTILITIES_EXPORT std::string GetFunctionFilename(
const std::string &name,
const std::string &name,
const std::string &variable,
const int pDomain = 0) const;
/// Returns the filename to be loaded for a given variable index.
LIB_UTILITIES_EXPORT std::string GetFunctionFilename(
const std::string &name,
const std::string &name,
const unsigned int &var,
const int pDomain = 0) const;
/// Returns the filename variable to be loaded for a given variable
......@@ -377,7 +377,7 @@ namespace Nektar
/// Returns the instance of AnalyticExpressionEvaluator specific to
/// this session.
LIB_UTILITIES_EXPORT AnalyticExpressionEvaluator&
LIB_UTILITIES_EXPORT AnalyticExpressionEvaluator&
GetExpressionEvaluator();
/* ------ TAGS ------ */
......@@ -386,7 +386,7 @@ namespace Nektar
const std::string& pName) const;
/// Sets a specified tag.
LIB_UTILITIES_EXPORT void SetTag(
const std::string& pName,
const std::string& pName,
const std::string& pValue);
/// Returns the value of a specified tag.
LIB_UTILITIES_EXPORT const std::string &GetTag(
......@@ -407,10 +407,10 @@ namespace Nektar
return m_cmdLineOptions.find(pName)->second.as<T>();
}
/// Registers a command-line argument with the session reader.
LIB_UTILITIES_EXPORT inline static std::string
LIB_UTILITIES_EXPORT inline static std::string
RegisterCmdLineArgument(
const std::string &pName,
const std::string &pShortName,
const std::string &pName,
const std::string &pShortName,
const std::string &pDescription);
/// Registers a command-line flag with the session reader.
LIB_UTILITIES_EXPORT inline static std::string
......@@ -462,7 +462,7 @@ namespace Nektar
/// Map of original boundary region ordering for parallel periodic
/// bcs.
BndRegionOrdering m_bndRegOrder;
/// String to enumeration map for Solver Info parameters.
/// String to enumeration map for Solver Info parameters.
LIB_UTILITIES_EXPORT static EnumMapList& GetSolverInfoEnums();
/// Default solver info options.
LIB_UTILITIES_EXPORT static SolverInfoMap& GetSolverInfoDefaults();
......@@ -473,7 +473,7 @@ namespace Nektar
/// Main constructor
LIB_UTILITIES_EXPORT SessionReader(
int argc,
int argc,
char *argv[]);
LIB_UTILITIES_EXPORT void InitSession();
......@@ -500,7 +500,7 @@ namespace Nektar
/// Loads the given XML document and instantiates an appropriate
/// communication object.
LIB_UTILITIES_EXPORT void CreateComm(
int &argc,
int &argc,
char* argv[]);
/// Partitions the mesh when running in parallel.
......@@ -525,6 +525,8 @@ namespace Nektar
LIB_UTILITIES_EXPORT void ReadFilters(TiXmlElement *filters);
/// Enforce parameters from command line arguments.
LIB_UTILITIES_EXPORT void CmdLineOverride();
/// Check values of solver info options are valid.
LIB_UTILITIES_EXPORT void VerifySolverInfo();
/// Parse a string in the form lhs = rhs.
LIB_UTILITIES_EXPORT void ParseEquals(
......@@ -617,7 +619,7 @@ namespace Nektar
* @param pString A valid value for the property.
* @param pEnumValue An enumeration value corresponding to this
* value.
*
*
* @return The value for the property provided by #pString.
*/
inline std::string SessionReader::RegisterEnumValue(
......@@ -641,7 +643,7 @@ namespace Nektar
* using this function. The property will take this value until it is
* overwritten by a value specified in the XML document, or specified
* as a command-line argument. Usage has the form:
*
*
* @code
* std::string GlobalLinSys::def
* = LibUtilities::SessionReader::RegisterDefaultSolverInfo(
......@@ -650,11 +652,11 @@ namespace Nektar
*
* @param pName The name of the property.
* @param pValue The default value of the property.
*
*
* @return The default value of the property provided by #pValue.
*/
inline std::string SessionReader::RegisterDefaultSolverInfo(
const std::string &pName,
const std::string &pName,
const std::string &pValue)
{
std::string vName = boost::to_upper_copy(pName);
......@@ -667,8 +669,8 @@ namespace Nektar
*
*/
inline std::string SessionReader::RegisterCmdLineArgument(
const std::string &pName,
const std::string &pShortName,
const std::string &pName,
const std::string &pShortName,
const std::string &pDescription)
{
ASSERTL0(!pName.empty(), "Empty name for cmdline argument.");
......
......@@ -55,23 +55,25 @@ namespace Nektar
typedef vector<map<int, int> > DofGraph;
MULTI_REGIONS_EXPORT
pair<int, StdRegions::Orientation> DeterminePeriodicEdgeOrientId(
int meshEdgeId,
StdRegions::Orientation edgeOrient,
const vector<PeriodicEntity> &periodicEdges);
MULTI_REGIONS_EXPORT
StdRegions::Orientation DeterminePeriodicFaceOrient(
StdRegions::Orientation faceOrient1,
StdRegions::Orientation faceOrient2);
/// Constructs mappings for the C0 scalar continuous Galerkin formulation.
class AssemblyMapCG: public AssemblyMap
{
typedef Array<OneD, const ExpListSharedPtr> BndCondExp;
typedef Array<OneD, const SpatialDomains::BoundaryConditionShPtr>
BndCond;
public:
/// Default constructor.
MULTI_REGIONS_EXPORT AssemblyMapCG(
......
......@@ -38,253 +38,285 @@
namespace Nektar
{
namespace SolverUtils
namespace SolverUtils
{
/**
* Constructor
*/
DriverArnoldi::DriverArnoldi(const LibUtilities::SessionReaderSharedPtr pSession)
: Driver(pSession)
{
m_session->LoadParameter("IO_InfoSteps", m_infosteps, 1);
};
/**
* Destructor
*/
DriverArnoldi::~DriverArnoldi()
{
};
/**
* Arnoldi driver initialisation
*/
void DriverArnoldi::v_InitObject(ostream &out)
{
Driver::v_InitObject(out);
m_session->MatchSolverInfo("SolverType",
"VelocityCorrectionScheme",
m_timeSteppingAlgorithm, false);
if (m_timeSteppingAlgorithm)
{
/**
*
*/
DriverArnoldi::DriverArnoldi(const LibUtilities::SessionReaderSharedPtr pSession)
: Driver(pSession)
{
m_session->LoadParameter("IO_InfoSteps", m_infosteps, 1);
};
m_period = m_session->GetParameter("TimeStep")
* m_session->GetParameter("NumSteps");
m_nfields = m_equ[0]->UpdateFields().num_elements() - 1;
}
else
{
m_period = 1.0;
m_nfields = m_equ[0]->UpdateFields().num_elements();
}
/**
*
*/
DriverArnoldi::~DriverArnoldi()
if(m_session->DefinesSolverInfo("ModeType") &&
(boost::iequals(m_session->GetSolverInfo("ModeType"),
"SingleMode")||
boost::iequals(m_session->GetSolverInfo("ModeType"),
"HalfMode")))
{
for(int i = 0; i < m_nfields; ++i)
{
};
m_equ[0]->UpdateFields()[i]->SetWaveSpace(true);
}
}
m_negatedOp = m_equ[0]->v_NegatedOp();
m_session->LoadParameter("kdim", m_kdim, 16);
m_session->LoadParameter("nvec", m_nvec, 2);
m_session->LoadParameter("nits", m_nits, 500);
m_session->LoadParameter("evtol", m_evtol, 1e-06);
m_session->LoadParameter("realShift", m_realShift, 0.0);
m_equ[0]->SetLambda(m_realShift);
/**
*
*/
void DriverArnoldi::v_InitObject(ostream &out)
m_session->LoadParameter("imagShift", m_imagShift, 0.0);
}
void DriverArnoldi::ArnoldiSummary(std::ostream &out)
{
if (m_comm->GetRank() == 0)
{
if(m_session->DefinesSolverInfo("ModeType") &&
boost::iequals(m_session->GetSolverInfo("ModeType"),
"SingleMode"))
{
Driver::v_InitObject(out);
m_session->MatchSolverInfo("SolverType","VelocityCorrectionScheme",m_timeSteppingAlgorithm, false);
if(m_timeSteppingAlgorithm)
{
m_period = m_session->GetParameter("TimeStep")* m_session->GetParameter("NumSteps");
m_nfields = m_equ[0]->UpdateFields().num_elements() - 1;
if(m_session->DefinesSolverInfo("ModeType") &&
(m_session->GetSolverInfo("ModeType")=="SingleMode"|| m_session->GetSolverInfo("ModeType")=="HalfMode") )
{
for(int i = 0; i < m_nfields; ++i)
{
m_equ[0]->UpdateFields()[i]->SetWaveSpace(true);
}
}
}
else
{
m_period = 1.0;
ASSERTL0(m_session->DefinesFunction("BodyForce"),"A BodyForce section needs to be defined for this solver type");
m_nfields = m_equ[0]->UpdateFields().num_elements();
}
m_session->LoadParameter("kdim", m_kdim, 16);
m_session->LoadParameter("nvec", m_nvec, 2);
m_session->LoadParameter("nits", m_nits, 500);
m_session->LoadParameter("evtol", m_evtol, 1e-06);
m_session->LoadParameter("realShift", m_realShift, 0.0);
m_equ[0]->SetLambda(m_realShift);
m_session->LoadParameter("imagShift", m_imagShift, 0.0);
out << "\tSingle Fourier mode : true " << endl;
ASSERTL0(m_session->DefinesSolverInfo("Homogeneous"),
"Expected a homogeneous expansion to be defined "
"with single mode");
}
else
{
out << "\tSingle Fourier mode : false " << endl;
}
if(m_session->DefinesSolverInfo("BetaZero"))
{
out << "\tBeta set to Zero : true (overrides LHom)"
<< endl;
}
else
{
out << "\tBeta set to Zero : false " << endl;
}
void DriverArnoldi::ArnoldiSummary(std::ostream &out)
if(m_timeSteppingAlgorithm)
{
if (m_comm->GetRank() == 0)
{
if(m_session->DefinesSolverInfo("SingleMode"))
{
out << "\tSingle Fourier mode : true " << endl;
ASSERTL0(m_session->DefinesSolverInfo("Homogeneous"),
"Expected a homogeneous expansion to be defined "
"with single mode");
}
else
{
out << "\tSingle Fourier mode : false " << endl;
}