Nektar++
Public Member Functions | Static Public Member Functions | Public Attributes | Static Public Attributes | Protected Member Functions | Static Protected Attributes | Private Member Functions | Private Attributes | Friends | List of all members
Nektar::CoupledLinearNS Class Reference

#include <CoupledLinearNS.h>

Inheritance diagram for Nektar::CoupledLinearNS:
[legend]

Public Member Functions

void SetUpCoupledMatrix (const NekDouble lambda=0.0, const Array< OneD, Array< OneD, NekDouble > > &Advfield=NullNekDoubleArrayOfArray, bool IsLinearNSEquation=true)
 
const SpatialDomains::ExpansionInfoMap & GenPressureExp (const SpatialDomains::ExpansionInfoMap &VelExp)
 
void Solve (void)
 
void SolveLinearNS (const Array< OneD, Array< OneD, NekDouble > > &forcing)
 
void SolveLinearNS (Array< OneD, Array< OneD, NekDouble > > &forcing, Array< OneD, MultiRegions::ExpListSharedPtr > &fields, MultiRegions::ExpListSharedPtr &pressure, const int HomogeneousMode=0)
 
void SolveUnsteadyStokesSystem (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time, const NekDouble a_iixDt)
 
void EvaluateAdvection (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
 
void SolveSteadyNavierStokes (void)
 
void Continuation (void)
 
void EvaluateNewtonRHS (Array< OneD, Array< OneD, NekDouble > > &Velocity, Array< OneD, Array< OneD, NekDouble > > &outarray)
 
void InfNorm (Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
 
void L2Norm (Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
 
void DefineForcingTerm (void)
 
- Public Member Functions inherited from Nektar::IncNavierStokes
 ~IncNavierStokes () override
 
void v_InitObject (bool DeclareField=true) override
 Initialisation object for EquationSystem. More...
 
int GetNConvectiveFields (void)
 
void AddForcing (const SolverUtils::ForcingSharedPtr &pForce)
 
void v_GetPressure (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure) override
 
void v_GetDensity (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &density) override
 
bool v_HasConstantDensity () override
 
void v_GetVelocity (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity) override
 
void v_SetMovingFrameVelocities (const Array< OneD, NekDouble > &vFrameVels) override
 
void v_GetMovingFrameVelocities (Array< OneD, NekDouble > &vFrameVels) override
 
void v_SetMovingFrameDisp (const Array< OneD, NekDouble > &vFrameDisp) override
 
void v_GetMovingFrameDisp (Array< OneD, NekDouble > &vFrameDisp) override
 
void v_SetMovingFrameProjectionMat (const bnu::matrix< NekDouble > &vProjMat) override
 
void v_GetMovingFrameProjectionMat (bnu::matrix< NekDouble > &vProjMat) override
 
void v_SetAeroForce (Array< OneD, NekDouble > forces) override
 
void v_GetAeroForce (Array< OneD, NekDouble > forces) override
 
bool DefinedForcing (const std::string &sForce)
 
- Public Member Functions inherited from Nektar::SolverUtils::AdvectionSystem
SOLVER_UTILS_EXPORT AdvectionSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 
SOLVER_UTILS_EXPORT ~AdvectionSystem () override
 
SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareField=true) override
 Initialisation object for EquationSystem. More...
 
SOLVER_UTILS_EXPORT AdvectionSharedPtr GetAdvObject ()
 Returns the advection object held by this instance. More...
 
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > GetElmtCFLVals (const bool FlagAcousticCFL=true)
 
SOLVER_UTILS_EXPORT NekDouble GetCFLEstimate (int &elmtid)
 
- Public Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
SOLVER_UTILS_EXPORT ~UnsteadySystem () override
 Destructor. More...
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep (const Array< OneD, const Array< OneD, NekDouble > > &inarray)
 Calculate the larger time-step mantaining the problem stable. More...
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void SetTimeStep (const NekDouble timestep)
 
SOLVER_UTILS_EXPORT void SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeSharedPtr & GetTimeIntegrationScheme ()
 Returns the time integration scheme. More...
 
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeOperators & GetTimeIntegrationSchemeOperators ()
 Returns the time integration scheme operators. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT void InitObject (bool DeclareField=true)
 Initialises the members of this object. More...
 
SOLVER_UTILS_EXPORT void DoInitialise (bool dumpInitialConditions=true)
 Perform any initialisation necessary before solving the problem. More...
 
SOLVER_UTILS_EXPORT void DoSolve ()
 Solve the problem. More...
 
SOLVER_UTILS_EXPORT void TransCoeffToPhys ()
 Transform from coefficient to physical space. More...
 
SOLVER_UTILS_EXPORT void TransPhysToCoeff ()
 Transform from physical to coefficient space. More...
 
SOLVER_UTILS_EXPORT void Output ()
 Perform output operations after solve. More...
 
SOLVER_UTILS_EXPORT std::string GetSessionName ()
 Get Session name. More...
 
template<class T >
std::shared_ptr< T > as ()
 
SOLVER_UTILS_EXPORT void ResetSessionName (std::string newname)
 Reset Session name. More...
 
SOLVER_UTILS_EXPORT LibUtilities::SessionReaderSharedPtr GetSession ()
 Get Session name. More...
 
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure ()
 Get pressure field if available. More...
 
SOLVER_UTILS_EXPORT void ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 
SOLVER_UTILS_EXPORT void PrintSummary (std::ostream &out)
 Print a summary of parameters and solver characteristics. More...
 
SOLVER_UTILS_EXPORT void SetLambda (NekDouble lambda)
 Set parameter m_lambda. More...
 
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction (std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
 Get a SessionFunction by name. More...
 
SOLVER_UTILS_EXPORT void SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 Initialise the data in the dependent fields. More...
 
SOLVER_UTILS_EXPORT void EvaluateExactSolution (int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 Evaluates an exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised=false)
 Compute the L2 error between fields and a given exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, bool Normalised=false)
 Compute the L2 error of the fields. More...
 
SOLVER_UTILS_EXPORT NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation. More...
 
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > ErrorExtraPoints (unsigned int field)
 Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf]. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n)
 Write checkpoint file of m_fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write checkpoint file of custom data fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow (const int n)
 Write base flow file of m_fields. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname)
 Write field data to the given filename. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write input fields to the given filename. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 Input field data from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFldToMultiDomains (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int ndomains)
 Input field data from the given file to multiple domains. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, std::vector< std::string > &fieldStr, Array< OneD, Array< OneD, NekDouble > > &coeffs)
 Output a field. Input field data into array from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, MultiRegions::ExpListSharedPtr &pField, std::string &pFieldName)
 Output a field. Input field data into ExpList from the given file. More...
 
SOLVER_UTILS_EXPORT void SessionSummary (SummaryList &vSummary)
 Write out a session summary. More...
 
SOLVER_UTILS_EXPORT Array< OneD, MultiRegions::ExpListSharedPtr > & UpdateFields ()
 
SOLVER_UTILS_EXPORT LibUtilities::FieldMetaDataMap & UpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated. More...
 
SOLVER_UTILS_EXPORT NekDouble GetTime ()
 Return final time. More...
 
SOLVER_UTILS_EXPORT int GetNcoeffs ()
 
SOLVER_UTILS_EXPORT int GetNcoeffs (const int eid)
 
SOLVER_UTILS_EXPORT int GetNumExpModes ()
 
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp ()
 
SOLVER_UTILS_EXPORT int GetNvariables ()
 
SOLVER_UTILS_EXPORT const std::string GetVariable (unsigned int i)
 
SOLVER_UTILS_EXPORT int GetTraceTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTraceNpoints ()
 
SOLVER_UTILS_EXPORT int GetExpSize ()
 
SOLVER_UTILS_EXPORT int GetPhys_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetCoeff_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTotPoints (int n)
 
SOLVER_UTILS_EXPORT int GetNpoints ()
 
SOLVER_UTILS_EXPORT int GetSteps ()
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void CopyFromPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void CopyToPhysField (const int i, const Array< OneD, const NekDouble > &input)
 
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > & UpdatePhysField (const int i)
 
SOLVER_UTILS_EXPORT void SetSteps (const int steps)
 
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
 
SOLVER_UTILS_EXPORT void FwdTransFields ()
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int GetCheckpointNumber ()
 
SOLVER_UTILS_EXPORT void SetCheckpointNumber (int num)
 
SOLVER_UTILS_EXPORT int GetCheckpointSteps ()
 
SOLVER_UTILS_EXPORT void SetCheckpointSteps (int num)
 
SOLVER_UTILS_EXPORT int GetInfoSteps ()
 
SOLVER_UTILS_EXPORT void SetInfoSteps (int num)
 
SOLVER_UTILS_EXPORT void SetIterationNumberPIT (int num)
 
SOLVER_UTILS_EXPORT void SetWindowNumberPIT (int num)
 
SOLVER_UTILS_EXPORT Array< OneD, const Array< OneD, NekDouble > > GetTraceNormals ()
 
SOLVER_UTILS_EXPORT void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetTimeStep (const NekDouble timestep)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 
SOLVER_UTILS_EXPORT bool NegatedOp ()
 Identify if operator is negated in DoSolve. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::FluidInterface
virtual ~FluidInterface ()=default
 
