Nektar++
Loading...
Searching...
No Matches
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.
 
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.
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep (const Array< OneD, const Array< OneD, NekDouble > > &inarray)
 Calculate the larger time-step mantaining the problem stable.
 
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.
 
SOLVER_UTILS_EXPORT LibUtilities::TimeIntegrationSchemeOperators & GetTimeIntegrationSchemeOperators ()
 Returns the time integration scheme operators.
 
- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor.
 
SOLVER_UTILS_EXPORT void InitObject (bool DeclareField=true)
 Initialises the members of this object.
 
SOLVER_UTILS_EXPORT void DoInitialise (bool dumpInitialConditions=true)
 Perform any initialisation necessary before solving the problem.
 
SOLVER_UTILS_EXPORT void DoSolve ()
 Solve the problem.
 
SOLVER_UTILS_EXPORT void TransCoeffToPhys ()
 Transform from coefficient to physical space.
 
SOLVER_UTILS_EXPORT void TransPhysToCoeff ()
 Transform from physical to coefficient space.
 
SOLVER_UTILS_EXPORT void Output ()
 Perform output operations after solve.
 
SOLVER_UTILS_EXPORT std::string GetSessionName ()
 Get Session name.
 
template<class T >
std::shared_ptr< T > as ()
 
SOLVER_UTILS_EXPORT void ResetSessionName (std::string newname)
 Reset Session name.
 
SOLVER_UTILS_EXPORT LibUtilities::SessionReaderSharedPtr GetSession ()
 Get Session name.
 
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure ()
 Get pressure field if available.
 
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.
 
SOLVER_UTILS_EXPORT void SetLambda (NekDouble lambda)
 Set parameter m_lambda.
 
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction (std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
 Get a SessionFunction by name.
 
SOLVER_UTILS_EXPORT void SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 Initialise the data in the dependent fields.
 
SOLVER_UTILS_EXPORT void EvaluateExactSolution (int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 Evaluates an exact solution.
 
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.
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, bool Normalised=false)
 Compute the L2 error of the fields.
 
SOLVER_UTILS_EXPORT NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation.
 
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.
 
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].
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n)
 Write checkpoint file of m_fields.
 
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.
 
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow (const int n)
 Write base flow file of m_fields.
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname)
 Write field data to the given filename.
 
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.
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 Input field data from the given file.
 
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.
 
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.
 
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.
 
SOLVER_UTILS_EXPORT void SessionSummary (SummaryList &vSummary)
 Write out a session summary.
 
SOLVER_UTILS_EXPORT Array< OneD, MultiRegions::ExpListSharedPtr > & UpdateFields ()
 
SOLVER_UTILS_EXPORT LibUtilities::FieldMetaDataMap & UpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated.
 
SOLVER_UTILS_EXPORT NekDouble GetTime ()
 Return final time.
 
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.
 
SOLVER_UTILS_EXPORT bool NegatedOp ()
 Identify if operator is negated in DoSolve.
 
- 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.
 
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)
 
SOLVER_UTILS_EXPORT void ResetMatricesNormal (Array< OneD, Array< OneD, NekDouble > > &traceNormals)
 
SOLVER_UTILS_EXPORT void UpdateNormalsFlag ()
 
const Array< OneD, const Array< OneD, NekDouble > > & GetGridVelocity ()
 
bool & GetUpdateNormalsFlag ()
 
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)
 
