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::ExpansionInfoMapGenPressureExp (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)
 
Array< OneD, int > & GetVelocity (void)
 
void AddForcing (const SolverUtils::ForcingSharedPtr &pForce)
 
virtual void GetPressure (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure)
 Extract array with pressure from physfield. More...
 
virtual void GetDensity (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &density)
 Extract array with density from physfield. More...
 
virtual bool HasConstantDensity ()
 
virtual void GetVelocity (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
 Extract array with velocity from physfield. More...
 
- 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, NekDoubleGetElmtCFLVals (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 ()
 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, NekDoubleErrorExtraPoints (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::FieldMetaDataMapUpdateFieldMetaDataMap ()
 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, Array< OneD, NekDouble > &output)
 
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 void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::FluidInterface
virtual ~FluidInterface ()=default
 

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, CoupledLocalToGlobalC0ContMapSharedPtrm_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 = SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction("CoupledLinearisedNS", CoupledLinearNS::create)
 Name of class. More...
 

Protected Member Functions

 CoupledLinearNS (const LibUtilities::SessionReaderSharedPtr &pSesssion, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 
virtual void v_InitObject ()
 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...
 
virtual MultiRegions::ExpListSharedPtr v_GetPressure ()
 
virtual Array< OneD, NekDoublev_GetMaxStdVelocity (const NekDouble SpeedSoundFactor)
 
virtual bool v_PreIntegrate (int step)
 
- Protected Member Functions inherited from Nektar::SolverUtils::AdvectionSystem
virtual SOLVER_UTILS_EXPORT bool v_PostIntegrate (int step)
 
- 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 void v_AppendOutput1D (Array< OneD, Array< OneD, NekDouble >> &solution1D)
 Print the solution at each solution point in a txt file. 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 UpdateTimeStepCheck ()
 
- 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)
 

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)
 Print a summary of time stepping parameters. More...
 
virtual void v_DoInitialise (void)
 Sets up initial conditions. More...
 
virtual void v_DoSolve (void)
 Solves an unsteady problem. More...
 
virtual bool v_NegatedOp (void)
 
virtual void v_TransCoeffToPhys (void)
 Virtual function for transformation to physical space. More...
 
virtual void v_TransPhysToCoeff (void)
 Virtual function for transformation to coefficient space. More...
 
virtual void v_Output (void)
 
virtual int v_GetForceDimension (void)
 

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, CoupledSolverMatricesm_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::ForcingSharedPtrm_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...
 
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_infosteps
 Number of time steps between outputting status information. More...
 
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, NekDoublem_magnitdEstimat
 estimate the magnitude of each conserved varibles More...
 
Array< OneD, NekDoublem_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
 
bool m_root
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader. More...
 
std::map< std::string, SolverUtils::SessionFunctionSharedPtrm_sessionFunctions
 Map of known SessionFunctions. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
Array< OneD, MultiRegions::ExpListSharedPtrm_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_spacedim
 Spatial dimension (>= expansion dim). More...
 
int m_expdim
 Expansion dimension. More...
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform. More...
 
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations. More...
 
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More...
 
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous. More...
 
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction. More...
 
Array< OneD, bool > m_checkIfSystemSingular
 Flag to indicate if the fields should be checked for singularity. More...
 
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
 Map to identify relevant solver info to dump in output fields. More...
 
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::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 58 of file CoupledLinearNS.cpp.

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

Member Function Documentation

◆ Continuation()

void Nektar::CoupledLinearNS::Continuation ( void  )

Definition at line 1645 of file CoupledLinearNS.cpp.

1646  {
1647  Array<OneD, Array<OneD, NekDouble> > u_N(m_velocity.size());
1648  Array<OneD, Array<OneD, NekDouble> > tmp_RHS(m_velocity.size());
1649  Array<OneD, Array<OneD, NekDouble> > RHS(m_velocity.size());
1650  Array<OneD, Array<OneD, NekDouble> > u_star(m_velocity.size());
1651 
1652  cout << "We apply the continuation method: " <<endl;
1653 
1654  for(int i = 0; i < m_velocity.size(); ++i)
1655  {
1656  u_N[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1657  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), u_N[i]);
1658 
1659  RHS[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1660  tmp_RHS[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1661 
1662  m_fields[m_velocity[i]]->PhysDeriv(i, u_N[i], tmp_RHS[i]);
1663  Vmath::Smul(tmp_RHS[i].size(), m_kinvis, tmp_RHS[i], 1, tmp_RHS[i], 1);
1664 
1665  bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1666  m_fields[m_velocity[i]]->SetWaveSpace(true);
1667  m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_RHS[i], RHS[i]);
1668  m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1669  }
1670 
1671  SetUpCoupledMatrix(0.0, u_N, true);
1672  SolveLinearNS(RHS);
1673 
1674  for(int i = 0; i < m_velocity.size(); ++i)
1675  {
1676  u_star[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1677  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), u_star[i]);
1678 
1679  //u_star(k+1) = u_N(k) + DeltaKinvis * u_star(k)
1680  Vmath::Smul(u_star[i].size(), m_kinvis, u_star[i], 1, u_star[i], 1);
1681  Vmath::Vadd(u_star[i].size(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1682 
1683  m_fields[m_velocity[i]]->FwdTrans(u_star[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1684  }
1685 
1687  }
void SolveLinearNS(const Array< OneD, Array< OneD, NekDouble > > &forcing)
void SetUpCoupledMatrix(const NekDouble lambda=0.0, const Array< OneD, Array< OneD, NekDouble > > &Advfield=NullNekDoubleArrayOfArray, bool IsLinearNSEquation=true)
NekDouble m_kinvis
Kinematic viscosity.
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetTotPoints()
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:322
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:225

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  {
97  p->InitObject();
98  return p;
99  }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.

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

◆ DefineForcingTerm()

void Nektar::CoupledLinearNS::DefineForcingTerm ( void  )

Definition at line 1517 of file CoupledLinearNS.cpp.

1518  {
1519  m_ForcingTerm = Array<OneD, Array<OneD, NekDouble> > (m_velocity.size());
1520  m_ForcingTerm_Coeffs = Array<OneD, Array<OneD, NekDouble> > (m_velocity.size());
1521 
1522  for(int i = 0; i < m_velocity.size(); ++i)
1523  {
1524  m_ForcingTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1525  m_ForcingTerm_Coeffs[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetNcoeffs(),0.0);
1526  }
1527 
1528  if(m_session->DefinesFunction("ForcingTerm"))
1529  {
1530  std::vector<std::string> fieldStr;
1531  for(int i = 0; i < m_velocity.size(); ++i)
1532  {
1533  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1534  }
1535  GetFunction( "ForcingTerm")->Evaluate(fieldStr, m_ForcingTerm);
1536  for(int i = 0; i < m_velocity.size(); ++i)
1537  {
1538  m_fields[m_velocity[i]]->FwdTrans_IterPerExp(m_ForcingTerm[i], m_ForcingTerm_Coeffs[i]);
1539  }
1540  }
1541  else
1542  {
1543  cout << "'ForcingTerm' section has not been defined in the input file => forcing=0" << endl;
1544  }
1545  }
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm_Coeffs
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetNcoeffs()
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.

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

1344  {
1345  // evaluate convection terms
1346  EvaluateAdvectionTerms(inarray,outarray,time);
1347 
1348  for (auto &x : m_forcing)
1349  {
1350  x->Apply(m_fields, outarray, outarray, time);
1351  }
1352  }
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)

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

1725  {
1726  Array<OneD, Array<OneD, NekDouble> > Eval_Adv(m_velocity.size());
1727  Array<OneD, Array<OneD, NekDouble> > tmp_DerVel(m_velocity.size());
1728  Array<OneD, Array<OneD, NekDouble> > AdvTerm(m_velocity.size());
1729  Array<OneD, Array<OneD, NekDouble> > ViscTerm(m_velocity.size());
1730  Array<OneD, Array<OneD, NekDouble> > Forc(m_velocity.size());
1731 
1732  for(int i = 0; i < m_velocity.size(); ++i)
1733  {
1734  Eval_Adv[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1735  tmp_DerVel[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1736 
1737  AdvTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1738  ViscTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1739  Forc[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1740  outarray[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1741 
1742  m_fields[m_velocity[i]]->PhysDeriv(i, Velocity[i], tmp_DerVel[i]);
1743 
1744  Vmath::Smul(tmp_DerVel[i].size(), m_kinvis, tmp_DerVel[i], 1, tmp_DerVel[i], 1);
1745  }
1746 
1747  EvaluateAdvectionTerms(Velocity, Eval_Adv, m_time);
1748 
1749  for(int i = 0; i < m_velocity.size(); ++i)
1750  {
1751  bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1752  m_fields[m_velocity[i]]->SetWaveSpace(true);
1753  m_fields[m_velocity[i]]->IProductWRTBase(Eval_Adv[i], AdvTerm[i]); //(w, (u.grad)u)
1754  m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_DerVel[i], ViscTerm[i]); //(grad w, grad u)
1755  m_fields[m_velocity[i]]->IProductWRTBase(m_ForcingTerm[i], Forc[i]); //(w, f)
1756  m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1757 
1758  Vmath::Vsub(outarray[i].size(), outarray[i], 1, AdvTerm[i], 1, outarray[i], 1);
1759  Vmath::Vsub(outarray[i].size(), outarray[i], 1, ViscTerm[i], 1, outarray[i], 1);
1760 
1761  Vmath::Vadd(outarray[i].size(), outarray[i], 1, Forc[i], 1, outarray[i], 1);
1762  }
1763  }
NekDouble m_time
Current time of simulation.
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:372

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

1767  {
1768  int i;
1772 
1773  int nummodes;
1774 
1775  for (auto &expMapIter : VelExp)
1776  {
1778 
1779  for(i = 0; i < expMapIter.second->m_basisKeyVector.size(); ++i)
1780  {
1781  LibUtilities::BasisKey B = expMapIter.second->m_basisKeyVector[i];
1782  nummodes = B.GetNumModes();
1783  ASSERTL0(nummodes > 3,"Velocity polynomial space not sufficiently high (>= 4)");
1784  // Should probably set to be an orthogonal basis.
1785  LibUtilities::BasisKey newB(B.GetBasisType(),nummodes-2,B.GetPointsKey());
1786  BasisVec.push_back(newB);
1787  }
1788 
1789  // Put new expansion into list.
1790  SpatialDomains::ExpansionInfoShPtr expansionElementShPtr =
1792  AllocateSharedPtr(expMapIter.second->m_geomShPtr, BasisVec);
1793  (*returnval)[expMapIter.first] = expansionElementShPtr;
1794  }
1795 
1796  // Save expansion into graph.
1797  m_graph->SetExpansionInfo("p",returnval);
1798 
1799  return *returnval;
1800  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
std::vector< BasisKey > BasisKeyVector
Name for a vector of BasisKeys.
std::shared_ptr< ExpansionInfoMap > ExpansionInfoMapShPtr
Definition: MeshGraph.h:145
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 1690 of file CoupledLinearNS.cpp.

1692  {
1693  for(int i = 0; i < m_velocity.size(); ++i)
1694  {
1695  outarray[i] = 0.0;
1696  for(int j = 0; j < inarray[i].size(); ++j)
1697  {
1698  if(inarray[i][j] > outarray[i])
1699  {
1700  outarray[i] = inarray[i][j];
1701  }
1702  }
1703  cout << "InfNorm["<<i<<"] = "<< outarray[i] <<endl;
1704  }
1705  }

References Nektar::IncNavierStokes::m_velocity.

◆ L2Norm()

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

Definition at line 1707 of file CoupledLinearNS.cpp.

1709  {
1710  for(int i = 0; i < m_velocity.size(); ++i)
1711  {
1712  outarray[i] = 0.0;
1713  for(int j = 0; j < inarray[i].size(); ++j)
1714  {
1715  outarray[i] += inarray[i][j]*inarray[i][j];
1716  }
1717  outarray[i]=sqrt(outarray[i]);
1718  cout << "L2Norm["<<i<<"] = "<< outarray[i] <<endl;
1719  }
1720  }
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267

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

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

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), 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 182 of file CoupledLinearNS.cpp.

183  {
184 
185  int nz;
186  if(m_singleMode)
187  {
188 
189  NekDouble lambda_imag;
190 
191  // load imaginary component of any potential shift
192  // Probably should be called from DriverArpack but not yet
193  // clear how to do this
194  m_session->LoadParameter("imagShift",lambda_imag,NekConstants::kNekUnsetDouble);
195  nz = 1;
196  m_mat = Array<OneD, CoupledSolverMatrices> (nz);
197 
198  ASSERTL1(m_npointsZ <=2,"Only expected a maxmimum of two planes in single mode linear NS solver");
199 
200  if(m_zeroMode)
201  {
202  SetUpCoupledMatrix(lambda,Advfield,IsLinearNSEquation,0,m_mat[0],m_locToGloMap[0],lambda_imag);
203  }
204  else
205  {
206  NekDouble beta = 2*M_PI/m_LhomZ;
207  NekDouble lam = lambda + m_kinvis*beta*beta;
208 
209  SetUpCoupledMatrix(lam,Advfield,IsLinearNSEquation,1,m_mat[0],m_locToGloMap[0],lambda_imag);
210  }
211  }
212  else
213  {
214  int n;
215  if(m_npointsZ > 1)
216  {
217  nz = m_npointsZ/2;
218  }
219  else
220  {
221  nz = 1;
222  }
223 
224  m_mat = Array<OneD, CoupledSolverMatrices> (nz);
225 
226  // mean mode or 2D mode.
227  SetUpCoupledMatrix(lambda,Advfield,IsLinearNSEquation,0,m_mat[0],m_locToGloMap[0]);
228 
229  for(n = 1; n < nz; ++n)
230  {
231  NekDouble beta = 2*M_PI*n/m_LhomZ;
232 
233  NekDouble lam = lambda + m_kinvis*beta*beta;
234 
235  SetUpCoupledMatrix(lam,Advfield,IsLinearNSEquation,n,m_mat[n],m_locToGloMap[n]);
236  }
237  }
238 
239  }
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
Array< OneD, CoupledLocalToGlobalC0ContMapSharedPtr > m_locToGloMap
Array< OneD, CoupledSolverMatrices > m_mat
int m_npointsZ
number of points in Z direction (if homogeneous)

References ASSERTL1, 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 1489 of file CoupledLinearNS.cpp.

1490  {
1491  const unsigned int ncmpt = m_velocity.size();
1492  Array<OneD, Array<OneD, NekDouble> > forcing_phys(ncmpt);
1493  Array<OneD, Array<OneD, NekDouble> > forcing (ncmpt);
1494 
1495  for(int i = 0; i < ncmpt; ++i)
1496  {
1497  forcing_phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[0]]->GetNpoints(), 0.0);
1498  forcing[i] = Array<OneD, NekDouble> (m_fields[m_velocity[0]]->GetNcoeffs(),0.0);
1499  }
1500 
1501  for (auto &x : m_forcing)
1502  {
1503  const NekDouble time = 0;
1504  x->Apply(m_fields, forcing_phys, forcing_phys, time);
1505  }
1506  for (unsigned int i = 0; i < ncmpt; ++i)
1507  {
1508  bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1509  m_fields[i]->SetWaveSpace(true);
1510  m_fields[i]->IProductWRTBase(forcing_phys[i], forcing[i]);
1511  m_fields[i]->SetWaveSpace(waveSpace);
1512  }
1513 
1514  SolveLinearNS(forcing);
1515  }
SOLVER_UTILS_EXPORT int GetNpoints()

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

1915  {
1916  int i,j,k,n,cnt,cnt1;
1917  int nbnd,nint,offset;
1918  int nvel = m_velocity.size();
1919  int nel = fields[0]->GetNumElmts();
1920  Array<OneD, unsigned int> bmap, imap;
1921 
1922  Array<OneD, NekDouble > f_bnd(m_mat[mode].m_BCinv->GetRows());
1923  NekVector< NekDouble > F_bnd(f_bnd.size(), f_bnd, eWrapper);
1924  Array<OneD, NekDouble > f_int(m_mat[mode].m_BCinv->GetColumns());
1925  NekVector< NekDouble > F_int(f_int.size(),f_int, eWrapper);
1926 
1927  int nz_loc;
1928  int nplanecoeffs = fields[m_velocity[0]]->GetNcoeffs();// this is fine since we pass the nplane coeff data.
1929 
1930  if(mode) // Homogeneous mode flag
1931  {
1932  nz_loc = 2;
1933  }
1934  else
1935  {
1936  if(m_singleMode)
1937  {
1938  nz_loc = 2;
1939  }
1940  else
1941  {
1942  nz_loc = 1;
1944  {
1945  Array<OneD, NekDouble> tmp;
1946  // Zero fields to set complex mode to zero;
1947  for(i = 0; i < fields.size(); ++i)
1948  {
1949  Vmath::Zero(nplanecoeffs,tmp = fields[i]->UpdateCoeffs()+nplanecoeffs,1);
1950  }
1951  Vmath::Zero(2*pressure->GetNcoeffs(),pressure->UpdateCoeffs(),1);
1952  }
1953  }
1954  }
1955 
1956  for(k = 0; k < nvel; ++k)
1957  {
1959  std::dynamic_pointer_cast<MultiRegions::ContField>(fields[k]);
1960 
1961  Array<OneD, NekDouble> sign = cfield->GetLocalToGlobalMap()->
1962  GetBndCondCoeffsToLocalCoeffsSign();
1963  const Array<OneD, const int> map= cfield->GetLocalToGlobalMap()->
1964  GetBndCondCoeffsToLocalCoeffsMap();
1965 
1966  // Add weak boundary conditions to forcing
1967  const Array<OneD, SpatialDomains::BoundaryConditionShPtr>
1968  bndConds = fields[k]->GetBndConditions();
1969  Array<OneD, const MultiRegions::ExpListSharedPtr> bndCondExp;
1970 
1972  {
1973  bndCondExp = m_fields[k]->GetPlane(2*mode)->GetBndCondExpansions();
1974  }
1975  else
1976  {
1977  bndCondExp = m_fields[k]->GetBndCondExpansions();
1978  }
1979 
1980  for(n = 0; n < nz_loc; ++n)
1981  {
1982  int bndcnt = 0;
1983  for(i = 0; i < bndCondExp.size(); ++i)
1984  {
1985  const Array<OneD, const NekDouble > bndcoeffs =
1986  bndCondExp[i]->GetCoeffs();
1987 
1988  cnt = 0;
1989  if(bndConds[i]->GetBoundaryConditionType() ==
1991  bndConds[i]->GetBoundaryConditionType() ==
1993  {
1994  if(m_locToGloMap[mode]->GetSignChange())
1995  {
1996  for(j = 0; j < (bndCondExp[i])->GetNcoeffs(); j++)
1997  {
1998  forcing[k][n*nplanecoeffs + map[bndcnt+j]] += sign[bndcnt+j] *
1999  bndcoeffs[j];
2000  }
2001  }
2002  else
2003  {
2004  for(j = 0; j < (bndCondExp[i])->GetNcoeffs(); j++)
2005  {
2006  forcing[k][n*nplanecoeffs + map[bndcnt+j]] += bndcoeffs[j];
2007  }
2008  }
2009  }
2010 
2011  bndcnt += bndCondExp[i]->GetNcoeffs();
2012  }
2013  }
2014  }
2015 
2016  Array<OneD, NekDouble > bnd (m_locToGloMap[mode]->GetNumLocalCoeffs(),0.0);
2017 
2018  // Construct f_bnd and f_int and fill in bnd from inarray
2019  // (with Dirichlet BCs imposed)
2020  int bndoffset = 0;
2021  cnt = cnt1 = 0;
2022  for(i = 0; i < nel; ++i) // loop over elements
2023  {
2024  fields[m_velocity[0]]->GetExp(i)->GetBoundaryMap(bmap);
2025  fields[m_velocity[0]]->GetExp(i)->GetInteriorMap(imap);
2026  nbnd = bmap.size();
2027  nint = imap.size();
2028  offset = fields[m_velocity[0]]->GetCoeff_Offset(i);
2029 
2030  for(j = 0; j < nvel; ++j) // loop over velocity fields
2031  {
2032  Array<OneD, NekDouble> incoeffs = fields[j]->UpdateCoeffs();
2033 
2034  for(n = 0; n < nz_loc; ++n)
2035  {
2036  for(k = 0; k < nbnd; ++k)
2037  {
2038  f_bnd[cnt+k] = forcing[j][n*nplanecoeffs +
2039  offset+bmap[k]];
2040 
2041  bnd[bndoffset + (n + j*nz_loc)*nbnd + k] =
2042  incoeffs[n*nplanecoeffs + offset + bmap[k]];
2043  }
2044  for(k = 0; k < nint; ++k)
2045  {
2046  f_int[cnt1+k] = forcing[j][n*nplanecoeffs +
2047  offset+imap[k]];
2048  }
2049 
2050  cnt += nbnd;
2051  cnt1 += nint;
2052  }
2053  }
2054  bndoffset += nvel*nz_loc*nbnd + nz_loc*(pressure->GetExp(i)->GetNcoeffs());
2055  }
2056 
2057  Array<OneD, NekDouble > f_p(m_mat[mode].m_D_int->GetRows());
2058  NekVector< NekDouble > F_p(f_p.size(),f_p,eWrapper);
2059  NekVector< NekDouble > F_p_tmp(m_mat[mode].m_Cinv->GetRows());
2060 
2061  // fbnd does not currently hold the pressure mean
2062  F_bnd = F_bnd - (*m_mat[mode].m_BCinv)*F_int;
2063  F_p_tmp = (*m_mat[mode].m_Cinv)*F_int;
2064  F_p = (*m_mat[mode].m_D_int) * F_p_tmp;
2065 
2066  // construct inner forcing
2067  Array<OneD, NekDouble > fh_bnd(m_locToGloMap[mode]->GetNumLocalCoeffs(),0.0);
2068 
2069  offset = cnt = 0;
2070  for(i = 0; i < nel; ++i)
2071  {
2072  nbnd = nz_loc*fields[0]->GetExp(i)->NumBndryCoeffs();
2073 
2074  for(j = 0; j < nvel; ++j)
2075  {
2076  for(k = 0; k < nbnd; ++k)
2077  {
2078  fh_bnd[offset + j*nbnd + k] =
2079  f_bnd[cnt+k];
2080  }
2081  cnt += nbnd;
2082  }
2083 
2084  nint = pressure->GetExp(i)->GetNcoeffs();
2085  offset += nvel*nbnd + nint*nz_loc;
2086  }
2087 
2088  offset = cnt1 = 0;
2089  for(i = 0; i < nel; ++i)
2090  {
2091  nbnd = nz_loc*fields[0]->GetExp(i)->NumBndryCoeffs();
2092  nint = pressure->GetExp(i)->GetNcoeffs();
2093 
2094  for(n = 0; n < nz_loc; ++n)
2095  {
2096  for(j = 0; j < nint; ++j)
2097  {
2098  fh_bnd[offset + nvel*nbnd + n*nint+j] = f_p[cnt1+j];
2099  }
2100  cnt1 += nint;
2101  }
2102  offset += nvel*nbnd + nz_loc*nint;
2103  }
2104  m_mat[mode].m_CoupledBndSys->Solve(fh_bnd,bnd,m_locToGloMap[mode]);
2105 
2106  // unpack pressure and velocity boundary systems.
2107  offset = cnt = 0;
2108  int totpcoeffs = pressure->GetNcoeffs();
2109  Array<OneD, NekDouble> p_coeffs = pressure->UpdateCoeffs();
2110  for(i = 0; i < nel; ++i)
2111  {
2112  nbnd = nz_loc*fields[0]->GetExp(i)->NumBndryCoeffs();
2113  nint = pressure->GetExp(i)->GetNcoeffs();
2114  for(j = 0; j < nvel; ++j)
2115  {
2116  for(k = 0; k < nbnd; ++k)
2117  {
2118  f_bnd[cnt+k] = bnd[offset + j*nbnd + k];
2119  }
2120  cnt += nbnd;
2121  }
2122  offset += nvel*nbnd + nint*nz_loc;
2123  }
2124 
2125  pressure->SetPhysState(false);
2126 
2127  offset = cnt = cnt1 = 0;
2128  for(i = 0; i < nel; ++i)
2129  {
2130  nint = pressure->GetExp(i)->GetNcoeffs();
2131  nbnd = fields[0]->GetExp(i)->NumBndryCoeffs();
2132  cnt1 = pressure->GetCoeff_Offset(i);
2133 
2134  for(n = 0; n < nz_loc; ++n)
2135  {
2136  for(j = 0; j < nint; ++j)
2137  {
2138  p_coeffs[n*totpcoeffs + cnt1+j] =
2139  f_p[cnt+j] = bnd[offset +
2140  nvel*nz_loc*nbnd +
2141  n*nint + j];
2142  }
2143  cnt += nint;
2144  }
2145  offset += (nvel*nbnd + nint)*nz_loc;
2146  }
2147 
2148  // Back solve first level of static condensation for interior
2149  // velocity space and store in F_int
2150  F_int = F_int + Transpose(*m_mat[mode].m_D_int)*F_p
2151  - Transpose(*m_mat[mode].m_Btilde)*F_bnd;
2152  F_int = (*m_mat[mode].m_Cinv)*F_int;
2153 
2154  // Unpack solution from Bnd and F_int to v_coeffs
2155  cnt = cnt1 = 0;
2156  for(i = 0; i < nel; ++i) // loop over elements
2157  {
2158  fields[0]->GetExp(i)->GetBoundaryMap(bmap);
2159  fields[0]->GetExp(i)->GetInteriorMap(imap);
2160  nbnd = bmap.size();
2161  nint = imap.size();
2162  offset = fields[0]->GetCoeff_Offset(i);
2163 
2164  for(j = 0; j < nvel; ++j) // loop over velocity fields
2165  {
2166  for(n = 0; n < nz_loc; ++n)
2167  {
2168  for(k = 0; k < nbnd; ++k)
2169  {
2170  fields[j]->SetCoeff(n*nplanecoeffs +
2171  offset+bmap[k],
2172  f_bnd[cnt+k]);
2173  }
2174 
2175  for(k = 0; k < nint; ++k)
2176  {
2177  fields[j]->SetCoeff(n*nplanecoeffs +
2178  offset+imap[k],
2179  f_int[cnt1+k]);
2180  }
2181  cnt += nbnd;
2182  cnt1 += nint;
2183  }
2184  }
2185  }
2186 
2187  for(j = 0; j < nvel; ++j)
2188  {
2189  fields[j]->SetPhysState(false);
2190  }
2191  }
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:292

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

1867  {
1868  int i,n;
1869  Array<OneD, MultiRegions::ExpListSharedPtr> vel_fields(m_velocity.size());
1870  Array<OneD, Array<OneD, NekDouble> > force(m_velocity.size());
1871 
1872  // Impose Dirichlet conditions on velocity fields
1873  for(i = 0; i < m_velocity.size(); ++i)
1874  {
1875  Vmath::Zero(m_fields[i]->GetNcoeffs(), m_fields[i]->UpdateCoeffs(),1);
1876  m_fields[i]->ImposeDirichletConditions(m_fields[i]->UpdateCoeffs());
1877  }
1878 
1880  {
1881  int ncoeffsplane = m_fields[m_velocity[0]]->GetPlane(0)->GetNcoeffs();
1882  for(n = 0; n < m_npointsZ/2; ++n)
1883  {
1884  // Get the a Fourier mode of velocity and forcing components.
1885  for(i = 0; i < m_velocity.size(); ++i)
1886  {
1887  vel_fields[i] = m_fields[m_velocity[i]]->GetPlane(2*n);
1888  // Note this needs to correlate with how we pass forcing
1889  force[i] = forcing[i] + 2*n*ncoeffsplane;
1890  }
1891 
1892  SolveLinearNS(force,vel_fields,m_pressure->GetPlane(2*n),n);
1893  }
1894  for(i = 0; i < m_velocity.size(); ++i)
1895  {
1896  m_fields[m_velocity[i]]->SetPhysState(false);
1897  }
1898  m_pressure->SetPhysState(false);
1899  }
1900  else
1901  {
1902  for(i = 0; i < m_velocity.size(); ++i)
1903  {
1904  vel_fields[i] = m_fields[m_velocity[i]];
1905  // Note this needs to correlate with how we pass forcing
1906  force[i] = forcing[i];
1907  }
1908  SolveLinearNS(force,vel_fields,m_pressure);
1909  }
1910  }

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

1548  {
1549  LibUtilities::Timer Newtontimer;
1550  Newtontimer.Start();
1551 
1552  Array<OneD, Array<OneD, NekDouble> > RHS_Coeffs(m_velocity.size());
1553  Array<OneD, Array<OneD, NekDouble> > RHS_Phys(m_velocity.size());
1554  Array<OneD, Array<OneD, NekDouble> > delta_velocity_Phys(m_velocity.size());
1555  Array<OneD, Array<OneD, NekDouble> >Velocity_Phys(m_velocity.size());
1556  Array<OneD, NekDouble > L2_norm(m_velocity.size(), 1.0);
1557  Array<OneD, NekDouble > Inf_norm(m_velocity.size(), 1.0);
1558 
1559  for(int i = 0; i < m_velocity.size(); ++i)
1560  {
1561  delta_velocity_Phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),1.0);
1562  Velocity_Phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1563  }
1564 
1565  m_counter=1;
1566 
1567  L2Norm(delta_velocity_Phys, L2_norm);
1568 
1569  //while(max(Inf_norm[0], Inf_norm[1]) > m_tol)
1570  while(max(L2_norm[0], L2_norm[1]) > m_tol)
1571  {
1572  if(m_counter == 1)
1573  //At the first Newton step, we use the solution of the
1574  //Stokes problem (if Restart=0 in input file) Or the
1575  //solution of the .rst file (if Restart=1 in input
1576  //file)
1577  {
1578  for(int i = 0; i < m_velocity.size(); ++i)
1579  {
1580  RHS_Coeffs[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetNcoeffs(),0.0);
1581  RHS_Phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1582  }
1583 
1584  for(int i = 0; i < m_velocity.size(); ++i)
1585  {
1586  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), Velocity_Phys[i]);
1587  }
1588 
1589  m_initialStep = true;
1590  EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1591  SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1592  SolveLinearNS(RHS_Coeffs);
1593  m_initialStep = false;
1594  }
1595  if(m_counter > 1)
1596  {
1597  EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1598 
1599  if(m_counter%m_MatrixSetUpStep == 0) //Setting Up the matrix is expensive. We do it at each "m_MatrixSetUpStep" step.
1600  {
1601  SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1602  }
1603  SolveLinearNS(RHS_Coeffs);
1604  }
1605 
1606  for(int i = 0; i < m_velocity.size(); ++i)
1607  {
1608  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(RHS_Coeffs[i], RHS_Phys[i]);
1609  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), delta_velocity_Phys[i]);
1610  }
1611 
1612  for(int i = 0; i < m_velocity.size(); ++i)
1613  {
1614  Vmath::Vadd(Velocity_Phys[i].size(),Velocity_Phys[i], 1, delta_velocity_Phys[i], 1,
1615  Velocity_Phys[i], 1);
1616  }
1617 
1618  //InfNorm(delta_velocity_Phys, Inf_norm);
1619  L2Norm(delta_velocity_Phys, L2_norm);
1620 
1621  if(max(Inf_norm[0], Inf_norm[1]) > 100)
1622  {
1623  cout<<"\nThe Newton method has failed at m_kinvis = "<<m_kinvis<<" (<=> Re = " << 1/m_kinvis << ")"<< endl;
1624  ASSERTL0(0, "The Newton method has failed... \n");
1625  }
1626 
1627 
1628  cout << "\n";
1629  m_counter++;
1630  }
1631 
1632  if (m_counter > 1) //We save u:=u+\delta u in u->Coeffs
1633  {
1634  for(int i = 0; i < m_velocity.size(); ++i)
1635  {
1636  m_fields[m_velocity[i]]->FwdTrans(Velocity_Phys[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1637  }
1638  }
1639 
1640  Newtontimer.Stop();
1641  cout<<"We have done "<< m_counter-1 << " iteration(s) in " << Newtontimer.TimePerTest(1)/60 << " minute(s). \n\n";
1642  }
void L2Norm(Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
void EvaluateNewtonRHS(Array< OneD, Array< OneD, NekDouble > > &Velocity, Array< OneD, Array< OneD, NekDouble > > &outarray)

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

1358  {
1359  int i;
1360  Array<OneD, Array< OneD, NekDouble> > F(m_nConvectiveFields);
1361  NekDouble lambda = 1.0/aii_Dt;
1362  static NekDouble lambda_store;
1363  Array <OneD, Array<OneD, NekDouble> > forcing(m_velocity.size());
1364  // Matrix solution
1365  if(fabs(lambda_store - lambda) > 1e-10)
1366  {
1367  SetUpCoupledMatrix(lambda);
1368  lambda_store = lambda;
1369  }
1370 
1371  // Forcing for advection solve
1372  for(i = 0; i < m_velocity.size(); ++i)
1373  {
1374  bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1375  m_fields[m_velocity[i]]->SetWaveSpace(true);
1376  m_fields[m_velocity[i]]->IProductWRTBase(inarray[i],m_fields[m_velocity[i]]->UpdateCoeffs());
1377  m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1378  Vmath::Smul(m_fields[m_velocity[i]]->GetNcoeffs(),lambda,m_fields[m_velocity[i]]->GetCoeffs(), 1,m_fields[m_velocity[i]]->UpdateCoeffs(),1);
1379  forcing[i] = m_fields[m_velocity[i]]->GetCoeffs();
1380  }
1381 
1382  SolveLinearNS(forcing);
1383 
1384  for(i = 0; i < m_velocity.size(); ++i)
1385  {
1386  m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),outarray[i]);
1387  }
1388  }
int m_nConvectiveFields
Number of fields to be convected;.

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  )
privatevirtual