SOLVER_UTILS_EXPORT void GetVelocity (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
 Extract array with velocity from physfield. More...
 
SOLVER_UTILS_EXPORT bool HasConstantDensity ()
 
SOLVER_UTILS_EXPORT void GetDensity (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &density)
 Extract array with density from physfield. More...
 
SOLVER_UTILS_EXPORT void GetPressure (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure)
 Extract array with pressure from physfield. More...
 
SOLVER_UTILS_EXPORT void SetMovingFrameVelocities (const Array< OneD, NekDouble > &vFrameVels)
 
SOLVER_UTILS_EXPORT void GetMovingFrameVelocities (Array< OneD, NekDouble > &vFrameVels)
 
SOLVER_UTILS_EXPORT void SetMovingFrameProjectionMat (const boost::numeric::ublas::matrix< NekDouble > &vProjMat)
 
SOLVER_UTILS_EXPORT void GetMovingFrameProjectionMat (boost::numeric::ublas::matrix< NekDouble > &vProjMat)
 
SOLVER_UTILS_EXPORT void SetMovingFrameDisp (const Array< OneD, NekDouble > &vFrameDisp)
 
SOLVER_UTILS_EXPORT void GetMovingFrameDisp (Array< OneD, NekDouble > &vFrameDisp)
 
SOLVER_UTILS_EXPORT void SetAeroForce (Array< OneD, NekDouble > forces)
 Set aerodynamic force and moment. More...
 
SOLVER_UTILS_EXPORT void GetAeroForce (Array< OneD, NekDouble > forces)
 Get aerodynamic force and moment. More...
 

Static Public Member Functions

static SolverUtils::EquationSystemSharedPtr create (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Creates an instance of this class. More...
 

Public Attributes

Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm
 
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm_Coeffs
 
Array< OneD, CoupledLocalToGlobalC0ContMapSharedPtr > m_locToGloMap
 

Static Public Attributes

static std::string className
 Name of class. More...
 
- Static Public Attributes inherited from Nektar::SolverUtils::UnsteadySystem
static std::string cmdSetStartTime
 
static std::string cmdSetStartChkNum
 

Protected Member Functions

 CoupledLinearNS (const LibUtilities::SessionReaderSharedPtr &pSesssion, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 
void v_InitObject (bool DeclareField=true) override
 Initialisation object for EquationSystem. More...
 
- Protected Member Functions inherited from Nektar::IncNavierStokes
 IncNavierStokes (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Constructor. More...
 
EquationType GetEquationType (void)
 
void EvaluateAdvectionTerms (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
 
void WriteModalEnergy (void)
 
void SetBoundaryConditions (NekDouble time)
 time dependent boundary conditions updating More...
 
void SetRadiationBoundaryForcing (int fieldid)
 Set Radiation forcing term. More...
 
void SetZeroNormalVelocity ()
 Set Normal Velocity Component to Zero. More...
 
void SetWomersleyBoundary (const int fldid, const int bndid)
 Set Womersley Profile if specified. More...
 
void SetUpWomersley (const int fldid, const int bndid, std::string womstr)
 Set Up Womersley details. More...
 
void SetMovingReferenceFrameBCs (const NekDouble &time)
 Set the moving reference frame boundary conditions. More...
 
void SetMRFWallBCs (const NekDouble &time)
 
void SetMRFDomainVelBCs (const NekDouble &time)
 
MultiRegions::ExpListSharedPtr v_GetPressure () override
 
void v_TransCoeffToPhys (void) override
 Virtual function for transformation to physical space. More...
 
void v_TransPhysToCoeff (void) override
 Virtual function for transformation to coefficient space. More...
 
virtual int v_GetForceDimension ()=0
 
Array< OneD, NekDouble > v_GetMaxStdVelocity (const NekDouble SpeedSoundFactor) override
 
bool v_PreIntegrate (int step) override
 
- Protected Member Functions inherited from Nektar::SolverUtils::AdvectionSystem
SOLVER_UTILS_EXPORT bool v_PostIntegrate (int step) override
 
virtual SOLVER_UTILS_EXPORT Array< OneD, NekDouble > v_GetMaxStdVelocity (const NekDouble SpeedSoundFactor=1.0)
 
- Protected Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
SOLVER_UTILS_EXPORT UnsteadySystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises UnsteadySystem class members. More...
 
SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareField=true) override
 Init object for UnsteadySystem class. More...
 
SOLVER_UTILS_EXPORT void v_DoSolve () override
 Solves an unsteady problem. More...
 
virtual SOLVER_UTILS_EXPORT void v_PrintStatusInformation (const int step, const NekDouble cpuTime)
 Print Status Information. More...
 
virtual SOLVER_UTILS_EXPORT void v_PrintSummaryStatistics (const NekDouble intTime)
 Print Summary Statistics. More...
 
SOLVER_UTILS_EXPORT void v_DoInitialise (bool dumpInitialConditions=true) override
 Sets up initial conditions. More...
 
SOLVER_UTILS_EXPORT void v_GenerateSummary (SummaryList &s) override
 Print a summary of time stepping parameters. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_GetTimeStep (const Array< OneD, const Array< OneD, NekDouble > > &inarray)
 Return the timestep to be used for the next step in the time-marching loop. More...
 
virtual SOLVER_UTILS_EXPORT bool v_PreIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_RequireFwdTrans ()
 
virtual SOLVER_UTILS_EXPORT void v_SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
virtual SOLVER_UTILS_EXPORT bool v_UpdateTimeStepCheck ()
 
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator ()
 Get the maximum timestep estimator for cfl control. More...
 
SOLVER_UTILS_EXPORT void CheckForRestartTime (NekDouble &time, int &nchk)
 
SOLVER_UTILS_EXPORT void SVVVarDiffCoeff (const Array< OneD, Array< OneD, NekDouble > > vel, StdRegions::VarCoeffMap &varCoeffMap)
 Evaluate the SVV diffusion coefficient according to Moura's paper where it should proportional to h time velocity. More...
 
SOLVER_UTILS_EXPORT void DoDummyProjection (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
 Perform dummy projection. More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises EquationSystem class members. More...
 
virtual SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareFeld=true)
 Initialisation object for EquationSystem. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoInitialise (bool dumpInitialConditions=true)
 Virtual function for initialisation implementation. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoSolve ()
 Virtual function for solve implementation. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Virtual function for the L_inf error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
 Virtual function for the L_2 error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT void v_TransCoeffToPhys ()
 Virtual function for transformation to physical space. More...
 
virtual SOLVER_UTILS_EXPORT void v_TransPhysToCoeff ()
 Virtual function for transformation to coefficient space. More...
 
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary (SummaryList &l)
 Virtual function for generating summary information. More...
 
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 
virtual SOLVER_UTILS_EXPORT void v_EvaluateExactSolution (unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 
virtual SOLVER_UTILS_EXPORT void v_Output (void)
 
virtual SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr v_GetPressure (void)
 
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp (void)
 Virtual function to identify if operator is negated in DoSolve. More...
 
virtual SOLVER_UTILS_EXPORT void v_ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 
- Protected Member Functions inherited from Nektar::SolverUtils::FluidInterface
virtual SOLVER_UTILS_EXPORT void v_GetVelocity (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)=0
 
virtual SOLVER_UTILS_EXPORT bool v_HasConstantDensity ()=0
 
virtual SOLVER_UTILS_EXPORT void v_GetDensity (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &density)=0
 
virtual SOLVER_UTILS_EXPORT void v_GetPressure (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure)=0
 
virtual SOLVER_UTILS_EXPORT void v_SetMovingFrameVelocities (const Array< OneD, NekDouble > &vFrameVels)
 
virtual SOLVER_UTILS_EXPORT void v_GetMovingFrameVelocities (Array< OneD, NekDouble > &vFrameVels)
 
virtual SOLVER_UTILS_EXPORT void v_SetMovingFrameProjectionMat (const boost::numeric::ublas::matrix< NekDouble > &vProjMat)
 
virtual SOLVER_UTILS_EXPORT void v_GetMovingFrameProjectionMat (boost::numeric::ublas::matrix< NekDouble > &vProjMat)
 
virtual SOLVER_UTILS_EXPORT void v_SetMovingFrameDisp (const Array< OneD, NekDouble > &vFrameDisp)
 
virtual SOLVER_UTILS_EXPORT void v_GetMovingFrameDisp (Array< OneD, NekDouble > &vFrameDisp)
 
virtual SOLVER_UTILS_EXPORT void v_SetAeroForce (Array< OneD, NekDouble > forces)
 
virtual SOLVER_UTILS_EXPORT void v_GetAeroForce (Array< OneD, NekDouble > forces)
 

Static Protected Attributes

static std::string solverTypeLookupId
 
- Static Protected Attributes inherited from Nektar::IncNavierStokes
static std::string eqTypeLookupIds []
 
- Static Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
static std::string equationSystemTypeLookupIds []
 
static std::string projectionTypeLookupIds []
 

Private Member Functions

void SetUpCoupledMatrix (const NekDouble lambda, const Array< OneD, Array< OneD, NekDouble > > &Advfield, bool IsLinearNSEquation, const int HomogeneousMode, CoupledSolverMatrices &mat, CoupledLocalToGlobalC0ContMapSharedPtr &locToGloMap, const NekDouble lambda_imag=NekConstants::kNekUnsetDouble)
 
void v_GenerateSummary (SolverUtils::SummaryList &s) override
 Print a summary of time stepping parameters. More...
 
void v_DoInitialise (bool dumpInitialConditions=true) override
 Sets up initial conditions. More...
 
void v_DoSolve (void) override
 Solves an unsteady problem. More...
 
bool v_NegatedOp (void) override
 
void v_TransCoeffToPhys (void) override
 Virtual function for transformation to physical space. More...
 
void v_TransPhysToCoeff (void) override
 Virtual function for transformation to coefficient space. More...
 
void v_Output (void) override
 
int v_GetForceDimension (void) override
 

Private Attributes

bool m_zeroMode
 Id to identify when single mode is mean mode (i.e. beta=0);. More...
 
int m_counter
 
bool m_initialStep
 
NekDouble m_tol
 
int m_maxIt
 
int m_Restart
 
int m_MatrixSetUpStep
 
NekDouble m_kinvisMin
 
NekDouble m_kinvisStep
 
NekDouble m_KinvisPercentage
 
Array< OneD, CoupledSolverMatrices > m_mat
 

Friends

class MemoryManager< CoupledLinearNS >
 

Additional Inherited Members

- Protected Types inherited from Nektar::SolverUtils::EquationSystem
enum  HomogeneousType { eHomogeneous1D , eHomogeneous2D , eHomogeneous3D , eNotHomogeneous }
 Parameter for homogeneous expansions. More...
 
- Protected Attributes inherited from Nektar::IncNavierStokes
ExtrapolateSharedPtr m_extrapolation
 
std::ofstream m_mdlFile
 modal energy file More...
 
bool m_SmoothAdvection
 bool to identify if advection term smoothing is requested More...
 
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
 Forcing terms. More...
 
int m_nConvectiveFields
 Number of fields to be convected;. More...
 
Array< OneD, int > m_velocity
 int which identifies which components of m_fields contains the velocity (u,v,w); More...
 
MultiRegions::ExpListSharedPtr m_pressure
 Pointer to field holding pressure field. More...
 
NekDouble m_kinvis
 Kinematic viscosity. More...
 
int m_energysteps
 dump energy to file at steps time More...
 
EquationType m_equationType
 equation type; More...
 
Array< OneD, Array< OneD, int > > m_fieldsBCToElmtID
 Mapping from BCs to Elmt IDs. More...
 
Array< OneD, Array< OneD, int > > m_fieldsBCToTraceID
 Mapping from BCs to Elmt Edge IDs. More...
 
Array< OneD, Array< OneD, NekDouble > > m_fieldsRadiationFactor
 RHS Factor for Radiation Condition. More...
 
int m_intSteps
 Number of time integration steps AND Order of extrapolation for pressure boundary conditions. More...
 
Array< OneD, NekDouble > m_pivotPoint
 
Array< OneD, NekDouble > m_aeroForces
 
std::map< int, std::map< int, WomersleyParamsSharedPtr > > m_womersleyParams
 Womersley parameters if required. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::AdvectionSystem
SolverUtils::AdvectionSharedPtr m_advObject
 Advection term. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::UnsteadySystem
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
 Wrapper to the time integration scheme. More...
 
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use. More...
 
Array< OneD, Array< OneD, NekDouble > > m_previousSolution
 Storage for previous solution for steady-state check. More...
 
std::vector< int > m_intVariables
 
NekDouble m_cflSafetyFactor
 CFL safety factor (comprise between 0 to 1). More...
 
NekDouble m_CFLGrowth
 CFL growth rate. More...
 
NekDouble m_CFLEnd
 Maximun cfl in cfl growth. More...
 
int m_abortSteps
 Number of steps between checks for abort conditions. More...
 
bool m_explicitDiffusion
 Indicates if explicit or implicit treatment of diffusion is used. More...
 
bool m_explicitAdvection
 Indicates if explicit or implicit treatment of advection is used. More...
 
bool m_explicitReaction
 Indicates if explicit or implicit treatment of reaction is used. More...
 
int m_steadyStateSteps
 Check for steady state at step interval. More...
 
NekDouble m_steadyStateTol
 Tolerance to which steady state should be evaluated at. More...
 
int m_filtersInfosteps
 Number of time steps between outputting filters information. More...
 
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters
 
bool m_homoInitialFwd
 Flag to determine if simulation should start in homogeneous forward transformed state. More...
 
std::ofstream m_errFile
 
NekDouble m_epsilon
 Diffusion coefficient. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
bool m_verbose
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader. More...
 
std::map< std::string, SolverUtils::SessionFunctionSharedPtr > m_sessionFunctions
 Map of known SessionFunctions. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
 Array holding all dependent variables. More...
 
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
 Pointer to boundary conditions object. More...
 
SpatialDomains::MeshGraphSharedPtr m_graph
 Pointer to graph defining mesh. More...
 
std::string m_sessionName
 Name of the session. More...
 
NekDouble m_time
 Current time of simulation. More...
 
int m_initialStep
 Number of the step where the simulation should begin. More...
 
NekDouble m_fintime
 Finish time of the simulation. More...
 
NekDouble m_timestep
 Time step size. More...
 
NekDouble m_lambda
 Lambda constant in real system if one required. More...
 
NekDouble m_checktime
 Time between checkpoints. More...
 
NekDouble m_lastCheckTime
 
NekDouble m_TimeIncrementFactor
 
int m_nchk
 Number of checkpoints written so far. More...
 
int m_steps
 Number of steps to take. More...
 
int m_checksteps
 Number of steps between checkpoints. More...
 
int m_infosteps
 Number of time steps between outputting status information. More...
 
int m_iterPIT = 0
 Number of parallel-in-time time iteration. More...
 
int m_windowPIT = 0
 Index of windows for parallel-in-time time iteration. More...
 
int m_spacedim
 Spatial dimension (>= expansion dim). More...
 
int m_expdim
 Expansion dimension. More...
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform. More...
 
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations. More...
 
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More...
 
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous. More...
 
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction. More...
 
Array< OneD, bool > m_checkIfSystemSingular
 Flag to indicate if the fields should be checked for singularity. More...
 
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
 Map to identify relevant solver info to dump in output fields. More...
 
Array< OneD, NekDouble > m_movingFrameVelsxyz
 Moving frame of reference velocities (u, v, w, omega_x, omega_y, omega_z, a_x, a_y, a_z, domega_x, domega_y, domega_z) More...
 
Array< OneD, NekDouble > m_movingFrameData
 Moving frame of reference angles with respect to the. More...
 
boost::numeric::ublas::matrix< NekDouble > m_movingFrameProjMat
 Projection matrix for transformation between inertial and moving. More...
 
int m_NumQuadPointsError
 Number of Quadrature points used to work out the error. More...
 
enum HomogeneousType m_HomogeneousType
 
NekDouble m_LhomX
 physical length in X direction (if homogeneous) More...
 
NekDouble m_LhomY
 physical length in Y direction (if homogeneous) More...
 
NekDouble m_LhomZ
 physical length in Z direction (if homogeneous) More...
 
int m_npointsX
 number of points in X direction (if homogeneous) More...
 
int m_npointsY
 number of points in Y direction (if homogeneous) More...
 
int m_npointsZ
 number of points in Z direction (if homogeneous) More...
 
int m_HomoDirec
 number of homogenous directions More...
 

Detailed Description

Set up expansion field for velocity and pressure, the local to global mapping arrays and the basic memory definitions for coupled matrix system

Definition at line 86 of file CoupledLinearNS.h.

Constructor & Destructor Documentation

◆ CoupledLinearNS()

Nektar::CoupledLinearNS::CoupledLinearNS ( const LibUtilities::SessionReaderSharedPtr &  pSesssion,
const SpatialDomains::MeshGraphSharedPtr &  pGraph 
)
protected

Definition at line 63 of file CoupledLinearNS.cpp.

66 : UnsteadySystem(pSession, pGraph), IncNavierStokes(pSession, pGraph),
67 m_zeroMode(false)
68{
69}
Nektar::CoupledLinearNS::m_zeroMode
bool m_zeroMode
Id to identify when single mode is mean mode (i.e. beta=0);.
Definition: CoupledLinearNS.h:174
Nektar::IncNavierStokes::IncNavierStokes
IncNavierStokes(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Constructor.
Definition: IncNavierStokes.cpp:83
Nektar::SolverUtils::UnsteadySystem::UnsteadySystem
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Initialises UnsteadySystem class members.
Definition: UnsteadySystem.cpp:74

Member Function Documentation

◆ Continuation()

void Nektar::CoupledLinearNS::Continuation ( void  )

Definition at line 1844 of file CoupledLinearNS.cpp.

1845{
1846 Array<OneD, Array<OneD, NekDouble>> u_N(m_velocity.size());
1847 Array<OneD, Array<OneD, NekDouble>> tmp_RHS(m_velocity.size());
1848 Array<OneD, Array<OneD, NekDouble>> RHS(m_velocity.size());
1849 Array<OneD, Array<OneD, NekDouble>> u_star(m_velocity.size());
1850
1851 cout << "We apply the continuation method: " << endl;
1852
1853 for (size_t i = 0; i < m_velocity.size(); ++i)
1854 {
1855 u_N[i] = Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetTotPoints(),
1856 0.0);
1857 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1858 u_N[i]);
1859
1860 RHS[i] = Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetTotPoints(),
1861 0.0);
1862 tmp_RHS[i] = Array<OneD, NekDouble>(
1863 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1864
1865 m_fields[m_velocity[i]]->PhysDeriv(i, u_N[i], tmp_RHS[i]);
1866 Vmath::Smul(tmp_RHS[i].size(), m_kinvis, tmp_RHS[i], 1, tmp_RHS[i], 1);
1867
1868 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1869 m_fields[m_velocity[i]]->SetWaveSpace(true);
1870 m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_RHS[i], RHS[i]);
1871 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1872 }
1873
1874 SetUpCoupledMatrix(0.0, u_N, true);
1875 SolveLinearNS(RHS);
1876
1877 for (size_t i = 0; i < m_velocity.size(); ++i)
1878 {
1879 u_star[i] = Array<OneD, NekDouble>(
1880 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1881 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1882 u_star[i]);
1883
1884 // u_star(k+1) = u_N(k) + DeltaKinvis * u_star(k)
1885 Vmath::Smul(u_star[i].size(), m_kinvis, u_star[i], 1, u_star[i], 1);
1886 Vmath::Vadd(u_star[i].size(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1887
1888 m_fields[m_velocity[i]]->FwdTrans(
1889 u_star[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1890 }
1891
1892 m_kinvis -= m_kinvis * m_KinvisPercentage / 100;
1893}
Nektar::CoupledLinearNS::SolveLinearNS
void SolveLinearNS(const Array< OneD, Array< OneD, NekDouble > > &forcing)
Definition: CoupledLinearNS.cpp:2085
Nektar::CoupledLinearNS::SetUpCoupledMatrix
void SetUpCoupledMatrix(const NekDouble lambda=0.0, const Array< OneD, Array< OneD, NekDouble > > &Advfield=NullNekDoubleArrayOfArray, bool IsLinearNSEquation=true)
Definition: CoupledLinearNS.cpp:202
Nektar::IncNavierStokes::m_kinvis
NekDouble m_kinvis
Kinematic viscosity.
Definition: IncNavierStokes.h:219
Nektar::IncNavierStokes::m_velocity
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
Definition: IncNavierStokes.h:214
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:359
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:764
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100

References Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::IncNavierStokes::m_kinvis, m_KinvisPercentage, Nektar::IncNavierStokes::m_velocity, SetUpCoupledMatrix(), Vmath::Smul(), SolveLinearNS(), and Vmath::Vadd().

Referenced by v_DoSolve().

◆ create()

static SolverUtils::EquationSystemSharedPtr Nektar::CoupledLinearNS::create ( const LibUtilities::SessionReaderSharedPtr &  pSession,
const SpatialDomains::MeshGraphSharedPtr &  pGraph 
)
inlinestatic

Creates an instance of this class.

Definition at line 92 of file CoupledLinearNS.h.

95 {
96 SolverUtils::EquationSystemSharedPtr p =
97 MemoryManager<CoupledLinearNS>::AllocateSharedPtr(pSession, pGraph);
98 p->InitObject();
99 return p;
100 }
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:166
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::SolverUtils::EquationSystemSharedPtr
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
Definition: EquationSystem.h:65

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

◆ DefineForcingTerm()

void Nektar::CoupledLinearNS::DefineForcingTerm ( void  )

Definition at line 1700 of file CoupledLinearNS.cpp.

1701{
1702 m_ForcingTerm = Array<OneD, Array<OneD, NekDouble>>(m_velocity.size());
1703 m_ForcingTerm_Coeffs =
1704 Array<OneD, Array<OneD, NekDouble>>(m_velocity.size());
1705
1706 for (size_t i = 0; i < m_velocity.size(); ++i)
1707 {
1708 m_ForcingTerm[i] = Array<OneD, NekDouble>(
1709 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1710 m_ForcingTerm_Coeffs[i] =
1711 Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetNcoeffs(), 0.0);
1712 }
1713
1714 if (m_session->DefinesFunction("ForcingTerm"))
1715 {
1716 std::vector<std::string> fieldStr;
1717 for (size_t i = 0; i < m_velocity.size(); ++i)
1718 {
1719 fieldStr.push_back(
1720 m_boundaryConditions->GetVariable(m_velocity[i]));
1721 }
1722 GetFunction("ForcingTerm")->Evaluate(fieldStr, m_ForcingTerm);
1723 for (size_t i = 0; i < m_velocity.size(); ++i)
1724 {
1725 m_fields[m_velocity[i]]->FwdTransLocalElmt(m_ForcingTerm[i],
1726 m_ForcingTerm_Coeffs[i]);
1727 }
1728 }
1729 else
1730 {
1731 cout << "'ForcingTerm' section has not been defined in the input file "
1732 "=> forcing=0"
1733 << endl;
1734 }
1735}
Nektar::CoupledLinearNS::m_ForcingTerm_Coeffs
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm_Coeffs
Definition: CoupledLinearNS.h:160
Nektar::CoupledLinearNS::m_ForcingTerm
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm
Definition: CoupledLinearNS.h:159
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:352
Nektar::SolverUtils::EquationSystem::GetNcoeffs
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:706
Nektar::SolverUtils::EquationSystem::m_boundaryConditions
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
Definition: EquationSystem.h:361
Nektar::SolverUtils::EquationSystem::GetFunction
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
Definition: EquationSystem.cpp:741

References Nektar::SolverUtils::EquationSystem::GetFunction(), Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::m_boundaryConditions, Nektar::SolverUtils::EquationSystem::m_fields, m_ForcingTerm, m_ForcingTerm_Coeffs, Nektar::SolverUtils::EquationSystem::m_session, and Nektar::IncNavierStokes::m_velocity.

Referenced by v_DoInitialise().

◆ EvaluateAdvection()

void Nektar::CoupledLinearNS::EvaluateAdvection ( const Array< OneD, const Array< OneD, NekDouble > > &  inarray,
Array< OneD, Array< OneD, NekDouble > > &  outarray,
const NekDouble  time 
)

Definition at line 1515 of file CoupledLinearNS.cpp.

1518{
1519 // evaluate convection terms
1520 EvaluateAdvectionTerms(inarray, outarray, time);
1521
1522 for (auto &x : m_forcing)
1523 {
1524 x->Apply(m_fields, outarray, outarray, time);
1525 }
1526}
Nektar::IncNavierStokes::m_forcing
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
Definition: IncNavierStokes.h:207
Nektar::IncNavierStokes::EvaluateAdvectionTerms
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Definition: IncNavierStokes.cpp:301

References Nektar::IncNavierStokes::EvaluateAdvectionTerms(), Nektar::SolverUtils::EquationSystem::m_fields, and Nektar::IncNavierStokes::m_forcing.

Referenced by v_DoInitialise().

◆ EvaluateNewtonRHS()

void Nektar::CoupledLinearNS::EvaluateNewtonRHS ( Array< OneD, Array< OneD, NekDouble > > &  Velocity,
Array< OneD, Array< OneD, NekDouble > > &  outarray 
)

Definition at line 1927 of file CoupledLinearNS.cpp.

1930{
1931 Array<OneD, Array<OneD, NekDouble>> Eval_Adv(m_velocity.size());
1932 Array<OneD, Array<OneD, NekDouble>> tmp_DerVel(m_velocity.size());
1933 Array<OneD, Array<OneD, NekDouble>> AdvTerm(m_velocity.size());
1934 Array<OneD, Array<OneD, NekDouble>> ViscTerm(m_velocity.size());
1935 Array<OneD, Array<OneD, NekDouble>> Forc(m_velocity.size());
1936
1937 for (size_t i = 0; i < m_velocity.size(); ++i)
1938 {
1939 Eval_Adv[i] = Array<OneD, NekDouble>(
1940 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1941 tmp_DerVel[i] = Array<OneD, NekDouble>(
1942 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1943
1944 AdvTerm[i] = Array<OneD, NekDouble>(
1945 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1946 ViscTerm[i] = Array<OneD, NekDouble>(
1947 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1948 Forc[i] = Array<OneD, NekDouble>(
1949 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1950 outarray[i] = Array<OneD, NekDouble>(
1951 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1952
1953 m_fields[m_velocity[i]]->PhysDeriv(i, Velocity[i], tmp_DerVel[i]);
1954
1955 Vmath::Smul(tmp_DerVel[i].size(), m_kinvis, tmp_DerVel[i], 1,
1956 tmp_DerVel[i], 1);
1957 }
1958
1959 EvaluateAdvectionTerms(Velocity, Eval_Adv, m_time);
1960
1961 for (size_t i = 0; i < m_velocity.size(); ++i)
1962 {
1963 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1964 m_fields[m_velocity[i]]->SetWaveSpace(true);
1965 m_fields[m_velocity[i]]->IProductWRTBase(Eval_Adv[i],
1966 AdvTerm[i]); //(w, (u.grad)u)
1967 m_fields[m_velocity[i]]->IProductWRTDerivBase(
1968 i, tmp_DerVel[i], ViscTerm[i]); //(grad w, grad u)
1969 m_fields[m_velocity[i]]->IProductWRTBase(m_ForcingTerm[i],
1970 Forc[i]); //(w, f)
1971 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1972
1973 Vmath::Vsub(outarray[i].size(), outarray[i], 1, AdvTerm[i], 1,
1974 outarray[i], 1);
1975 Vmath::Vsub(outarray[i].size(), outarray[i], 1, ViscTerm[i], 1,
1976 outarray[i], 1);
1977
1978 Vmath::Vadd(outarray[i].size(), outarray[i], 1, Forc[i], 1, outarray[i],
1979 1);
1980 }
1981}
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:367
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.hpp:220

References Nektar::IncNavierStokes::EvaluateAdvectionTerms(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_ForcingTerm, Nektar::IncNavierStokes::m_kinvis, Nektar::SolverUtils::EquationSystem::m_time, Nektar::IncNavierStokes::m_velocity, Vmath::Smul(), Vmath::Vadd(), and Vmath::Vsub().

Referenced by SolveSteadyNavierStokes().

◆ GenPressureExp()

const SpatialDomains::ExpansionInfoMap & Nektar::CoupledLinearNS::GenPressureExp ( const SpatialDomains::ExpansionInfoMap &  VelExp)

Definition at line 1983 of file CoupledLinearNS.cpp.

1985{
1986 SpatialDomains::ExpansionInfoMapShPtr returnval;
1987 returnval =
1988 MemoryManager<SpatialDomains::ExpansionInfoMap>::AllocateSharedPtr();
1989
1990 int nummodes;
1991
1992 for (auto &expMapIter : VelExp)
1993 {
1994 LibUtilities::BasisKeyVector BasisVec;
1995
1996 for (size_t i = 0; i < expMapIter.second->m_basisKeyVector.size(); ++i)
1997 {
1998 LibUtilities::BasisKey B = expMapIter.second->m_basisKeyVector[i];
1999 nummodes = B.GetNumModes();
2000 ASSERTL0(nummodes > 3,
2001 "Velocity polynomial space not sufficiently high (>= 4)");
2002 // Should probably set to be an orthogonal basis.
2003 LibUtilities::BasisKey newB(B.GetBasisType(), nummodes - 2,
2004 B.GetPointsKey());
2005 BasisVec.push_back(newB);
2006 }
2007
2008 // Put new expansion into list.
2009 SpatialDomains::ExpansionInfoShPtr expansionElementShPtr =
2010 MemoryManager<SpatialDomains::ExpansionInfo>::AllocateSharedPtr(
2011 expMapIter.second->m_geomShPtr, BasisVec);
2012 (*returnval)[expMapIter.first] = expansionElementShPtr;
2013 }
2014
2015 // Save expansion into graph.
2016 m_graph->SetExpansionInfo("p", returnval);
2017
2018 return *returnval;
2019}
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
Nektar::SolverUtils::EquationSystem::m_graph
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
Definition: EquationSystem.h:363
Nektar::LibUtilities::BasisKeyVector
std::vector< BasisKey > BasisKeyVector
Name for a vector of BasisKeys.
Definition: FoundationsFwd.hpp:68
Nektar::SpatialDomains::ExpansionInfoMapShPtr
std::shared_ptr< ExpansionInfoMap > ExpansionInfoMapShPtr
Definition: MeshGraph.h:143
Nektar::SpatialDomains::ExpansionInfoShPtr
std::shared_ptr< ExpansionInfo > ExpansionInfoShPtr
Definition: MeshGraph.h:140

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::BasisKey::GetBasisType(), Nektar::LibUtilities::BasisKey::GetNumModes(), Nektar::LibUtilities::BasisKey::GetPointsKey(), and Nektar::SolverUtils::EquationSystem::m_graph.

Referenced by v_InitObject().

◆ InfNorm()

void Nektar::CoupledLinearNS::InfNorm ( Array< OneD, Array< OneD, NekDouble > > &  inarray,
Array< OneD, NekDouble > &  outarray 
)

Definition at line 1895 of file CoupledLinearNS.cpp.

1897{
1898 for (size_t i = 0; i < m_velocity.size(); ++i)
1899 {
1900 outarray[i] = 0.0;
1901 for (size_t j = 0; j < inarray[i].size(); ++j)
1902 {
1903 if (inarray[i][j] > outarray[i])
1904 {
1905 outarray[i] = inarray[i][j];
1906 }
1907 }
1908 cout << "InfNorm[" << i << "] = " << outarray[i] << endl;
1909 }
1910}

References Nektar::IncNavierStokes::m_velocity.

◆ L2Norm()

void Nektar::CoupledLinearNS::L2Norm ( Array< OneD, Array< OneD, NekDouble > > &  inarray,
Array< OneD, NekDouble > &  outarray 
)

Definition at line 1912 of file CoupledLinearNS.cpp.

1914{
1915 for (size_t i = 0; i < m_velocity.size(); ++i)
1916 {
1917 outarray[i] = 0.0;
1918 for (size_t j = 0; j < inarray[i].size(); ++j)
1919 {
1920 outarray[i] += inarray[i][j] * inarray[i][j];
1921 }
1922 outarray[i] = sqrt(outarray[i]);
1923 cout << "L2Norm[" << i << "] = " << outarray[i] << endl;
1924 }
1925}
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294

References Nektar::IncNavierStokes::m_velocity, and tinysimd::sqrt().

Referenced by SolveSteadyNavierStokes().

◆ SetUpCoupledMatrix() [1/2]

void Nektar::CoupledLinearNS::SetUpCoupledMatrix ( const NekDouble  lambda,
const Array< OneD, Array< OneD, NekDouble > > &  Advfield,
bool  IsLinearNSEquation,
const int  HomogeneousMode,
CoupledSolverMatrices &  mat,
CoupledLocalToGlobalC0ContMapSharedPtr &  locToGloMap,
const NekDouble  lambda_imag = NekConstants::kNekUnsetDouble 
)
private

Generate the linearised Navier Stokes solver based on the static condensation of the interior velocity space and pressure modes. This call also allows for a Fourier mode to be specified, however if HomogeneneousMode= 0 then can be used for a standared 2D and hopefully 3D (in the future).

Set up a coupled linearised Naviers-Stokes solve in the following manner:

Consider the linearised Navier-Stokes element matrix denoted as

\[ \left [ \begin{array}{ccc} A & -D_{bnd}^T & B\\ -D_{bnd}& 0 & -D_{int}\\ \tilde{B}^T & -D_{int}^T & C \end{array} \right ] \left [ \begin{array}{c} {\bf v}_{bnd} \\ p\\ {\bf v}_{int} \end{array} \right ] = \left [ \begin{array}{c} {\bf f}_{bnd} \\ 0\\ {\bf f}_{int} \end{array} \right ] \]

where \({\bf v}_{bnd}, {\bf f}_{bnd}\) denote the degrees of freedom of the elemental velocities on the boundary of the element, \({\bf v}_{int}, {\bf f}_{int}\) denote the degrees of freedom of the elemental velocities on the interior of the element and \(p\) is the piecewise continuous pressure. The matrices have the interpretation

\[ A[n,m] = (\nabla \phi^b_n, \nu \nabla \phi^b_m) + (\phi^b_n,{\bf U \cdot \nabla} \phi^b_m) + (\phi^b_n \nabla^T {\bf U} \phi^b_m), \]

\[ B[n,m] = (\nabla \phi^b_n, \nu \nabla \phi^i_m) + (\phi^b_n,{\bf U \cdot \nabla} \phi^i_m) + (\phi^b_n \nabla^T {\bf U} \phi^i_m),\]

\[\tilde{B}^T[n,m] = (\nabla \phi^i_n, \nu \nabla \phi^b_m) + (\phi^i_n,{\bf U \cdot \nabla} \phi^b_m) + (\phi^i_n \nabla^T {\bf U} \phi^b_m) \]

\[ C[n,m] = (\nabla \phi^i_n, \nu \nabla \phi^i_m) + (\phi^i_n,{\bf U \cdot \nabla} \phi^i_m) + (\phi^i_n \nabla^T {\bf U} \phi^i_m),\]

\[ D_{bnd}[n,m] = (\psi_m,\nabla \phi^b),\]

\[ D_{int}[n,m] = (\psi_m,\nabla \phi^i) \]

where \(\psi\) is the space of pressures typically at order \(P-2\) and \(\phi\) is the velocity vector space of polynomial order \(P\). (Note the above definitions for the transpose terms shoudl be interpreted with care since we have used a tensor differential for the \(\nabla^T \) terms)

Note \(B = \tilde{B}^T\) if just considering the stokes operator and then \(C\) is also symmetric. Also note that \(A,B\) and \(C\) are block diagonal in the Oseen equation when \(\nabla^T {\bf U}\) are zero.

Since \(C\) is invertible we can premultiply the governing elemental equation by the following matrix:

\[ \left [ \begin{array}{ccc} I & 0 & -BC^{-1}\\ 0& I & D_{int}C^{-1}\\ 0 & 0 & I \end{array} \right ] \left \{ \left [ \begin{array}{ccc} A & -D_{bnd}^T & B\\ -D_{bnd}& 0 & -D_{int}\\ \tilde{B}^T & -D_{int}^T & C \end{array} \right ] \left [ \begin{array}{c} {\bf v}_{bnd} \\ p\\ {\bf v_{int}} \end{array} \right ] = \left [ \begin{array}{c} {\bf f}_{bnd} \\ 0\\ {\bf f_{int}} \end{array} \right ] \right \} \]

which if we multiply out the matrix equation we obtain

\[ \left [ \begin{array}{ccc} A - B C^{-1}\tilde{B}^T & -D_{bnd}^T +B C^{-1} D_{int}^T& 0\\ -D_{bnd}+D_{int} C^{-1} \tilde{B}^T & -D_{int} C^{-1} D_{int}^T & 0\\ \tilde{B}^T & -D_{int}^T & C \end{array} \right ] \left [ \begin{array}{c} {\bf v}_{bnd} \\ p\\ {\bf v_{int}} \end{array} \right ] = \left [ \begin{array}{c} {\bf f}_{bnd} -B C^{-1} {\bf f}_{int}\\ f_p = D_{int} C^{-1} {\bf f}_{int}\\ {\bf f_{int}} \end{array} \right ] \]

In the above equation the \({\bf v}_{int}\) degrees of freedom are decoupled and so we need to solve for the \({\bf v}_{bnd}, p\) degrees of freedom. The final step is to perform a second level of static condensation but where we will lump the mean pressure mode (or a pressure degree of freedom containing a mean component) with the velocity boundary degrees of freedom. To do we define \({\bf b} = [{\bf v}_{bnd}, p_0]\) where \(p_0\) is the mean pressure mode and \(\hat{p}\) to be the remainder of the pressure space. We now repartition the top \(2 \times 2\) block of matrices of previous matrix equation as

\[ \left [ \begin{array}{cc} \hat{A} & \hat{B}\\ \hat{C} & \hat{D} \end{array} \right ] \left [ \begin{array}{c} {\bf b} \\ \hat{p} \end{array} \right ] = \left [ \begin{array}{c} \hat{\bf f}_{bnd} \\ \hat{f}_p \end{array} \right ] \label{eqn.linNS_fac2} \]

where

\[ \hat{A}_{ij} = \left [ \begin{array}{cc} A - B C^{-1}\tilde{B}^T & [-D_{bnd}^T +B C^{-1} D_{int}^T]_{i0}\\ {[}-D_{bnd}+D_{int} C^{-1} \tilde{B}^T]_{0j} & -[D_{int} C^{-1} D_{int}^T ]_{00} \end{array} \right ] \]

\[ \hat{B}_{ij} = \left [ \begin{array}{c} [-D_{bnd}^T +B C^{-1} D_{int}^T]_{i,j+1} \\ {[} -D_{int} C^{-1} D^T_{int} ]_{0j}\end{array} \right ] \]

\[ \hat{C}_{ij} = \left [\begin{array}{cc} -D_{bnd} + D_{int} C^{-1} \tilde{B}^T, & {[} -D_{int} C^{-1} D^T_{int} ]_{i+1,0}\end{array} \right ] \]

\[ \hat{D}_{ij} = \left [\begin{array}{c} {[} -D_{int} C^{-1} D^T_{int} ]_{i+1,j+1}\end{array} \right ] \]

and

\[ fh\_{bnd} = \left [ \begin{array}{c} {\bf f}_{bnd} -B C^{-1} {\bf f}_{int}\\ {[}D_{int} C^{-1} {\bf f}_{int}]_0 \end{array}\right ] \hspace{1cm} [fh\_p_{i} = \left [ \begin{array}{c} {[}D_{int} C^{-1} {\bf f}_{iint}]_{i+1} \end{array}\right ] \]

Since the \(\hat{D}\) is decoupled and invertible we can now statically condense the previous matrix equationto decouple \({\bf b}\) from \(\hat{p}\) by solving the following system

\[ \left [ \begin{array}{cc} \hat{A} - \hat{B} \hat{D}^{-1} \hat{C} & 0 \\ \hat{C} & \hat{D} \end{array} \right ] \left [ \begin{array}{c} {\bf b} \\ \hat{p} \end{array} \right ] = \left [ \begin{array}{c} \hat{\bf f}_{bnd} - \hat{B} \hat{D}^{-1} \hat{f}_p\\ \hat{f}_p \end{array} \right ] \]

The matrix \(\hat{A} - \hat{B} \hat{D}^{-1} \hat{C}\) has to be globally assembled and solved iteratively or directly. One we obtain the solution to \({\bf b}\) we can use the second row of equation fourth matrix equation to solve for \(\hat{p}\) and finally the last row of equation second matrix equation to solve for \({\bf v}_{int}\).

Definition at line 404 of file CoupledLinearNS.cpp.

410{
411 size_t n, i, j, k;
412 size_t nel = m_fields[m_velocity[0]]->GetNumElmts();
413 size_t nvel = m_velocity.size();
414
415 // if Advfield is defined can assume it is an Oseen or LinearNS equation
416 bool AddAdvectionTerms =
417 (Advfield == NullNekDoubleArrayOfArray) ? false : true;
418 // bool AddAdvectionTerms = true; // Temporary debugging trip
419
420 // call velocity space Static condensation and fill block
421 // matrices. Need to set this up differently for Oseen and
422 // Lin NS. Ideally should make block diagonal for Stokes and
423 // Oseen problems.
424 DNekScalMatSharedPtr loc_mat;
425 StdRegions::StdExpansionSharedPtr locExp;
426 NekDouble one = 1.0;
427 size_t nint, nbndry;
428 size_t rows, cols;
429 NekDouble zero = 0.0;
430 Array<OneD, unsigned int> bmap, imap;
431
432 Array<OneD, unsigned int> nsize_bndry(nel);
433 Array<OneD, unsigned int> nsize_bndry_p1(nel);
434 Array<OneD, unsigned int> nsize_int(nel);
435 Array<OneD, unsigned int> nsize_p(nel);
436 Array<OneD, unsigned int> nsize_p_m1(nel);
437
438 size_t nz_loc;
439
440 if (HomogeneousMode) // Homogeneous mode flag
441 {
442 nz_loc = 2;
443 }
444 else
445 {
446 if (m_singleMode)
447 {
448 nz_loc = 2;
449 }
450 else
451 {
452 nz_loc = 1;
453 }
454 }
455
456 // Set up block matrix sizes -
457 for (n = 0; n < nel; ++n)
458 {
459 nsize_bndry[n] = nvel *
460 m_fields[m_velocity[0]]->GetExp(n)->NumBndryCoeffs() *
461 nz_loc;
462 nsize_bndry_p1[n] = nsize_bndry[n] + nz_loc;
463 nsize_int[n] =
464 (nvel * m_fields[m_velocity[0]]->GetExp(n)->GetNcoeffs() * nz_loc -
465 nsize_bndry[n]);
466 nsize_p[n] = m_pressure->GetExp(n)->GetNcoeffs() * nz_loc;
467 nsize_p_m1[n] = nsize_p[n] - nz_loc;
468 }
469
470 MatrixStorage blkmatStorage = eDIAGONAL;
471 DNekScalBlkMatSharedPtr pAh =
472 MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
473 nsize_bndry_p1, nsize_bndry_p1, blkmatStorage);
474 mat.m_BCinv = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
475 nsize_bndry, nsize_int, blkmatStorage);
476 mat.m_Btilde = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
477 nsize_bndry, nsize_int, blkmatStorage);
478 mat.m_Cinv = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
479 nsize_int, nsize_int, blkmatStorage);
480
481 mat.m_D_bnd = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
482 nsize_p, nsize_bndry, blkmatStorage);
483
484 mat.m_D_int = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
485 nsize_p, nsize_int, blkmatStorage);
486
487 // Final level static condensation matrices.
488 DNekScalBlkMatSharedPtr pBh =
489 MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
490 nsize_bndry_p1, nsize_p_m1, blkmatStorage);
491 DNekScalBlkMatSharedPtr pCh =
492 MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
493 nsize_p_m1, nsize_bndry_p1, blkmatStorage);
494 DNekScalBlkMatSharedPtr pDh =
495 MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nsize_p_m1, nsize_p_m1,
496 blkmatStorage);
497
498 LibUtilities::Timer timer;
499 timer.Start();
500
501 for (n = 0; n < nel; ++n)
502 {
503 nbndry = nsize_bndry[n];
504 nint = nsize_int[n];
505 k = nsize_bndry_p1[n];
506 DNekMatSharedPtr Ah =
507 MemoryManager<DNekMat>::AllocateSharedPtr(k, k, zero);
508 Array<OneD, NekDouble> Ah_data = Ah->GetPtr();
509 int AhRows = k;
510 DNekMatSharedPtr B =
511 MemoryManager<DNekMat>::AllocateSharedPtr(nbndry, nint, zero);
512 Array<OneD, NekDouble> B_data = B->GetPtr();
513 DNekMatSharedPtr C =
514 MemoryManager<DNekMat>::AllocateSharedPtr(nbndry, nint, zero);
515 Array<OneD, NekDouble> C_data = C->GetPtr();
516 DNekMatSharedPtr D =
517 MemoryManager<DNekMat>::AllocateSharedPtr(nint, nint, zero);
518 Array<OneD, NekDouble> D_data = D->GetPtr();
519
520 DNekMatSharedPtr Dbnd = MemoryManager<DNekMat>::AllocateSharedPtr(
521 nsize_p[n], nsize_bndry[n], zero);
522 DNekMatSharedPtr Dint = MemoryManager<DNekMat>::AllocateSharedPtr(
523 nsize_p[n], nsize_int[n], zero);
524
525 locExp = m_fields[m_velocity[0]]->GetExp(n);
526 locExp->GetBoundaryMap(bmap);
527 locExp->GetInteriorMap(imap);
528 StdRegions::ConstFactorMap factors;
529 factors[StdRegions::eFactorLambda] = lambda / m_kinvis;
530 LocalRegions::MatrixKey helmkey(
531 StdRegions::eHelmholtz, locExp->DetShapeType(), *locExp, factors);
532
533 size_t ncoeffs = m_fields[m_velocity[0]]->GetExp(n)->GetNcoeffs();
534 size_t nphys = m_fields[m_velocity[0]]->GetExp(n)->GetTotPoints();
535 size_t nbmap = bmap.size();
536 size_t nimap = imap.size();
537
538 Array<OneD, NekDouble> coeffs(ncoeffs);
539 Array<OneD, NekDouble> phys(nphys);
540 size_t psize = m_pressure->GetExp(n)->GetNcoeffs();
541 size_t pqsize = m_pressure->GetExp(n)->GetTotPoints();
542
543 Array<OneD, NekDouble> deriv(pqsize);
544 Array<OneD, NekDouble> pcoeffs(psize);
545
546 if (AddAdvectionTerms == false) // use static condensed managed matrices
547 {
548 // construct velocity matrices using statically
549 // condensed elemental matrices and then construct
550 // pressure matrix systems
551 DNekScalBlkMatSharedPtr CondMat;
552 CondMat = locExp->GetLocStaticCondMatrix(helmkey);
553
554 for (k = 0; k < nvel * nz_loc; ++k)
555 {
556 DNekScalMat &SubBlock = *CondMat->GetBlock(0, 0);
557 rows = SubBlock.GetRows();
558 cols = SubBlock.GetColumns();
559 for (i = 0; i < rows; ++i)
560 {
561 for (j = 0; j < cols; ++j)
562 {
563 (*Ah)(i + k * rows, j + k * cols) =
564 m_kinvis * SubBlock(i, j);
565 }
566 }
567 }
568
569 for (k = 0; k < nvel * nz_loc; ++k)
570 {
571 DNekScalMat &SubBlock = *CondMat->GetBlock(0, 1);
572 DNekScalMat &SubBlock1 = *CondMat->GetBlock(1, 0);
573 rows = SubBlock.GetRows();
574 cols = SubBlock.GetColumns();
575 for (i = 0; i < rows; ++i)
576 {
577 for (j = 0; j < cols; ++j)
578 {
579 (*B)(i + k * rows, j + k * cols) = SubBlock(i, j);
580 (*C)(i + k * rows, j + k * cols) =
581 m_kinvis * SubBlock1(j, i);
582 }
583 }
584 }
585
586 for (k = 0; k < nvel * nz_loc; ++k)
587 {
588 DNekScalMat &SubBlock = *CondMat->GetBlock(1, 1);
589 NekDouble inv_kinvis = 1.0 / m_kinvis;
590 rows = SubBlock.GetRows();
591 cols = SubBlock.GetColumns();
592 for (i = 0; i < rows; ++i)
593 {
594 for (j = 0; j < cols; ++j)
595 {
596 (*D)(i + k * rows, j + k * cols) =
597 inv_kinvis * SubBlock(i, j);
598 }
599 }
600 }
601
602 // Loop over pressure space and construct boundary block matrices.
603 for (i = 0; i < bmap.size(); ++i)
604 {
605 // Fill element with mode
606 Vmath::Zero(ncoeffs, coeffs, 1);
607 coeffs[bmap[i]] = 1.0;
608 m_fields[m_velocity[0]]->GetExp(n)->BwdTrans(coeffs, phys);
609
610 // Differentiation & Inner product wrt base.
611 for (j = 0; j < nvel; ++j)
612 {
613 if ((nz_loc == 2) && (j == 2)) // handle d/dz derivative
614 {
615 NekDouble beta = -2 * M_PI * HomogeneousMode / m_LhomZ;
616
617 Vmath::Smul(
618 m_fields[m_velocity[0]]->GetExp(n)->GetTotPoints(),
619 beta, phys, 1, deriv, 1);
620
621 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
622
623 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
624 Dbnd->GetRawPtr() +
625 ((nz_loc * j + 1) * bmap.size() + i) *
626 nsize_p[n],
627 1);
628
629 Vmath::Neg(psize, pcoeffs, 1);
630 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
631 Dbnd->GetRawPtr() +
632 ((nz_loc * j) * bmap.size() + i) *
633 nsize_p[n] +
634 psize,
635 1);
636 }
637 else
638 {
639 if (j < 2) // required for mean mode of homogeneous
640 // expansion
641 {
642 locExp->PhysDeriv(MultiRegions::DirCartesianMap[j],
643 phys, deriv);
644 m_pressure->GetExp(n)->IProductWRTBase(deriv,
645 pcoeffs);
646 // copy into column major storage.
647 for (k = 0; k < nz_loc; ++k)
648 {
649 Blas::Dcopy(
650 psize, &(pcoeffs)[0], 1,
651 Dbnd->GetRawPtr() +
652 ((nz_loc * j + k) * bmap.size() + i) *
653 nsize_p[n] +
654 k * psize,
655 1);
656 }
657 }
658 }
659 }
660 }
661
662 for (i = 0; i < imap.size(); ++i)
663 {
664 // Fill element with mode
665 Vmath::Zero(ncoeffs, coeffs, 1);
666 coeffs[imap[i]] = 1.0;
667 m_fields[m_velocity[0]]->GetExp(n)->BwdTrans(coeffs, phys);
668
669 // Differentiation & Inner product wrt base.
670 for (j = 0; j < nvel; ++j)
671 {
672 if ((nz_loc == 2) && (j == 2)) // handle d/dz derivative
673 {
674 NekDouble beta = -2 * M_PI * HomogeneousMode / m_LhomZ;
675
676 Vmath::Smul(
677 m_fields[m_velocity[0]]->GetExp(n)->GetTotPoints(),
678 beta, phys, 1, deriv, 1);
679
680 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
681
682 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
683 Dint->GetRawPtr() +
684 ((nz_loc * j + 1) * imap.size() + i) *
685 nsize_p[n],
686 1);
687
688 Vmath::Neg(psize, pcoeffs, 1);
689 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
690 Dint->GetRawPtr() +
691 ((nz_loc * j) * imap.size() + i) *
692 nsize_p[n] +
693 psize,
694 1);
695 }
696 else
697 {
698 if (j < 2) // required for mean mode of homogeneous
699 // expansion
700 {
701 // m_fields[m_velocity[0]]->GetExp(n)->PhysDeriv(j,phys,
702 // deriv);
703 locExp->PhysDeriv(MultiRegions::DirCartesianMap[j],
704 phys, deriv);
705
706 m_pressure->GetExp(n)->IProductWRTBase(deriv,
707 pcoeffs);
708
709 // copy into column major storage.
710 for (k = 0; k < nz_loc; ++k)
711 {
712 Blas::Dcopy(
713 psize, &(pcoeffs)[0], 1,
714 Dint->GetRawPtr() +
715 ((nz_loc * j + k) * imap.size() + i) *
716 nsize_p[n] +
717 k * psize,
718 1);
719 }
720 }
721 }
722 }
723 }
724 }
725 else
726 {
727 // construct velocity matrices and pressure systems at
728 // the same time resusing differential of velocity
729 // space
730
731 DNekScalMat &HelmMat =
732 *(locExp->as<LocalRegions::Expansion>()->GetLocMatrix(helmkey));
733 DNekScalMatSharedPtr MassMat;
734
735 Array<OneD, const NekDouble> HelmMat_data =
736 HelmMat.GetOwnedMatrix()->GetPtr();
737 NekDouble HelmMatScale = HelmMat.Scale();
738 int HelmMatRows = HelmMat.GetRows();
739
740 if ((lambda_imag != NekConstants::kNekUnsetDouble) && (nz_loc == 2))
741 {
742 LocalRegions::MatrixKey masskey(
743 StdRegions::eMass, locExp->DetShapeType(), *locExp);
744 MassMat = locExp->as<LocalRegions::Expansion>()->GetLocMatrix(
745 masskey);
746 }
747
748 Array<OneD, NekDouble> Advtmp;
749 Array<OneD, Array<OneD, NekDouble>> AdvDeriv(nvel * nvel);
750 // Use ExpList phys array for temporaary storage
751 Array<OneD, NekDouble> tmpphys = m_fields[0]->UpdatePhys();
752 int phys_offset = m_fields[m_velocity[0]]->GetPhys_Offset(n);
753 size_t nv;
754 int npoints = locExp->GetTotPoints();
755
756 // Calculate derivative of base flow
757 if (IsLinearNSEquation)
758 {
759 size_t nv1;
760 size_t cnt = 0;
761 AdvDeriv[0] = Array<OneD, NekDouble>(nvel * nvel * npoints);
762 for (nv = 0; nv < nvel; ++nv)
763 {
764 for (nv1 = 0; nv1 < nvel; ++nv1)
765 {
766 if (cnt < nvel * nvel - 1)
767 {
768 AdvDeriv[cnt + 1] = AdvDeriv[cnt] + npoints;
769 ++cnt;
770 }
771
772 if ((nv1 == 2) && (m_HomogeneousType == eHomogeneous1D))
773 {
774 Vmath::Zero(npoints, AdvDeriv[nv * nvel + nv1],
775 1); // dU/dz = 0
776 }
777 else
778 {
779 locExp->PhysDeriv(
780 MultiRegions::DirCartesianMap[nv1],
781 Advfield[nv] + phys_offset,
782 AdvDeriv[nv * nvel + nv1]);
783 }
784 }
785 }
786 }
787
788 for (i = 0; i < nbmap; ++i)
789 {
790
791 // Fill element with mode
792 Vmath::Zero(ncoeffs, coeffs, 1);
793 coeffs[bmap[i]] = 1.0;
794 locExp->BwdTrans(coeffs, phys);
795
796 for (k = 0; k < nvel * nz_loc; ++k)
797 {
798 for (j = 0; j < nbmap; ++j)
799 {
800 // Ah_data[i+k*nbmap +
801 // (j+k*nbmap)*AhRows] +=
802 // m_kinvis*HelmMat(bmap[i],bmap[j]);
803 Ah_data[i + k * nbmap + (j + k * nbmap) * AhRows] +=
804 m_kinvis * HelmMatScale *
805 HelmMat_data[bmap[i] + HelmMatRows * bmap[j]];
806 }
807
808 for (j = 0; j < nimap; ++j)
809 {
810 B_data[i + k * nbmap + (j + k * nimap) * nbndry] +=
811 m_kinvis * HelmMatScale *
812 HelmMat_data[bmap[i] + HelmMatRows * imap[j]];
813 }
814 }
815
816 if ((lambda_imag != NekConstants::kNekUnsetDouble) &&
817 (nz_loc == 2))
818 {
819 for (k = 0; k < nvel; ++k)
820 {
821 for (j = 0; j < nbmap; ++j)
822 {
823 Ah_data[i + 2 * k * nbmap +
824 (j + (2 * k + 1) * nbmap) * AhRows] -=
825 lambda_imag * (*MassMat)(bmap[i], bmap[j]);
826 }
827
828 for (j = 0; j < nbmap; ++j)
829 {
830 Ah_data[i + (2 * k + 1) * nbmap +
831 (j + 2 * k * nbmap) * AhRows] +=
832 lambda_imag * (*MassMat)(bmap[i], bmap[j]);
833 }
834
835 for (j = 0; j < nimap; ++j)
836 {
837 B_data[i + 2 * k * nbmap +
838 (j + (2 * k + 1) * nimap) * nbndry] -=
839 lambda_imag * (*MassMat)(bmap[i], imap[j]);
840 }
841
842 for (j = 0; j < nimap; ++j)
843 {
844 B_data[i + (2 * k + 1) * nbmap +
845 (j + 2 * k * nimap) * nbndry] +=
846 lambda_imag * (*MassMat)(bmap[i], imap[j]);
847 }
848 }
849 }
850
851 for (k = 0; k < nvel; ++k)
852 {
853 if ((nz_loc == 2) && (k == 2)) // handle d/dz derivative
854 {
855 NekDouble beta = -2 * M_PI * HomogeneousMode / m_LhomZ;
856
857 // Real Component
858 Vmath::Smul(npoints, beta, phys, 1, deriv, 1);
859
860 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
861 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
862 Dbnd->GetRawPtr() +
863 ((nz_loc * k + 1) * bmap.size() + i) *
864 nsize_p[n],
865 1);
866
867 // Imaginary Component
868 Vmath::Neg(psize, pcoeffs, 1);
869 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
870 Dbnd->GetRawPtr() +
871 ((nz_loc * k) * bmap.size() + i) *
872 nsize_p[n] +
873 psize,
874 1);
875
876 // now do advection terms
877 Vmath::Vmul(npoints, Advtmp = Advfield[k] + phys_offset,
878 1, deriv, 1, tmpphys, 1);
879
880 locExp->IProductWRTBase(tmpphys, coeffs);
881
882 // real contribution
883 for (nv = 0; nv < nvel; ++nv)
884 {
885 for (j = 0; j < nbmap; ++j)
886 {
887 Ah_data[j + 2 * nv * nbmap +
888 (i + (2 * nv + 1) * nbmap) * AhRows] +=
889 coeffs[bmap[j]];
890 }
891
892 for (j = 0; j < nimap; ++j)
893 {
894 C_data[i + (2 * nv + 1) * nbmap +
895 (j + 2 * nv * nimap) * nbndry] +=
896 coeffs[imap[j]];
897 }
898 }
899
900 Vmath::Neg(ncoeffs, coeffs, 1);
901 // imaginary contribution
902 for (nv = 0; nv < nvel; ++nv)
903 {
904 for (j = 0; j < nbmap; ++j)
905 {
906 Ah_data[j + (2 * nv + 1) * nbmap +
907 (i + 2 * nv * nbmap) * AhRows] +=
908 coeffs[bmap[j]];
909 }
910
911 for (j = 0; j < nimap; ++j)
912 {
913 C_data[i + 2 * nv * nbmap +
914 (j + (2 * nv + 1) * nimap) * nbndry] +=
915 coeffs[imap[j]];
916 }
917 }
918 }
919 else
920 {
921 if (k < 2)
922 {
923 locExp->PhysDeriv(MultiRegions::DirCartesianMap[k],
924 phys, deriv);
925 Vmath::Vmul(npoints,
926 Advtmp = Advfield[k] + phys_offset, 1,
927 deriv, 1, tmpphys, 1);
928 locExp->IProductWRTBase(tmpphys, coeffs);
929
930 for (nv = 0; nv < nvel * nz_loc; ++nv)
931 {
932 for (j = 0; j < nbmap; ++j)
933 {
934 Ah_data[j + nv * nbmap +
935 (i + nv * nbmap) * AhRows] +=
936 coeffs[bmap[j]];
937 }
938
939 for (j = 0; j < nimap; ++j)
940 {
941 C_data[i + nv * nbmap +
942 (j + nv * nimap) * nbndry] +=
943 coeffs[imap[j]];
944 }
945 }
946
947 // copy into column major storage.
948 m_pressure->GetExp(n)->IProductWRTBase(deriv,
949 pcoeffs);
950 for (j = 0; j < nz_loc; ++j)
951 {
952 Blas::Dcopy(
953 psize, &(pcoeffs)[0], 1,
954 Dbnd->GetRawPtr() +
955 ((nz_loc * k + j) * bmap.size() + i) *
956 nsize_p[n] +
957 j * psize,
958 1);
959 }
960 }
961 }
962
963 if (IsLinearNSEquation)
964 {
965 for (nv = 0; nv < nvel; ++nv)
966 {
967 // u' . Grad U terms
968 Vmath::Vmul(npoints, phys, 1,
969 AdvDeriv[k * nvel + nv], 1, tmpphys, 1);
970 locExp->IProductWRTBase(tmpphys, coeffs);
971
972 for (size_t n1 = 0; n1 < nz_loc; ++n1)
973 {
974 for (j = 0; j < nbmap; ++j)
975 {
976 Ah_data[j + (k * nz_loc + n1) * nbmap +
977 (i + (nv * nz_loc + n1) * nbmap) *
978 AhRows] += coeffs[bmap[j]];
979 }
980
981 for (j = 0; j < nimap; ++j)
982 {
983 C_data[i + (nv * nz_loc + n1) * nbmap +
984 (j + (k * nz_loc + n1) * nimap) *
985 nbndry] += coeffs[imap[j]];
986 }
987 }
988 }
989 }
990 }
991 }
992
993 for (i = 0; i < nimap; ++i)
994 {
995 // Fill element with mode
996 Vmath::Zero(ncoeffs, coeffs, 1);
997 coeffs[imap[i]] = 1.0;
998 locExp->BwdTrans(coeffs, phys);
999
1000 for (k = 0; k < nvel * nz_loc; ++k)
1001 {
1002 for (j = 0; j < nbmap; ++j) // C set up as transpose
1003 {
1004 C_data[j + k * nbmap + (i + k * nimap) * nbndry] +=
1005 m_kinvis * HelmMatScale *
1006 HelmMat_data[imap[i] + HelmMatRows * bmap[j]];
1007 }
1008
1009 for (j = 0; j < nimap; ++j)
1010 {
1011 D_data[i + k * nimap + (j + k * nimap) * nint] +=
1012 m_kinvis * HelmMatScale *
1013 HelmMat_data[imap[i] + HelmMatRows * imap[j]];
1014 }
1015 }
1016
1017 if ((lambda_imag != NekConstants::kNekUnsetDouble) &&
1018 (nz_loc == 2))
1019 {
1020 for (k = 0; k < nvel; ++k)
1021 {
1022 for (j = 0; j < nbmap; ++j) // C set up as transpose
1023 {
1024 C_data[j + 2 * k * nbmap +
1025 (i + (2 * k + 1) * nimap) * nbndry] +=
1026 lambda_imag * (*MassMat)(bmap[j], imap[i]);
1027 }
1028
1029 for (j = 0; j < nbmap; ++j) // C set up as transpose
1030 {
1031 C_data[j + (2 * k + 1) * nbmap +
1032 (i + 2 * k * nimap) * nbndry] -=
1033 lambda_imag * (*MassMat)(bmap[j], imap[i]);
1034 }
1035
1036 for (j = 0; j < nimap; ++j)
1037 {
1038 D_data[i + 2 * k * nimap +
1039 (j + (2 * k + 1) * nimap) * nint] -=
1040 lambda_imag * (*MassMat)(imap[i], imap[j]);
1041 }
1042
1043 for (j = 0; j < nimap; ++j)
1044 {
1045 D_data[i + (2 * k + 1) * nimap +
1046 (j + 2 * k * nimap) * nint] +=
1047 lambda_imag * (*MassMat)(imap[i], imap[j]);
1048 }
1049 }
1050 }
1051
1052 for (k = 0; k < nvel; ++k)
1053 {
1054 if ((nz_loc == 2) && (k == 2)) // handle d/dz derivative
1055 {
1056 NekDouble beta = -2 * M_PI * HomogeneousMode / m_LhomZ;
1057
1058 // Real Component
1059 Vmath::Smul(npoints, beta, phys, 1, deriv, 1);
1060 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
1061 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
1062 Dint->GetRawPtr() +
1063 ((nz_loc * k + 1) * imap.size() + i) *
1064 nsize_p[n],
1065 1);
1066 // Imaginary Component
1067 Vmath::Neg(psize, pcoeffs, 1);
1068 Blas::Dcopy(psize, &(pcoeffs)[0], 1,
1069 Dint->GetRawPtr() +
1070 ((nz_loc * k) * imap.size() + i) *
1071 nsize_p[n] +
1072 psize,
1073 1);
1074
1075 // Advfield[k] *d/dx_k to all velocity
1076 // components on diagonal
1077 Vmath::Vmul(npoints, Advtmp = Advfield[k] + phys_offset,
1078 1, deriv, 1, tmpphys, 1);
1079 locExp->IProductWRTBase(tmpphys, coeffs);
1080
1081 // Real Components
1082 for (nv = 0; nv < nvel; ++nv)
1083 {
1084 for (j = 0; j < nbmap; ++j)
1085 {
1086 B_data[j + 2 * nv * nbmap +
1087 (i + (2 * nv + 1) * nimap) * nbndry] +=
1088 coeffs[bmap[j]];
1089 }
1090
1091 for (j = 0; j < nimap; ++j)
1092 {
1093 D_data[j + 2 * nv * nimap +
1094 (i + (2 * nv + 1) * nimap) * nint] +=
1095 coeffs[imap[j]];
1096 }
1097 }
1098 Vmath::Neg(ncoeffs, coeffs, 1);
1099 // Imaginary
1100 for (nv = 0; nv < nvel; ++nv)
1101 {
1102 for (j = 0; j < nbmap; ++j)
1103 {
1104 B_data[j + (2 * nv + 1) * nbmap +
1105 (i + 2 * nv * nimap) * nbndry] +=
1106 coeffs[bmap[j]];
1107 }
1108
1109 for (j = 0; j < nimap; ++j)
1110 {
1111 D_data[j + (2 * nv + 1) * nimap +
1112 (i + 2 * nv * nimap) * nint] +=
1113 coeffs[imap[j]];
1114 }
1115 }
1116 }
1117 else
1118 {
1119 if (k < 2)
1120 {
1121 // Differentiation & Inner product wrt base.
1122 // locExp->PhysDeriv(k,phys, deriv);
1123 locExp->PhysDeriv(MultiRegions::DirCartesianMap[k],
1124 phys, deriv);
1125 Vmath::Vmul(npoints,
1126 Advtmp = Advfield[k] + phys_offset, 1,
1127 deriv, 1, tmpphys, 1);
1128 locExp->IProductWRTBase(tmpphys, coeffs);
1129
1130 for (nv = 0; nv < nvel * nz_loc; ++nv)
1131 {
1132 for (j = 0; j < nbmap; ++j)
1133 {
1134 B_data[j + nv * nbmap +
1135 (i + nv * nimap) * nbndry] +=
1136 coeffs[bmap[j]];
1137 }
1138
1139 for (j = 0; j < nimap; ++j)
1140 {
1141 D_data[j + nv * nimap +
1142 (i + nv * nimap) * nint] +=
1143 coeffs[imap[j]];
1144 }
1145 }
1146 // copy into column major storage.
1147 m_pressure->GetExp(n)->IProductWRTBase(deriv,
1148 pcoeffs);
1149 for (j = 0; j < nz_loc; ++j)
1150 {
1151 Blas::Dcopy(
1152 psize, &(pcoeffs)[0], 1,
1153 Dint->GetRawPtr() +
1154 ((nz_loc * k + j) * imap.size() + i) *
1155 nsize_p[n] +
1156 j * psize,
1157 1);
1158 }
1159 }
1160 }
1161
1162 if (IsLinearNSEquation)
1163 {
1164 size_t n1;
1165 for (nv = 0; nv < nvel; ++nv)
1166 {
1167 // u'.Grad U terms
1168 Vmath::Vmul(npoints, phys, 1,
1169 AdvDeriv[k * nvel + nv], 1, tmpphys, 1);
1170 locExp->IProductWRTBase(tmpphys, coeffs);
1171
1172 for (n1 = 0; n1 < nz_loc; ++n1)
1173 {
1174 for (j = 0; j < nbmap; ++j)
1175 {
1176 B_data[j + (k * nz_loc + n1) * nbmap +
1177 (i + (nv * nz_loc + n1) * nimap) *
1178 nbndry] += coeffs[bmap[j]];
1179 }
1180
1181 for (j = 0; j < nimap; ++j)
1182 {
1183 D_data[j + (k * nz_loc + n1) * nimap +
1184 (i + (nv * nz_loc + n1) * nimap) *
1185 nint] += coeffs[imap[j]];
1186 }
1187 }
1188 }
1189 }
1190 }
1191 }
1192
1193 D->Invert();
1194 (*B) = (*B) * (*D);
1195
1196 // perform (*Ah) = (*Ah) - (*B)*(*C) but since size of
1197 // Ah is larger than (*B)*(*C) easier to call blas
1198 // directly
1199 Blas::Dgemm('N', 'T', B->GetRows(), C->GetRows(), B->GetColumns(),
1200 -1.0, B->GetRawPtr(), B->GetRows(), C->GetRawPtr(),
1201 C->GetRows(), 1.0, Ah->GetRawPtr(), Ah->GetRows());
1202 }
1203
1204 mat.m_BCinv->SetBlock(
1205 n, n,
1206 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, B));
1207 mat.m_Btilde->SetBlock(
1208 n, n,
1209 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, C));
1210 mat.m_Cinv->SetBlock(
1211 n, n,
1212 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, D));
1213 mat.m_D_bnd->SetBlock(
1214 n, n,
1215 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, Dbnd));
1216 mat.m_D_int->SetBlock(
1217 n, n,
1218 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, Dint));
1219
1220 // Do matrix manipulations and get final set of block matries
1221 // reset boundary to put mean mode into boundary system.
1222
1223 DNekMatSharedPtr Cinv, BCinv, Btilde;
1224 DNekMat DintCinvDTint, BCinvDTint_m_DTbnd, DintCinvBTtilde_m_Dbnd;
1225
1226 Cinv = D;
1227 BCinv = B;
1228 Btilde = C;
1229
1230 DintCinvDTint = (*Dint) * (*Cinv) * Transpose(*Dint);
1231 BCinvDTint_m_DTbnd = (*BCinv) * Transpose(*Dint) - Transpose(*Dbnd);
1232
1233 // This could be transpose of BCinvDint in some cases
1234 DintCinvBTtilde_m_Dbnd =
1235 (*Dint) * (*Cinv) * Transpose(*Btilde) - (*Dbnd);
1236
1237 // Set up final set of matrices.
1238 DNekMatSharedPtr Bh = MemoryManager<DNekMat>::AllocateSharedPtr(
1239 nsize_bndry_p1[n], nsize_p_m1[n], zero);
1240 DNekMatSharedPtr Ch = MemoryManager<DNekMat>::AllocateSharedPtr(
1241 nsize_p_m1[n], nsize_bndry_p1[n], zero);
1242 DNekMatSharedPtr Dh = MemoryManager<DNekMat>::AllocateSharedPtr(
1243 nsize_p_m1[n], nsize_p_m1[n], zero);
1244 Array<OneD, NekDouble> Dh_data = Dh->GetPtr();
1245
1246 // Copy matrices into final structures.
1247 size_t n1, n2;
1248 for (n1 = 0; n1 < nz_loc; ++n1)
1249 {
1250 for (i = 0; i < psize - 1; ++i)
1251 {
1252 for (n2 = 0; n2 < nz_loc; ++n2)
1253 {
1254 for (j = 0; j < psize - 1; ++j)
1255 {
1256 //(*Dh)(i+n1*(psize-1),j+n2*(psize-1)) =
1257 //-DintCinvDTint(i+1+n1*psize,j+1+n2*psize);
1258 Dh_data[(i + n1 * (psize - 1)) +
1259 (j + n2 * (psize - 1)) * Dh->GetRows()] =
1260 -DintCinvDTint(i + 1 + n1 * psize,
1261 j + 1 + n2 * psize);
1262 }
1263 }
1264 }
1265 }
1266
1267 for (n1 = 0; n1 < nz_loc; ++n1)
1268 {
1269 for (i = 0; i < nsize_bndry_p1[n] - nz_loc; ++i)
1270 {
1271 (*Ah)(i, nsize_bndry_p1[n] - nz_loc + n1) =
1272 BCinvDTint_m_DTbnd(i, n1 * psize);
1273 (*Ah)(nsize_bndry_p1[n] - nz_loc + n1, i) =
1274 DintCinvBTtilde_m_Dbnd(n1 * psize, i);
1275 }
1276 }
1277
1278 for (n1 = 0; n1 < nz_loc; ++n1)
1279 {
1280 for (n2 = 0; n2 < nz_loc; ++n2)
1281 {
1282 (*Ah)(nsize_bndry_p1[n] - nz_loc + n1,
1283 nsize_bndry_p1[n] - nz_loc + n2) =
1284 -DintCinvDTint(n1 * psize, n2 * psize);
1285 }
1286 }
1287
1288 for (n1 = 0; n1 < nz_loc; ++n1)
1289 {
1290 for (j = 0; j < psize - 1; ++j)
1291 {
1292 for (i = 0; i < nsize_bndry_p1[n] - nz_loc; ++i)
1293 {
1294 (*Bh)(i, j + n1 * (psize - 1)) =
1295 BCinvDTint_m_DTbnd(i, j + 1 + n1 * psize);
1296 (*Ch)(j + n1 * (psize - 1), i) =
1297 DintCinvBTtilde_m_Dbnd(j + 1 + n1 * psize, i);
1298 }
1299 }
1300 }
1301
1302 for (n1 = 0; n1 < nz_loc; ++n1)
1303 {
1304 for (n2 = 0; n2 < nz_loc; ++n2)
1305 {
1306 for (j = 0; j < psize - 1; ++j)
1307 {
1308 (*Bh)(nsize_bndry_p1[n] - nz_loc + n1,
1309 j + n2 * (psize - 1)) =
1310 -DintCinvDTint(n1 * psize, j + 1 + n2 * psize);
1311 (*Ch)(j + n2 * (psize - 1),
1312 nsize_bndry_p1[n] - nz_loc + n1) =
1313 -DintCinvDTint(j + 1 + n2 * psize, n1 * psize);
1314 }
1315 }
1316 }
1317
1318 // Do static condensation
1319 Dh->Invert();
1320 (*Bh) = (*Bh) * (*Dh);
1321 //(*Ah) = (*Ah) - (*Bh)*(*Ch);
1322 Blas::Dgemm('N', 'N', Bh->GetRows(), Ch->GetColumns(), Bh->GetColumns(),
1323 -1.0, Bh->GetRawPtr(), Bh->GetRows(), Ch->GetRawPtr(),
1324 Ch->GetRows(), 1.0, Ah->GetRawPtr(), Ah->GetRows());
1325
1326 // Set matrices for later inversion. Probably do not need to be
1327 // attached to class
1328 pAh->SetBlock(
1329 n, n,
1330 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, Ah));
1331 pBh->SetBlock(
1332 n, n,
1333 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, Bh));
1334 pCh->SetBlock(
1335 n, n,
1336 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, Ch));
1337 pDh->SetBlock(
1338 n, n,
1339 loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, Dh));
1340 }
1341 timer.Stop();
1342 cout << "Matrix Setup Costs: " << timer.TimePerTest(1) << endl;
1343
1344 timer.Start();
1345 // Set up global coupled boundary solver.
1346 // This is a key to define the solution matrix type
1347 // currently we are giving it a argument of eLInearAdvectionReaction
1348 // since this then makes the matrix storage of type eFull
1349 MultiRegions::GlobalLinSysKey key(StdRegions::eLinearAdvectionReaction,
1350 locToGloMap);
1351 mat.m_CoupledBndSys =
1352 MemoryManager<MultiRegions::GlobalLinSysDirectStaticCond>::
1353 AllocateSharedPtr(key, m_fields[0], pAh, pBh, pCh, pDh,
1354 locToGloMap);
1355 mat.m_CoupledBndSys->Initialise(locToGloMap);
1356}
Nektar::IncNavierStokes::m_pressure
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
Definition: IncNavierStokes.h:217
Nektar::SolverUtils::EquationSystem::m_singleMode
bool m_singleMode
Flag to determine if single homogeneous mode is used.
Definition: EquationSystem.h:400
Nektar::SolverUtils::EquationSystem::m_HomogeneousType
enum HomogeneousType m_HomogeneousType
Definition: EquationSystem.h:453
Nektar::SolverUtils::EquationSystem::m_LhomZ
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
Definition: EquationSystem.h:457
Blas::Dcopy
static void Dcopy(const int &n, const double *x, const int &incx, double *y, const int &incy)
BLAS level 1: Copy x to y.
Definition: Blas.hpp:128
Blas::Dgemm
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition: Blas.hpp:383
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:59
Nektar::MultiRegions::DirCartesianMap
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
Nektar::NekConstants::kNekUnsetDouble
static const NekDouble kNekUnsetDouble
Definition: NektarUnivConsts.hpp:44
Nektar::StdRegions::StdExpansionSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1768
Nektar::StdRegions::ConstFactorMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:402
Nektar::VarcoeffHashingTest::factors
StdRegions::ConstFactorMap factors
Definition: TestVarcoeffHashing.cpp:51
Nektar::DNekScalMat
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
Definition: NekMatrixFwd.hpp:65
Nektar::DNekScalMatSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:66
Nektar::DNekMat
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:50
Nektar::NullNekDoubleArrayOfArray
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
Definition: SharedArray.hpp:889
Nektar::DNekScalBlkMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:235
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::LibUtilities::beta, Blas::Dcopy(), Blas::Dgemm(), Nektar::MultiRegions::DirCartesianMap, Nektar::eDIAGONAL, Nektar::StdRegions::eFactorLambda, Nektar::StdRegions::eHelmholtz, Nektar::SolverUtils::EquationSystem::eHomogeneous1D, Nektar::StdRegions::eLinearAdvectionReaction, Nektar::StdRegions::eMass, Nektar::VarcoeffHashingTest::factors, Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::NekConstants::kNekUnsetDouble, Nektar::coupledSolverMatrices::m_BCinv, Nektar::coupledSolverMatrices::m_Btilde, Nektar::coupledSolverMatrices::m_Cinv, Nektar::coupledSolverMatrices::m_CoupledBndSys, Nektar::coupledSolverMatrices::m_D_bnd, Nektar::coupledSolverMatrices::m_D_int, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_HomogeneousType, Nektar::IncNavierStokes::m_kinvis, Nektar::SolverUtils::EquationSystem::m_LhomZ, Nektar::IncNavierStokes::m_pressure, Nektar::SolverUtils::EquationSystem::m_singleMode, Nektar::IncNavierStokes::m_velocity, Vmath::Neg(), Nektar::NullNekDoubleArrayOfArray, Vmath::Smul(), Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), Nektar::LibUtilities::Timer::TimePerTest(), Nektar::Transpose(), Vmath::Vmul(), and Vmath::Zero().

