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
int GetNConvectiveFields (void)
 
void AddForcing (const SolverUtils::ForcingSharedPtr &pForce)
 
bool DefinedForcing (const std::string &sForce)
 
- Public Member Functions inherited from Nektar::SolverUtils::AdvectionSystem
SOLVER_UTILS_EXPORT AdvectionSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 
SOLVER_UTILS_EXPORT ~AdvectionSystem () override=default
 
SOLVER_UTILS_EXPORT AdvectionSharedPtr GetAdvObject ()
 Returns the advection object held by this instance. More...
 
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > GetElmtCFLVals (const bool FlagAcousticCFL=true)
 
SOLVER_UTILS_EXPORT NekDouble GetCFLEstimate (int &elmtid)
 
- Public Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
SOLVER_UTILS_EXPORT ~UnsteadySystem () override=default
 Destructor. More...
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep (const Array< OneD, const Array< OneD, NekDouble > > &inarray)
 Calculate the larger time-step mantaining the problem stable. More...
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void SetTimeStep (const NekDouble timestep)
 
SOLVER_UTILS_EXPORT void SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeSharedPtr & GetTimeIntegrationScheme ()
 Returns the time integration scheme. More...
 
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeOperators & GetTimeIntegrationSchemeOperators ()
 Returns the time integration scheme operators. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT void InitObject (bool DeclareField=true)
 Initialises the members of this object. More...
 
SOLVER_UTILS_EXPORT void DoInitialise (bool dumpInitialConditions=true)
 Perform any initialisation necessary before solving the problem. More...
 
SOLVER_UTILS_EXPORT void DoSolve ()
 Solve the problem. More...
 
SOLVER_UTILS_EXPORT void TransCoeffToPhys ()
 Transform from coefficient to physical space. More...
 
SOLVER_UTILS_EXPORT void TransPhysToCoeff ()
 Transform from physical to coefficient space. More...
 
SOLVER_UTILS_EXPORT void Output ()
 Perform output operations after solve. More...
 
SOLVER_UTILS_EXPORT std::string GetSessionName ()
 Get Session name. More...
 
template<class T >
std::shared_ptr< T > as ()
 
SOLVER_UTILS_EXPORT void ResetSessionName (std::string newname)
 Reset Session name. More...
 
SOLVER_UTILS_EXPORT LibUtilities::SessionReaderSharedPtr GetSession ()
 Get Session name. More...
 
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure ()
 Get pressure field if available. More...
 
SOLVER_UTILS_EXPORT void ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 
SOLVER_UTILS_EXPORT void PrintSummary (std::ostream &out)
 Print a summary of parameters and solver characteristics. More...
 
SOLVER_UTILS_EXPORT void SetLambda (NekDouble lambda)
 Set parameter m_lambda. More...
 
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction (std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
 Get a SessionFunction by name. More...
 
SOLVER_UTILS_EXPORT void SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 Initialise the data in the dependent fields. More...
 
SOLVER_UTILS_EXPORT void EvaluateExactSolution (int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 Evaluates an exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised=false)
 Compute the L2 error between fields and a given exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, bool Normalised=false)
 Compute the L2 error of the fields. More...
 
SOLVER_UTILS_EXPORT NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation. More...
 
SOLVER_UTILS_EXPORT NekDouble H1Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised=false)
 Compute the H1 error between fields and a given exact solution. More...
 
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > ErrorExtraPoints (unsigned int field)
 Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf]. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n)
 Write checkpoint file of m_fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write checkpoint file of custom data fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow (const int n)
 Write base flow file of m_fields. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname)
 Write field data to the given filename. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write input fields to the given filename. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 Input field data from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFldToMultiDomains (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int ndomains)
 Input field data from the given file to multiple domains. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, std::vector< std::string > &fieldStr, Array< OneD, Array< OneD, NekDouble > > &coeffs)
 Output a field. Input field data into array from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, MultiRegions::ExpListSharedPtr &pField, std::string &pFieldName)
 Output a field. Input field data into ExpList from the given file. More...
 
SOLVER_UTILS_EXPORT void SessionSummary (SummaryList &vSummary)
 Write out a session summary. More...
 
SOLVER_UTILS_EXPORT Array< OneD, MultiRegions::ExpListSharedPtr > & UpdateFields ()
 
