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
virtual ~IncNavierStokes ()
 
virtual void v_InitObject (bool DeclareField=true) override
 Init object for UnsteadySystem class. More...
 
int GetNConvectiveFields (void)
 
void AddForcing (const SolverUtils::ForcingSharedPtr &pForce)
 
virtual void v_GetPressure (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure) override
 
virtual void v_GetDensity (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &density) override
 
virtual bool v_HasConstantDensity () override
 
virtual void v_GetVelocity (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity) override
 
virtual void v_SetMovingFrameVelocities (const Array< OneD, NekDouble > &vFrameVels) override
 
virtual void v_GetMovingFrameVelocities (Array< OneD, NekDouble > &vFrameVels) override
 
virtual void v_SetMovingFrameAngles (const Array< OneD, NekDouble > &vFrameTheta) override
 
virtual void v_GetMovingFrameAngles (Array< OneD, NekDouble > &vFrameTheta) override
 
virtual void v_SetMovingFrameProjectionMat (const bnu::matrix< NekDouble > &vProjMat) override
 
virtual void v_GetMovingFrameProjectionMat (bnu::matrix< NekDouble > &vProjMat) override
 
bool DefinedForcing (const std::string &sForce)
 
void GetPivotPoint (Array< OneD, NekDouble > &vPivotPoint)
 
- Public Member Functions inherited from Nektar::SolverUtils::AdvectionSystem
SOLVER_UTILS_EXPORT AdvectionSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 
virtual SOLVER_UTILS_EXPORT ~AdvectionSystem ()
 
virtual SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareField=true) override
 Init object for UnsteadySystem class. 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
virtual SOLVER_UTILS_EXPORT ~UnsteadySystem ()
 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 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 GetFinalTime ()
 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 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...
 
SOLVER_UTILS_EXPORT bool ParallelInTime ()
 Check if solver use Parallel-in-Time. 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 SetMovingFrameAngles (const Array< OneD, NekDouble > &vFrameTheta)
 
SOLVER_UTILS_EXPORT void GetMovingFrameAngles (Array< OneD, NekDouble > &vFrameTheta)
 

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)
 
virtual void v_InitObject (bool DeclareField=true) override
 Init object for UnsteadySystem class. 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)
 
virtual MultiRegions::ExpListSharedPtr v_GetPressure () override
 
virtual void v_TransCoeffToPhys (void) override
 Virtual function for transformation to physical space. More...
 
virtual void v_TransPhysToCoeff (void) override
 Virtual function for transformation to coefficient space. More...
 
virtual int v_GetForceDimension ()=0
 
virtual Array< OneD, NekDouble > v_GetMaxStdVelocity (const NekDouble SpeedSoundFactor) override
 
virtual bool v_PreIntegrate (int step) override
 
- Protected Member Functions inherited from Nektar::SolverUtils::AdvectionSystem
virtual 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...
 
virtual SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareField=true) override
 Init object for UnsteadySystem class. More...
 
virtual 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...
 
virtual SOLVER_UTILS_EXPORT void v_DoInitialise (bool dumpInitialConditions=true) override
 Sets up initial conditions. More...
 
virtual 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_SetMovingFrameAngles (const Array< OneD, NekDouble > &vFrameTheta)
 
virtual SOLVER_UTILS_EXPORT void v_GetMovingFrameAngles (Array< OneD, NekDouble > &vFrameTheta)
 

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)
 
virtual void v_GenerateSummary (SolverUtils::SummaryList &s) override
 Print a summary of time stepping parameters. More...
 
virtual void v_DoInitialise (bool dumpInitialConditions=true) override
 Sets up initial conditions. More...
 
virtual void v_DoSolve (void) override
 Solves an unsteady problem. More...
 
virtual bool v_NegatedOp (void) override
 
virtual void v_TransCoeffToPhys (void) override
 Virtual function for transformation to physical space. More...
 
virtual void v_TransPhysToCoeff (void) override
 Virtual function for transformation to coefficient space. More...
 
virtual void v_Output (void) override
 
virtual 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
 
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. More...
 
Array< OneD, NekDouble > m_movingFrameTheta
 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:76

Member Function Documentation

◆ Continuation()

void Nektar::CoupledLinearNS::Continuation ( void  )

Definition at line 1846 of file CoupledLinearNS.cpp.