SOLVER_UTILS_EXPORT void ExtraFldOutputGrid (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.
 
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.
 
SOLVER_UTILS_EXPORT void GetPressure (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure)
 Extract array with pressure from physfield.
 
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.
 
SOLVER_UTILS_EXPORT void GetAeroForce (Array< OneD, NekDouble > forces)
 Get aerodynamic force and moment.
 

Static Public Member Functions

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

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.
 
- 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.
 
- Protected Member Functions inherited from Nektar::IncNavierStokes
 IncNavierStokes (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Constructor.
 
 ~IncNavierStokes () override=default
 
void v_InitObject (bool DeclareField=true) override
 Initialisation object for EquationSystem.
 
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
 
void SetRadiationBoundaryForcing (int fieldid)
 Set Radiation forcing term.
 
void SetZeroNormalVelocity ()
 Set Normal Velocity Component to Zero.
 
void SetWomersleyBoundary (const int fldid, const int bndid)
 Set Womersley Profile if specified.
 
void SetUpWomersley (const int fldid, const int bndid, std::string womstr)
 Set Up Womersley details.
 
MultiRegions::ExpListSharedPtr v_GetPressure () override
 
void v_TransCoeffToPhys (void) override
 Virtual function for transformation to physical space.
 
void v_TransPhysToCoeff (void) override
 Virtual function for transformation to coefficient space.
 
Array< OneD, NekDouble > v_GetMaxStdVelocity (const NekDouble SpeedSoundFactor) override
 
bool v_PreIntegrate (int step) override
 
- Protected Member Functions inherited from Nektar::SolverUtils::AdvectionSystem
SOLVER_UTILS_EXPORT bool v_PostIntegrate (int step) override
 
- Protected Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
SOLVER_UTILS_EXPORT UnsteadySystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises UnsteadySystem class members.
 
SOLVER_UTILS_EXPORT void v_InitObject (bool DeclareField=true) override
 Init object for UnsteadySystem class.
 
SOLVER_UTILS_EXPORT void v_DoSolve () override
 Solves an unsteady problem.
 
virtual SOLVER_UTILS_EXPORT void v_PrintStatusInformation (const int step, const NekDouble cpuTime)
 Print Status Information.
 
virtual SOLVER_UTILS_EXPORT void v_PrintSummaryStatistics (const NekDouble intTime)
 Print Summary Statistics.
 
SOLVER_UTILS_EXPORT void v_DoInitialise (bool dumpInitialConditions=true) override
 Sets up initial conditions.
 
SOLVER_UTILS_EXPORT void v_GenerateSummary (SummaryList &s) override
 Print a summary of time stepping parameters.
 
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.
 
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.
 
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.
 
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.
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises EquationSystem class members.
 
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.
 
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.
 
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.
 
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_ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 

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

1829{
1830 Array<OneD, Array<OneD, NekDouble>> u_N(m_velocity.size());
1831 Array<OneD, Array<OneD, NekDouble>> tmp_RHS(m_velocity.size());
1832 Array<OneD, Array<OneD, NekDouble>> RHS(m_velocity.size());
1833 Array<OneD, Array<OneD, NekDouble>> u_star(m_velocity.size());
1834
1835 std::cout << "We apply the continuation method: " << std::endl;
1836
1837 for (size_t i = 0; i < m_velocity.size(); ++i)
1838 {
1839 u_N[i] = Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetTotPoints(),
1840 0.0);
1841 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1842 u_N[i]);
1843
1844 RHS[i] = Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetTotPoints(),
1845 0.0);
1846 tmp_RHS[i] = Array<OneD, NekDouble>(
1847 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1848
1849 m_fields[m_velocity[i]]->PhysDeriv(i, u_N[i], tmp_RHS[i]);
1850 Vmath::Smul(tmp_RHS[i].size(), m_kinvis, tmp_RHS[i], 1, tmp_RHS[i], 1);
1851
1852 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1853 m_fields[m_velocity[i]]->SetWaveSpace(true);
1854 m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_RHS[i], RHS[i]);
1855 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1856 }
1857
1858 SetUpCoupledMatrix(0.0, u_N, true);
1859 SolveLinearNS(RHS);
1860
1861 for (size_t i = 0; i < m_velocity.size(); ++i)
1862 {
1863 u_star[i] = Array<OneD, NekDouble>(
1864 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1865 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1866 u_star[i]);
1867
1868 // u_star(k+1) = u_N(k) + DeltaKinvis * u_star(k)
1869 Vmath::Smul(u_star[i].size(), m_kinvis, u_star[i], 1, u_star[i], 1);
1870 Vmath::Vadd(u_star[i].size(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1871
1872 m_fields[m_velocity[i]]->FwdTrans(
1873 u_star[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1874 }
1875
1876 m_kinvis -= m_kinvis * m_KinvisPercentage / 100;
1877}
Nektar::CoupledLinearNS::SolveLinearNS
void SolveLinearNS(const Array< OneD, Array< OneD, NekDouble > > &forcing)
Definition CoupledLinearNS.cpp:2069
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
Nektar::SolverUtils::EquationSystemSharedPtr
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
Definition EquationSystem.h:65
Nektar::UnitTests::p
std::vector< double > p(NPUPPER)

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

◆ DefineForcingTerm()

void Nektar::CoupledLinearNS::DefineForcingTerm ( void  )

Definition at line 1682 of file CoupledLinearNS.cpp.

1683{
1684 m_ForcingTerm = Array<OneD, Array<OneD, NekDouble>>(m_velocity.size());
1685 m_ForcingTerm_Coeffs =
1686 Array<OneD, Array<OneD, NekDouble>>(m_velocity.size());
1687
1688 for (size_t i = 0; i < m_velocity.size(); ++i)
1689 {
1690 m_ForcingTerm[i] = Array<OneD, NekDouble>(
1691 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1692 m_ForcingTerm_Coeffs[i] =
1693 Array<OneD, NekDouble>(m_fields[m_velocity[i]]->GetNcoeffs(), 0.0);
1694 }
1695
1696 if (m_session->DefinesFunction("ForcingTerm"))
1697 {
1698 std::vector<std::string> fieldStr;
1699 for (size_t i = 0; i < m_velocity.size(); ++i)
1700 {
1701 fieldStr.push_back(
1702 m_boundaryConditions->GetVariable(m_velocity[i]));
1703 }
1704 GetFunction("ForcingTerm")->Evaluate(fieldStr, m_ForcingTerm);
1705 for (size_t i = 0; i < m_velocity.size(); ++i)
1706 {
1707 m_fields[m_velocity[i]]->FwdTransLocalElmt(m_ForcingTerm[i],
1708 m_ForcingTerm_Coeffs[i]);
1709 }
1710 }
1711 else
1712 {
1713 std::cout
1714 << "'ForcingTerm' section has not been defined in the input file "
1715 "=> forcing=0"
1716 << std::endl;
1717 }
1718}
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:748

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

1496{
1497 // evaluate convection terms
1498 EvaluateAdvectionTerms(inarray, outarray, time);
1499
1500 for (auto &x : m_forcing)
1501 {
1502 x->Apply(m_fields, outarray, outarray, time);
1503 }
1504}
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 1911 of file CoupledLinearNS.cpp.

1914{
1915 Array<OneD, Array<OneD, NekDouble>> Eval_Adv(m_velocity.size());
1916 Array<OneD, Array<OneD, NekDouble>> tmp_DerVel(m_velocity.size());
1917 Array<OneD, Array<OneD, NekDouble>> AdvTerm(m_velocity.size());
1918 Array<OneD, Array<OneD, NekDouble>> ViscTerm(m_velocity.size());
1919 Array<OneD, Array<OneD, NekDouble>> Forc(m_velocity.size());
1920
1921 for (size_t i = 0; i < m_velocity.size(); ++i)
1922 {
1923 Eval_Adv[i] = Array<OneD, NekDouble>(
1924 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1925 tmp_DerVel[i] = Array<OneD, NekDouble>(
1926 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1927
1928 AdvTerm[i] = Array<OneD, NekDouble>(
1929 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1930 ViscTerm[i] = Array<OneD, NekDouble>(
1931 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1932 Forc[i] = Array<OneD, NekDouble>(
1933 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1934 outarray[i] = Array<OneD, NekDouble>(
1935 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1936
1937 m_fields[m_velocity[i]]->PhysDeriv(i, Velocity[i], tmp_DerVel[i]);
1938
1939 Vmath::Smul(tmp_DerVel[i].size(), m_kinvis, tmp_DerVel[i], 1,
1940 tmp_DerVel[i], 1);
1941 }
1942
1943 EvaluateAdvectionTerms(Velocity, Eval_Adv, m_time);
1944
1945 for (size_t i = 0; i < m_velocity.size(); ++i)
1946 {
1947 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1948 m_fields[m_velocity[i]]->SetWaveSpace(true);
1949 m_fields[m_velocity[i]]->IProductWRTBase(Eval_Adv[i],
1950 AdvTerm[i]); //(w, (u.grad)u)
1951 m_fields[m_velocity[i]]->IProductWRTDerivBase(
1952 i, tmp_DerVel[i], ViscTerm[i]); //(grad w, grad u)
1953 m_fields[m_velocity[i]]->IProductWRTBase(m_ForcingTerm[i],
1954 Forc[i]); //(w, f)
1955 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1956
1957 Vmath::Vsub(outarray[i].size(), outarray[i], 1, AdvTerm[i], 1,
1958 outarray[i], 1);
1959 Vmath::Vsub(outarray[i].size(), outarray[i], 1, ViscTerm[i], 1,
1960 outarray[i], 1);
1961
1962 Vmath::Vadd(outarray[i].size(), outarray[i], 1, Forc[i], 1, outarray[i],
1963 1);
1964 }
1965}
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 1967 of file CoupledLinearNS.cpp.

1969{
1970 SpatialDomains::ExpansionInfoMapShPtr returnval;
1971 returnval =
1972 MemoryManager<SpatialDomains::ExpansionInfoMap>::AllocateSharedPtr();
1973
1974 int nummodes;
1975
1976 for (auto &expMapIter : VelExp)
1977 {
1978 LibUtilities::BasisKeyVector BasisVec;
1979
1980 for (size_t i = 0; i < expMapIter.second->m_basisKeyVector.size(); ++i)
1981 {
1982 LibUtilities::BasisKey B = expMapIter.second->m_basisKeyVector[i];
1983 nummodes = B.GetNumModes();
1984 ASSERTL0(nummodes > 3,
1985 "Velocity polynomial space not sufficiently high (>= 4)");
1986 // Should probably set to be an orthogonal basis.
1987 LibUtilities::BasisKey newB(B.GetBasisType(), nummodes - 2,
1988 B.GetPointsKey());
1989 BasisVec.push_back(newB);
1990 }
1991
1992 // Put new expansion into list.
1993 SpatialDomains::ExpansionInfoShPtr expansionElementShPtr =
1994 MemoryManager<SpatialDomains::ExpansionInfo>::AllocateSharedPtr(
1995 expMapIter.second->m_geomPtr, BasisVec);
1996 (*returnval)[expMapIter.first] = expansionElementShPtr;
1997 }
1998
1999 // Save expansion into graph.
2000 m_graph->SetExpansionInfo("p", returnval);
2001
2002 return *returnval;
2003}
ASSERTL0
#define ASSERTL0(condition, msg)
Definition ErrorUtil.hpp:216
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:186
Nektar::SpatialDomains::ExpansionInfoShPtr
std::shared_ptr< ExpansionInfo > ExpansionInfoShPtr
Definition MeshGraph.h:183

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

1881{
1882 for (size_t i = 0; i < m_velocity.size(); ++i)
1883 {
1884 outarray[i] = 0.0;
1885 for (size_t j = 0; j < inarray[i].size(); ++j)
1886 {
1887 if (inarray[i][j] > outarray[i])
1888 {
1889 outarray[i] = inarray[i][j];
1890 }
1891 }
1892 std::cout << "InfNorm[" << i << "] = " << outarray[i] << std::endl;
1893 }
1894}

References Nektar::IncNavierStokes::m_velocity.

◆ L2Norm()

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

Definition at line 1896 of file CoupledLinearNS.cpp.

1898{
1899 for (size_t i = 0; i < m_velocity.size(); ++i)
1900 {
1901 outarray[i] = 0.0;
1902 for (size_t j = 0; j < inarray[i].size(); ++j)
1903 {
1904 outarray[i] += inarray[i][j] * inarray[i][j];
1905 }
1906 outarray[i] = std::sqrt(outarray[i]);
1907 std::cout << "L2Norm[" << i << "] = " << outarray[i] << std::endl;
1908 }
1909}

References Nektar::IncNavierStokes::m_velocity.

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_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:1783
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:929
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
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::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:250
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 1652 of file CoupledLinearNS.cpp.

1653{
1654 const size_t ncmpt = m_velocity.size();
1655 Array<OneD, Array<OneD, NekDouble>> forcing_phys(ncmpt);
1656 Array<OneD, Array<OneD, NekDouble>> forcing(ncmpt);
1657
1658 for (size_t i = 0; i < ncmpt; ++i)
1659 {
1660 forcing_phys[i] =
1661 Array<OneD, NekDouble>(m_fields[m_velocity[0]]->GetNpoints(), 0.0);
1662 forcing[i] =
1663 Array<OneD, NekDouble>(m_fields[m_velocity[0]]->GetNcoeffs(), 0.0);
1664 }
1665
1666 for (auto &x : m_forcing)
1667 {
1668 const NekDouble time = 0;
1669 x->Apply(m_fields, forcing_phys, forcing_phys, time);
1670 }
1671 for (size_t i = 0; i < ncmpt; ++i)
1672 {
1673 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1674 m_fields[i]->SetWaveSpace(true);
1675 m_fields[i]->IProductWRTBase(forcing_phys[i], forcing[i]);
1676 m_fields[i]->SetWaveSpace(waveSpace);
1677 }
1678
1679 SolveLinearNS(forcing);
1680}
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 2116 of file CoupledLinearNS.cpp.

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

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

2071{
2072 size_t i;
2073 Array<OneD, MultiRegions::ExpListSharedPtr> vel_fields(m_velocity.size());
2074 Array<OneD, Array<OneD, NekDouble>> force(m_velocity.size());
2075
2076 // Impose Dirichlet conditions on velocity fields
2077 for (i = 0; i < m_velocity.size(); ++i)
2078 {
2079 Vmath::Zero(m_fields[i]->GetNcoeffs(), m_fields[i]->UpdateCoeffs(), 1);
2080 m_fields[i]->ImposeDirichletConditions(m_fields[i]->UpdateCoeffs());
2081 }
2082
2083 if (m_HomogeneousType == eHomogeneous1D)
2084 {
2085 int ncoeffsplane = m_fields[m_velocity[0]]->GetPlane(0)->GetNcoeffs();
2086 for (int n = 0; n < m_npointsZ / 2; ++n)
2087 {
2088 // Get the a Fourier mode of velocity and forcing components.
2089 for (i = 0; i < m_velocity.size(); ++i)
2090 {
2091 vel_fields[i] = m_fields[m_velocity[i]]->GetPlane(2 * n);
2092 // Note this needs to correlate with how we pass forcing
2093 force[i] = forcing[i] + 2 * n * ncoeffsplane;
2094 }
2095
2096 SolveLinearNS(force, vel_fields, m_pressure->GetPlane(2 * n), n);
2097 }
2098 for (i = 0; i < m_velocity.size(); ++i)
2099 {
2100 m_fields[m_velocity[i]]->SetPhysState(false);
2101 }
2102 m_pressure->SetPhysState(false);
2103 }
2104 else
2105 {
2106 for (i = 0; i < m_velocity.size(); ++i)
2107 {
2108 vel_fields[i] = m_fields[m_velocity[i]];
2109 // Note this needs to correlate with how we pass forcing
2110 force[i] = forcing[i];
2111 }
2112 SolveLinearNS(force, vel_fields, m_pressure);
2113 }
2114}

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

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

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

1510{
1511 size_t i;
1512 Array<OneD, Array<OneD, NekDouble>> F(m_nConvectiveFields);
1513 NekDouble lambda = 1.0 / aii_Dt;
1514 static NekDouble lambda_store;
1515 Array<OneD, Array<OneD, NekDouble>> forcing(m_velocity.size());
1516 // Matrix solution
1517 if (fabs(lambda_store - lambda) > 1e-10)
1518 {
1519 SetUpCoupledMatrix(lambda);
1520 lambda_store = lambda;
1521 }
1522
1523 // Forcing for advection solve
1524 for (i = 0; i < m_velocity.size(); ++i)
1525 {
1526 bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1527 m_fields[m_velocity[i]]->SetWaveSpace(true);
1528 m_fields[m_velocity[i]]->IProductWRTBase(
1529 inarray[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1530 m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1531 Vmath::Smul(m_fields[m_velocity[i]]->GetNcoeffs(), lambda,
1532 m_fields[m_velocity[i]]->GetCoeffs(), 1,
1533 m_fields[m_velocity[i]]->UpdateCoeffs(), 1);
1534 forcing[i] = m_fields[m_velocity[i]]->GetCoeffs();
1535 }
1536
1537 SolveLinearNS(forcing);
1538
1539 for (i = 0; i < m_velocity.size(); ++i)
1540 {
1541 m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
1542 outarray[i]);
1543 }
1544}
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

Virtual function for initialisation implementation.

By default, nothing needs initialising at the EquationSystem level.

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1361 of file CoupledLinearNS.cpp.

1362{
1363 switch (m_equationType)
1364 {
1365 case eUnsteadyStokes:
1366 case eUnsteadyNavierStokes:
1367 {
1368 // Could defind this from IncNavierStokes class?
1369 m_ode.DefineOdeRhs(&CoupledLinearNS::EvaluateAdvection, this);
1370
1371 m_ode.DefineImplicitSolve(
1372 &CoupledLinearNS::SolveUnsteadyStokesSystem, this);
1373
1374 // Set initial condition using time t=0
1375
1376 SetInitialConditions(0.0, dumpInitialConditions);
1377 break;
1378 }
1379 case eSteadyStokes:
1380 SetUpCoupledMatrix(0.0);
1381 break;
1382 case eSteadyOseen:
1383 {
1384 Array<OneD, Array<OneD, NekDouble>> AdvField(m_velocity.size());
1385 for (size_t i = 0; i < m_velocity.size(); ++i)
1386 {
1387 AdvField[i] = Array<OneD, NekDouble>(
1388 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1389 }
1390
1391 ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1392 "Advection Velocity section must be defined in "
1393 "session file.");
1394
1395 std::vector<std::string> fieldStr;
1396 for (size_t i = 0; i < m_velocity.size(); ++i)
1397 {
1398 fieldStr.push_back(
1399 m_boundaryConditions->GetVariable(m_velocity[i]));
1400 }
1401 GetFunction("AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1402
1403 SetUpCoupledMatrix(0.0, AdvField, false);
1404 }
1405 break;
1406 case eSteadyNavierStokes:
1407 {
1408 m_session->LoadParameter("KinvisMin", m_kinvisMin);
1409 m_session->LoadParameter("KinvisPercentage", m_KinvisPercentage);
1410 m_session->LoadParameter("Tolerence", m_tol);
1411 m_session->LoadParameter("MaxIteration", m_maxIt);
1412 m_session->LoadParameter("MatrixSetUpStep", m_MatrixSetUpStep);
1413 m_session->LoadParameter("Restart", m_Restart);
1414
1415 DefineForcingTerm();
1416
1417 if (m_Restart == 1)
1418 {
1419 ASSERTL0(m_session->DefinesFunction("Restart"),
1420 "Restart section must be defined in session file.");
1421
1422 Array<OneD, Array<OneD, NekDouble>> Restart(m_velocity.size());
1423 for (size_t i = 0; i < m_velocity.size(); ++i)
1424 {
1425 Restart[i] = Array<OneD, NekDouble>(
1426 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1427 }
1428 std::vector<std::string> fieldStr;
1429 for (size_t i = 0; i < m_velocity.size(); ++i)
1430 {
1431 fieldStr.push_back(
1432 m_boundaryConditions->GetVariable(m_velocity[i]));
1433 }
1434 GetFunction("Restart")->Evaluate(fieldStr, Restart);
1435
1436 for (size_t i = 0; i < m_velocity.size(); ++i)
1437 {
1438 m_fields[m_velocity[i]]->FwdTransLocalElmt(
1439 Restart[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1440 }
1441 std::cout << "Saving the RESTART file for m_kinvis = "
1442 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1443 << std::endl;
1444 }
1445 else // We solve the Stokes Problem
1446 {
1447
1448 SetUpCoupledMatrix(0.0);
1449 m_initialStep = true;
1450 m_counter = 1;
1451 // SolveLinearNS(m_ForcingTerm_Coeffs);
1452 Solve();
1453 m_initialStep = false;
1454 std::cout << "Saving the Stokes Flow for m_kinvis = "
1455 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1456 << std::endl;
1457 }
1458 }
1459 break;
1460 case eSteadyLinearisedNS:
1461 {
1462 SetInitialConditions(0.0, dumpInitialConditions);
1463
1464 Array<OneD, Array<OneD, NekDouble>> AdvField(m_velocity.size());
1465 for (size_t i = 0; i < m_velocity.size(); ++i)
1466 {
1467 AdvField[i] = Array<OneD, NekDouble>(
1468 m_fields[m_velocity[i]]->GetTotPoints(), 0.0);
1469 }
1470
1471 ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1472 "Advection Velocity section must be defined in "
1473 "session file.");
1474
1475 std::vector<std::string> fieldStr;
1476 for (size_t i = 0; i < m_velocity.size(); ++i)
1477 {
1478 fieldStr.push_back(
1479 m_boundaryConditions->GetVariable(m_velocity[i]));
1480 }
1481 GetFunction("AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1482
1483 SetUpCoupledMatrix(m_lambda, AdvField, true);
1484 }
1485 break;
1486 case eNoEquationType:
1487 default:
1488 ASSERTL0(false,
1489 "Unknown or undefined equation type for CoupledLinearNS");
1490 }
1491}
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:1506
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:1493
Nektar::IncNavierStokes::m_equationType
EquationType m_equationType
equation type;
Definition IncNavierStokes.h:188
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::eUnsteadyNavierStokes
@ eUnsteadyNavierStokes
Definition IncNavierStokes.h:75
Nektar::eSteadyLinearisedNS
@ eSteadyLinearisedNS
Definition IncNavierStokes.h:72

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

Virtual function for solve implementation.

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1568 of file CoupledLinearNS.cpp.

1569{
1570 switch (m_equationType)
1571 {
1572 case eUnsteadyStokes:
1573 case eUnsteadyNavierStokes:
1574 // AdvanceInTime(m_steps);
1575 UnsteadySystem::v_DoSolve();
1576 break;
1577 case eSteadyStokes:
1578 case eSteadyOseen:
1579 case eSteadyLinearisedNS:
1580 {
1581 Solve();
1582 break;
1583 }
1584 case eSteadyNavierStokes:
1585 {
1586 LibUtilities::Timer Generaltimer;
1587 Generaltimer.Start();
1588
1589 int Check(0);
1590
1591 // Saving the init datas
1592 Checkpoint_Output(Check);
1593 Check++;
1594
1595 std::cout
1596 << "We execute INITIALLY SolveSteadyNavierStokes for m_kinvis "
1597 "= "
1598 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1599 << std::endl;
1600 SolveSteadyNavierStokes();
1601
1602 while (m_kinvis > m_kinvisMin)
1603 {
1604 if (Check == 1)
1605 {
1606 std::cout
1607 << "We execute SolveSteadyNavierStokes for m_kinvis = "
1608 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1609 << std::endl;
1610 SolveSteadyNavierStokes();
1611 Checkpoint_Output(Check);
1612 Check++;
1613 }
1614
1615 Continuation();
1616
1617 if (m_kinvis > m_kinvisMin)
1618 {
1619 std::cout
1620 << "We execute SolveSteadyNavierStokes for m_kinvis = "
1621 << m_kinvis << " (<=> Re = " << 1 / m_kinvis << ")"
1622 << std::endl;
1623 SolveSteadyNavierStokes();
1624 Checkpoint_Output(Check);
1625 Check++;
1626 }
1627 }
1628
1629 Generaltimer.Stop();
1630 std::cout << "\nThe total calculation time is : "
1631 << Generaltimer.TimePerTest(1) / 60 << " minute(s). \n\n";
1632
1633 break;
1634 }
1635 case eNoEquationType:
1636 default:
1637 ASSERTL0(false,
1638 "Unknown or undefined equation type for CoupledLinearNS");
1639 }
1640}
Nektar::SolverUtils::EquationSystem::Checkpoint_Output
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
Definition EquationSystem.cpp:1260
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 &  l)
overrideprivatevirtual

Virtual function for generating summary information.

Reimplemented from Nektar::SolverUtils::EquationSystem.

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

2464{
2465 return m_session->GetVariables().size();
2466}

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:1967
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:142
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:184
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 1647 of file CoupledLinearNS.cpp.

1648{
1649 return true;
1650}

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

2404{
2405 std::vector<Array<OneD, NekDouble>> fieldcoeffs(m_fields.size() + 1);
2406 std::vector<std::string> variables(m_fields.size() + 1);
2407 size_t i;
2408
2409 for (i = 0; i < m_fields.size(); ++i)
2410 {
2411 fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
2412 variables[i] = m_boundaryConditions->GetVariable(i);
2413 }
2414
2415 fieldcoeffs[i] = Array<OneD, NekDouble>(m_fields[0]->GetNcoeffs());
2416 // project pressure field to velocity space
2417 if (m_singleMode == true)
2418 {
2419 Array<OneD, NekDouble> tmpfieldcoeffs(m_fields[0]->GetNcoeffs() / 2);
2420 m_pressure->GetPlane(0)->BwdTrans(
2421 m_pressure->GetPlane(0)->GetCoeffs(),
2422 m_pressure->GetPlane(0)->UpdatePhys());
2423 m_pressure->GetPlane(1)->BwdTrans(
2424 m_pressure->GetPlane(1)->GetCoeffs(),
2425 m_pressure->GetPlane(1)->UpdatePhys());
2426 m_fields[0]->GetPlane(0)->FwdTransLocalElmt(
2427 m_pressure->GetPlane(0)->GetPhys(), fieldcoeffs[i]);
2428 m_fields[0]->GetPlane(1)->FwdTransLocalElmt(
2429 m_pressure->GetPlane(1)->GetPhys(), tmpfieldcoeffs);
2430 for (int e = 0; e < m_fields[0]->GetNcoeffs() / 2; e++)
2431 {
2432 fieldcoeffs[i][e + m_fields[0]->GetNcoeffs() / 2] =
2433 tmpfieldcoeffs[e];
2434 }
2435 }
2436 else
2437 {
2438 m_pressure->BwdTrans(m_pressure->GetCoeffs(), m_pressure->UpdatePhys());
2439 m_fields[0]->FwdTransLocalElmt(m_pressure->GetPhys(), fieldcoeffs[i]);
2440 }
2441 variables[i] = "p";
2442
2443 if (!m_comm->IsParallelInTime())
2444 {
2445 WriteFld(m_sessionName + ".fld", m_fields[0], fieldcoeffs, variables);
2446 }
2447 else
2448 {
2449 std::string newdir = m_sessionName + ".pit";
2450 if (!fs::is_directory(newdir))
2451 {
2452 fs::create_directory(newdir);
2453 }
2454 WriteFld(
2455 newdir + "/" + m_sessionName + "_" +
2456 std::to_string(m_windowPIT * m_comm->GetTimeComm()->GetSize() +
2457 m_comm->GetTimeComm()->GetRank() + 1) +
2458 ".fld",
2459 m_fields[0], fieldcoeffs, variables);
2460 }
2461}
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:1322
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 1546 of file CoupledLinearNS.cpp.

1547{
1548 size_t nfields = m_fields.size();
1549 for (size_t k = 0; k < nfields; ++k)
1550 {
1551 // Backward Transformation in physical space for time evolution
1552 m_fields[k]->BwdTrans(m_fields[k]->GetCoeffs(),
1553 m_fields[k]->UpdatePhys());
1554 }
1555}

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

1558{
1559 size_t nfields = m_fields.size();
1560 for (size_t k = 0; k < nfields; ++k)
1561 {
1562 // Forward Transformation in physical space for time evolution
1563 m_fields[k]->FwdTransLocalElmt(m_fields[k]->GetPhys(),
1564 m_fields[k]->UpdateCoeffs());
1565 }
1566}

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

Friends And Related Symbol 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:196
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 Sun Nov 9 2025 20:17:02 for Nektar++ by doxygen 1.9.8