Sets up initial conditions.

Sets the initial conditions.

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 1205 of file CoupledLinearNS.cpp.

1206  {
1207  switch(m_equationType)
1208  {
1209  case eUnsteadyStokes:
1210  case eUnsteadyNavierStokes:
1211  {
1212  // LibUtilities::TimeIntegrationMethod intMethod;
1213  // std::string TimeIntStr = m_session->GetSolverInfo("TIMEINTEGRATIONMETHOD");
1214  // int i;
1215  // for(i = 0; i < (int) LibUtilities::SIZE_TimeIntegrationMethod; ++i)
1216  // {
1217  // if(boost::iequals(LibUtilities::TimeIntegrationMethodMap[i],TimeIntStr))
1218  // {
1219  // intMethod = (LibUtilities::TimeIntegrationMethod)i;
1220  // break;
1221  // }
1222  // }
1223  //
1224  // ASSERTL0(i != (int) LibUtilities::SIZE_TimeIntegrationMethod, "Invalid time integration type.");
1225  //
1226  // m_integrationScheme = LibUtilities::GetTimeIntegrationWrapperFactory().CreateInstance(LibUtilities::TimeIntegrationMethodMap[intMethod]);
1227 
1228  // Could defind this from IncNavierStokes class?
1230 
1232 
1233  // Set initial condition using time t=0
1234 
1235  SetInitialConditions(0.0);
1236 
1237  }
1238  case eSteadyStokes:
1239  SetUpCoupledMatrix(0.0);
1240  break;
1241  case eSteadyOseen:
1242  {
1243  Array<OneD, Array<OneD, NekDouble> > AdvField(m_velocity.size());
1244  for(int i = 0; i < m_velocity.size(); ++i)
1245  {
1246  AdvField[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1247  }
1248 
1249  ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1250  "Advection Velocity section must be defined in "
1251  "session file.");
1252 
1253  std::vector<std::string> fieldStr;
1254  for(int i = 0; i < m_velocity.size(); ++i)
1255  {
1256  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1257  }
1258  GetFunction("AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1259 
1260  SetUpCoupledMatrix(0.0,AdvField,false);
1261  }
1262  break;
1263  case eSteadyNavierStokes:
1264  {
1265  m_session->LoadParameter("KinvisMin", m_kinvisMin);
1266  m_session->LoadParameter("KinvisPercentage", m_KinvisPercentage);
1267  m_session->LoadParameter("Tolerence", m_tol);
1268  m_session->LoadParameter("MaxIteration", m_maxIt);
1269  m_session->LoadParameter("MatrixSetUpStep", m_MatrixSetUpStep);
1270  m_session->LoadParameter("Restart", m_Restart);
1271 
1272 
1274 
1275  if (m_Restart == 1)
1276  {
1277  ASSERTL0(m_session->DefinesFunction("Restart"),
1278  "Restart section must be defined in session file.");
1279 
1280  Array<OneD, Array<OneD, NekDouble> > Restart(m_velocity.size());
1281  for(int i = 0; i < m_velocity.size(); ++i)
1282  {
1283  Restart[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1284  }
1285  std::vector<std::string> fieldStr;
1286  for(int i = 0; i < m_velocity.size(); ++i)
1287  {
1288  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1289  }
1290  GetFunction( "Restart")->Evaluate(fieldStr, Restart);
1291 
1292  for(int i = 0; i < m_velocity.size(); ++i)
1293  {
1294  m_fields[m_velocity[i]]->FwdTrans_IterPerExp(Restart[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1295  }
1296  cout << "Saving the RESTART file for m_kinvis = "<< m_kinvis << " (<=> Re = " << 1/m_kinvis << ")" <<endl;
1297  }
1298  else //We solve the Stokes Problem
1299  {
1300 
1301  SetUpCoupledMatrix(0.0);
1302  m_initialStep = true;
1303  m_counter=1;
1304  //SolveLinearNS(m_ForcingTerm_Coeffs);
1305  Solve();
1306  m_initialStep = false;
1307  cout << "Saving the Stokes Flow for m_kinvis = "<< m_kinvis << " (<=> Re = " << 1/m_kinvis << ")" <<endl;
1308  }
1309  }
1310  break;
1311  case eSteadyLinearisedNS:
1312  {
1313  SetInitialConditions(0.0);
1314 
1315  Array<OneD, Array<OneD, NekDouble> > AdvField(m_velocity.size());
1316  for(int i = 0; i < m_velocity.size(); ++i)
1317  {
1318  AdvField[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1319  }
1320 
1321  ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1322  "Advection Velocity section must be defined in "
1323  "session file.");
1324 
1325  std::vector<std::string> fieldStr;
1326  for(int i = 0; i < m_velocity.size(); ++i)
1327  {
1328  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1329  }
1330  GetFunction("AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1331 
1332  SetUpCoupledMatrix(m_lambda,AdvField,true);
1333  }
1334  break;
1335  case eNoEquationType:
1336  default:
1337  ASSERTL0(false,"Unknown or undefined equation type for CoupledLinearNS");
1338  }
1339  }
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)
EquationType m_equationType
equation type;
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
NekDouble m_lambda
Lambda constant in real system if one required.
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
@ eSteadyNavierStokes
@ eUnsteadyStokes
@ eUnsteadyNavierStokes
@ eSteadyLinearisedNS
@ eNoEquationType

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  )
privatevirtual

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

1416  {
1417  switch(m_equationType)
1418  {
1419  case eUnsteadyStokes:
1420  case eUnsteadyNavierStokes:
1421  //AdvanceInTime(m_steps);
1423  break;
1424  case eSteadyStokes:
1425  case eSteadyOseen:
1426  case eSteadyLinearisedNS:
1427  {
1428  Solve();
1429  break;
1430  }
1431  case eSteadyNavierStokes:
1432  {
1433  LibUtilities::Timer Generaltimer;
1434  Generaltimer.Start();
1435 
1436  int Check(0);
1437 
1438  //Saving the init datas
1439  Checkpoint_Output(Check);
1440  Check++;
1441 
1442  cout<<"We execute INITIALLY SolveSteadyNavierStokes for m_kinvis = "<<m_kinvis<<" (<=> Re = "<<1/m_kinvis<<")"<<endl;
1444 
1445  while(m_kinvis > m_kinvisMin)
1446  {
1447  if (Check == 1)
1448  {
1449  cout<<"We execute SolveSteadyNavierStokes for m_kinvis = "<<m_kinvis<<" (<=> Re = "<<1/m_kinvis<<")"<<endl;
1451  Checkpoint_Output(Check);
1452  Check++;
1453  }
1454 
1455  Continuation();
1456 
1457  if (m_kinvis > m_kinvisMin)
1458  {
1459  cout<<"We execute SolveSteadyNavierStokes for m_kinvis = "<<m_kinvis<<" (<=> Re = "<<1/m_kinvis<<")"<<endl;
1461  Checkpoint_Output(Check);
1462  Check++;
1463  }
1464  }
1465 
1466 
1467  Generaltimer.Stop();
1468  cout<<"\nThe total calculation time is : " << Generaltimer.TimePerTest(1)/60 << " minute(s). \n\n";
1469 
1470  break;
1471  }
1472  case eNoEquationType:
1473  default:
1474  ASSERTL0(false,"Unknown or undefined equation type for CoupledLinearNS");
1475  }
1476  }
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
virtual SOLVER_UTILS_EXPORT void v_DoSolve()
Solves an unsteady problem.

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)
privatevirtual

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

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

References Nektar::SolverUtils::AddSummaryItem().

◆ v_GetForceDimension()

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

Implements Nektar::IncNavierStokes.

Definition at line 2231 of file CoupledLinearNS.cpp.

2232  {
2233  return m_session->GetVariables().size();
2234  }

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

◆ v_InitObject()

void Nektar::CoupledLinearNS::v_InitObject ( )
protectedvirtual

Init object for UnsteadySystem class.

Initialization object for UnsteadySystem class.

Reimplemented from Nektar::IncNavierStokes.

Definition at line 67 of file CoupledLinearNS.cpp.

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

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  )
privatevirtual

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

1485  {
1486  return true;
1487  }

◆ v_Output()

void Nektar::CoupledLinearNS::v_Output ( void  )
privatevirtual

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

2194  {
2195  std::vector<Array<OneD, NekDouble> > fieldcoeffs(m_fields.size()+1);
2196  std::vector<std::string> variables(m_fields.size()+1);
2197  int i;
2198 
2199  for(i = 0; i < m_fields.size(); ++i)
2200  {
2201  fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
2202  variables[i] = m_boundaryConditions->GetVariable(i);
2203  }
2204 
2205  fieldcoeffs[i] = Array<OneD, NekDouble>(m_fields[0]->GetNcoeffs());
2206  // project pressure field to velocity space
2207  if(m_singleMode==true)
2208  {
2209  Array<OneD, NekDouble > tmpfieldcoeffs (m_fields[0]->GetNcoeffs()/2);
2210  m_pressure->GetPlane(0)->BwdTrans_IterPerExp(m_pressure->GetPlane(0)->GetCoeffs(), m_pressure->GetPlane(0)->UpdatePhys());
2211  m_pressure->GetPlane(1)->BwdTrans_IterPerExp(m_pressure->GetPlane(1)->GetCoeffs(), m_pressure->GetPlane(1)->UpdatePhys());
2212  m_fields[0]->GetPlane(0)->FwdTrans_IterPerExp(m_pressure->GetPlane(0)->GetPhys(),fieldcoeffs[i]);
2213  m_fields[0]->GetPlane(1)->FwdTrans_IterPerExp(m_pressure->GetPlane(1)->GetPhys(),tmpfieldcoeffs);
2214  for(int e=0; e<m_fields[0]->GetNcoeffs()/2; e++)
2215  {
2216  fieldcoeffs[i][e+m_fields[0]->GetNcoeffs()/2] = tmpfieldcoeffs[e];
2217  }
2218  }
2219  else
2220  {
2221  m_pressure->BwdTrans_IterPerExp(m_pressure->GetCoeffs(),m_pressure->UpdatePhys());
2222  m_fields[0]->FwdTrans_IterPerExp(m_pressure->GetPhys(),fieldcoeffs[i]);
2223  }
2224  variables[i] = "p";
2225 
2226  std::string outname = m_sessionName + ".fld";
2227 
2228  WriteFld(outname,m_fields[0],fieldcoeffs,variables);
2229  }
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
std::string m_sessionName
Name of the session.

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

◆ v_TransCoeffToPhys()

void Nektar::CoupledLinearNS::v_TransCoeffToPhys ( void  )
privatevirtual

Virtual function for transformation to physical space.

Reimplemented from Nektar::IncNavierStokes.

Definition at line 1391 of file CoupledLinearNS.cpp.

1392  {
1393  int nfields = m_fields.size();
1394  for (int k=0 ; k < nfields; ++k)
1395  {
1396  //Backward Transformation in physical space for time evolution
1397  m_fields[k]->BwdTrans_IterPerExp(m_fields[k]->GetCoeffs(),
1398  m_fields[k]->UpdatePhys());
1399  }
1400 
1401  }

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

◆ v_TransPhysToCoeff()

void Nektar::CoupledLinearNS::v_TransPhysToCoeff ( void  )
privatevirtual

Virtual function for transformation to coefficient space.

Reimplemented from Nektar::IncNavierStokes.

Definition at line 1403 of file CoupledLinearNS.cpp.

1404  {
1405  int nfields = m_fields.size();
1406  for (int k=0 ; k < nfields; ++k)
1407  {
1408  //Forward Transformation in physical space for time evolution
1409  m_fields[k]->FwdTrans_IterPerExp(m_fields[k]->GetPhys(),
1410  m_fields[k]->UpdateCoeffs());
1411 
1412  }
1413  }

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 = SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction("CoupledLinearisedNS", CoupledLinearNS::create)
static

Name of class.

Definition at line 101 of file CoupledLinearNS.h.

◆ m_counter

int Nektar::CoupledLinearNS::m_counter
private

Definition at line 172 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

◆ m_ForcingTerm

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

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

Referenced by DefineForcingTerm().

◆ m_initialStep

bool Nektar::CoupledLinearNS::m_initialStep
private

Definition at line 173 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

◆ m_kinvisMin

NekDouble Nektar::CoupledLinearNS::m_kinvisMin
private

Definition at line 178 of file CoupledLinearNS.h.

Referenced by v_DoInitialise(), and v_DoSolve().

◆ m_KinvisPercentage

NekDouble Nektar::CoupledLinearNS::m_KinvisPercentage
private

Definition at line 180 of file CoupledLinearNS.h.

Referenced by Continuation(), and v_DoInitialise().

◆ m_kinvisStep

NekDouble Nektar::CoupledLinearNS::m_kinvisStep
private

Definition at line 179 of file CoupledLinearNS.h.

◆ m_locToGloMap

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

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

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

◆ m_maxIt

int Nektar::CoupledLinearNS::m_maxIt
private

Definition at line 175 of file CoupledLinearNS.h.

Referenced by v_DoInitialise().

◆ m_Restart

int Nektar::CoupledLinearNS::m_Restart
private

Definition at line 176 of file CoupledLinearNS.h.

Referenced by v_DoInitialise().

◆ m_tol

NekDouble Nektar::CoupledLinearNS::m_tol
private

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

Referenced by SetUpCoupledMatrix(), and v_InitObject().