SOLVER_UTILS_EXPORT LibUtilities::FieldMetaDataMap & UpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated. More...
 
SOLVER_UTILS_EXPORT NekDouble GetTime ()
 Return final time. More...
 
SOLVER_UTILS_EXPORT int GetNcoeffs ()
 
SOLVER_UTILS_EXPORT int GetNcoeffs (const int eid)
 
SOLVER_UTILS_EXPORT int GetNumExpModes ()
 
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp ()
 
SOLVER_UTILS_EXPORT int GetNvariables ()
 
SOLVER_UTILS_EXPORT const std::string GetVariable (unsigned int i)
 
SOLVER_UTILS_EXPORT int GetTraceTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTraceNpoints ()
 
SOLVER_UTILS_EXPORT int GetExpSize ()
 
SOLVER_UTILS_EXPORT int GetPhys_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetCoeff_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTotPoints (int n)
 
SOLVER_UTILS_EXPORT int GetNpoints ()
 
SOLVER_UTILS_EXPORT int GetSteps ()
 
SOLVER_UTILS_EXPORT void SetSteps (const int steps)
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void CopyFromPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void CopyToPhysField (const int i, const Array< OneD, const NekDouble > &input)
 
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > & UpdatePhysField (const int i)
 
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
 
SOLVER_UTILS_EXPORT void FwdTransFields ()
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int GetCheckpointNumber ()
 
SOLVER_UTILS_EXPORT void SetCheckpointNumber (int num)
 
SOLVER_UTILS_EXPORT int GetCheckpointSteps ()
 
SOLVER_UTILS_EXPORT void SetCheckpointSteps (int num)
 
SOLVER_UTILS_EXPORT int GetInfoSteps ()
 
SOLVER_UTILS_EXPORT void SetInfoSteps (int num)
 
SOLVER_UTILS_EXPORT void SetIterationNumberPIT (int num)
 
SOLVER_UTILS_EXPORT void SetWindowNumberPIT (int num)
 
SOLVER_UTILS_EXPORT Array< OneD, const Array< OneD, NekDouble > > GetTraceNormals ()
 
SOLVER_UTILS_EXPORT void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetTimeStep (const NekDouble timestep)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 
SOLVER_UTILS_EXPORT bool NegatedOp ()
 Identify if operator is negated in DoSolve. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::ALEHelper
virtual ~ALEHelper ()=default
 