◆ SetUpCoupledMatrix() [2/2]

void Nektar::CoupledLinearNS::SetUpCoupledMatrix ( const NekDouble  lambda = 0.0,
const Array< OneD, Array< OneD, NekDouble > > &  Advfield = NullNekDoubleArrayOfArray,
bool  IsLinearNSEquation = true 
)

Generate the linearised Navier Stokes solver based on the static condensation of the interior velocity space and pressure modes.

Set up a coupled linearised Naviers-Stokes solve. Primarily a wrapper fuction around the full mode by mode version of SetUpCoupledMatrix.

Definition at line 202 of file CoupledLinearNS.cpp.

205{
206
207 int nz;
208 if (m_singleMode)
209 {
210
211 NekDouble lambda_imag;
212
213 // load imaginary component of any potential shift
214 // Probably should be called from DriverArpack but not yet
215 // clear how to do this
216 m_session->LoadParameter("imagShift", lambda_imag,
217 NekConstants::kNekUnsetDouble);
218 nz = 1;
219 m_mat = Array<OneD, CoupledSolverMatrices>(nz);
220
221 ASSERTL1(m_npointsZ <= 2, "Only expected a maxmimum of two planes in "
222 "single mode linear NS solver");
223
224 if (m_zeroMode)
225 {
226 SetUpCoupledMatrix(lambda, Advfield, IsLinearNSEquation, 0,
227 m_mat[0], m_locToGloMap[0], lambda_imag);
228 }
229 else
230 {
231 NekDouble beta = 2 * M_PI / m_LhomZ;
232 NekDouble lam = lambda + m_kinvis * beta * beta;
233
234 SetUpCoupledMatrix(lam, Advfield, IsLinearNSEquation, 1, m_mat[0],
235 m_locToGloMap[0], lambda_imag);
236 }
237 }
238 else
239 {
240 int n;
241 if (m_npointsZ > 1)
242 {
243 nz = m_npointsZ / 2;
244 }
245 else
246 {
247 nz = 1;
248 }
249
250 m_mat = Array<OneD, CoupledSolverMatrices>(nz);
251
252 // mean mode or 2D mode.
253 SetUpCoupledMatrix(lambda, Advfield, IsLinearNSEquation, 0, m_mat[0],
254 m_locToGloMap[0]);
255
256 for (n = 1; n < nz; ++n)
257 {
258 NekDouble beta = 2 * M_PI * n / m_LhomZ;
259
260 NekDouble lam = lambda + m_kinvis * beta * beta;
261
262 SetUpCoupledMatrix(lam, Advfield, IsLinearNSEquation, n, m_mat[n],
263 m_locToGloMap[n]);
264 }
265 }
266}
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
Nektar::CoupledLinearNS::m_locToGloMap
Array< OneD, CoupledLocalToGlobalC0ContMapSharedPtr > m_locToGloMap
Definition: CoupledLinearNS.h:162
Nektar::CoupledLinearNS::m_mat
Array< OneD, CoupledSolverMatrices > m_mat
Definition: CoupledLinearNS.h:187
Nektar::SolverUtils::EquationSystem::m_npointsZ
int m_npointsZ
number of points in Z direction (if homogeneous)
Definition: EquationSystem.h:461