1847{
1848 Array<OneD, Array<OneD, NekDouble>> u_N(m_velocity.size());
1849 Array<OneD, Array<OneD, NekDouble>> tmp_RHS(m_velocity.size());
1850 Array<OneD, Array<OneD, NekDouble>> RHS(m_velocity.size());
1851 Array<OneD, Array<OneD, NekDouble>> u_star(m_velocity.size());
1852
1853 cout << "We apply the continuation method: " << endl;
1854
1855 for (size_t i = 0; i < m_velocity.size(); ++i)
1856 {
1857 u_N[i] = Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetTotPoints(),
1858 0.0);
1859 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1860 u_N[i]);
1861
1862 RHS[i] = Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetTotPoints(),
1863 0.0);
1864 tmp_RHS[i] = Array<OneD, NekDouble>(
1865 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1866
1867 m_fields[m_velocity[i]]->PhysDeriv(i, u_N[i], tmp_RHS[i]);
1868 Vmath::Smul(tmp_RHS[i].size(), m_kinvis, tmp_RHS[i], 1, tmp_RHS[i], 1);
1869
1870 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1871 m_fields[m_velocity[i]]->SetWaveSpace(true);
1872 m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_RHS[i], RHS[i]);
1873 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1874 }
1875
1876 SetUpCoupledMatrix(0.0, u_N, true);
1877 SolveLinearNS(RHS);
1878
1879 for (size_t i = 0; i < m_velocity.size(); ++i)
1880 {
1881 u_star[i] = Array<OneD, NekDouble>(
1882 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1883 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1884 u_star[i]);
1885
1886 // u_star(k+1) = u_N(k) + DeltaKinvis * u_star(k)
1887 Vmath::Smul(u_star[i].size(), m_kinvis, u_star[i], 1, u_star[i], 1);
1888 Vmath::Vadd(u_star[i].size(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1889
1890 m_fields[m_velocity[i]]->FwdTrans(
1891 u_star[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1892 }
1893
1894 m_kinvis -= m_kinvis * m_KinvisPercentage / 100;
1895}
Nektar::CoupledLinearNS::SolveLinearNS
void SolveLinearNS(const Array< OneD, Array< OneD, NekDouble > > &forcing)
Definition: CoupledLinearNS.cpp:2087
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:363
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:765
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.cpp:354
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.cpp:245

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:168
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:66

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

◆ DefineForcingTerm()

void Nektar::CoupledLinearNS::DefineForcingTerm ( void  )

Definition at line 1702 of file CoupledLinearNS.cpp.

1703{
1704 m_ForcingTerm = Array<OneD, Array<OneD, NekDouble>>(m_velocity.size());
1705 m_ForcingTerm_Coeffs =
1706 Array<OneD, Array<OneD, NekDouble>>(m_velocity.size());
1707
1708 for (size_t i = 0; i < m_velocity.size(); ++i)
1709 {
1710 m_ForcingTerm[i] = Array<OneD, NekDouble>(
1711 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1712 m_ForcingTerm_Coeffs[i] =
1713 Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetNcoeffs(), 0.0);
1714 }
1715
1716 if (m_session->DefinesFunction("ForcingTerm"))
1717 {
1718 std::vector<std::string> fieldStr;
1719 for (size_t i = 0; i < m_velocity.size(); ++i)
1720 {
1721 fieldStr.push_back(
1722 m_boundaryConditions->GetVariable(m_velocity[i]));
1723 }
1724 GetFunction("ForcingTerm")->Evaluate(fieldStr, m_ForcingTerm);
1725 for (size_t i = 0; i < m_velocity.size(); ++i)
1726 {
1727 m_fields[m_velocity[i]]->FwdTransLocalElmt(m_ForcingTerm[i],
1728 m_ForcingTerm_Coeffs[i]);
1729 }
1730 }
1731 else
1732 {
1733 cout << "'ForcingTerm' section has not been defined in the input file "
1734 "=> forcing=0"
1735 << endl;
1736 }
1737}
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:356
Nektar::SolverUtils::EquationSystem::GetNcoeffs
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:707
Nektar::SolverUtils::EquationSystem::m_boundaryConditions
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
Definition: EquationSystem.h:365
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:737

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:399

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 1929 of file CoupledLinearNS.cpp.

1932{
1933 Array<OneD, Array<OneD, NekDouble>> Eval_Adv(m_velocity.size());
1934 Array<OneD, Array<OneD, NekDouble>> tmp_DerVel(m_velocity.size());
1935 Array<OneD, Array<OneD, NekDouble>> AdvTerm(m_velocity.size());
1936 Array<OneD, Array<OneD, NekDouble>> ViscTerm(m_velocity.size());
1937 Array<OneD, Array<OneD, NekDouble>> Forc(m_velocity.size());
1938
1939 for (size_t i = 0; i < m_velocity.size(); ++i)
1940 {
1941 Eval_Adv[i] = Array<OneD, NekDouble>(
1942 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1943 tmp_DerVel[i] = Array<OneD, NekDouble>(
1944 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1945
1946 AdvTerm[i] = Array<OneD, NekDouble>(
1947 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1948 ViscTerm[i] = Array<OneD, NekDouble>(
1949 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1950 Forc[i] = Array<OneD, NekDouble>(
1951 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1952 outarray[i] = Array<OneD, NekDouble>(
1953 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1954
1955 m_fields[m_velocity[i]]->PhysDeriv(i, Velocity[i], tmp_DerVel[i]);
1956
1957 Vmath::Smul(tmp_DerVel[i].size(), m_kinvis, tmp_DerVel[i], 1,
1958 tmp_DerVel[i], 1);
1959 }
1960
1961 EvaluateAdvectionTerms(Velocity, Eval_Adv, m_time);
1962
1963 for (size_t i = 0; i < m_velocity.size(); ++i)
1964 {
1965 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1966 m_fields[m_velocity[i]]->SetWaveSpace(true);
1967 m_fields[m_velocity[i]]->IProductWRTBase(Eval_Adv[i],
1968 AdvTerm[i]); //(w, (u.grad)u)
1969 m_fields[m_velocity[i]]->IProductWRTDerivBase(
1970 i, tmp_DerVel[i], ViscTerm[i]); //(grad w, grad u)
1971 m_fields[m_velocity[i]]->IProductWRTBase(m_ForcingTerm[i],
1972 Forc[i]); //(w, f)
1973 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1974
1975 Vmath::Vsub(outarray[i].size(), outarray[i], 1, AdvTerm[i], 1,
1976 outarray[i], 1);
1977 Vmath::Vsub(outarray[i].size(), outarray[i], 1, ViscTerm[i], 1,
1978 outarray[i], 1);
1979
1980 Vmath::Vadd(outarray[i].size(), outarray[i], 1, Forc[i], 1, outarray[i],
1981 1);
1982 }
1983}
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:371
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.cpp:414

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 1985 of file CoupledLinearNS.cpp.

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

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 1897 of file CoupledLinearNS.cpp.

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

References Nektar::IncNavierStokes::m_velocity.

◆ L2Norm()

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

Definition at line 1914 of file CoupledLinearNS.cpp.

1916{
1917 for (size_t i = 0; i < m_velocity.size(); ++i)
1918 {
1919 outarray[i] = 0.0;
1920 for (size_t j = 0; j < inarray[i].size(); ++j)
1921 {
1922 outarray[i] += inarray[i][j] * inarray[i][j];
1923 }
1924 outarray[i] = sqrt(outarray[i]);
1925 cout << "L2Norm[" << i << "] = " << outarray[i] << endl;
1926 }
1927}
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:404
Nektar::SolverUtils::EquationSystem::m_HomogeneousType
enum HomogeneousType m_HomogeneousType
Definition: EquationSystem.h:454
Nektar::SolverUtils::EquationSystem::m_LhomZ
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
Definition: EquationSystem.h:458
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:130
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:385
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:61
Nektar::MultiRegions::DirCartesianMap
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:90
Nektar::NekConstants::kNekUnsetDouble
static const NekDouble kNekUnsetDouble
Definition: NektarUnivConsts.hpp:46
Nektar::StdRegions::StdExpansionSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1775
Nektar::StdRegions::ConstFactorMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:408
Nektar::VarcoeffHashingTest::factors
StdRegions::ConstFactorMap factors
Definition: TestVarcoeffHashing.cpp:52
Nektar::DNekScalMat
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
Definition: NekMatrixFwd.hpp:67
Nektar::DNekScalMatSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:68
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.cpp:207
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:513
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:487

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:249
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:462

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 1672 of file CoupledLinearNS.cpp.

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

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 2134 of file CoupledLinearNS.cpp.

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

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 2087 of file CoupledLinearNS.cpp.

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

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 1739 of file CoupledLinearNS.cpp.

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

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 boost::ignore_unused(time);
1534
1535 size_t i;
1536 Array<OneD, Array<OneD, NekDouble>> F(m_nConvectiveFields);
1537 NekDouble lambda = 1.0 / aii_Dt;
1538 static NekDouble lambda_store;
1539 Array<OneD, Array<OneD, NekDouble>> forcing(m_velocity.size());
1540 // Matrix solution
1541 if (fabs(lambda_store - lambda) > 1e-10)
1542 {
1543 SetUpCoupledMatrix(lambda);
1544 lambda_store = lambda;
1545 }
1546
1547 // Forcing for advection solve
1548 for (i = 0; i < m_velocity.size(); ++i)
1549 {
1550 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1551 m_fields[m_velocity[i]]->SetWaveSpace(true);
1552 m_fields[m_velocity[i]]->IProductWRTBase(
1553 inarray[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1554 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1555 Vmath::Smul(m_fields[m_velocity[i]]->GetNcoeffs(), lambda,
1556 m_fields[m_velocity[i]]->GetCoeffs(), 1,
1557 m_fields[m_velocity[i]]->UpdateCoeffs(), 1);
1558 forcing[i] = m_fields[m_velocity[i]]->GetCoeffs();
1559 }
1560
1561 SolveLinearNS(forcing);
1562
1563 for (i = 0; i < m_velocity.size(); ++i)
1564 {
1565 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1566 outarray[i]);
1567 }
1568}
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:88
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineImplicitSolve
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationSchemeOperators.h:104
Nektar::SolverUtils::EquationSystem::m_lambda
NekDouble m_lambda
Lambda constant in real system if one required.
Definition: EquationSystem.h:380
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:675
Nektar::SolverUtils::UnsteadySystem::m_ode
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:79
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 1592 of file CoupledLinearNS.cpp.

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

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:49

References Nektar::SolverUtils::AddSummaryItem().

◆ v_GetForceDimension()

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

Implements Nektar::IncNavierStokes.

Definition at line 2481 of file CoupledLinearNS.cpp.

2482{
2483 return m_session->GetVariables().size();
2484}

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

◆ v_InitObject()

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

Init object for UnsteadySystem class.

Initialization object for UnsteadySystem class.

Reimplemented from Nektar::IncNavierStokes.

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:1985
Nektar::IncNavierStokes::v_InitObject
virtual void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
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:144
Nektar::SolverUtils::AdvectionSystem::m_advObject
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
Definition: AdvectionSystem.h:72
Nektar::SolverUtils::EquationSystem::m_useFFT
bool m_useFFT
Flag to determine if FFT is used for homogeneous transform.
Definition: EquationSystem.h:410
Nektar::SolverUtils::EquationSystem::m_homogen_dealiasing
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
Definition: EquationSystem.h:415
Nektar::SpatialDomains::ExpansionInfoMap
std::map< int, ExpansionInfoShPtr > ExpansionInfoMap
Definition: MeshGraph.h:143
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, and Nektar::IncNavierStokes::v_InitObject().

◆ 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 1667 of file CoupledLinearNS.cpp.

1668{
1669 return true;
1670}

◆ 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 2421 of file CoupledLinearNS.cpp.

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

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, Nektar::SolverUtils::EquationSystem::ParallelInTime(), and Nektar::SolverUtils::EquationSystem::WriteFld().

◆ v_TransCoeffToPhys()

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

Virtual function for transformation to physical space.

Reimplemented from Nektar::IncNavierStokes.

Definition at line 1570 of file CoupledLinearNS.cpp.

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

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::IncNavierStokes.

Definition at line 1581 of file CoupledLinearNS.cpp.

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

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:198
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:103

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:583
Nektar::eCoupledLinearisedNS
@ eCoupledLinearisedNS
Definition: IncNavierStokes.h:57

Definition at line 170 of file CoupledLinearNS.h.

Generated on Wed Jan 31 2024 08:26:52 for Nektar++ by doxygen 1.9.4