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::ExpansionMapGenPressureExp (const SpatialDomains::ExpansionMap &VelExp)
 
void Solve (void)
 
void SolveLinearNS (const Array< OneD, Array< OneD, NekDouble > > &forcing)
 
void SolveLinearNS (const 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 (void)
 
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...
 
- 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...
 

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

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)
 
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 ()
 
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 ()
 
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...
 
- 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::TimeIntegrationWrapperSharedPtr m_intScheme
 Wrapper to the time integration scheme. More...
 
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use. More...
 
LibUtilities::TimeIntegrationSolutionSharedPtr m_intSoln
 
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...
 
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...
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
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_lambda
 Lambda constant in real system if one required. More...
 
NekDouble m_checktime
 Time between checkpoints. More...
 
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 87 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  }
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.
bool m_zeroMode
Id to identify when single mode is mean mode (i.e. beta=0);.

Member Function Documentation

◆ Continuation()

void Nektar::CoupledLinearNS::Continuation ( void  )

Definition at line 1660 of file CoupledLinearNS.cpp.

References 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().

1661  {
1662  Array<OneD, Array<OneD, NekDouble> > u_N(m_velocity.num_elements());
1663  Array<OneD, Array<OneD, NekDouble> > tmp_RHS(m_velocity.num_elements());
1664  Array<OneD, Array<OneD, NekDouble> > RHS(m_velocity.num_elements());
1665  Array<OneD, Array<OneD, NekDouble> > u_star(m_velocity.num_elements());
1666 
1667  cout << "We apply the continuation method: " <<endl;
1668 
1669  for(int i = 0; i < m_velocity.num_elements(); ++i)
1670  {
1671  u_N[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1672  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), u_N[i]);
1673 
1674  RHS[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1675  tmp_RHS[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1676 
1677  m_fields[m_velocity[i]]->PhysDeriv(i, u_N[i], tmp_RHS[i]);
1678  Vmath::Smul(tmp_RHS[i].num_elements(), m_kinvis, tmp_RHS[i], 1, tmp_RHS[i], 1);
1679 
1680  bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1681  m_fields[m_velocity[i]]->SetWaveSpace(true);
1682  m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_RHS[i], RHS[i]);
1683  m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1684  }
1685 
1686  SetUpCoupledMatrix(0.0, u_N, true);
1687  SolveLinearNS(RHS);
1688 
1689  for(int i = 0; i < m_velocity.num_elements(); ++i)
1690  {
1691  u_star[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1692  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), u_star[i]);
1693 
1694  //u_star(k+1) = u_N(k) + DeltaKinvis * u_star(k)
1695  Vmath::Smul(u_star[i].num_elements(), m_kinvis, u_star[i], 1, u_star[i], 1);
1696  Vmath::Vadd(u_star[i].num_elements(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1697 
1698  m_fields[m_velocity[i]]->FwdTrans(u_star[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1699  }
1700 
1702  }
NekDouble m_kinvis
Kinematic viscosity.
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
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)
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
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:302

◆ 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 93 of file CoupledLinearNS.h.

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

96  {
98  p->InitObject();
99  return p;
100  }
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.

◆ DefineForcingTerm()

void Nektar::CoupledLinearNS::DefineForcingTerm ( void  )

Definition at line 1532 of file CoupledLinearNS.cpp.

References Nektar::SolverUtils::EquationSystem::GetFunction(), 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().

1533  {
1534  m_ForcingTerm = Array<OneD, Array<OneD, NekDouble> > (m_velocity.num_elements());
1535  m_ForcingTerm_Coeffs = Array<OneD, Array<OneD, NekDouble> > (m_velocity.num_elements());
1536 
1537  for(int i = 0; i < m_velocity.num_elements(); ++i)
1538  {
1539  m_ForcingTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1540  m_ForcingTerm_Coeffs[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetNcoeffs(),0.0);
1541  }
1542 
1543  if(m_session->DefinesFunction("ForcingTerm"))
1544  {
1545  std::vector<std::string> fieldStr;
1546  for(int i = 0; i < m_velocity.num_elements(); ++i)
1547  {
1548  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1549  }
1550  GetFunction( "ForcingTerm")->Evaluate(fieldStr, m_ForcingTerm);
1551  for(int i = 0; i < m_velocity.num_elements(); ++i)
1552  {
1553  m_fields[m_velocity[i]]->FwdTrans_IterPerExp(m_ForcingTerm[i], m_ForcingTerm_Coeffs[i]);
1554  }
1555  }
1556  else
1557  {
1558  cout << "'ForcingTerm' section has not been defined in the input file => forcing=0" << endl;
1559  }
1560  }
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm_Coeffs
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.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.

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

References Nektar::IncNavierStokes::EvaluateAdvectionTerms(), Nektar::SolverUtils::EquationSystem::m_fields, and Nektar::IncNavierStokes::m_forcing.

Referenced by v_DoInitialise().

1356  {
1357  // evaluate convection terms
1358  EvaluateAdvectionTerms(inarray,outarray);
1359 
1360  for (auto &x : m_forcing)
1361  {
1362  x->Apply(m_fields, outarray, outarray, time);
1363  }
1364  }
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.

◆ EvaluateNewtonRHS()

void Nektar::CoupledLinearNS::EvaluateNewtonRHS ( Array< OneD, Array< OneD, NekDouble > > &  Velocity,
Array< OneD, Array< OneD, NekDouble > > &  outarray 
)

Definition at line 1738 of file CoupledLinearNS.cpp.

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

Referenced by SolveSteadyNavierStokes().

1740  {
1741  Array<OneD, Array<OneD, NekDouble> > Eval_Adv(m_velocity.num_elements());
1742  Array<OneD, Array<OneD, NekDouble> > tmp_DerVel(m_velocity.num_elements());
1743  Array<OneD, Array<OneD, NekDouble> > AdvTerm(m_velocity.num_elements());
1744  Array<OneD, Array<OneD, NekDouble> > ViscTerm(m_velocity.num_elements());
1745  Array<OneD, Array<OneD, NekDouble> > Forc(m_velocity.num_elements());
1746 
1747  for(int i = 0; i < m_velocity.num_elements(); ++i)
1748  {
1749  Eval_Adv[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1750  tmp_DerVel[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1751 
1752  AdvTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1753  ViscTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1754  Forc[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1755  outarray[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1756 
1757  m_fields[m_velocity[i]]->PhysDeriv(i, Velocity[i], tmp_DerVel[i]);
1758 
1759  Vmath::Smul(tmp_DerVel[i].num_elements(), m_kinvis, tmp_DerVel[i], 1, tmp_DerVel[i], 1);
1760  }
1761 
1762  EvaluateAdvectionTerms(Velocity, Eval_Adv);
1763 
1764  for(int i = 0; i < m_velocity.num_elements(); ++i)
1765  {
1766  bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1767  m_fields[m_velocity[i]]->SetWaveSpace(true);
1768  m_fields[m_velocity[i]]->IProductWRTBase(Eval_Adv[i], AdvTerm[i]); //(w, (u.grad)u)
1769  m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_DerVel[i], ViscTerm[i]); //(grad w, grad u)
1770  m_fields[m_velocity[i]]->IProductWRTBase(m_ForcingTerm[i], Forc[i]); //(w, f)
1771  m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1772 
1773  Vmath::Vsub(outarray[i].num_elements(), outarray[i], 1, AdvTerm[i], 1, outarray[i], 1);
1774  Vmath::Vsub(outarray[i].num_elements(), outarray[i], 1, ViscTerm[i], 1, outarray[i], 1);
1775 
1776  Vmath::Vadd(outarray[i].num_elements(), outarray[i], 1, Forc[i], 1, outarray[i], 1);
1777  }
1778  }
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, Array< OneD, NekDouble > > m_ForcingTerm
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
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:346
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
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:302

◆ GenPressureExp()

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

Definition at line 1782 of file CoupledLinearNS.cpp.

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().

1783  {
1784  int i;
1786 
1788 
1789  int nummodes;
1790 
1791  for (auto &expMapIter : VelExp)
1792  {
1794 
1795  for(i = 0; i < expMapIter.second->m_basisKeyVector.size(); ++i)
1796  {
1797  LibUtilities::BasisKey B = expMapIter.second->m_basisKeyVector[i];
1798  nummodes = B.GetNumModes();
1799  ASSERTL0(nummodes > 3,"Velocity polynomial space not sufficiently high (>= 4)");
1800  // Should probably set to be an orthogonal basis.
1801  LibUtilities::BasisKey newB(B.GetBasisType(),nummodes-2,B.GetPointsKey());
1802  BasisVec.push_back(newB);
1803  }
1804 
1805  // Put new expansion into list.
1806  SpatialDomains::ExpansionShPtr expansionElementShPtr =
1807  MemoryManager<SpatialDomains::Expansion>::AllocateSharedPtr(expMapIter.second->m_geomShPtr, BasisVec);
1808  (*returnval)[expMapIter.first] = expansionElementShPtr;
1809  }
1810 
1811  // Save expansion into graph.
1812  m_graph->SetExpansions("p",returnval);
1813 
1814  return *returnval;
1815  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
std::vector< BasisKey > BasisKeyVector
Name for a vector of BasisKeys.
std::shared_ptr< Expansion > ExpansionShPtr
Definition: MeshGraph.h:151
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
std::shared_ptr< ExpansionMap > ExpansionMapShPtr
Definition: MeshGraph.h:154

◆ InfNorm()

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

Definition at line 1705 of file CoupledLinearNS.cpp.

References Nektar::IncNavierStokes::m_velocity.

1707  {
1708  for(int i = 0; i < m_velocity.num_elements(); ++i)
1709  {
1710  outarray[i] = 0.0;
1711  for(int j = 0; j < inarray[i].num_elements(); ++j)
1712  {
1713  if(inarray[i][j] > outarray[i])
1714  {
1715  outarray[i] = inarray[i][j];
1716  }
1717  }
1718  cout << "InfNorm["<<i<<"] = "<< outarray[i] <<endl;
1719  }
1720  }
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);

◆ L2Norm()

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

Definition at line 1722 of file CoupledLinearNS.cpp.

References Nektar::IncNavierStokes::m_velocity.

Referenced by SolveSteadyNavierStokes().

1724  {
1725  for(int i = 0; i < m_velocity.num_elements(); ++i)
1726  {
1727  outarray[i] = 0.0;
1728  for(int j = 0; j < inarray[i].num_elements(); ++j)
1729  {
1730  outarray[i] += inarray[i][j]*inarray[i][j];
1731  }
1732  outarray[i]=sqrt(outarray[i]);
1733  cout << "L2Norm["<<i<<"] = "<< outarray[i] <<endl;
1734  }
1735  }
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);

◆ SetUpCoupledMatrix() [1/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 181 of file CoupledLinearNS.cpp.

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().

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

◆ SetUpCoupledMatrix() [2/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 377 of file CoupledLinearNS.cpp.

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().

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

◆ Solve()

void Nektar::CoupledLinearNS::Solve ( void  )

Definition at line 1504 of file CoupledLinearNS.cpp.

References Nektar::SolverUtils::EquationSystem::m_fields, Nektar::IncNavierStokes::m_forcing, Nektar::IncNavierStokes::m_velocity, and SolveLinearNS().

Referenced by v_DoInitialise(), and v_DoSolve().

1505  {
1506  const unsigned int ncmpt = m_velocity.num_elements();
1507  Array<OneD, Array<OneD, NekDouble> > forcing_phys(ncmpt);
1508  Array<OneD, Array<OneD, NekDouble> > forcing (ncmpt);
1509 
1510  for(int i = 0; i < ncmpt; ++i)
1511  {
1512  forcing_phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[0]]->GetNpoints(), 0.0);
1513  forcing[i] = Array<OneD, NekDouble> (m_fields[m_velocity[0]]->GetNcoeffs(),0.0);
1514  }
1515 
1516  for (auto &x : m_forcing)
1517  {
1518  const NekDouble time = 0;
1519  x->Apply(m_fields, forcing_phys, forcing_phys, time);
1520  }
1521  for (unsigned int i = 0; i < ncmpt; ++i)
1522  {
1523  bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1524  m_fields[i]->SetWaveSpace(true);
1525  m_fields[i]->IProductWRTBase(forcing_phys[i], forcing[i]);
1526  m_fields[i]->SetWaveSpace(waveSpace);
1527  }
1528 
1529  SolveLinearNS(forcing);
1530  }
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
void SolveLinearNS(const Array< OneD, Array< OneD, NekDouble > > &forcing)
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
double NekDouble
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.

◆ SolveLinearNS() [1/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 1881 of file CoupledLinearNS.cpp.

References Nektar::SolverUtils::EquationSystem::eHomogeneous1D, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_HomogeneousType, Nektar::SolverUtils::EquationSystem::m_npointsZ, Nektar::IncNavierStokes::m_pressure, and Nektar::IncNavierStokes::m_velocity.

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

1882  {
1883  int i,n;
1884  Array<OneD, MultiRegions::ExpListSharedPtr> vel_fields(m_velocity.num_elements());
1885  Array<OneD, Array<OneD, NekDouble> > force(m_velocity.num_elements());
1886 
1888  {
1889  int ncoeffsplane = m_fields[m_velocity[0]]->GetPlane(0)->GetNcoeffs();
1890  for(n = 0; n < m_npointsZ/2; ++n)
1891  {
1892  // Get the a Fourier mode of velocity and forcing components.
1893  for(i = 0; i < m_velocity.num_elements(); ++i)
1894  {
1895  vel_fields[i] = m_fields[m_velocity[i]]->GetPlane(2*n);
1896  // Note this needs to correlate with how we pass forcing
1897  force[i] = forcing[i] + 2*n*ncoeffsplane;
1898  }
1899 
1900  SolveLinearNS(force,vel_fields,m_pressure->GetPlane(2*n),n);
1901  }
1902  for(i = 0; i < m_velocity.num_elements(); ++i)
1903  {
1904  m_fields[m_velocity[i]]->SetPhysState(false);
1905  }
1906  m_pressure->SetPhysState(false);
1907  }
1908  else
1909  {
1910  for(i = 0; i < m_velocity.num_elements(); ++i)
1911  {
1912  vel_fields[i] = m_fields[m_velocity[i]];
1913  // Note this needs to correlate with how we pass forcing
1914  force[i] = forcing[i];
1915  }
1916  SolveLinearNS(force,vel_fields,m_pressure);
1917  }
1918  }
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
int m_npointsZ
number of points in Z direction (if homogeneous)
void SolveLinearNS(const Array< OneD, Array< OneD, NekDouble > > &forcing)
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
enum HomogeneousType m_HomogeneousType

◆ SolveLinearNS() [2/2]

void Nektar::CoupledLinearNS::SolveLinearNS ( const Array< OneD, Array< OneD, NekDouble > > &  forcing,
Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
MultiRegions::ExpListSharedPtr pressure,
const int  HomogeneousMode = 0 
)

Definition at line 1920 of file CoupledLinearNS.cpp.

References Nektar::SpatialDomains::eDirichlet, Nektar::SolverUtils::EquationSystem::eHomogeneous1D, Nektar::SpatialDomains::ePeriodic, Nektar::eSteadyNavierStokes, Nektar::eWrapper, Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::IncNavierStokes::m_equationType, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_HomogeneousType, m_initialStep, m_locToGloMap, m_mat, Nektar::SolverUtils::EquationSystem::m_singleMode, Nektar::IncNavierStokes::m_velocity, Nektar::Transpose(), and Vmath::Zero().

1921  {
1922  int i,j,k,n,eid,cnt,cnt1;
1923  int nbnd,nint,offset;
1924  int nvel = m_velocity.num_elements();
1925  int nel = fields[0]->GetNumElmts();
1926  Array<OneD, unsigned int> bmap, imap;
1927 
1928  Array<OneD, NekDouble > f_bnd(m_mat[mode].m_BCinv->GetRows());
1929  NekVector< NekDouble > F_bnd(f_bnd.num_elements(), f_bnd, eWrapper);
1930  Array<OneD, NekDouble > f_int(m_mat[mode].m_BCinv->GetColumns());
1931  NekVector< NekDouble > F_int(f_int.num_elements(),f_int, eWrapper);
1932 
1933  int nz_loc;
1934  int nplanecoeffs = fields[m_velocity[0]]->GetNcoeffs();// this is fine since we pass the nplane coeff data.
1935 
1936  if(mode) // Homogeneous mode flag
1937  {
1938  nz_loc = 2;
1939  }
1940  else
1941  {
1942  if(m_singleMode)
1943  {
1944  nz_loc = 2;
1945  }
1946  else
1947  {
1948  nz_loc = 1;
1950  {
1951  // Zero fields to set complex mode to zero;
1952  for(i = 0; i < fields.num_elements(); ++i)
1953  {
1954  Vmath::Zero(2*fields[i]->GetNcoeffs(),fields[i]->UpdateCoeffs(),1);
1955  }
1956  Vmath::Zero(2*pressure->GetNcoeffs(),pressure->UpdateCoeffs(),1);
1957  }
1958  }
1959  }
1960 
1961  // Assemble f_bnd and f_int
1962  cnt = cnt1 = 0;
1963  for(i = 0; i < nel; ++i) // loop over elements
1964  {
1965  eid = i;
1966  fields[m_velocity[0]]->GetExp(eid)->GetBoundaryMap(bmap);
1967  fields[m_velocity[0]]->GetExp(eid)->GetInteriorMap(imap);
1968  nbnd = bmap.num_elements();
1969  nint = imap.num_elements();
1970  offset = fields[m_velocity[0]]->GetCoeff_Offset(eid);
1971 
1972  for(j = 0; j < nvel; ++j) // loop over velocity fields
1973  {
1974  for(n = 0; n < nz_loc; ++n)
1975  {
1976  for(k = 0; k < nbnd; ++k)
1977  {
1978  f_bnd[cnt+k] = forcing[j][n*nplanecoeffs +
1979  offset+bmap[k]];
1980  }
1981  for(k = 0; k < nint; ++k)
1982  {
1983  f_int[cnt1+k] = forcing[j][n*nplanecoeffs +
1984  offset+imap[k]];
1985  }
1986  cnt += nbnd;
1987  cnt1 += nint;
1988  }
1989  }
1990  }
1991 
1992  Array<OneD, NekDouble > f_p(m_mat[mode].m_D_int->GetRows());
1993  NekVector< NekDouble > F_p(f_p.num_elements(),f_p,eWrapper);
1994  NekVector< NekDouble > F_p_tmp(m_mat[mode].m_Cinv->GetRows());
1995 
1996  // fbnd does not currently hold the pressure mean
1997  F_bnd = F_bnd - (*m_mat[mode].m_BCinv)*F_int;
1998  F_p_tmp = (*m_mat[mode].m_Cinv)*F_int;
1999  F_p = (*m_mat[mode].m_D_int) * F_p_tmp;
2000 
2001  // construct inner forcing
2002  Array<OneD, NekDouble > bnd (m_locToGloMap[mode]->GetNumGlobalCoeffs(),0.0);
2003  Array<OneD, NekDouble > fh_bnd(m_locToGloMap[mode]->GetNumGlobalCoeffs(),0.0);
2004 
2005  const Array<OneD,const int>& loctoglomap
2006  = m_locToGloMap[mode]->GetLocalToGlobalMap();
2007  const Array<OneD,const NekDouble>& loctoglosign
2008  = m_locToGloMap[mode]->GetLocalToGlobalSign();
2009 
2010  offset = cnt = 0;
2011  for(i = 0; i < nel; ++i)
2012  {
2013  eid = i;
2014  nbnd = nz_loc*fields[0]->GetExp(eid)->NumBndryCoeffs();
2015 
2016  for(j = 0; j < nvel; ++j)
2017  {
2018  for(k = 0; k < nbnd; ++k)
2019  {
2020  fh_bnd[loctoglomap[offset+j*nbnd+k]] +=
2021  loctoglosign[offset+j*nbnd+k]*f_bnd[cnt+k];
2022  }
2023  cnt += nbnd;
2024  }
2025 
2026  nint = pressure->GetExp(eid)->GetNcoeffs();
2027  offset += nvel*nbnd + nint*nz_loc;
2028  }
2029 
2030  offset = cnt1 = 0;
2031  for(i = 0; i < nel; ++i)
2032  {
2033  eid = i;
2034  nbnd = nz_loc*fields[0]->GetExp(eid)->NumBndryCoeffs();
2035  nint = pressure->GetExp(eid)->GetNcoeffs();
2036 
2037  for(n = 0; n < nz_loc; ++n)
2038  {
2039  for(j = 0; j < nint; ++j)
2040  {
2041  fh_bnd[loctoglomap[offset + nvel*nbnd + n*nint+j]] = f_p[cnt1+j];
2042  }
2043  cnt1 += nint;
2044  }
2045  offset += nvel*nbnd + nz_loc*nint;
2046  }
2047 
2048  // Set Weak BC into f_bnd and Dirichlet Dofs in bnd
2049  const Array<OneD,const int>& bndmap
2050  = m_locToGloMap[mode]->GetBndCondCoeffsToGlobalCoeffsMap();
2051 
2052  // Forcing function with weak boundary conditions and
2053  // Dirichlet conditions
2054  int bndcnt=0;
2055 
2056  for(k = 0; k < nvel; ++k)
2057  {
2058  const Array<OneD, SpatialDomains::BoundaryConditionShPtr> bndConds = fields[k]->GetBndConditions();
2059  Array<OneD, const MultiRegions::ExpListSharedPtr> bndCondExp;
2061  {
2062  bndCondExp = m_fields[k]->GetPlane(2*mode)->GetBndCondExpansions();
2063  }
2064  else
2065  {
2066  bndCondExp = m_fields[k]->GetBndCondExpansions();
2067  }
2068 
2069  for(i = 0; i < bndCondExp.num_elements(); ++i)
2070  {
2071  const Array<OneD, const NekDouble > bndCondCoeffs = bndCondExp[i]->GetCoeffs();
2072  cnt = 0;
2073  for(n = 0; n < nz_loc; ++n)
2074  {
2075  if(bndConds[i]->GetBoundaryConditionType()
2077  {
2078  for(j = 0; j < (bndCondExp[i])->GetNcoeffs(); j++)
2079  {
2081  {
2082  //This condition set all the Dirichlet BC at 0 after
2083  //the initial step of the Newton method
2084  bnd[bndmap[bndcnt++]] = 0;
2085  }
2086  else
2087  {
2088  bnd[bndmap[bndcnt++]] = bndCondCoeffs[cnt++];
2089  }
2090  }
2091  }
2092  else if (bndConds[i]->GetBoundaryConditionType()
2094  {
2095  for(j = 0; j < (bndCondExp[i])->GetNcoeffs(); j++)
2096  {
2097  fh_bnd[bndmap[bndcnt++]]
2098  += bndCondCoeffs[cnt++];
2099  }
2100  }
2101  }
2102  }
2103  }
2104 
2105  m_mat[mode].m_CoupledBndSys->Solve(fh_bnd,bnd,m_locToGloMap[mode]);
2106 
2107  // unpack pressure and velocity boundary systems.
2108  offset = cnt = 0;
2109  int totpcoeffs = pressure->GetNcoeffs();
2110  Array<OneD, NekDouble> p_coeffs = pressure->UpdateCoeffs();
2111  for(i = 0; i < nel; ++i)
2112  {
2113  eid = i;
2114  nbnd = nz_loc*fields[0]->GetExp(eid)->NumBndryCoeffs();
2115  nint = pressure->GetExp(eid)->GetNcoeffs();
2116 
2117  for(j = 0; j < nvel; ++j)
2118  {
2119  for(k = 0; k < nbnd; ++k)
2120  {
2121  f_bnd[cnt+k] = loctoglosign[offset+j*nbnd+k]*bnd[loctoglomap[offset + j*nbnd + k]];
2122  }
2123  cnt += nbnd;
2124  }
2125  offset += nvel*nbnd + nint*nz_loc;
2126  }
2127 
2128  pressure->SetPhysState(false);
2129 
2130  offset = cnt = cnt1 = 0;
2131  for(i = 0; i < nel; ++i)
2132  {
2133  eid = i;
2134  nint = pressure->GetExp(eid)->GetNcoeffs();
2135  nbnd = fields[0]->GetExp(eid)->NumBndryCoeffs();
2136  cnt1 = pressure->GetCoeff_Offset(eid);
2137 
2138  for(n = 0; n < nz_loc; ++n)
2139  {
2140  for(j = 0; j < nint; ++j)
2141  {
2142  p_coeffs[n*totpcoeffs + cnt1+j] =
2143  f_p[cnt+j] = bnd[loctoglomap[offset +
2144  (nvel*nz_loc)*nbnd +
2145  n*nint + j]];
2146  }
2147  cnt += nint;
2148  }
2149  offset += (nvel*nbnd + nint)*nz_loc;
2150  }
2151 
2152  // Back solve first level of static condensation for interior
2153  // velocity space and store in F_int
2154  F_int = F_int + Transpose(*m_mat[mode].m_D_int)*F_p
2155  - Transpose(*m_mat[mode].m_Btilde)*F_bnd;
2156  F_int = (*m_mat[mode].m_Cinv)*F_int;
2157 
2158  // Unpack solution from Bnd and F_int to v_coeffs
2159  cnt = cnt1 = 0;
2160  for(i = 0; i < nel; ++i) // loop over elements
2161  {
2162  eid = i;
2163  fields[0]->GetExp(eid)->GetBoundaryMap(bmap);
2164  fields[0]->GetExp(eid)->GetInteriorMap(imap);
2165  nbnd = bmap.num_elements();
2166  nint = imap.num_elements();
2167  offset = fields[0]->GetCoeff_Offset(eid);
2168 
2169  for(j = 0; j < nvel; ++j) // loop over velocity fields
2170  {
2171  for(n = 0; n < nz_loc; ++n)
2172  {
2173  for(k = 0; k < nbnd; ++k)
2174  {
2175  fields[j]->SetCoeff(n*nplanecoeffs +
2176  offset+bmap[k],
2177  f_bnd[cnt+k]);
2178  }
2179 
2180  for(k = 0; k < nint; ++k)
2181  {
2182  fields[j]->SetCoeff(n*nplanecoeffs +
2183  offset+imap[k],
2184  f_int[cnt1+k]);
2185  }
2186  cnt += nbnd;
2187  cnt1 += nint;
2188  }
2189  }
2190  }
2191 
2192  for(j = 0; j < nvel; ++j)
2193  {
2194  fields[j]->SetPhysState(false);
2195  }
2196  }
EquationType m_equationType
equation type;
bool m_singleMode
Flag to determine if single homogeneous mode is used.
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
Array< OneD, CoupledSolverMatrices > m_mat
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Array< OneD, CoupledLocalToGlobalC0ContMapSharedPtr > m_locToGloMap
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNcoeffs()
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:376
enum HomogeneousType m_HomogeneousType

◆ SolveSteadyNavierStokes()

void Nektar::CoupledLinearNS::SolveSteadyNavierStokes ( void  )

Definition at line 1562 of file CoupledLinearNS.cpp.

References ASSERTL0, EvaluateNewtonRHS(), 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().

1563  {
1564  LibUtilities::Timer Newtontimer;
1565  Newtontimer.Start();
1566 
1567  Array<OneD, Array<OneD, NekDouble> > RHS_Coeffs(m_velocity.num_elements());
1568  Array<OneD, Array<OneD, NekDouble> > RHS_Phys(m_velocity.num_elements());
1569  Array<OneD, Array<OneD, NekDouble> > delta_velocity_Phys(m_velocity.num_elements());
1570  Array<OneD, Array<OneD, NekDouble> >Velocity_Phys(m_velocity.num_elements());
1571  Array<OneD, NekDouble > L2_norm(m_velocity.num_elements(), 1.0);
1572  Array<OneD, NekDouble > Inf_norm(m_velocity.num_elements(), 1.0);
1573 
1574  for(int i = 0; i < m_velocity.num_elements(); ++i)
1575  {
1576  delta_velocity_Phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),1.0);
1577  Velocity_Phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1578  }
1579 
1580  m_counter=1;
1581 
1582  L2Norm(delta_velocity_Phys, L2_norm);
1583 
1584  //while(max(Inf_norm[0], Inf_norm[1]) > m_tol)
1585  while(max(L2_norm[0], L2_norm[1]) > m_tol)
1586  {
1587  if(m_counter == 1)
1588  //At the first Newton step, we use the solution of the
1589  //Stokes problem (if Restart=0 in input file) Or the
1590  //solution of the .rst file (if Restart=1 in input
1591  //file)
1592  {
1593  for(int i = 0; i < m_velocity.num_elements(); ++i)
1594  {
1595  RHS_Coeffs[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetNcoeffs(),0.0);
1596  RHS_Phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1597  }
1598 
1599  for(int i = 0; i < m_velocity.num_elements(); ++i)
1600  {
1601  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), Velocity_Phys[i]);
1602  }
1603 
1604  m_initialStep = true;
1605  EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1606  SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1607  SolveLinearNS(RHS_Coeffs);
1608  m_initialStep = false;
1609  }
1610  if(m_counter > 1)
1611  {
1612  EvaluateNewtonRHS(Velocity_Phys, RHS_Coeffs);
1613 
1614  if(m_counter%m_MatrixSetUpStep == 0) //Setting Up the matrix is expensive. We do it at each "m_MatrixSetUpStep" step.
1615  {
1616  SetUpCoupledMatrix(0.0, Velocity_Phys, true);
1617  }
1618  SolveLinearNS(RHS_Coeffs);
1619  }
1620 
1621  for(int i = 0; i < m_velocity.num_elements(); ++i)
1622  {
1623  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(RHS_Coeffs[i], RHS_Phys[i]);
1624  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), delta_velocity_Phys[i]);
1625  }
1626 
1627  for(int i = 0; i < m_velocity.num_elements(); ++i)
1628  {
1629  Vmath::Vadd(Velocity_Phys[i].num_elements(),Velocity_Phys[i], 1, delta_velocity_Phys[i], 1,
1630  Velocity_Phys[i], 1);
1631  }
1632 
1633  //InfNorm(delta_velocity_Phys, Inf_norm);
1634  L2Norm(delta_velocity_Phys, L2_norm);
1635 
1636  if(max(Inf_norm[0], Inf_norm[1]) > 100)
1637  {
1638  cout<<"\nThe Newton method has failed at m_kinvis = "<<m_kinvis<<" (<=> Re = " << 1/m_kinvis << ")"<< endl;
1639  ASSERTL0(0, "The Newton method has failed... \n");
1640  }
1641 
1642 
1643  cout << "\n";
1644  m_counter++;
1645  }
1646 
1647  if (m_counter > 1) //We save u:=u+\delta u in u->Coeffs
1648  {
1649  for(int i = 0; i < m_velocity.num_elements(); ++i)
1650  {
1651  m_fields[m_velocity[i]]->FwdTrans(Velocity_Phys[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1652  }
1653  }
1654 
1655  Newtontimer.Stop();
1656  cout<<"We have done "<< m_counter-1 << " iteration(s) in " << Newtontimer.TimePerTest(1)/60 << " minute(s). \n\n";
1657  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
NekDouble m_kinvis
Kinematic viscosity.
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
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)
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)
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
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:302

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

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().

1370  {
1371  int i;
1372  Array<OneD, Array< OneD, NekDouble> > F(m_nConvectiveFields);
1373  NekDouble lambda = 1.0/aii_Dt;
1374  static NekDouble lambda_store;
1375  Array <OneD, Array<OneD, NekDouble> > forcing(m_velocity.num_elements());
1376  // Matrix solution
1377  if(fabs(lambda_store - lambda) > 1e-10)
1378  {
1379  cout << "Setting up Stokes matrix problem [.";
1380  fflush(stdout);
1381  SetUpCoupledMatrix(lambda);
1382  cout << "]" << endl;
1383  lambda_store = lambda;
1384  }
1385 
1386  // Forcing for advection solve
1387  for(i = 0; i < m_velocity.num_elements(); ++i)
1388  {
1389  bool waveSpace = m_fields[m_velocity[i]]->GetWaveSpace();
1390  m_fields[m_velocity[i]]->SetWaveSpace(true);
1391  m_fields[m_velocity[i]]->IProductWRTBase(inarray[i],m_fields[m_velocity[i]]->UpdateCoeffs());
1392  m_fields[m_velocity[i]]->SetWaveSpace(waveSpace);
1393  Vmath::Smul(m_fields[m_velocity[i]]->GetNcoeffs(),lambda,m_fields[m_velocity[i]]->GetCoeffs(), 1,m_fields[m_velocity[i]]->UpdateCoeffs(),1);
1394  forcing[i] = m_fields[m_velocity[i]]->GetCoeffs();
1395  }
1396 
1397  SolveLinearNS(forcing);
1398 
1399  for(i = 0; i < m_velocity.num_elements(); ++i)
1400  {
1401  m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),outarray[i]);
1402  }
1403  }
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
int m_nConvectiveFields
Number of fields to be convected;.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
void SolveLinearNS(const Array< OneD, Array< OneD, NekDouble > > &forcing)
double NekDouble
void SetUpCoupledMatrix(const NekDouble lambda=0.0, const Array< OneD, Array< OneD, NekDouble > > &Advfield=NullNekDoubleArrayofArray, bool IsLinearNSEquation=true)
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNcoeffs()

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

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::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().