virtual SOLVER_UTILS_EXPORT void v_ALEInitObject (int spaceDim, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
 
SOLVER_UTILS_EXPORT void InitObject (int spaceDim, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
 
virtual SOLVER_UTILS_EXPORT void v_UpdateGridVelocity (const NekDouble &time)
 
virtual SOLVER_UTILS_EXPORT void v_ALEPreMultiplyMass (Array< OneD, Array< OneD, NekDouble > > &fields)
 
SOLVER_UTILS_EXPORT void ALEDoElmtInvMass (Array< OneD, Array< OneD, NekDouble > > &traceNormals, Array< OneD, Array< OneD, NekDouble > > &fields, NekDouble time)
 Update m_fields with u^n by multiplying by inverse mass matrix. That's then used in e.g. checkpoint output and L^2 error calculation. More...
 
SOLVER_UTILS_EXPORT void ALEDoElmtInvMassBwdTrans (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
 
SOLVER_UTILS_EXPORT void MoveMesh (const NekDouble &time, Array< OneD, Array< OneD, NekDouble > > &traceNormals)
 
const Array< OneD, const Array< OneD, NekDouble > > & GetGridVelocity ()
 
SOLVER_UTILS_EXPORT const Array< OneD, const Array< OneD, NekDouble > > & GetGridVelocityTrace ()
 
SOLVER_UTILS_EXPORT void ExtraFldOutputGridVelocity (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 
- 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, const int step)
 
SOLVER_UTILS_EXPORT bool GetMovingFrameVelocities (Array< OneD, NekDouble > &vFrameVels, const int step)
 
SOLVER_UTILS_EXPORT void SetMovingFrameDisp (const Array< OneD, NekDouble > &vFrameDisp, const int step)
 
SOLVER_UTILS_EXPORT void SetMovingFramePivot (const Array< OneD, NekDouble > &vFramePivot)
 
SOLVER_UTILS_EXPORT bool GetMovingFrameDisp (Array< OneD, NekDouble > &vFrameDisp, const int step)
 
SOLVER_UTILS_EXPORT void SetAeroForce (Array< OneD, NekDouble > forces)
 Set aerodynamic force and moment. More...
 
SOLVER_UTILS_EXPORT void GetAeroForce (Array< OneD, NekDouble > forces)
 Get aerodynamic force and moment. More...
 

Static Public Member Functions

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

Public Attributes

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

Static Public Attributes

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

Protected Member Functions

 CoupledLinearNS (const LibUtilities::SessionReaderSharedPtr &pSesssion, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 
void v_InitObject (bool DeclareField=true) override
 Initialisation object for EquationSystem. More...
 
- Protected Member Functions inherited from Nektar::IncNavierStokes
 IncNavierStokes (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Constructor. More...
 
 ~IncNavierStokes () override=default
 
void v_InitObject (bool DeclareField=true) override
 Initialisation object for EquationSystem. More...
 
void v_GetPressure (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure) override
 
void v_GetDensity (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &density) override
 
bool v_HasConstantDensity () override
 
void v_GetVelocity (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity) override
 
void v_SetMovingFrameVelocities (const Array< OneD, NekDouble > &vFrameVels, const int step) override
 
bool v_GetMovingFrameVelocities (Array< OneD, NekDouble > &vFrameVels, const int step) override
 
void v_SetMovingFrameDisp (const Array< OneD, NekDouble > &vFrameDisp, const int step) override
 
void v_SetMovingFramePivot (const Array< OneD, NekDouble > &vFramePivot) override
 
bool v_GetMovingFrameDisp (Array< OneD, NekDouble > &vFrameDisp, const int step) override
 
void v_SetAeroForce (Array< OneD, NekDouble > forces) override
 
void v_GetAeroForce (Array< OneD, NekDouble > forces) override
 
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...
 
MultiRegions::ExpListSharedPtr v_GetPressure () override
 
void v_TransCoeffToPhys (void) override
 Virtual function for transformation to physical space. More...
 
void v_TransPhysToCoeff (void) override
 Virtual function for transformation to coefficient space. More...
 
virtual int v_GetForceDimension ()=0
 
Array< OneD, NekDouble > v_GetMaxStdVelocity (const NekDouble SpeedSoundFactor) override
 
bool v_PreIntegrate (int step) override
 
- Protected Member Functions inherited from Nektar::SolverUtils::AdvectionSystem
SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareField=true) override
 Initialisation object for EquationSystem. More...
 
SOLVER_UTILS_EXPORT bool v_PostIntegrate (int step) override
 
virtual SOLVER_UTILS_EXPORT Array< OneD, NekDouble > v_GetMaxStdVelocity (const NekDouble SpeedSoundFactor=1.0)
 
- Protected Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
SOLVER_UTILS_EXPORT UnsteadySystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises UnsteadySystem class members. More...
 
SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareField=true) override
 Init object for UnsteadySystem class. More...
 
SOLVER_UTILS_EXPORT void v_DoSolve () override
 Solves an unsteady problem. More...
 
virtual SOLVER_UTILS_EXPORT void v_PrintStatusInformation (const int step, const NekDouble cpuTime)
 Print Status Information. More...
 
virtual SOLVER_UTILS_EXPORT void v_PrintSummaryStatistics (const NekDouble intTime)
 Print Summary Statistics. More...
 
SOLVER_UTILS_EXPORT void v_DoInitialise (bool dumpInitialConditions=true) override
 Sets up initial conditions. More...
 
SOLVER_UTILS_EXPORT void v_GenerateSummary (SummaryList &s) override
 Print a summary of time stepping parameters. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_GetTimeStep (const Array< OneD, const Array< OneD, NekDouble > > &inarray)
 Return the timestep to be used for the next step in the time-marching loop. More...
 
virtual SOLVER_UTILS_EXPORT bool v_PreIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_RequireFwdTrans ()
 
virtual SOLVER_UTILS_EXPORT void v_SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
virtual SOLVER_UTILS_EXPORT bool v_UpdateTimeStepCheck ()
 
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator ()
 Get the maximum timestep estimator for cfl control. More...
 
SOLVER_UTILS_EXPORT void CheckForRestartTime (NekDouble &time, int &nchk)
 
SOLVER_UTILS_EXPORT void SVVVarDiffCoeff (const Array< OneD, Array< OneD, NekDouble > > vel, StdRegions::VarCoeffMap &varCoeffMap)
 Evaluate the SVV diffusion coefficient according to Moura's paper where it should proportional to h time velocity. More...
 
SOLVER_UTILS_EXPORT void DoDummyProjection (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
 Perform dummy projection. More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises EquationSystem class members. More...
 
virtual SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareFeld=true)
 Initialisation object for EquationSystem. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoInitialise (bool dumpInitialConditions=true)
 Virtual function for initialisation implementation. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoSolve ()
 Virtual function for solve implementation. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Virtual function for the L_inf error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
 Virtual function for the L_2 error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_H1Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
 Virtual function for the H_1 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)
 
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, const int step)
 
virtual SOLVER_UTILS_EXPORT bool v_GetMovingFrameVelocities (Array< OneD, NekDouble > &vFrameVels, const int step)
 
virtual SOLVER_UTILS_EXPORT void v_SetMovingFrameDisp (const Array< OneD, NekDouble > &vFrameDisp, const int step)
 
virtual SOLVER_UTILS_EXPORT void v_SetMovingFramePivot (const Array< OneD, NekDouble > &vFramePivot)
 
virtual SOLVER_UTILS_EXPORT bool v_GetMovingFrameDisp (Array< OneD, NekDouble > &vFrameDisp, const int step)
 
virtual SOLVER_UTILS_EXPORT void v_SetAeroForce (Array< OneD, NekDouble > forces)
 
virtual SOLVER_UTILS_EXPORT void v_GetAeroForce (Array< OneD, NekDouble > forces)
 

Static Protected Attributes

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

Private Member Functions

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

Private Attributes

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

Friends

class MemoryManager< CoupledLinearNS >
 

Additional Inherited Members

- Protected Types inherited from Nektar::SolverUtils::EquationSystem
enum  HomogeneousType { eHomogeneous1D , eHomogeneous2D , eHomogeneous3D , eNotHomogeneous }
 Parameter for homogeneous expansions. More...
 
- Protected Attributes inherited from Nektar::IncNavierStokes
ExtrapolateSharedPtr m_extrapolation
 
IncBoundaryConditionsSharedPtr m_IncNavierStokesBCs
 
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
 pivot point for moving reference frame More...
 
Array< OneD, NekDouble > m_aeroForces
 
std::map< int, std::map< int, WomersleyParamsSharedPtr > > m_womersleyParams
 Womersley parameters if required. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::AdvectionSystem
SolverUtils::AdvectionSharedPtr m_advObject
 Advection term. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::UnsteadySystem
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
 Wrapper to the time integration scheme. More...
 
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use. More...
 
Array< OneD, Array< OneD, NekDouble > > m_previousSolution
 Storage for previous solution for steady-state check. More...
 
std::vector< int > m_intVariables
 
NekDouble m_cflSafetyFactor
 CFL safety factor (comprise between 0 to 1). More...
 
NekDouble m_CFLGrowth
 CFL growth rate. More...
 
NekDouble m_CFLEnd
 Maximun cfl in cfl growth. More...
 
int m_abortSteps
 Number of steps between checks for abort conditions. More...
 
bool m_explicitDiffusion
 Indicates if explicit or implicit treatment of diffusion is used. More...
 
bool m_explicitAdvection
 Indicates if explicit or implicit treatment of advection is used. More...
 
bool m_explicitReaction
 Indicates if explicit or implicit treatment of reaction is used. More...
 
int m_steadyStateSteps
 Check for steady state at step interval. More...
 
NekDouble m_steadyStateTol
 Tolerance to which steady state should be evaluated at. More...
 
int m_filtersInfosteps
 Number of time steps between outputting filters information. More...
 
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters
 
bool m_homoInitialFwd
 Flag to determine if simulation should start in homogeneous forward transformed state. More...
 
std::ofstream m_errFile
 
NekDouble m_epsilon
 Diffusion coefficient. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
bool m_verbose
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader. More...
 
std::map< std::string, SolverUtils::SessionFunctionSharedPtr > m_sessionFunctions
 Map of known SessionFunctions. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
 Array holding all dependent variables. More...
 
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
 Pointer to boundary conditions object. More...
 
SpatialDomains::MeshGraphSharedPtr m_graph
 Pointer to graph defining mesh. More...
 
std::string m_sessionName
 Name of the session. More...
 
NekDouble m_time
 Current time of simulation. More...
 
int m_initialStep
 Number of the step where the simulation should begin. More...
 
NekDouble m_fintime
 Finish time of the simulation. More...
 
NekDouble m_timestep
 Time step size. More...
 
NekDouble m_lambda
 Lambda constant in real system if one required. More...
 
NekDouble m_checktime
 Time between checkpoints. More...
 
NekDouble m_lastCheckTime
 
NekDouble m_TimeIncrementFactor
 
int m_nchk
 Number of checkpoints written so far. More...
 
int m_steps
 Number of steps to take. More...
 
int m_checksteps
 Number of steps between checkpoints. More...
 
int m_infosteps
 Number of time steps between outputting status information. More...
 
int m_iterPIT = 0
 Number of parallel-in-time time iteration. More...
 
int m_windowPIT = 0
 Index of windows for parallel-in-time time iteration. More...
 
int m_spacedim
 Spatial dimension (>= expansion dim). More...
 
int m_expdim
 Expansion dimension. More...
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform. More...
 
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations. More...
 
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More...
 
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous. More...
 
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction. More...
 
Array< OneD, bool > m_checkIfSystemSingular
 Flag to indicate if the fields should be checked for singularity. More...
 
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
 Map to identify relevant solver info to dump in output fields. More...
 
Array< OneD, NekDouble > m_movingFrameData
 Moving reference frame status in the inertial frame X, Y, Z, Theta_x, Theta_y, Theta_z, U, V, W, Omega_x, Omega_y, Omega_z, A_x, A_y, A_z, DOmega_x, DOmega_y, DOmega_z, pivot_x, pivot_y, pivot_z. More...
 
std::vector< std::string > m_strFrameData
 variable name in m_movingFrameData 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...
 
- Protected Attributes inherited from Nektar::SolverUtils::ALEHelper
Array< OneD, MultiRegions::ExpListSharedPtr > m_fieldsALE
 
Array< OneD, Array< OneD, NekDouble > > m_gridVelocity
 
Array< OneD, Array< OneD, NekDouble > > m_gridVelocityTrace
 
std::vector< ALEBaseShPtr > m_ALEs
 
bool m_ALESolver = false
 
bool m_ImplicitALESolver = false
 
NekDouble m_prevStageTime = 0.0
 
int m_spaceDim
 

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

64 : UnsteadySystem(pSession, pGraph), IncNavierStokes(pSession, pGraph),
65 m_zeroMode(false)
66{
67}
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:81
Nektar::SolverUtils::UnsteadySystem::UnsteadySystem
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Initialises UnsteadySystem class members.
Definition: UnsteadySystem.cpp:74

Member Function Documentation

◆ Continuation()

void Nektar::CoupledLinearNS::Continuation ( void  )

Definition at line 1850 of file CoupledLinearNS.cpp.

1851{
1852 Array<OneD, Array<OneD, NekDouble>> u_N(m_velocity.size());
1853 Array<OneD, Array<OneD, NekDouble>> tmp_RHS(m_velocity.size());
1854 Array<OneD, Array<OneD, NekDouble>> RHS(m_velocity.size());
1855 Array<OneD, Array<OneD, NekDouble>> u_star(m_velocity.size());
1856
1857 std::cout << "We apply the continuation method: " << std::endl;
1858
1859 for (size_t i = 0; i < m_velocity.size(); ++i)
1860 {
1861 u_N[i] = Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetTotPoints(),
1862 0.0);
1863 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1864 u_N[i]);
1865
1866 RHS[i] = Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetTotPoints(),
1867 0.0);
1868 tmp_RHS[i] = Array<OneD, NekDouble>(
1869 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1870
1871 m_fields[m_velocity[i]]->PhysDeriv(i, u_N[i], tmp_RHS[i]);
1872 Vmath::Smul(tmp_RHS[i].size(), m_kinvis, tmp_RHS[i], 1, tmp_RHS[i], 1);
1873
1874 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1875 m_fields[m_velocity[i]]->SetWaveSpace(true);
1876 m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_RHS[i], RHS[i]);
1877 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1878 }
1879
1880 SetUpCoupledMatrix(0.0, u_N, true);
1881 SolveLinearNS(RHS);
1882
1883 for (size_t i = 0; i < m_velocity.size(); ++i)
1884 {
1885 u_star[i] = Array<OneD, NekDouble>(
1886 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1887 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1888 u_star[i]);
1889
1890 // u_star(k+1) = u_N(k) + DeltaKinvis * u_star(k)
1891 Vmath::Smul(u_star[i].size(), m_kinvis, u_star[i], 1, u_star[i], 1);
1892 Vmath::Vadd(u_star[i].size(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1893
1894 m_fields[m_velocity[i]]->FwdTrans(
1895 u_star[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1896 }
1897
1898 m_kinvis -= m_kinvis * m_KinvisPercentage / 100;
1899}
Nektar::CoupledLinearNS::SolveLinearNS
void SolveLinearNS(const Array< OneD, Array< OneD, NekDouble > > &forcing)
Definition: CoupledLinearNS.cpp:2091
Nektar::CoupledLinearNS::SetUpCoupledMatrix
void SetUpCoupledMatrix(const NekDouble lambda=0.0, const Array< OneD, Array< OneD, NekDouble > > &Advfield=NullNekDoubleArrayOfArray, bool IsLinearNSEquation=true)
Definition: CoupledLinearNS.cpp:200
Nektar::IncNavierStokes::m_kinvis
NekDouble m_kinvis
Kinematic viscosity.
Definition: IncNavierStokes.h:183
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:178
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:467
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:341
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100

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

Referenced by v_DoSolve().

◆ create()

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

Creates an instance of this class.

Definition at line 92 of file CoupledLinearNS.h.

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

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

◆ DefineForcingTerm()

void Nektar::CoupledLinearNS::DefineForcingTerm ( void  )

Definition at line 1704 of file CoupledLinearNS.cpp.

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

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

Referenced by v_DoInitialise().

◆ EvaluateAdvection()

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

Definition at line 1515 of file CoupledLinearNS.cpp.

1518{
1519 // evaluate convection terms
1520 EvaluateAdvectionTerms(inarray, outarray, time);
1521
1522 for (auto &x : m_forcing)
1523 {
1524 x->Apply(m_fields, outarray, outarray, time);
1525 }
1526}
Nektar::IncNavierStokes::m_forcing
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
Definition: IncNavierStokes.h:171
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:298

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

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

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

Referenced by SolveSteadyNavierStokes().

◆ GenPressureExp()

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

Definition at line 1989 of file CoupledLinearNS.cpp.

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

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

Referenced by v_InitObject().

◆ InfNorm()

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

Definition at line 1901 of file CoupledLinearNS.cpp.

1903{
1904 for (size_t i = 0; i < m_velocity.size(); ++i)
1905 {
1906 outarray[i] = 0.0;
1907 for (size_t j = 0; j < inarray[i].size(); ++j)
1908 {
1909 if (inarray[i][j] > outarray[i])
1910 {
1911 outarray[i] = inarray[i][j];
1912 }
1913 }
1914 std::cout << "InfNorm[" << i << "] = " << outarray[i] << std::endl;
1915 }
1916}

References Nektar::IncNavierStokes::m_velocity.

◆ L2Norm()

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

Definition at line 1918 of file CoupledLinearNS.cpp.

1920{
1921 for (size_t i = 0; i < m_velocity.size(); ++i)
1922 {
1923 outarray[i] = 0.0;
1924 for (size_t j = 0; j < inarray[i].size(); ++j)
1925 {
1926 outarray[i] += inarray[i][j] * inarray[i][j];
1927 }
1928 outarray[i] = std::sqrt(outarray[i]);
1929 std::cout << "L2Norm[" << i << "] = " << outarray[i] << std::endl;
1930 }
1931}
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285

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

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

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

◆ SetUpCoupledMatrix() [2/2]

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

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

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

Definition at line 200 of file CoupledLinearNS.cpp.

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

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

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

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

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

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

◆ SolveLinearNS() [2/2]

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

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

Parameters
forcingA list of forcing functions for each velocity component

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

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

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

This vector is further manipulated into

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

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

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

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

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

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

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

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

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

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

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

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

Definition at line 2091 of file CoupledLinearNS.cpp.

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

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

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

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

Referenced by v_DoSolve().

◆ SolveUnsteadyStokesSystem()

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

Definition at line 1528 of file CoupledLinearNS.cpp.

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

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

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

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

◆ v_DoSolve()

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

Solves an unsteady problem.

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

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 1590 of file CoupledLinearNS.cpp.

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

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

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

References Nektar::SolverUtils::AddSummaryItem().

◆ v_GetForceDimension()

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

Implements Nektar::IncNavierStokes.

Definition at line 2485 of file CoupledLinearNS.cpp.

2486{
2487 return m_session->GetVariables().size();
2488}

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

◆ v_InitObject()

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

Initialisation object for EquationSystem.

Continuous field

Setting up the normals

Setting up the normals

Reimplemented from Nektar::SolverUtils::AdvectionSystem.

Definition at line 69 of file CoupledLinearNS.cpp.

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

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

1670{
1671 return true;
1672}

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

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

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

◆ v_TransCoeffToPhys()

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

Virtual function for transformation to physical space.

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1568 of file CoupledLinearNS.cpp.

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

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

◆ v_TransPhysToCoeff()

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

Virtual function for transformation to coefficient space.

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1579 of file CoupledLinearNS.cpp.

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

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

Friends And Related Function Documentation

◆ MemoryManager< CoupledLinearNS >

friend class MemoryManager< CoupledLinearNS >
friend

Definition at line 83 of file CoupledLinearNS.h.

Member Data Documentation

◆ className

std::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: BasicUtils/NekFactory.hpp:197
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:99

Name of class.

Definition at line 102 of file CoupledLinearNS.h.

◆ m_counter

int Nektar::CoupledLinearNS::m_counter
private

Definition at line 176 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

◆ m_ForcingTerm

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

Definition at line 159 of file CoupledLinearNS.h.

Referenced by DefineForcingTerm(), and EvaluateNewtonRHS().

◆ m_ForcingTerm_Coeffs

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

Definition at line 160 of file CoupledLinearNS.h.

Referenced by DefineForcingTerm().

◆ m_initialStep

bool Nektar::CoupledLinearNS::m_initialStep
private

Definition at line 177 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

◆ m_kinvisMin

NekDouble Nektar::CoupledLinearNS::m_kinvisMin
private

Definition at line 183 of file CoupledLinearNS.h.

Referenced by v_DoInitialise(), and v_DoSolve().

◆ m_KinvisPercentage

NekDouble Nektar::CoupledLinearNS::m_KinvisPercentage
private

Definition at line 185 of file CoupledLinearNS.h.

Referenced by Continuation(), and v_DoInitialise().

◆ m_kinvisStep

NekDouble Nektar::CoupledLinearNS::m_kinvisStep
private

Definition at line 184 of file CoupledLinearNS.h.

◆ m_locToGloMap

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

Definition at line 162 of file CoupledLinearNS.h.

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

◆ m_mat

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

Definition at line 187 of file CoupledLinearNS.h.

Referenced by SetUpCoupledMatrix(), and SolveLinearNS().

◆ m_MatrixSetUpStep

int Nektar::CoupledLinearNS::m_MatrixSetUpStep
private

Definition at line 182 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

◆ m_maxIt

int Nektar::CoupledLinearNS::m_maxIt
private

Definition at line 179 of file CoupledLinearNS.h.

Referenced by v_DoInitialise().

◆ m_Restart

int Nektar::CoupledLinearNS::m_Restart
private

Definition at line 180 of file CoupledLinearNS.h.

Referenced by v_DoInitialise().

◆ m_tol

NekDouble Nektar::CoupledLinearNS::m_tol
private

Definition at line 178 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

◆ m_zeroMode

bool Nektar::CoupledLinearNS::m_zeroMode
private

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

Definition at line 174 of file CoupledLinearNS.h.

Referenced by SetUpCoupledMatrix(), and v_InitObject().

◆ solverTypeLookupId

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

Definition at line 170 of file CoupledLinearNS.h.

Generated on Wed May 28 2025 19:03:43 for Nektar++ by doxygen 1.9.4