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

#include <CoupledLinearNS.h>

Inheritance diagram for Nektar::CoupledLinearNS:
[legend]

Public Member Functions

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

Static Public Member Functions

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

Public Attributes

Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm
 
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm_Coeffs
 
Array< OneD, CoupledLocalToGlobalC0ContMapSharedPtr > m_locToGloMap
 
- Public Attributes inherited from Nektar::SolverUtils::UnsteadySystem
NekDouble m_cflSafetyFactor
 CFL safety factor (comprise between 0 to 1). More...
 
NekDouble m_cflNonAcoustic
 
NekDouble m_CFLGrowth
 CFL growth rate. More...
 
NekDouble m_CFLEnd
 maximun cfl in cfl growth More...
 

Static Public Attributes

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

Protected Member Functions

 CoupledLinearNS (const LibUtilities::SessionReaderSharedPtr &pSesssion, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 
virtual void v_InitObject (bool DeclareField=true) override
 Init object for UnsteadySystem class. More...
 
- Protected Member Functions inherited from Nektar::IncNavierStokes
 IncNavierStokes (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Constructor. More...
 
EquationType GetEquationType (void)
 
void EvaluateAdvectionTerms (const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
 
void WriteModalEnergy (void)
 
void SetBoundaryConditions (NekDouble time)
 time dependent boundary conditions updating More...
 
void SetRadiationBoundaryForcing (int fieldid)
 Set Radiation forcing term. More...
 
void SetZeroNormalVelocity ()
 Set Normal Velocity Component to Zero. More...
 
void SetWomersleyBoundary (const int fldid, const int bndid)
 Set Womersley Profile if specified. More...
 
void SetUpWomersley (const int fldid, const int bndid, std::string womstr)
 Set Up Womersley details. More...
 
void SetMovingReferenceFrameBCs (const NekDouble &time)
 Set the moving reference frame boundary conditions. More...
 
void SetMRFWallBCs (const NekDouble &time)
 
void SetMRFDomainVelBCs (const NekDouble &time)
 
virtual MultiRegions::ExpListSharedPtr v_GetPressure () override
 
virtual Array< OneD, NekDouble > v_GetMaxStdVelocity (const NekDouble SpeedSoundFactor) override
 
virtual bool v_PreIntegrate (int step) override
 
- Protected Member Functions inherited from Nektar::SolverUtils::AdvectionSystem
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate (int step) override
 
- Protected Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
SOLVER_UTILS_EXPORT UnsteadySystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises UnsteadySystem class members. More...
 
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator ()
 Get the maximum timestep estimator for cfl control. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_GetTimeStep (const Array< OneD, const Array< OneD, NekDouble >> &inarray)
 Return the timestep to be used for the next step in the time-marching loop. More...
 
virtual SOLVER_UTILS_EXPORT bool v_RequireFwdTrans ()
 
virtual SOLVER_UTILS_EXPORT void v_SteadyStateResidual (int step, Array< OneD, NekDouble > &L2)
 
SOLVER_UTILS_EXPORT void CheckForRestartTime (NekDouble &time, int &nchk)
 
SOLVER_UTILS_EXPORT void SVVVarDiffCoeff (const Array< OneD, Array< OneD, NekDouble >> vel, StdRegions::VarCoeffMap &varCoeffMap)
 Evaluate the SVV diffusion coefficient according to Moura's paper where it should proportional to h time velocity. More...
 
virtual SOLVER_UTILS_EXPORT bool v_UpdateTimeStepCheck ()
 
SOLVER_UTILS_EXPORT void DoDummyProjection (const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
 Perform dummy projection. More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Initialises EquationSystem class members. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Virtual function for the L_inf error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
 Virtual function for the L_2 error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT void v_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)
 
- Protected Member Functions inherited from Nektar::SolverUtils::FluidInterface
virtual SOLVER_UTILS_EXPORT void v_SetMovingFrameProjectionMat (const boost::numeric::ublas::matrix< NekDouble > &vProjMat)
 
virtual SOLVER_UTILS_EXPORT void v_GetMovingFrameProjectionMat (boost::numeric::ublas::matrix< NekDouble > &vProjMat)
 

Private Member Functions

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

Private Attributes

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

Friends

class MemoryManager< CoupledLinearNS >
 

Additional Inherited Members

- Protected Types inherited from Nektar::SolverUtils::EquationSystem
enum  HomogeneousType { eHomogeneous1D , eHomogeneous2D , eHomogeneous3D , eNotHomogeneous }
 Parameter for homogeneous expansions. More...
 
- Protected Attributes inherited from Nektar::IncNavierStokes
ExtrapolateSharedPtr m_extrapolation
 
std::ofstream m_mdlFile
 modal energy file More...
 
bool m_SmoothAdvection
 bool to identify if advection term smoothing is requested More...
 
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
 Forcing terms. More...
 
int m_nConvectiveFields
 Number of fields to be convected;. More...
 
Array< OneD, int > m_velocity
 int which identifies which components of m_fields contains the velocity (u,v,w); More...
 
MultiRegions::ExpListSharedPtr m_pressure
 Pointer to field holding pressure field. More...
 
NekDouble m_kinvis
 Kinematic viscosity. More...
 
int m_energysteps
 dump energy to file at steps time More...
 
EquationType m_equationType
 equation type; More...
 
Array< OneD, Array< OneD, int > > m_fieldsBCToElmtID
 Mapping from BCs to Elmt IDs. More...
 
Array< OneD, Array< OneD, int > > m_fieldsBCToTraceID
 Mapping from BCs to Elmt Edge IDs. More...
 
Array< OneD, Array< OneD, NekDouble > > m_fieldsRadiationFactor
 RHS Factor for Radiation Condition. More...
 
int m_intSteps
 Number of time integration steps AND Order of extrapolation for pressure boundary conditions. More...
 
Array< OneD, NekDouble > m_pivotPoint
 
std::map< int, std::map< int, WomersleyParamsSharedPtr > > m_womersleyParams
 Womersley parameters if required. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::AdvectionSystem
SolverUtils::AdvectionSharedPtr m_advObject
 Advection term. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::UnsteadySystem
int m_abortSteps
 Number of steps between checks for abort conditions. More...
 
int m_filtersInfosteps
 Number of time steps between outputting filters information. More...
 
int m_nanSteps
 
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
 Wrapper to the time integration scheme. More...
 
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use. More...
 
NekDouble m_epsilon
 
bool m_explicitDiffusion
 Indicates if explicit or implicit treatment of diffusion is used. More...
 
bool m_explicitAdvection
 Indicates if explicit or implicit treatment of advection is used. More...
 
bool m_explicitReaction
 Indicates if explicit or implicit treatment of reaction is used. More...
 
bool m_homoInitialFwd
 Flag to determine if simulation should start in homogeneous forward transformed state. More...
 
NekDouble m_steadyStateTol
 Tolerance to which steady state should be evaluated at. More...
 
int m_steadyStateSteps
 Check for steady state at step interval. More...
 
NekDouble m_steadyStateRes = 1.0
 
NekDouble m_steadyStateRes0 = 1.0
 
Array< OneD, Array< OneD, NekDouble > > m_previousSolution
 Storage for previous solution for steady-state check. More...
 
std::ofstream m_errFile
 
std::vector< int > m_intVariables
 
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters
 
NekDouble m_filterTimeWarning
 Number of time steps between outputting status information. More...
 
NekDouble m_TimeIntegLambda = 0.0
 coefff of spacial derivatives(rhs or m_F in GLM) in calculating the residual of the whole equation(used in unsteady time integrations) More...
 
bool m_flagImplicitItsStatistics
 
bool m_flagImplicitSolver = false
 
Array< OneD, NekDouble > m_magnitdEstimat
 estimate the magnitude of each conserved varibles More...
 
Array< OneD, NekDouble > m_locTimeStep
 local time step(notice only for jfnk other see m_cflSafetyFactor) More...
 
NekDouble m_inArrayNorm = -1.0
 
int m_TotLinItePerStep = 0
 
int m_StagesPerStep = 1
 
bool m_flagUpdatePreconMat
 
int m_maxLinItePerNewton
 
int m_TotNewtonIts = 0
 
int m_TotLinIts = 0
 
int m_TotImpStages = 0
 
bool m_CalcPhysicalAV = true
 flag to update artificial viscosity More...
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
bool m_verbose
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader. More...
 
std::map< std::string, SolverUtils::SessionFunctionSharedPtr > m_sessionFunctions
 Map of known SessionFunctions. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
 Array holding all dependent variables. More...
 
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
 Pointer to boundary conditions object. More...
 
SpatialDomains::MeshGraphSharedPtr m_graph
 Pointer to graph defining mesh. More...
 
std::string m_sessionName
 Name of the session. More...
 
NekDouble m_time
 Current time of simulation. More...
 
int m_initialStep
 Number of the step where the simulation should begin. More...
 
NekDouble m_fintime
 Finish time of the simulation. More...
 
NekDouble m_timestep
 Time step size. More...
 
NekDouble m_timestepMax = -1.0
 Time step size. More...
 
NekDouble m_lambda
 Lambda constant in real system if one required. More...
 
NekDouble m_checktime
 Time between checkpoints. More...
 
NekDouble m_lastCheckTime
 
NekDouble m_TimeIncrementFactor
 
int m_nchk
 Number of checkpoints written so far. More...
 
int m_steps
 Number of steps to take. More...
 
int m_checksteps
 Number of steps between checkpoints. More...
 
int m_infosteps
 Number of time steps between outputting status information. More...
 
int m_pararealIter
 Number of parareal time iteration. More...
 
int m_spacedim
 Spatial dimension (>= expansion dim). More...
 
int m_expdim
 Expansion dimension. More...
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform. More...
 
bool m_useInitialCondition
 Flag to determine if IC are used. More...
 
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations. More...
 
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More...
 
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous. More...
 
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction. More...
 
Array< OneD, bool > m_checkIfSystemSingular
 Flag to indicate if the fields should be checked for singularity. More...
 
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
 Map to identify relevant solver info to dump in output fields. More...
 
Array< OneD, NekDouble > m_movingFrameVelsxyz
 Moving frame of reference velocities. More...
 
Array< OneD, NekDouble > m_movingFrameTheta
 Moving frame of reference angles with respect to the. More...
 
boost::numeric::ublas::matrix< NekDouble > m_movingFrameProjMat
 Projection matrix for transformation between inertial and moving. More...
 
int m_NumQuadPointsError
 Number of Quadrature points used to work out the error. More...
 
enum HomogeneousType m_HomogeneousType
 
NekDouble m_LhomX
 physical length in X direction (if homogeneous) More...
 
NekDouble m_LhomY
 physical length in Y direction (if homogeneous) More...
 
NekDouble m_LhomZ
 physical length in Z direction (if homogeneous) More...
 
int m_npointsX
 number of points in X direction (if homogeneous) More...
 
int m_npointsY
 number of points in Y direction (if homogeneous) More...
 
int m_npointsZ
 number of points in Z direction (if homogeneous) More...
 
int m_HomoDirec
 number of homogenous directions More...
 
- Static Protected Attributes inherited from Nektar::IncNavierStokes
static std::string eqTypeLookupIds []
 
- Static Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
static std::string equationSystemTypeLookupIds []
 

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

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

Member Function Documentation

◆ Continuation()

void Nektar::CoupledLinearNS::Continuation ( void  )

Definition at line 1842 of file CoupledLinearNS.cpp.

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

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

Referenced by v_DoSolve().

◆ create()

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

Creates an instance of this class.

Definition at line 92 of file CoupledLinearNS.h.

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

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

◆ DefineForcingTerm()

void Nektar::CoupledLinearNS::DefineForcingTerm ( void  )

Definition at line 1698 of file CoupledLinearNS.cpp.

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

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

1514 {
1515  // evaluate convection terms
1516  EvaluateAdvectionTerms(inarray, outarray, time);
1517 
1518  for (auto &x : m_forcing)
1519  {
1520  x->Apply(m_fields, outarray, outarray, time);
1521  }
1522 }
Nektar::IncNavierStokes::EvaluateAdvectionTerms
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Definition: IncNavierStokes.cpp:399
Nektar::IncNavierStokes::m_forcing
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
Definition: IncNavierStokes.h:194

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

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

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

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

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

Referenced by v_InitObject().

◆ InfNorm()

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

Definition at line 1893 of file CoupledLinearNS.cpp.

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

References Nektar::IncNavierStokes::m_velocity.

◆ L2Norm()

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

Definition at line 1910 of file CoupledLinearNS.cpp.

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

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

Referenced by SolveSteadyNavierStokes().

◆ SetUpCoupledMatrix() [1/2]

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

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

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

Consider the linearised Navier-Stokes element matrix denoted as

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

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

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

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

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

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

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

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

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

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

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

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

which if we multiply out the matrix equation we obtain

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

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

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

where

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

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

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

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

and

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

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

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

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

Definition at line 400 of file CoupledLinearNS.cpp.

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

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

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

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, and m_zeroMode.

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

◆ Solve()

void Nektar::CoupledLinearNS::Solve ( void  )

Definition at line 1668 of file CoupledLinearNS.cpp.

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

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

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

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

◆ SolveLinearNS() [2/2]

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

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

Parameters
forcingA list of forcing functions for each velocity component

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

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

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

This vector is further manipulated into

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

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

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

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

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

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

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

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

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

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

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

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

Definition at line 2083 of file CoupledLinearNS.cpp.

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

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, and Vmath::Zero().

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

◆ SolveSteadyNavierStokes()

void Nektar::CoupledLinearNS::SolveSteadyNavierStokes ( void  )

Definition at line 1735 of file CoupledLinearNS.cpp.

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

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

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

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 ( void  )
overrideprivatevirtual

Sets up initial conditions.

Sets the initial conditions.

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 1359 of file CoupledLinearNS.cpp.

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

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

◆ v_DoSolve()

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

Solves an unsteady problem.

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

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 1588 of file CoupledLinearNS.cpp.

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

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

◆ v_GenerateSummary()

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

Print a summary of time stepping parameters.

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

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 1354 of file CoupledLinearNS.cpp.

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

References Nektar::SolverUtils::AddSummaryItem().

◆ v_GetForceDimension()

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

Implements Nektar::IncNavierStokes.

Definition at line 2476 of file CoupledLinearNS.cpp.

2477 {
2478  return m_session->GetVariables().size();
2479 }

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

◆ v_InitObject()

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

Init object for UnsteadySystem class.

Initialization object for UnsteadySystem class.

Reimplemented from Nektar::IncNavierStokes.

Definition at line 67 of file CoupledLinearNS.cpp.

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

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

◆ v_NegatedOp()

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

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1663 of file CoupledLinearNS.cpp.

1664 {
1665  return true;
1666 }

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

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

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

◆ v_TransCoeffToPhys()

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

Virtual function for transformation to physical space.

Reimplemented from Nektar::IncNavierStokes.

Definition at line 1566 of file CoupledLinearNS.cpp.

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

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

◆ v_TransPhysToCoeff()

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

Virtual function for transformation to coefficient space.

Reimplemented from Nektar::IncNavierStokes.

Definition at line 1577 of file CoupledLinearNS.cpp.

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

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

Friends And Related Function Documentation

◆ MemoryManager< CoupledLinearNS >

friend class MemoryManager< CoupledLinearNS >
friend

Definition at line 83 of file CoupledLinearNS.h.

Member Data Documentation

◆ className

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

Name of class.

Definition at line 102 of file CoupledLinearNS.h.

◆ m_counter

int Nektar::CoupledLinearNS::m_counter
private

Definition at line 174 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 175 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

◆ m_kinvisMin

NekDouble Nektar::CoupledLinearNS::m_kinvisMin
private

Definition at line 181 of file CoupledLinearNS.h.

Referenced by v_DoInitialise(), and v_DoSolve().

◆ m_KinvisPercentage

NekDouble Nektar::CoupledLinearNS::m_KinvisPercentage
private

Definition at line 183 of file CoupledLinearNS.h.

Referenced by Continuation(), and v_DoInitialise().

◆ m_kinvisStep

NekDouble Nektar::CoupledLinearNS::m_kinvisStep
private

Definition at line 182 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 185 of file CoupledLinearNS.h.

Referenced by SetUpCoupledMatrix(), and SolveLinearNS().

◆ m_MatrixSetUpStep

int Nektar::CoupledLinearNS::m_MatrixSetUpStep
private

Definition at line 180 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

◆ m_maxIt

int Nektar::CoupledLinearNS::m_maxIt
private

Definition at line 177 of file CoupledLinearNS.h.

Referenced by v_DoInitialise().

◆ m_Restart

int Nektar::CoupledLinearNS::m_Restart
private

Definition at line 178 of file CoupledLinearNS.h.

Referenced by v_DoInitialise().

◆ m_tol

NekDouble Nektar::CoupledLinearNS::m_tol
private

Definition at line 176 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 172 of file CoupledLinearNS.h.

Referenced by SetUpCoupledMatrix(), and v_InitObject().

Generated on Sat Mar 25 2023 12:35:56 for Nektar++ by doxygen 1.9.1