1209  {
1210  switch(m_equationType)
1211  {
1212  case eUnsteadyStokes:
1213  case eUnsteadyNavierStokes:
1214  {
1215 
1216 // LibUtilities::TimeIntegrationMethod intMethod;
1217 // std::string TimeIntStr = m_session->GetSolverInfo("TIMEINTEGRATIONMETHOD");
1218 // int i;
1219 // for(i = 0; i < (int) LibUtilities::SIZE_TimeIntegrationMethod; ++i)
1220 // {
1221 // if(boost::iequals(LibUtilities::TimeIntegrationMethodMap[i],TimeIntStr))
1222 // {
1223 // intMethod = (LibUtilities::TimeIntegrationMethod)i;
1224 // break;
1225 // }
1226 // }
1227 //
1228 // ASSERTL0(i != (int) LibUtilities::SIZE_TimeIntegrationMethod, "Invalid time integration type.");
1229 //
1230 // m_integrationScheme = LibUtilities::GetTimeIntegrationWrapperFactory().CreateInstance(LibUtilities::TimeIntegrationMethodMap[intMethod]);
1231 
1232  // Could defind this from IncNavierStokes class?
1234 
1236 
1237  // Set initial condition using time t=0
1238 
1239  SetInitialConditions(0.0);
1240 
1241  }
1242  case eSteadyStokes:
1243  SetUpCoupledMatrix(0.0);
1244  break;
1245  case eSteadyOseen:
1246  {
1247  Array<OneD, Array<OneD, NekDouble> > AdvField(m_velocity.num_elements());
1248  for(int i = 0; i < m_velocity.num_elements(); ++i)
1249  {
1250  AdvField[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1251  }
1252 
1253  ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1254  "Advection Velocity section must be defined in "
1255  "session file.");
1256 
1257  std::vector<std::string> fieldStr;
1258  for(int i = 0; i < m_velocity.num_elements(); ++i)
1259  {
1260  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1261  }
1262  GetFunction("AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1263 
1264  SetUpCoupledMatrix(0.0,AdvField,false);
1265  }
1266  break;
1267  case eSteadyNavierStokes:
1268  {
1269  m_session->LoadParameter("KinvisMin", m_kinvisMin);
1270  m_session->LoadParameter("KinvisPercentage", m_KinvisPercentage);
1271  m_session->LoadParameter("Tolerence", m_tol);
1272  m_session->LoadParameter("MaxIteration", m_maxIt);
1273  m_session->LoadParameter("MatrixSetUpStep", m_MatrixSetUpStep);
1274  m_session->LoadParameter("Restart", m_Restart);
1275 
1276 
1278 
1279  if (m_Restart == 1)
1280  {
1281  ASSERTL0(m_session->DefinesFunction("Restart"),
1282  "Restart section must be defined in session file.");
1283 
1284  Array<OneD, Array<OneD, NekDouble> > Restart(m_velocity.num_elements());
1285  for(int i = 0; i < m_velocity.num_elements(); ++i)
1286  {
1287  Restart[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1288  }
1289  std::vector<std::string> fieldStr;
1290  for(int i = 0; i < m_velocity.num_elements(); ++i)
1291  {
1292  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1293  }
1294  GetFunction( "Restart")->Evaluate(fieldStr, Restart);
1295 
1296  for(int i = 0; i < m_velocity.num_elements(); ++i)
1297  {
1298  m_fields[m_velocity[i]]->FwdTrans_IterPerExp(Restart[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1299  }
1300  cout << "Saving the RESTART file for m_kinvis = "<< m_kinvis << " (<=> Re = " << 1/m_kinvis << ")" <<endl;
1301  }
1302  else //We solve the Stokes Problem
1303  {
1304 
1305  /*Array<OneD, Array<OneD, NekDouble> >ZERO(m_velocity.num_elements());
1306  *
1307  * for(int i = 0; i < m_velocity.num_elements(); ++i)
1308  * {
1309  * ZERO[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1310  * m_fields[m_velocity[i]]->FwdTrans(ZERO[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1311  }*/
1312 
1313  SetUpCoupledMatrix(0.0);
1314  m_initialStep = true;
1315  m_counter=1;
1316  //SolveLinearNS(m_ForcingTerm_Coeffs);
1317  Solve();
1318  m_initialStep = false;
1319  cout << "Saving the Stokes Flow for m_kinvis = "<< m_kinvis << " (<=> Re = " << 1/m_kinvis << ")" <<endl;
1320  }
1321  }
1322  break;
1323  case eSteadyLinearisedNS:
1324  {
1325  SetInitialConditions(0.0);
1326 
1327  Array<OneD, Array<OneD, NekDouble> > AdvField(m_velocity.num_elements());
1328  for(int i = 0; i < m_velocity.num_elements(); ++i)
1329  {
1330  AdvField[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1331  }
1332 
1333  ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1334  "Advection Velocity section must be defined in "
1335  "session file.");
1336 
1337  std::vector<std::string> fieldStr;
1338  for(int i = 0; i < m_velocity.num_elements(); ++i)
1339  {
1340  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1341  }
1342  GetFunction("AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1343 
1344  SetUpCoupledMatrix(m_lambda,AdvField,true);
1345  }
1346  break;
1347  case eNoEquationType:
1348  default:
1349  ASSERTL0(false,"Unknown or undefined equation type for CoupledLinearNS");
1350  }
1351  }
EquationType m_equationType
equation type;
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
NekDouble m_kinvis
Kinematic viscosity.
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_lambda
Lambda constant in real system if one required.
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
void SetUpCoupledMatrix(const NekDouble lambda=0.0, const Array< OneD, Array< OneD, NekDouble > > &Advfield=NullNekDoubleArrayofArray, bool IsLinearNSEquation=true)
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
void EvaluateAdvection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
void SolveUnsteadyStokesSystem(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time, const NekDouble a_iixDt)

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

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().

1431  {
1432  switch(m_equationType)
1433  {
1434  case eUnsteadyStokes:
1435  case eUnsteadyNavierStokes:
1436  //AdvanceInTime(m_steps);
1438  break;
1439  case eSteadyStokes:
1440  case eSteadyOseen:
1441  case eSteadyLinearisedNS:
1442  {
1443  Solve();
1444  break;
1445  }
1446  case eSteadyNavierStokes:
1447  {
1448  LibUtilities::Timer Generaltimer;
1449  Generaltimer.Start();
1450 
1451  int Check(0);
1452 
1453  //Saving the init datas
1454  Checkpoint_Output(Check);
1455  Check++;
1456 
1457  cout<<"We execute INITIALLY SolveSteadyNavierStokes for m_kinvis = "<<m_kinvis<<" (<=> Re = "<<1/m_kinvis<<")"<<endl;
1459 
1460  while(m_kinvis > m_kinvisMin)
1461  {
1462  if (Check == 1)
1463  {
1464  cout<<"We execute SolveSteadyNavierStokes for m_kinvis = "<<m_kinvis<<" (<=> Re = "<<1/m_kinvis<<")"<<endl;
1466  Checkpoint_Output(Check);
1467  Check++;
1468  }
1469 
1470  Continuation();
1471 
1472  if (m_kinvis > m_kinvisMin)
1473  {
1474  cout<<"We execute SolveSteadyNavierStokes for m_kinvis = "<<m_kinvis<<" (<=> Re = "<<1/m_kinvis<<")"<<endl;
1476  Checkpoint_Output(Check);
1477  Check++;
1478  }
1479  }
1480 
1481 
1482  Generaltimer.Stop();
1483  cout<<"\nThe total calculation time is : " << Generaltimer.TimePerTest(1)/60 << " minute(s). \n\n";
1484 
1485  break;
1486  }
1487  case eNoEquationType:
1488  default:
1489  ASSERTL0(false,"Unknown or undefined equation type for CoupledLinearNS");
1490  }
1491  }
virtual SOLVER_UTILS_EXPORT void v_DoSolve()
Solves an unsteady problem.
EquationType m_equationType
equation type;
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
NekDouble m_kinvis
Kinematic viscosity.
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.

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

References Nektar::SolverUtils::AddSummaryItem().

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

◆ v_GetForceDimension()

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

Implements Nektar::IncNavierStokes.

Definition at line 2236 of file CoupledLinearNS.cpp.

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

2237  {
2238  return m_session->GetVariables().size();
2239  }
LibUtilities::SessionReaderSharedPtr m_session
The session reader.

◆ 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.

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().

68  {
70 
71  int i;
72  int expdim = m_graph->GetMeshDimension();
73 
74  // Get Expansion list for orthogonal expansion at p-2
75  const SpatialDomains::ExpansionMap &pressure_exp = GenPressureExp(m_graph->GetExpansions("u"));
76 
77  m_nConvectiveFields = m_fields.num_elements();
78  if(boost::iequals(m_boundaryConditions->GetVariable(m_nConvectiveFields-1), "p"))
79  {
80  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");
81  }
82  // Decide how to declare explist for pressure.
83  if(expdim == 2)
84  {
85  int nz;
86 
88  {
89  ASSERTL0(m_fields.num_elements() > 2,"Expect to have three at least three components of velocity variables");
90  LibUtilities::BasisKey Homo1DKey = m_fields[0]->GetHomogeneousBasis()->GetBasisKey();
91 
93 
94  ASSERTL1(m_npointsZ%2==0,"Non binary number of planes have been specified");
95  nz = m_npointsZ/2;
96 
97  }
98  else
99  {
100  //m_pressure2 = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(m_session, pressure_exp);
102  nz = 1;
103  }
104 
105  Array<OneD, MultiRegions::ExpListSharedPtr> velocity(m_velocity.num_elements());
106  for(i =0 ; i < m_velocity.num_elements(); ++i)
107  {
108  velocity[i] = m_fields[m_velocity[i]];
109  }
110 
111  // Set up Array of mappings
112  m_locToGloMap = Array<OneD, CoupledLocalToGlobalC0ContMapSharedPtr> (nz);
113 
114  if(m_singleMode)
115  {
116 
117  ASSERTL0(nz <=2 ,"For single mode calculation can only have nz <= 2");
118  if(m_session->DefinesSolverInfo("BetaZero"))
119  {
120  m_zeroMode = true;
121  }
122  int nz_loc = 2;
124  }
125  else
126  {
127  // base mode
128  int nz_loc = 1;
130 
131  if(nz > 1)
132  {
133  nz_loc = 2;
134  // Assume all higher modes have the same boundary conditions and re-use mapping
136  // Note high order modes cannot be singular
137  for(i = 2; i < nz; ++i)
138  {
140  }
141  }
142  }
143  }
144  else if (expdim == 3)
145  {
146  //m_pressure = MemoryManager<MultiRegions::ExpList3D>::AllocateSharedPtr(pressure_exp);
147  ASSERTL0(false,"Setup mapping aray");
148  }
149  else
150  {
151  ASSERTL0(false,"Exp dimension not recognised");
152  }
153 
154  // creation of the extrapolation object
156  {
157  std::string vExtrapolation = "Standard";
158 
159  if (m_session->DefinesSolverInfo("Extrapolation"))
160  {
161  vExtrapolation = m_session->GetSolverInfo("Extrapolation");
162  }
163 
165  vExtrapolation,
166  m_session,
167  m_fields,
168  m_pressure,
169  m_velocity,
170  m_advObject);
171  }
172 
173  }
EquationType m_equationType
equation type;
bool m_singleMode
Flag to determine if single homogeneous mode is used.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:49
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
const SpatialDomains::ExpansionMap & GenPressureExp(const SpatialDomains::ExpansionMap &VelExp)
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
ExtrapolateSharedPtr m_extrapolation
bool m_useFFT
Flag to determine if FFT is used for homogeneous transform.
int m_npointsZ
number of points in Z direction (if homogeneous)
int m_nConvectiveFields
Number of fields to be convected;.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
Array< OneD, CoupledLocalToGlobalC0ContMapSharedPtr > m_locToGloMap
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
bool m_zeroMode
Id to identify when single mode is mean mode (i.e. beta=0);.
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250
virtual void v_InitObject()
Init object for UnsteadySystem class.
enum HomogeneousType m_HomogeneousType
std::map< int, ExpansionShPtr > ExpansionMap
Definition: MeshGraph.h:152

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

1500  {
1501  return true;
1502  }

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

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().

2199  {
2200  std::vector<Array<OneD, NekDouble> > fieldcoeffs(m_fields.num_elements()+1);
2201  std::vector<std::string> variables(m_fields.num_elements()+1);
2202  int i;
2203 
2204  for(i = 0; i < m_fields.num_elements(); ++i)
2205  {
2206  fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
2207  variables[i] = m_boundaryConditions->GetVariable(i);
2208  }
2209 
2210  fieldcoeffs[i] = Array<OneD, NekDouble>(m_fields[0]->GetNcoeffs());
2211  // project pressure field to velocity space
2212  if(m_singleMode==true)
2213  {
2214  Array<OneD, NekDouble > tmpfieldcoeffs (m_fields[0]->GetNcoeffs()/2);
2215  m_pressure->GetPlane(0)->BwdTrans_IterPerExp(m_pressure->GetPlane(0)->GetCoeffs(), m_pressure->GetPlane(0)->UpdatePhys());
2216  m_pressure->GetPlane(1)->BwdTrans_IterPerExp(m_pressure->GetPlane(1)->GetCoeffs(), m_pressure->GetPlane(1)->UpdatePhys());
2217  m_fields[0]->GetPlane(0)->FwdTrans_IterPerExp(m_pressure->GetPlane(0)->GetPhys(),fieldcoeffs[i]);
2218  m_fields[0]->GetPlane(1)->FwdTrans_IterPerExp(m_pressure->GetPlane(1)->GetPhys(),tmpfieldcoeffs);
2219  for(int e=0; e<m_fields[0]->GetNcoeffs()/2; e++)
2220  {
2221  fieldcoeffs[i][e+m_fields[0]->GetNcoeffs()/2] = tmpfieldcoeffs[e];
2222  }
2223  }
2224  else
2225  {
2226  m_pressure->BwdTrans_IterPerExp(m_pressure->GetCoeffs(),m_pressure->UpdatePhys());
2227  m_fields[0]->FwdTrans_IterPerExp(m_pressure->GetPhys(),fieldcoeffs[i]);
2228  }
2229  variables[i] = "p";
2230 
2231  std::string outname = m_sessionName + ".fld";
2232 
2233  WriteFld(outname,m_fields[0],fieldcoeffs,variables);
2234  }
bool m_singleMode
Flag to determine if single homogeneous mode is used.
std::string m_sessionName
Name of the session.
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNcoeffs()
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.

◆ v_TransCoeffToPhys()

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

Virtual function for transformation to physical space.

Reimplemented from Nektar::IncNavierStokes.

Definition at line 1406 of file CoupledLinearNS.cpp.

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

1407  {
1408  int nfields = m_fields.num_elements();
1409  for (int k=0 ; k < nfields; ++k)
1410  {
1411  //Backward Transformation in physical space for time evolution
1412  m_fields[k]->BwdTrans_IterPerExp(m_fields[k]->GetCoeffs(),
1413  m_fields[k]->UpdatePhys());
1414  }
1415 
1416  }
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.

◆ v_TransPhysToCoeff()

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

Virtual function for transformation to coefficient space.

Reimplemented from Nektar::IncNavierStokes.

Definition at line 1418 of file CoupledLinearNS.cpp.

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

1419  {
1420  int nfields = m_fields.num_elements();
1421  for (int k=0 ; k < nfields; ++k)
1422  {
1423  //Forward Transformation in physical space for time evolution
1424  m_fields[k]->FwdTrans_IterPerExp(m_fields[k]->GetPhys(),
1425  m_fields[k]->UpdateCoeffs());
1426 
1427  }
1428  }
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.

Friends And Related Function Documentation

◆ MemoryManager< CoupledLinearNS >

friend class MemoryManager< CoupledLinearNS >
friend

Definition at line 90 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 102 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 SolveLinearNS(), 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().