References ASSERTL1, Nektar::LibUtilities::beta, Nektar::NekConstants::kNekUnsetDouble, Nektar::IncNavierStokes::m_kinvis, Nektar::SolverUtils::EquationSystem::m_LhomZ, m_locToGloMap, m_mat, Nektar::SolverUtils::EquationSystem::m_npointsZ, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_singleMode, m_zeroMode, and SetUpCoupledMatrix().

Referenced by Continuation(), SetUpCoupledMatrix(), SolveSteadyNavierStokes(), SolveUnsteadyStokesSystem(), and v_DoInitialise().

◆ Solve()

void Nektar::CoupledLinearNS::Solve ( void  )

Definition at line 1670 of file CoupledLinearNS.cpp.

1671{
1672 const size_t ncmpt = m_velocity.size();
1673 Array<OneD, Array<OneD, NekDouble>> forcing_phys(ncmpt);
1674 Array<OneD, Array<OneD, NekDouble>> forcing(ncmpt);
1675
1676 for (size_t i = 0; i < ncmpt; ++i)
1677 {
1678 forcing_phys[i] =
1679 Array<OneD, NekDouble>(m_fields[m_velocity[0]]->GetNpoints(), 0.0);
1680 forcing[i] =
1681 Array<OneD, NekDouble>(m_fields[m_velocity[0]]->GetNcoeffs(), 0.0);
1682 }
1683
1684 for (auto &x : m_forcing)
1685 {
1686 const NekDouble time = 0;
1687 x->Apply(m_fields, forcing_phys, forcing_phys, time);
1688 }
1689 for (size_t i = 0; i < ncmpt; ++i)
1690 {
1691 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1692 m_fields[i]->SetWaveSpace(true);
1693 m_fields[i]->IProductWRTBase(forcing_phys[i], forcing[i]);
1694 m_fields[i]->SetWaveSpace(waveSpace);
1695 }
1696
1697 SolveLinearNS(forcing);
1698}
Nektar::SolverUtils::EquationSystem::GetNpoints
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:774

References Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::GetNpoints(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::IncNavierStokes::m_forcing, Nektar::IncNavierStokes::m_velocity, and SolveLinearNS().

Referenced by v_DoInitialise(), and v_DoSolve().

◆ SolveLinearNS() [1/2]

void Nektar::CoupledLinearNS::SolveLinearNS ( Array< OneD, Array< OneD, NekDouble > > &  forcing,
Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
MultiRegions::ExpListSharedPtr &  pressure,
const int  HomogeneousMode = 0 
)

Definition at line 2132 of file CoupledLinearNS.cpp.

2136{
2137 size_t i, j, k, n, cnt, cnt1;
2138 size_t nbnd, nint, offset;
2139 size_t nvel = m_velocity.size();
2140 size_t nel = fields[0]->GetNumElmts();
2141 Array<OneD, unsigned int> bmap, imap;
2142
2143 Array<OneD, NekDouble> f_bnd(m_mat[mode].m_BCinv->GetRows());
2144 NekVector<NekDouble> F_bnd(f_bnd.size(), f_bnd, eWrapper);
2145 Array<OneD, NekDouble> f_int(m_mat[mode].m_BCinv->GetColumns());
2146 NekVector<NekDouble> F_int(f_int.size(), f_int, eWrapper);
2147
2148 size_t nz_loc;
2149 int nplanecoeffs =
2150 fields[m_velocity[0]]
2151 ->GetNcoeffs(); // this is fine since we pass the nplane coeff data.
2152
2153 if (mode) // Homogeneous mode flag
2154 {
2155 nz_loc = 2;
2156 }
2157 else
2158 {
2159 if (m_singleMode)
2160 {
2161 nz_loc = 2;
2162 }
2163 else
2164 {
2165 nz_loc = 1;
2166 if (m_HomogeneousType == eHomogeneous1D)
2167 {
2168 Array<OneD, NekDouble> tmp;
2169 // Zero fields to set complex mode to zero;
2170 for (i = 0; i < fields.size(); ++i)
2171 {
2172 Vmath::Zero(nplanecoeffs,
2173 tmp = fields[i]->UpdateCoeffs() + nplanecoeffs,
2174 1);
2175 }
2176 Vmath::Zero(2 * pressure->GetNcoeffs(),
2177 pressure->UpdateCoeffs(), 1);
2178 }
2179 }
2180 }
2181
2182 for (k = 0; k < nvel; ++k)
2183 {
2184 MultiRegions::ContFieldSharedPtr cfield =
2185 std::dynamic_pointer_cast<MultiRegions::ContField>(fields[k]);
2186
2187 Array<OneD, NekDouble> sign =
2188 cfield->GetLocalToGlobalMap()->GetBndCondCoeffsToLocalCoeffsSign();
2189 const Array<OneD, const int> map =
2190 cfield->GetLocalToGlobalMap()->GetBndCondCoeffsToLocalCoeffsMap();
2191
2192 // Add weak boundary conditions to forcing
2193 const Array<OneD, SpatialDomains::BoundaryConditionShPtr> bndConds =
2194 fields[k]->GetBndConditions();
2195 Array<OneD, const MultiRegions::ExpListSharedPtr> bndCondExp;
2196
2197 if (m_HomogeneousType == eHomogeneous1D)
2198 {
2199 bndCondExp =
2200 m_fields[k]->GetPlane(2 * mode)->GetBndCondExpansions();
2201 }
2202 else
2203 {
2204 bndCondExp = m_fields[k]->GetBndCondExpansions();
2205 }
2206
2207 for (n = 0; n < nz_loc; ++n)
2208 {
2209 int bndcnt = 0;
2210 for (i = 0; i < bndCondExp.size(); ++i)
2211 {
2212 const Array<OneD, const NekDouble> bndcoeffs =
2213 bndCondExp[i]->GetCoeffs();
2214 size_t nCoeffs = (bndCondExp[i])->GetNcoeffs();
2215
2216 cnt = 0;
2217 if (bndConds[i]->GetBoundaryConditionType() ==
2218 SpatialDomains::eNeumann ||
2219 bndConds[i]->GetBoundaryConditionType() ==
2220 SpatialDomains::eRobin)
2221 {
2222 if (m_locToGloMap[mode]->GetSignChange())
2223 {
2224 for (j = 0; j < nCoeffs; j++)
2225 {
2226 forcing[k][n * nplanecoeffs + map[bndcnt + j]] +=
2227 sign[bndcnt + j] * bndcoeffs[j];
2228 }
2229 }
2230 else
2231 {
2232 for (j = 0; j < nCoeffs; j++)
2233 {
2234 forcing[k][n * nplanecoeffs + map[bndcnt + j]] +=
2235 bndcoeffs[j];
2236 }
2237 }
2238 }
2239
2240 bndcnt += bndCondExp[i]->GetNcoeffs();
2241 }
2242 }
2243 }
2244
2245 Array<OneD, NekDouble> bnd(m_locToGloMap[mode]->GetNumLocalCoeffs(), 0.0);
2246
2247 // Construct f_bnd and f_int and fill in bnd from inarray
2248 // (with Dirichlet BCs imposed)
2249 int bndoffset = 0;
2250 cnt = cnt1 = 0;
2251 for (i = 0; i < nel; ++i) // loop over elements
2252 {
2253 fields[m_velocity[0]]->GetExp(i)->GetBoundaryMap(bmap);
2254 fields[m_velocity[0]]->GetExp(i)->GetInteriorMap(imap);
2255 nbnd = bmap.size();
2256 nint = imap.size();
2257 offset = fields[m_velocity[0]]->GetCoeff_Offset(i);
2258
2259 for (j = 0; j < nvel; ++j) // loop over velocity fields
2260 {
2261 Array<OneD, NekDouble> incoeffs = fields[j]->UpdateCoeffs();
2262
2263 for (n = 0; n < nz_loc; ++n)
2264 {
2265 for (k = 0; k < nbnd; ++k)
2266 {
2267 f_bnd[cnt + k] =
2268 forcing[j][n * nplanecoeffs + offset + bmap[k]];
2269
2270 bnd[bndoffset + (n + j * nz_loc) * nbnd + k] =
2271 incoeffs[n * nplanecoeffs + offset + bmap[k]];
2272 }
2273 for (k = 0; k < nint; ++k)
2274 {
2275 f_int[cnt1 + k] =
2276 forcing[j][n * nplanecoeffs + offset + imap[k]];
2277 }
2278
2279 cnt += nbnd;
2280 cnt1 += nint;
2281 }
2282 }
2283 bndoffset +=
2284 nvel * nz_loc * nbnd + nz_loc * (pressure->GetExp(i)->GetNcoeffs());
2285 }
2286
2287 Array<OneD, NekDouble> f_p(m_mat[mode].m_D_int->GetRows());
2288 NekVector<NekDouble> F_p(f_p.size(), f_p, eWrapper);
2289 NekVector<NekDouble> F_p_tmp(m_mat[mode].m_Cinv->GetRows());
2290
2291 // fbnd does not currently hold the pressure mean
2292 F_bnd = F_bnd - (*m_mat[mode].m_BCinv) * F_int;
2293 F_p_tmp = (*m_mat[mode].m_Cinv) * F_int;
2294 F_p = (*m_mat[mode].m_D_int) * F_p_tmp;
2295
2296 // construct inner forcing
2297 Array<OneD, NekDouble> fh_bnd(m_locToGloMap[mode]->GetNumLocalCoeffs(),
2298 0.0);
2299
2300 offset = cnt = 0;
2301 for (i = 0; i < nel; ++i)
2302 {
2303 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2304
2305 for (j = 0; j < nvel; ++j)
2306 {
2307 for (k = 0; k < nbnd; ++k)
2308 {
2309 fh_bnd[offset + j * nbnd + k] = f_bnd[cnt + k];
2310 }
2311 cnt += nbnd;
2312 }
2313
2314 nint = pressure->GetExp(i)->GetNcoeffs();
2315 offset += nvel * nbnd + nint * nz_loc;
2316 }
2317
2318 offset = cnt1 = 0;
2319 for (i = 0; i < nel; ++i)
2320 {
2321 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2322 nint = pressure->GetExp(i)->GetNcoeffs();
2323
2324 for (n = 0; n < nz_loc; ++n)
2325 {
2326 for (j = 0; j < nint; ++j)
2327 {
2328 fh_bnd[offset + nvel * nbnd + n * nint + j] = f_p[cnt1 + j];
2329 }
2330 cnt1 += nint;
2331 }
2332 offset += nvel * nbnd + nz_loc * nint;
2333 }
2334 m_mat[mode].m_CoupledBndSys->Solve(fh_bnd, bnd, m_locToGloMap[mode]);
2335
2336 // unpack pressure and velocity boundary systems.
2337 offset = cnt = 0;
2338 int totpcoeffs = pressure->GetNcoeffs();
2339 Array<OneD, NekDouble> p_coeffs = pressure->UpdateCoeffs();
2340 for (i = 0; i < nel; ++i)
2341 {
2342 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2343 nint = pressure->GetExp(i)->GetNcoeffs();
2344 for (j = 0; j < nvel; ++j)
2345 {
2346 for (k = 0; k < nbnd; ++k)
2347 {
2348 f_bnd[cnt + k] = bnd[offset + j * nbnd + k];
2349 }
2350 cnt += nbnd;
2351 }
2352 offset += nvel * nbnd + nint * nz_loc;
2353 }
2354
2355 pressure->SetPhysState(false);
2356
2357 offset = cnt = cnt1 = 0;
2358 for (i = 0; i < nel; ++i)
2359 {
2360 nint = pressure->GetExp(i)->GetNcoeffs();
2361 nbnd = fields[0]->GetExp(i)->NumBndryCoeffs();
2362 cnt1 = pressure->GetCoeff_Offset(i);
2363
2364 for (n = 0; n < nz_loc; ++n)
2365 {
2366 for (j = 0; j < nint; ++j)
2367 {
2368 p_coeffs[n * totpcoeffs + cnt1 + j] = f_p[cnt + j] =
2369 bnd[offset + nvel * nz_loc * nbnd + n * nint + j];
2370 }
2371 cnt += nint;
2372 }
2373 offset += (nvel * nbnd + nint) * nz_loc;
2374 }
2375
2376 // Back solve first level of static condensation for interior
2377 // velocity space and store in F_int
2378 F_int = F_int + Transpose(*m_mat[mode].m_D_int) * F_p -
2379 Transpose(*m_mat[mode].m_Btilde) * F_bnd;
2380 F_int = (*m_mat[mode].m_Cinv) * F_int;
2381
2382 // Unpack solution from Bnd and F_int to v_coeffs
2383 cnt = cnt1 = 0;
2384 for (i = 0; i < nel; ++i) // loop over elements
2385 {
2386 fields[0]->GetExp(i)->GetBoundaryMap(bmap);
2387 fields[0]->GetExp(i)->GetInteriorMap(imap);
2388 nbnd = bmap.size();
2389 nint = imap.size();
2390 offset = fields[0]->GetCoeff_Offset(i);
2391
2392 for (j = 0; j < nvel; ++j) // loop over velocity fields
2393 {
2394 for (n = 0; n < nz_loc; ++n)
2395 {
2396 for (k = 0; k < nbnd; ++k)
2397 {
2398 fields[j]->SetCoeff(n * nplanecoeffs + offset + bmap[k],
2399 f_bnd[cnt + k]);
2400 }
2401
2402 for (k = 0; k < nint; ++k)
2403 {
2404 fields[j]->SetCoeff(n * nplanecoeffs + offset + imap[k],
2405 f_int[cnt1 + k]);
2406 }
2407 cnt += nbnd;
2408 cnt1 += nint;
2409 }
2410 }
2411 }
2412
2413 for (j = 0; j < nvel; ++j)
2414 {
2415 fields[j]->SetPhysState(false);
2416 }
2417}
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:47
CG_Iterations.pressure
pressure
Definition: CG_Iterations.py:229
Nektar::MultiRegions::ContFieldSharedPtr
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:268

References Nektar::SolverUtils::EquationSystem::eHomogeneous1D, Nektar::SpatialDomains::eNeumann, Nektar::SpatialDomains::eRobin, Nektar::eWrapper, Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_HomogeneousType, m_locToGloMap, m_mat, Nektar::SolverUtils::EquationSystem::m_singleMode, Nektar::IncNavierStokes::m_velocity, CG_Iterations::pressure, sign, Nektar::Transpose(), and Vmath::Zero().

◆ SolveLinearNS() [2/2]

void Nektar::CoupledLinearNS::SolveLinearNS ( const Array< OneD, Array< OneD, NekDouble > > &  forcing)

Solve the coupled linear Navier-Stokes solve using matrix systems set up at construction. The solution is stored in m_velocity and m_pressure.

Parameters
forcingA list of forcing functions for each velocity component

The routine involves two levels of static condensations. Initially we require a statically condensed forcing function which requires the following manipulation

\[ {F\_bnd} = {\bf f}_{bnd} -m\_B \,m\_Cinv\, {\bf f}_{int}, \hspace{1cm} F\_p = m\_D\_{int}\, m\_Cinv\, {\bf f}_{int} \]

Where \({\bf f}_{bnd}\) denote the forcing degrees of freedom of the elemental velocities on the boundary of the element, \({\bf f}_{int}\) denote the forcing degrees of freedom of the elemental velocities on the interior of the element. (see detailed description for more details).

This vector is further manipulated into

\[ Fh\_{bnd} = \left [ \begin{array}{c} f\_{bnd} -m\_B \, m\_Cinv\, {\bf f}_{int}\\ \left [m\_D\_{int} \, m\_Cinv \,{\bf f}_{int} \right]_0 \end{array}\right ] \hspace{1cm} [Fh\_p]_{i} = \begin{array}{c} [m\_D\_{int} \, m\_Cinv \, {\bf f}_{int}]_{i+1} \end{array} \]

where \(-{[}m\_D\_{int}^T\, m\_Cinv \,{\bf f}_{int}]_0\) which is corresponds to the mean mode of the pressure degrees of freedom and is now added to the boundary system and the remainder of the block becomes the interior forcing for the inner static condensation (see detailed description for more details) which is set up in a GlobalLinSysDirectStaticCond class.

Finally we perform the final maniplation of the forcing to using hte

\[ Fh\_{bnd} = Fh\_{bnd} - m\_Bh \,m\_Chinv \, Fh\_p \]

We can now call the solver to the global coupled boundary system (through the call to #m_CoupledBndSys->Solve) to obtain the velocity boundary solution as the mean pressure solution, i.e.

\[ {\cal A}^T(\hat{A} - \hat{C}^T \hat{D}^{-1} \hat{B} ){\cal A} \, Bnd = Fh\_{bnd} \]

Once we know the solution to the above the rest of the pressure modes are recoverable thorugh

\[ Ph = m\_Dhinv\, (Bnd - m\_Ch^T \, Fh_{bnd}) \]

We can now unpack \( Fh\_{bnd} \) (last elemental mode) and \( Ph \) into m_pressure and \( F_p\) and \( Fh\_{bnd}\) into a closed pack list of boundary velocoity degrees of freedom stored in \( F\_bnd\).

Finally the interior velocity degrees of freedom are then obtained through the relationship

\[ F\_{int} = m\_Cinv\ ( F\_{int} + m\_D\_{int}^T\, F\_p - m\_Btilde^T\, Bnd) \]

We then unpack the solution back to the MultiRegion structures m_velocity and m_pressure

Definition at line 2085 of file CoupledLinearNS.cpp.

2087{
2088 size_t i;
2089 Array<OneD, MultiRegions::ExpListSharedPtr> vel_fields(m_velocity.size());
2090 Array<OneD, Array<OneD, NekDouble>> force(m_velocity.size());
2091
2092 // Impose Dirichlet conditions on velocity fields
2093 for (i = 0; i < m_velocity.size(); ++i)
2094 {
2095 Vmath::Zero(m_fields[i]->GetNcoeffs(), m_fields[i]->UpdateCoeffs(), 1);
2096 m_fields[i]->ImposeDirichletConditions(m_fields[i]->UpdateCoeffs());
2097 }
2098
2099 if (m_HomogeneousType == eHomogeneous1D)
2100 {
2101 int ncoeffsplane = m_fields[m_velocity[0]]->GetPlane(0)->GetNcoeffs();
2102 for (int n = 0; n < m_npointsZ / 2; ++n)
2103 {
2104 // Get the a Fourier mode of velocity and forcing components.
2105 for (i = 0; i < m_velocity.size(); ++i)
2106 {
2107 vel_fields[i] = m_fields[m_velocity[i]]->GetPlane(2 * n);
2108 // Note this needs to correlate with how we pass forcing
2109 force[i] = forcing[i] + 2 * n * ncoeffsplane;
2110 }
2111
2112 SolveLinearNS(force, vel_fields, m_pressure->GetPlane(2 * n), n);
2113 }
2114 for (i = 0; i < m_velocity.size(); ++i)
2115 {
2116 m_fields[m_velocity[i]]->SetPhysState(false);
2117 }
2118 m_pressure->SetPhysState(false);
2119 }
2120 else
2121 {
2122 for (i = 0; i < m_velocity.size(); ++i)
2123 {
2124 vel_fields[i] = m_fields[m_velocity[i]];
2125 // Note this needs to correlate with how we pass forcing
2126 force[i] = forcing[i];
2127 }
2128 SolveLinearNS(force, vel_fields, m_pressure);
2129 }
2130}

References Nektar::SolverUtils::EquationSystem::eHomogeneous1D, Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_HomogeneousType, Nektar::SolverUtils::EquationSystem::m_npointsZ, Nektar::IncNavierStokes::m_pressure, Nektar::IncNavierStokes::m_velocity, SolveLinearNS(), and Vmath::Zero().

Referenced by Continuation(), Solve(), SolveLinearNS(), SolveSteadyNavierStokes(), and SolveUnsteadyStokesSystem().

◆ SolveSteadyNavierStokes()

void Nektar::CoupledLinearNS::SolveSteadyNavierStokes ( void  )

Definition at line 1737 of file CoupledLinearNS.cpp.

1738{
1739 LibUtilities::Timer Newtontimer;
1740 Newtontimer.Start();
1741
1742 Array<OneD, Array<OneD, NekDouble>> RHS_Coeffs(m_velocity.size());
1743 Array<OneD, Array<OneD, NekDouble>> RHS_Phys(m_velocity.size());
1744 Array<OneD, Array<OneD, NekDouble>> delta_velocity_Phys(m_velocity.size());
1745 Array<OneD, Array<OneD, NekDouble>> Velocity_Phys(m_velocity.size());
1746 Array<OneD, NekDouble> L2_norm(m_velocity.size(), 1.0);
1747 Array<OneD, NekDouble> Inf_norm(m_velocity.size(), 1.0);
1748
1749 for (size_t i = 0; i < m_velocity.size(); ++i)
1750 {
1751 delta_velocity_Phys[i] = Array<OneD, NekDouble>(
1752 m_fields[m_velocity[i]]->GetTotPoints(), 1.0);
1753 Velocity_Phys[i] = Array<OneD, NekDouble>(
1754 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1755 }
1756
1757 m_counter = 1;
1758
1759 L2Norm(delta_velocity_Phys, L2_norm);
1760
1761 // while(max(Inf_norm[0], Inf_norm[1]) > m_tol)
1762 while (max(L2_norm[0], L2_norm[1]) > m_tol)
1763 {
1764 if (m_counter == 1)
1765 // At the first Newton step, we use the solution of the
1766 // Stokes problem (if Restart=0 in input file) Or the
1767 // solution of the .rst file (if Restart=1 in input
1768 // file)
1769 {
1770 for (size_t i = 0; i < m_velocity.size(); ++i)
1771 {
1772 RHS_Coeffs[i] = Array<OneD, NekDouble>(
1773 m_fields[m_velocity[i]]->GetNcoeffs(), 0.0);
1774 RHS_Phys[i] = Array<OneD, NekDouble>(
1775 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1776 }
1777
1778 for (size_t i = 0; i < m_velocity.size(); ++i)
1779 {
1780 m_fields[m_velocity[i]]->BwdTrans(
1781 m_fields[m_velocity[i]]->GetCoeffs(), Velocity_Phys[i]);
1782 }
1783
1784 m_initialStep = true;
1785 EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1786 SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1787 SolveLinearNS(RHS_Coeffs);
1788 m_initialStep = false;
1789 }
1790 if (m_counter > 1)
1791 {
1792 EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1793
1794 if (m_counter % m_MatrixSetUpStep ==
1795 0) // Setting Up the matrix is expensive. We do it at each
1796 // "m_MatrixSetUpStep" step.
1797 {
1798 SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1799 }
1800 SolveLinearNS(RHS_Coeffs);
1801 }
1802
1803 for (size_t i = 0; i < m_velocity.size(); ++i)
1804 {
1805 m_fields[m_velocity[i]]->BwdTrans(RHS_Coeffs[i], RHS_Phys[i]);
1806 m_fields[m_velocity[i]]->BwdTrans(
1807 m_fields[m_velocity[i]]->GetCoeffs(), delta_velocity_Phys[i]);
1808 }
1809
1810 for (size_t i = 0; i < m_velocity.size(); ++i)
1811 {
1812 Vmath::Vadd(Velocity_Phys[i].size(), Velocity_Phys[i], 1,
1813 delta_velocity_Phys[i], 1, Velocity_Phys[i], 1);
1814 }
1815
1816 // InfNorm(delta_velocity_Phys, Inf_norm);
1817 L2Norm(delta_velocity_Phys, L2_norm);
1818
1819 if (max(Inf_norm[0], Inf_norm[1]) > 100)
1820 {
1821 cout << "\nThe Newton method has failed at m_kinvis = " << m_kinvis
1822 << " (<=> Re = " << 1 / m_kinvis << ")" << endl;
1823 ASSERTL0(0, "The Newton method has failed... \n");
1824 }
1825
1826 cout << "\n";
1827 m_counter++;
1828 }
1829
1830 if (m_counter > 1) // We save u:=u+\delta u in u->Coeffs
1831 {
1832 for (size_t i = 0; i < m_velocity.size(); ++i)
1833 {
1834 m_fields[m_velocity[i]]->FwdTrans(
1835 Velocity_Phys[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1836 }
1837 }
1838
1839 Newtontimer.Stop();
1840 cout << "We have done " << m_counter - 1 << " iteration(s) in "
1841 << Newtontimer.TimePerTest(1) / 60 << " minute(s). \n\n";
1842}
Nektar::CoupledLinearNS::L2Norm
void L2Norm(Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
Definition: CoupledLinearNS.cpp:1912
Nektar::CoupledLinearNS::EvaluateNewtonRHS
void EvaluateNewtonRHS(Array< OneD, Array< OneD, NekDouble > > &Velocity, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: CoupledLinearNS.cpp:1927

References ASSERTL0, EvaluateNewtonRHS(), Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), L2Norm(), m_counter, Nektar::SolverUtils::EquationSystem::m_fields, m_initialStep, Nektar::IncNavierStokes::m_kinvis, m_MatrixSetUpStep, m_tol, Nektar::IncNavierStokes::m_velocity, SetUpCoupledMatrix(), SolveLinearNS(), Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), Nektar::LibUtilities::Timer::TimePerTest(), and Vmath::Vadd().

Referenced by v_DoSolve().

◆ SolveUnsteadyStokesSystem()

void Nektar::CoupledLinearNS::SolveUnsteadyStokesSystem ( const Array< OneD, const Array< OneD, NekDouble > > &  inarray,
Array< OneD, Array< OneD, NekDouble > > &  outarray,
const NekDouble  time,
const NekDouble  a_iixDt 
)

Definition at line 1528 of file CoupledLinearNS.cpp.

1532{
1533 size_t i;
1534 Array<OneD, Array<OneD, NekDouble>> F(m_nConvectiveFields);
1535 NekDouble lambda = 1.0 / aii_Dt;
1536 static NekDouble lambda_store;
1537 Array<OneD, Array<OneD, NekDouble>> forcing(m_velocity.size());
1538 // Matrix solution
1539 if (fabs(lambda_store - lambda) > 1e-10)
1540 {
1541 SetUpCoupledMatrix(lambda);
1542 lambda_store = lambda;
1543 }
1544
1545 // Forcing for advection solve
1546 for (i = 0; i < m_velocity.size(); ++i)
1547 {
1548 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1549 m_fields[m_velocity[i]]->SetWaveSpace(true);
1550 m_fields[m_velocity[i]]->IProductWRTBase(
1551 inarray[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1552 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1553 Vmath::Smul(m_fields[m_velocity[i]]->GetNcoeffs(), lambda,
1554 m_fields[m_velocity[i]]->GetCoeffs(), 1,
1555 m_fields[m_velocity[i]]->UpdateCoeffs(), 1);
1556 forcing[i] = m_fields[m_velocity[i]]->GetCoeffs();
1557 }
1558
1559 SolveLinearNS(forcing);
1560
1561 for (i = 0; i < m_velocity.size(); ++i)
1562 {
1563 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1564 outarray[i]);
1565 }
1566}
Nektar::IncNavierStokes::m_nConvectiveFields
int m_nConvectiveFields
Number of fields to be convected;.
Definition: IncNavierStokes.h:210

References Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::IncNavierStokes::m_nConvectiveFields, Nektar::IncNavierStokes::m_velocity, SetUpCoupledMatrix(), Vmath::Smul(), and SolveLinearNS().

Referenced by v_DoInitialise().

◆ v_DoInitialise()

void Nektar::CoupledLinearNS::v_DoInitialise ( bool  dumpInitialConditions = true)
overrideprivatevirtual

Sets up initial conditions.

Sets the initial conditions.

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 1363 of file CoupledLinearNS.cpp.

1364{
1365 switch (m_equationType)
1366 {
1367 case eUnsteadyStokes:
1368 case eUnsteadyNavierStokes:
1369 {
1370 // LibUtilities::TimeIntegrationMethod intMethod;
1371 // std::string TimeIntStr =
1372 // m_session->GetSolverInfo("TIMEINTEGRATIONMETHOD");
1373 // int i;
1374 // for(i = 0; i < (int)
1375 // LibUtilities::SIZE_TimeIntegrationMethod; ++i)
1376 // {
1377 // if(boost::iequals(LibUtilities::TimeIntegrationMethodMap[i],TimeIntStr))
1378 // {
1379 // intMethod =
1380 // (LibUtilities::TimeIntegrationMethod)i;
1381 // break;
1382 // }
1383 // }
1384 //
1385 // ASSERTL0(i != (int)
1386 // LibUtilities::SIZE_TimeIntegrationMethod, "Invalid
1387 // time integration type.");
1388 //
1389 // m_integrationScheme =
1390 // LibUtilities::GetTimeIntegrationWrapperFactory().CreateInstance(LibUtilities::TimeIntegrationMethodMap[intMethod]);
1391
1392 // Could defind this from IncNavierStokes class?
1393 m_ode.DefineOdeRhs(&CoupledLinearNS::EvaluateAdvection, this);
1394
1395 m_ode.DefineImplicitSolve(
1396 &CoupledLinearNS::SolveUnsteadyStokesSystem, this);
1397
1398 // Set initial condition using time t=0
1399
1400 SetInitialConditions(0.0, dumpInitialConditions);
1401 break;
1402 }
1403 case eSteadyStokes:
1404 SetUpCoupledMatrix(0.0);
1405 break;
1406 case eSteadyOseen:
1407 {
1408 Array<OneD, Array<OneD, NekDouble>> AdvField(m_velocity.size());
1409 for (size_t i = 0; i < m_velocity.size(); ++i)
1410 {
1411 AdvField[i] = Array<OneD, NekDouble>(
1412 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1413 }
1414
1415 ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1416 "Advection Velocity section must be defined in "
1417 "session file.");
1418
1419 std::vector<std::string> fieldStr;
1420 for (size_t i = 0; i < m_velocity.size(); ++i)
1421 {
1422 fieldStr.push_back(
1423 m_boundaryConditions->GetVariable(m_velocity[i]));
1424 }
1425 GetFunction("AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1426
1427 SetUpCoupledMatrix(0.0, AdvField, false);
1428 }
1429 break;
1430 case eSteadyNavierStokes:
1431 {
1432 m_session->LoadParameter("KinvisMin", m_kinvisMin);
1433 m_session->LoadParameter("KinvisPercentage", m_KinvisPercentage);
1434 m_session->LoadParameter("Tolerence", m_tol);
1435 m_session->LoadParameter("MaxIteration", m_maxIt);
1436 m_session->LoadParameter("MatrixSetUpStep", m_MatrixSetUpStep);
1437 m_session->LoadParameter("Restart", m_Restart);
1438
1439 DefineForcingTerm();
1440
1441 if (m_Restart == 1)
1442 {
1443 ASSERTL0(m_session->DefinesFunction("Restart"),
1444 "Restart section must be defined in session file.");
1445
1446 Array<OneD, Array<OneD, NekDouble>> Restart(m_velocity.size());
1447 for (size_t i = 0; i < m_velocity.size(); ++i)
1448 {
1449 Restart[i] = Array<OneD, NekDouble>(
1450 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1451 }
1452 std::vector<std::string> fieldStr;
1453 for (size_t i = 0; i < m_velocity.size(); ++i)
1454 {
1455 fieldStr.push_back(
1456 m_boundaryConditions->GetVariable(m_velocity[i]));
1457 }
1458 GetFunction("Restart")->Evaluate(fieldStr, Restart);
1459
1460 for (size_t i = 0; i < m_velocity.size(); ++i)
1461 {
1462 m_fields[m_velocity[i]]->FwdTransLocalElmt(
1463 Restart[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1464 }
1465 cout << "Saving the RESTART file for m_kinvis = " << m_kinvis
1466 << " (<=> Re = " << 1 / m_kinvis << ")" << endl;
1467 }
1468 else // We solve the Stokes Problem
1469 {
1470
1471 SetUpCoupledMatrix(0.0);
1472 m_initialStep = true;
1473 m_counter = 1;
1474 // SolveLinearNS(m_ForcingTerm_Coeffs);
1475 Solve();
1476 m_initialStep = false;
1477 cout << "Saving the Stokes Flow for m_kinvis = " << m_kinvis
1478 << " (<=> Re = " << 1 / m_kinvis << ")" << endl;
1479 }
1480 }
1481 break;
1482 case eSteadyLinearisedNS:
1483 {
1484 SetInitialConditions(0.0, dumpInitialConditions);
1485
1486 Array<OneD, Array<OneD, NekDouble>> AdvField(m_velocity.size());
1487 for (size_t i = 0; i < m_velocity.size(); ++i)
1488 {
1489 AdvField[i] = Array<OneD, NekDouble>(
1490 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1491 }
1492
1493 ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1494 "Advection Velocity section must be defined in "
1495 "session file.");
1496
1497 std::vector<std::string> fieldStr;
1498 for (size_t i = 0; i < m_velocity.size(); ++i)
1499 {
1500 fieldStr.push_back(
1501 m_boundaryConditions->GetVariable(m_velocity[i]));
1502 }
1503 GetFunction("AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1504
1505 SetUpCoupledMatrix(m_lambda, AdvField, true);
1506 }
1507 break;
1508 case eNoEquationType:
1509 default:
1510 ASSERTL0(false,
1511 "Unknown or undefined equation type for CoupledLinearNS");
1512 }
1513}
Nektar::CoupledLinearNS::SolveUnsteadyStokesSystem
void SolveUnsteadyStokesSystem(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time, const NekDouble a_iixDt)
Definition: CoupledLinearNS.cpp:1528
Nektar::CoupledLinearNS::EvaluateAdvection
void EvaluateAdvection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Definition: CoupledLinearNS.cpp:1515
Nektar::IncNavierStokes::m_equationType
EquationType m_equationType
equation type;
Definition: IncNavierStokes.h:224
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineOdeRhs
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationSchemeOperators.h:84
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineImplicitSolve
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationSchemeOperators.h:100
Nektar::SolverUtils::EquationSystem::m_lambda
NekDouble m_lambda
Lambda constant in real system if one required.
Definition: EquationSystem.h:376
Nektar::SolverUtils::EquationSystem::SetInitialConditions
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
Definition: EquationSystem.h:674
Nektar::SolverUtils::UnsteadySystem::m_ode
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:88
Nektar::eSteadyNavierStokes
@ eSteadyNavierStokes
Definition: IncNavierStokes.h:76
Nektar::eUnsteadyStokes
@ eUnsteadyStokes
Definition: IncNavierStokes.h:73
Nektar::eUnsteadyNavierStokes
@ eUnsteadyNavierStokes
Definition: IncNavierStokes.h:75
Nektar::eSteadyLinearisedNS
@ eSteadyLinearisedNS
Definition: IncNavierStokes.h:72
Nektar::eNoEquationType
@ eNoEquationType
Definition: IncNavierStokes.h:69

References ASSERTL0, DefineForcingTerm(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineImplicitSolve(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineOdeRhs(), Nektar::eNoEquationType, Nektar::eSteadyLinearisedNS, Nektar::eSteadyNavierStokes, Nektar::eSteadyOseen, Nektar::eSteadyStokes, Nektar::eUnsteadyNavierStokes, Nektar::eUnsteadyStokes, EvaluateAdvection(), Nektar::SolverUtils::EquationSystem::GetFunction(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::m_boundaryConditions, m_counter, Nektar::IncNavierStokes::m_equationType, Nektar::SolverUtils::EquationSystem::m_fields, m_initialStep, Nektar::IncNavierStokes::m_kinvis, m_kinvisMin, m_KinvisPercentage, Nektar::SolverUtils::EquationSystem::m_lambda, m_MatrixSetUpStep, m_maxIt, Nektar::SolverUtils::UnsteadySystem::m_ode, m_Restart, Nektar::SolverUtils::EquationSystem::m_session, m_tol, Nektar::IncNavierStokes::m_velocity, Nektar::SolverUtils::EquationSystem::SetInitialConditions(), SetUpCoupledMatrix(), Solve(), and SolveUnsteadyStokesSystem().

◆ v_DoSolve()

void Nektar::CoupledLinearNS::v_DoSolve ( void  )
overrideprivatevirtual

Solves an unsteady problem.

Initialises the time integration scheme (as specified in the session file), and perform the time integration.

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 1590 of file CoupledLinearNS.cpp.

1591{
1592 switch (m_equationType)
1593 {
1594 case eUnsteadyStokes:
1595 case eUnsteadyNavierStokes:
1596 // AdvanceInTime(m_steps);
1597 UnsteadySystem::v_DoSolve();
1598 break;
1599 case eSteadyStokes:
1600 case eSteadyOseen:
1601 case eSteadyLinearisedNS:
1602 {
1603 Solve();
1604 break;
1605 }
1606 case eSteadyNavierStokes:
1607 {
1608 LibUtilities::Timer Generaltimer;
1609 Generaltimer.Start();
1610
1611 int Check(0);
1612
1613 // Saving the init datas
1614 Checkpoint_Output(Check);
1615 Check++;
1616
1617 cout << "We execute INITIALLY SolveSteadyNavierStokes for m_kinvis "
1618 "= "
1619 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")" << endl;
1620 SolveSteadyNavierStokes();
1621
1622 while (m_kinvis > m_kinvisMin)
1623 {
1624 if (Check == 1)
1625 {
1626 cout << "We execute SolveSteadyNavierStokes for m_kinvis = "
1627 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1628 << endl;
1629 SolveSteadyNavierStokes();
1630 Checkpoint_Output(Check);
1631 Check++;
1632 }
1633
1634 Continuation();
1635
1636 if (m_kinvis > m_kinvisMin)
1637 {
1638 cout << "We execute SolveSteadyNavierStokes for m_kinvis = "
1639 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1640 << endl;
1641 SolveSteadyNavierStokes();
1642 Checkpoint_Output(Check);
1643 Check++;
1644 }
1645 }
1646
1647 Generaltimer.Stop();
1648 cout << "\nThe total calculation time is : "
1649 << Generaltimer.TimePerTest(1) / 60 << " minute(s). \n\n";
1650
1651 break;
1652 }
1653 case eNoEquationType:
1654 default:
1655 ASSERTL0(false,
1656 "Unknown or undefined equation type for CoupledLinearNS");
1657 }
1658}
Nektar::SolverUtils::EquationSystem::Checkpoint_Output
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
Definition: EquationSystem.cpp:1176
Nektar::SolverUtils::UnsteadySystem::v_DoSolve
SOLVER_UTILS_EXPORT void v_DoSolve() override
Solves an unsteady problem.
Definition: UnsteadySystem.cpp:210

References ASSERTL0, Nektar::SolverUtils::EquationSystem::Checkpoint_Output(), Continuation(), Nektar::eNoEquationType, Nektar::eSteadyLinearisedNS, Nektar::eSteadyNavierStokes, Nektar::eSteadyOseen, Nektar::eSteadyStokes, Nektar::eUnsteadyNavierStokes, Nektar::eUnsteadyStokes, Nektar::IncNavierStokes::m_equationType, Nektar::IncNavierStokes::m_kinvis, m_kinvisMin, Solve(), SolveSteadyNavierStokes(), Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), Nektar::LibUtilities::Timer::TimePerTest(), and Nektar::SolverUtils::UnsteadySystem::v_DoSolve().

◆ v_GenerateSummary()

void Nektar::CoupledLinearNS::v_GenerateSummary ( SolverUtils::SummaryList &  s)
overrideprivatevirtual

Print a summary of time stepping parameters.

Prints a summary with some information regards the time-stepping.

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 1358 of file CoupledLinearNS.cpp.

1359{
1360 SolverUtils::AddSummaryItem(s, "Solver Type", "Coupled Linearised NS");
1361}
Nektar::SolverUtils::AddSummaryItem
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:47

References Nektar::SolverUtils::AddSummaryItem().

◆ v_GetForceDimension()

int Nektar::CoupledLinearNS::v_GetForceDimension ( void  )
overrideprivatevirtual

Implements Nektar::IncNavierStokes.

Definition at line 2479 of file CoupledLinearNS.cpp.

2480{
2481 return m_session->GetVariables().size();
2482}

References Nektar::SolverUtils::EquationSystem::m_session.

◆ v_InitObject()

void Nektar::CoupledLinearNS::v_InitObject ( bool  DeclareFeld = true)
overrideprotectedvirtual

Initialisation object for EquationSystem.

Continuous field

Setting up the normals

Setting up the normals

Reimplemented from Nektar::SolverUtils::AdvectionSystem.

Definition at line 71 of file CoupledLinearNS.cpp.

72{
73 IncNavierStokes::v_InitObject(DeclareField);
74
75 size_t i;
76 int expdim = m_graph->GetMeshDimension();
77
78 // Get Expansion list for orthogonal expansion at p-2
79 const SpatialDomains::ExpansionInfoMap &pressure_exp =
80 GenPressureExp(m_graph->GetExpansionInfo("u"));
81
82 m_nConvectiveFields = m_fields.size();
83 if (boost::iequals(
84 m_boundaryConditions->GetVariable(m_nConvectiveFields - 1), "p"))
85 {
86 ASSERTL0(false, "Last field is defined as pressure but this is not "
87 "suitable for this solver, please remove this field as "
88 "it is implicitly defined");
89 }
90 // Decide how to declare explist for pressure.
91 if (expdim == 2)
92 {
93 size_t nz;
94
95 if (m_HomogeneousType == eHomogeneous1D)
96 {
97 ASSERTL0(m_fields.size() > 2, "Expect to have three at least three "
98 "components of velocity variables");
99 LibUtilities::BasisKey Homo1DKey =
100 m_fields[0]->GetHomogeneousBasis()->GetBasisKey();
101
102 m_pressure = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::
103 AllocateSharedPtr(m_session, Homo1DKey, m_LhomZ, m_useFFT,
104 m_homogen_dealiasing, pressure_exp);
105
106 ASSERTL1(m_npointsZ % 2 == 0,
107 "Non binary number of planes have been specified");
108 nz = m_npointsZ / 2;
109 }
110 else
111 {
112 m_pressure =
113 MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
114 m_session, pressure_exp);
115 nz = 1;
116 }
117
118 Array<OneD, MultiRegions::ExpListSharedPtr> velocity(m_velocity.size());
119 for (i = 0; i < m_velocity.size(); ++i)
120 {
121 velocity[i] = m_fields[m_velocity[i]];
122 }
123
124 // Set up Array of mappings
125 m_locToGloMap = Array<OneD, CoupledLocalToGlobalC0ContMapSharedPtr>(nz);
126
127 if (m_singleMode)
128 {
129
130 ASSERTL0(nz <= 2,
131 "For single mode calculation can only have nz <= 2");
132 if (m_session->DefinesSolverInfo("BetaZero"))
133 {
134 m_zeroMode = true;
135 }
136 int nz_loc = 2;
137 m_locToGloMap[0] =
138 MemoryManager<CoupledLocalToGlobalC0ContMap>::AllocateSharedPtr(
139 m_session, m_graph, m_boundaryConditions, velocity,
140 m_pressure, nz_loc, m_zeroMode);
141 }
142 else
143 {
144 // base mode
145 int nz_loc = 1;
146 m_locToGloMap[0] =
147 MemoryManager<CoupledLocalToGlobalC0ContMap>::AllocateSharedPtr(
148 m_session, m_graph, m_boundaryConditions, velocity,
149 m_pressure, nz_loc);
150
151 if (nz > 1)
152 {
153 nz_loc = 2;
154 // Assume all higher modes have the same boundary conditions and
155 // re-use mapping
156 m_locToGloMap[1] =
157 MemoryManager<CoupledLocalToGlobalC0ContMap>::
158 AllocateSharedPtr(
159 m_session, m_graph, m_boundaryConditions, velocity,
160 m_pressure->GetPlane(2), nz_loc, false);
161 // Note high order modes cannot be singular
162 for (i = 2; i < nz; ++i)
163 {
164 m_locToGloMap[i] = m_locToGloMap[1];
165 }
166 }
167 }
168 }
169 else if (expdim == 3)
170 {
171 // m_pressure =
172 // MemoryManager<MultiRegions::ExpList3D>::AllocateSharedPtr(pressure_exp);
173 ASSERTL0(false, "Setup mapping aray");
174 }
175 else
176 {
177 ASSERTL0(false, "Exp dimension not recognised");
178 }
179
180 // creation of the extrapolation object
181 if (m_equationType == eUnsteadyNavierStokes)
182 {
183 std::string vExtrapolation = "Standard";
184
185 if (m_session->DefinesSolverInfo("Extrapolation"))
186 {
187 vExtrapolation = m_session->GetSolverInfo("Extrapolation");
188 }
189
190 m_extrapolation = GetExtrapolateFactory().CreateInstance(
191 vExtrapolation, m_session, m_fields, m_pressure, m_velocity,
192 m_advObject);
193 }
194}
Nektar::CoupledLinearNS::GenPressureExp
const SpatialDomains::ExpansionInfoMap & GenPressureExp(const SpatialDomains::ExpansionInfoMap &VelExp)
Definition: CoupledLinearNS.cpp:1983
Nektar::IncNavierStokes::v_InitObject
void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
Definition: IncNavierStokes.cpp:91
Nektar::IncNavierStokes::m_extrapolation
ExtrapolateSharedPtr m_extrapolation
Definition: IncNavierStokes.h:198
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:143
Nektar::SolverUtils::AdvectionSystem::m_advObject
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
Definition: AdvectionSystem.h:69
Nektar::SolverUtils::EquationSystem::m_useFFT
bool m_useFFT
Flag to determine if FFT is used for homogeneous transform.
Definition: EquationSystem.h:406
Nektar::SolverUtils::EquationSystem::m_homogen_dealiasing
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
Definition: EquationSystem.h:411
Nektar::MovementTests::velocity
const std::vector< NekDouble > velocity
Definition: TestMovement.cpp:57
Nektar::SpatialDomains::ExpansionInfoMap
std::map< int, ExpansionInfoShPtr > ExpansionInfoMap
Definition: MeshGraph.h:141
Nektar::GetExtrapolateFactory
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:48

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, ASSERTL1, Nektar::LibUtilities::NekFactory< tKey, tBase, tParam >::CreateInstance(), Nektar::SolverUtils::EquationSystem::eHomogeneous1D, Nektar::eUnsteadyNavierStokes, GenPressureExp(), Nektar::GetExtrapolateFactory(), Nektar::SolverUtils::AdvectionSystem::m_advObject, Nektar::SolverUtils::EquationSystem::m_boundaryConditions, Nektar::IncNavierStokes::m_equationType, Nektar::IncNavierStokes::m_extrapolation, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_graph, Nektar::SolverUtils::EquationSystem::m_homogen_dealiasing, Nektar::SolverUtils::EquationSystem::m_HomogeneousType, Nektar::SolverUtils::EquationSystem::m_LhomZ, m_locToGloMap, Nektar::IncNavierStokes::m_nConvectiveFields, Nektar::SolverUtils::EquationSystem::m_npointsZ, Nektar::IncNavierStokes::m_pressure, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_singleMode, Nektar::SolverUtils::EquationSystem::m_useFFT, Nektar::IncNavierStokes::m_velocity, m_zeroMode, Nektar::IncNavierStokes::v_InitObject(), and Nektar::MovementTests::velocity.

◆ v_NegatedOp()

bool Nektar::CoupledLinearNS::v_NegatedOp ( void  )
overrideprivatevirtual

Virtual function to define if operator in DoSolve is negated with regard to the strong form. This is currently only used in Arnoldi solves. For Coupledd solver this is true since Stokes operator is set up as a LHS rather than RHS operation

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1665 of file CoupledLinearNS.cpp.

1666{
1667 return true;
1668}

◆ v_Output()

void Nektar::CoupledLinearNS::v_Output ( void  )
overrideprivatevirtual

Write the field data to file. The file is named according to the session name with the extension .fld appended.

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 2419 of file CoupledLinearNS.cpp.

2420{
2421 std::vector<Array<OneD, NekDouble>> fieldcoeffs(m_fields.size() + 1);
2422 std::vector<std::string> variables(m_fields.size() + 1);
2423 size_t i;
2424
2425 for (i = 0; i < m_fields.size(); ++i)
2426 {
2427 fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
2428 variables[i] = m_boundaryConditions->GetVariable(i);
2429 }
2430
2431 fieldcoeffs[i] = Array<OneD, NekDouble>(m_fields[0]->GetNcoeffs());
2432 // project pressure field to velocity space
2433 if (m_singleMode == true)
2434 {
2435 Array<OneD, NekDouble> tmpfieldcoeffs(m_fields[0]->GetNcoeffs() / 2);
2436 m_pressure->GetPlane(0)->BwdTrans(
2437 m_pressure->GetPlane(0)->GetCoeffs(),
2438 m_pressure->GetPlane(0)->UpdatePhys());
2439 m_pressure->GetPlane(1)->BwdTrans(
2440 m_pressure->GetPlane(1)->GetCoeffs(),
2441 m_pressure->GetPlane(1)->UpdatePhys());
2442 m_fields[0]->GetPlane(0)->FwdTransLocalElmt(
2443 m_pressure->GetPlane(0)->GetPhys(), fieldcoeffs[i]);
2444 m_fields[0]->GetPlane(1)->FwdTransLocalElmt(
2445 m_pressure->GetPlane(1)->GetPhys(), tmpfieldcoeffs);
2446 for (int e = 0; e < m_fields[0]->GetNcoeffs() / 2; e++)
2447 {
2448 fieldcoeffs[i][e + m_fields[0]->GetNcoeffs() / 2] =
2449 tmpfieldcoeffs[e];
2450 }
2451 }
2452 else
2453 {
2454 m_pressure->BwdTrans(m_pressure->GetCoeffs(), m_pressure->UpdatePhys());
2455 m_fields[0]->FwdTransLocalElmt(m_pressure->GetPhys(), fieldcoeffs[i]);
2456 }
2457 variables[i] = "p";
2458
2459 if (!m_comm->IsParallelInTime())
2460 {
2461 WriteFld(m_sessionName + ".fld", m_fields[0], fieldcoeffs, variables);
2462 }
2463 else
2464 {
2465 std::string newdir = m_sessionName + ".pit";
2466 if (!fs::is_directory(newdir))
2467 {
2468 fs::create_directory(newdir);
2469 }
2470 WriteFld(
2471 newdir + "/" + m_sessionName + "_" +
2472 std::to_string(m_windowPIT * m_comm->GetTimeComm()->GetSize() +
2473 m_comm->GetTimeComm()->GetRank() + 1) +
2474 ".fld",
2475 m_fields[0], fieldcoeffs, variables);
2476 }
2477}
Nektar::SolverUtils::EquationSystem::m_comm
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: EquationSystem.h:349
Nektar::SolverUtils::EquationSystem::m_windowPIT
int m_windowPIT
Index of windows for parallel-in-time time iteration.
Definition: EquationSystem.h:394
Nektar::SolverUtils::EquationSystem::WriteFld
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition: EquationSystem.cpp:1238
Nektar::SolverUtils::EquationSystem::m_sessionName
std::string m_sessionName
Name of the session.
Definition: EquationSystem.h:365

References Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::m_boundaryConditions, Nektar::SolverUtils::EquationSystem::m_comm, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::IncNavierStokes::m_pressure, Nektar::SolverUtils::EquationSystem::m_sessionName, Nektar::SolverUtils::EquationSystem::m_singleMode, Nektar::SolverUtils::EquationSystem::m_windowPIT, and Nektar::SolverUtils::EquationSystem::WriteFld().

◆ v_TransCoeffToPhys()

void Nektar::CoupledLinearNS::v_TransCoeffToPhys ( void  )
overrideprivatevirtual

Virtual function for transformation to physical space.

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1568 of file CoupledLinearNS.cpp.

1569{
1570 size_t nfields = m_fields.size();
1571 for (size_t k = 0; k < nfields; ++k)
1572 {
1573 // Backward Transformation in physical space for time evolution
1574 m_fields[k]->BwdTrans(m_fields[k]->GetCoeffs(),
1575 m_fields[k]->UpdatePhys());
1576 }
1577}

References Nektar::SolverUtils::EquationSystem::m_fields.

◆ v_TransPhysToCoeff()

void Nektar::CoupledLinearNS::v_TransPhysToCoeff ( void  )
overrideprivatevirtual

Virtual function for transformation to coefficient space.

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1579 of file CoupledLinearNS.cpp.

1580{
1581 size_t nfields = m_fields.size();
1582 for (size_t k = 0; k < nfields; ++k)
1583 {
1584 // Forward Transformation in physical space for time evolution
1585 m_fields[k]->FwdTransLocalElmt(m_fields[k]->GetPhys(),
1586 m_fields[k]->UpdateCoeffs());
1587 }
1588}

References Nektar::SolverUtils::EquationSystem::m_fields.

Friends And Related Function Documentation

◆ MemoryManager< CoupledLinearNS >

friend class MemoryManager< CoupledLinearNS >
friend

Definition at line 83 of file CoupledLinearNS.h.

Member Data Documentation

◆ className

string Nektar::CoupledLinearNS::className
static
Initial value:
=
SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
"CoupledLinearisedNS", CoupledLinearNS::create)
Nektar::CoupledLinearNS::create
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: CoupledLinearNS.h:92
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:197
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:99

Name of class.

Definition at line 102 of file CoupledLinearNS.h.

◆ m_counter

int Nektar::CoupledLinearNS::m_counter
private

Definition at line 176 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

◆ m_ForcingTerm

Array<OneD, Array<OneD, NekDouble> > Nektar::CoupledLinearNS::m_ForcingTerm

Definition at line 159 of file CoupledLinearNS.h.

Referenced by DefineForcingTerm(), and EvaluateNewtonRHS().

◆ m_ForcingTerm_Coeffs

Array<OneD, Array<OneD, NekDouble> > Nektar::CoupledLinearNS::m_ForcingTerm_Coeffs

Definition at line 160 of file CoupledLinearNS.h.

Referenced by DefineForcingTerm().

◆ m_initialStep

bool Nektar::CoupledLinearNS::m_initialStep
private

Definition at line 177 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

◆ m_kinvisMin

NekDouble Nektar::CoupledLinearNS::m_kinvisMin
private

Definition at line 183 of file CoupledLinearNS.h.

Referenced by v_DoInitialise(), and v_DoSolve().

◆ m_KinvisPercentage

NekDouble Nektar::CoupledLinearNS::m_KinvisPercentage
private

Definition at line 185 of file CoupledLinearNS.h.

Referenced by Continuation(), and v_DoInitialise().

◆ m_kinvisStep

NekDouble Nektar::CoupledLinearNS::m_kinvisStep
private

Definition at line 184 of file CoupledLinearNS.h.

◆ m_locToGloMap

Array<OneD, CoupledLocalToGlobalC0ContMapSharedPtr> Nektar::CoupledLinearNS::m_locToGloMap

Definition at line 162 of file CoupledLinearNS.h.

Referenced by SetUpCoupledMatrix(), SolveLinearNS(), and v_InitObject().

◆ m_mat

Array<OneD, CoupledSolverMatrices> Nektar::CoupledLinearNS::m_mat
private

Definition at line 187 of file CoupledLinearNS.h.

Referenced by SetUpCoupledMatrix(), and SolveLinearNS().

◆ m_MatrixSetUpStep

int Nektar::CoupledLinearNS::m_MatrixSetUpStep
private

Definition at line 182 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

◆ m_maxIt

int Nektar::CoupledLinearNS::m_maxIt
private

Definition at line 179 of file CoupledLinearNS.h.

Referenced by v_DoInitialise().

◆ m_Restart

int Nektar::CoupledLinearNS::m_Restart
private

Definition at line 180 of file CoupledLinearNS.h.

Referenced by v_DoInitialise().

◆ m_tol

NekDouble Nektar::CoupledLinearNS::m_tol
private

Definition at line 178 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

◆ m_zeroMode

bool Nektar::CoupledLinearNS::m_zeroMode
private

Id to identify when single mode is mean mode (i.e. beta=0);.

Definition at line 174 of file CoupledLinearNS.h.

Referenced by SetUpCoupledMatrix(), and v_InitObject().

◆ solverTypeLookupId

string Nektar::CoupledLinearNS::solverTypeLookupId
staticprotected
Initial value:
=
LibUtilities::SessionReader::RegisterEnumValue(
"SolverType", "CoupledLinearisedNS", eCoupledLinearisedNS)
Nektar::LibUtilities::SessionReader::RegisterEnumValue
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
Definition: SessionReader.h:559
Nektar::eCoupledLinearisedNS
@ eCoupledLinearisedNS
Definition: IncNavierStokes.h:57

Definition at line 170 of file CoupledLinearNS.h.

Generated on Wed Jan 31 2024 09:36:00 for Nektar++ by doxygen 1.9.4