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:
Inheritance graph
[legend]
Collaboration diagram for Nektar::CoupledLinearNS:
Collaboration graph
[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 ()
 
virtual void v_GetFluxVector (const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
 
virtual void v_NumericalFlux (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
 
int GetNConvectiveFields (void)
 
Array< OneD, int > & GetVelocity (void)
 
Array< OneD, NekDoubleGetElmtCFLVals (void)
 
NekDouble GetCFLEstimate (int &elmtid)
 
void AddForcing (const SolverUtils::ForcingSharedPtr &pForce)
 
- Public Member Functions inherited from Nektar::SolverUtils::AdvectionSystem
SOLVER_UTILS_EXPORT AdvectionSystem (const LibUtilities::SessionReaderSharedPtr &pSession)
 
virtual SOLVER_UTILS_EXPORT ~AdvectionSystem ()
 
AdvectionSharedPtr GetAdvObject ()
 Returns the advection object held by this instance. More...
 
- 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 >
boost::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 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 void EvaluateFunction (Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
 Evaluates a function as specified in the session file. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (std::vector< std::string > pFieldNames, Array< OneD, Array< OneD, NekDouble > > &pFields, const std::string &pName, const int domain=0)
 Populate given fields with the function from session. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (std::vector< std::string > pFieldNames, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const std::string &pName, const int domain=0)
 Populate given fields with the function from session. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, const NekDouble &pTime=0.0, const int domain=0)
 
SOLVER_UTILS_EXPORT std::string DescribeFunction (std::string pFieldName, const std::string &pFunctionName, const int domain)
 Provide a description of a function for a given field name. More...
 
SOLVER_UTILS_EXPORT void InitialiseBaseFlow (Array< OneD, Array< OneD, NekDouble > > &base)
 Perform initialisation of the base flow. 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 WeakAdvectionGreensDivergenceForm (const Array< OneD, Array< OneD, NekDouble > > &F, Array< OneD, NekDouble > &outarray)
 Compute the inner product $ (\nabla \phi \cdot F) $. More...
 
SOLVER_UTILS_EXPORT void WeakAdvectionDivergenceForm (const Array< OneD, Array< OneD, NekDouble > > &F, Array< OneD, NekDouble > &outarray)
 Compute the inner product $ (\phi, \nabla \cdot F) $. More...
 
SOLVER_UTILS_EXPORT void WeakAdvectionNonConservativeForm (const Array< OneD, Array< OneD, NekDouble > > &V, const Array< OneD, const NekDouble > &u, Array< OneD, NekDouble > &outarray, bool UseContCoeffs=false)
 Compute the inner product $ (\phi, V\cdot \nabla u) $. More...
 
f SOLVER_UTILS_EXPORT void AdvectionNonConservativeForm (const Array< OneD, Array< OneD, NekDouble > > &V, const Array< OneD, const NekDouble > &u, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
 Compute the non-conservative advection. More...
 
SOLVER_UTILS_EXPORT void WeakDGAdvection (const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false, int nvariables=0)
 Calculate the weak discontinuous Galerkin advection. More...
 
SOLVER_UTILS_EXPORT void WeakDGDiffusion (const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false)
 Calculate weak DG Diffusion in the LDG form. 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 ScanForHistoryPoints ()
 Builds map of which element holds each history point. More...
 
SOLVER_UTILS_EXPORT void WriteHistoryData (std::ostream &out)
 Probe each history point and write to 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 GetNumElmVelocity ()
 
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 SetStepsToOne ()
 
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
 
SOLVER_UTILS_EXPORT void FwdTransFields ()
 
SOLVER_UTILS_EXPORT void GetFluxVector (const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
 
SOLVER_UTILS_EXPORT void GetFluxVector (const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &fluxX, Array< OneD, Array< OneD, NekDouble > > &fluxY)
 
SOLVER_UTILS_EXPORT void GetFluxVector (const int i, const int j, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
 
SOLVER_UTILS_EXPORT void NumericalFlux (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
 
SOLVER_UTILS_EXPORT void NumericalFlux (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxX, Array< OneD, Array< OneD, NekDouble > > &numfluxY)
 
SOLVER_UTILS_EXPORT void NumFluxforScalar (const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
 
SOLVER_UTILS_EXPORT void NumFluxforVector (const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int NoCaseStringCompare (const string &s1, const string &s2)
 Perform a case-insensitive string comparison. More...
 

Static Public Member Functions

static SolverUtils::EquationSystemSharedPtr create (const LibUtilities::SessionReaderSharedPtr &pSession)
 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)
 
virtual void v_InitObject ()
 Init object for UnsteadySystem class. More...
 
- Protected Member Functions inherited from Nektar::IncNavierStokes
 IncNavierStokes (const LibUtilities::SessionReaderSharedPtr &pSession)
 Constructor. More...
 
EquationType GetEquationType (void)
 
void EvaluateAdvectionTerms (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
 
void WriteModalEnergy (void)
 
void SetBoundaryConditions (NekDouble time)
 time dependent boundary conditions updating More...
 
void SetRadiationBoundaryForcing (int fieldid)
 Set Radiation forcing term. More...
 
bool CalcSteadyState (void)
 evaluate steady state More...
 
virtual MultiRegions::ExpListSharedPtr v_GetPressure ()
 
virtual bool v_PreIntegrate (int step)
 
virtual bool v_PostIntegrate (int step)
 
- Protected Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
SOLVER_UTILS_EXPORT UnsteadySystem (const LibUtilities::SessionReaderSharedPtr &pSession)
 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 void v_NumericalFlux (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxX, Array< OneD, Array< OneD, NekDouble > > &numfluxY)
 
virtual SOLVER_UTILS_EXPORT void v_NumFluxforScalar (const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
 
virtual SOLVER_UTILS_EXPORT void v_NumFluxforVector (const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
 
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_SteadyStateCheck (int step)
 
SOLVER_UTILS_EXPORT void CheckForRestartTime (NekDouble &time)
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession)
 Initialises EquationSystem class members. More...
 
int nocase_cmp (const string &s1, const string &s2)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. 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)
 
SOLVER_UTILS_EXPORT void SetUpBaseFields (SpatialDomains::MeshGraphSharedPtr &mesh)
 
SOLVER_UTILS_EXPORT void ImportFldBase (std::string pInfile, SpatialDomains::MeshGraphSharedPtr pGraph)
 
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 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_singleMode
 Identify if a single mode is required for stability analysis. More...
 
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_subSteppingScheme
 bool to identify if using a substepping scheme More...
 
bool m_SmoothAdvection
 bool to identify if advection term smoothing is requested More...
 
LibUtilities::TimeIntegrationWrapperSharedPtr m_subStepIntegrationScheme
 
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...
 
int m_cflsteps
 dump cfl estimate More...
 
int m_steadyStateSteps
 Check for steady state at step interval. More...
 
NekDouble m_steadyStateTol
 Tolerance to which steady state should be evaluated at. More...
 
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...
 
- 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...
 
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...
 
std::vector< int > m_intVariables
 
std::vector< FilterSharedPtrm_filters
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
map< std::string, Array< OneD, Array< OneD, float > > > m_interpWeights
 Map of the interpolation weights for a specific filename. More...
 
map< std::string, Array< OneD, Array< OneD, unsigned int > > > m_interpInds
 Map of the interpolation indices for a specific filename. More...
 
Array< OneD, MultiRegions::ExpListSharedPtrm_fields
 Array holding all dependent variables. More...
 
Array< OneD, MultiRegions::ExpListSharedPtrm_base
 Base fields. More...
 
Array< OneD, MultiRegions::ExpListSharedPtrm_derivedfields
 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...
 
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_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, Array< OneD, Array< OneD, NekDouble > > > m_gradtan
 1 x nvariable x nq More...
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tanbasis
 2 x m_spacedim x nq 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...
 
int m_NumMode
 Mode to use in case of single mode analysis. More...
 

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

Constructor & Destructor Documentation

Nektar::CoupledLinearNS::CoupledLinearNS ( const LibUtilities::SessionReaderSharedPtr pSesssion)
protected

Definition at line 57 of file CoupledLinearNS.cpp.

57  :
58  UnsteadySystem(pSession),
59  IncNavierStokes(pSession),
60  m_singleMode(false),
61  m_zeroMode(false)
62  {
63  }
bool m_singleMode
Identify if a single mode is required for stability analysis.
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr &pSession)
Initialises UnsteadySystem class members.
bool m_zeroMode
Id to identify when single mode is mean mode (i.e. beta=0);.
IncNavierStokes(const LibUtilities::SessionReaderSharedPtr &pSession)
Constructor.

Member Function Documentation

void Nektar::CoupledLinearNS::Continuation ( void  )

Definition at line 1646 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().

1647  {
1648  Array<OneD, Array<OneD, NekDouble> > u_N(m_velocity.num_elements());
1649  Array<OneD, Array<OneD, NekDouble> > tmp_RHS(m_velocity.num_elements());
1650  Array<OneD, Array<OneD, NekDouble> > RHS(m_velocity.num_elements());
1651  Array<OneD, Array<OneD, NekDouble> > u_star(m_velocity.num_elements());
1652 
1653  cout << "We apply the continuation method: " <<endl;
1654 
1655  for(int i = 0; i < m_velocity.num_elements(); ++i)
1656  {
1657  u_N[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1658  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), u_N[i]);
1659 
1660  RHS[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1661  tmp_RHS[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1662 
1663  m_fields[m_velocity[i]]->PhysDeriv(i, u_N[i], tmp_RHS[i]);
1664  Vmath::Smul(tmp_RHS[i].num_elements(), m_kinvis, tmp_RHS[i], 1, tmp_RHS[i], 1);
1665 
1666  m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_RHS[i], RHS[i]);
1667  }
1668 
1669  SetUpCoupledMatrix(0.0, u_N, true);
1670  SolveLinearNS(RHS);
1671 
1672  for(int i = 0; i < m_velocity.num_elements(); ++i)
1673  {
1674  u_star[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1675  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), u_star[i]);
1676 
1677  //u_star(k+1) = u_N(k) + DeltaKinvis * u_star(k)
1678  Vmath::Smul(u_star[i].num_elements(), m_kinvis, u_star[i], 1, u_star[i], 1);
1679  Vmath::Vadd(u_star[i].num_elements(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1680 
1681  m_fields[m_velocity[i]]->FwdTrans(u_star[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1682  }
1683 
1685  }
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:199
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:285
static SolverUtils::EquationSystemSharedPtr Nektar::CoupledLinearNS::create ( const LibUtilities::SessionReaderSharedPtr pSession)
inlinestatic

Creates an instance of this class.

Definition at line 95 of file CoupledLinearNS.h.

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

97  {
99  p->InitObject();
100  return p;
101  }
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
boost::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
void Nektar::CoupledLinearNS::DefineForcingTerm ( void  )

Definition at line 1521 of file CoupledLinearNS.cpp.

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

1522  {
1525 
1526  for(int i = 0; i < m_velocity.num_elements(); ++i)
1527  {
1528  m_ForcingTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1529  m_ForcingTerm_Coeffs[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetNcoeffs(),0.0);
1530  }
1531 
1532  if(m_session->DefinesFunction("ForcingTerm"))
1533  {
1534  std::vector<std::string> fieldStr;
1535  for(int i = 0; i < m_velocity.num_elements(); ++i)
1536  {
1537  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1538  }
1539  EvaluateFunction(fieldStr, m_ForcingTerm, "ForcingTerm");
1540  for(int i = 0; i < m_velocity.num_elements(); ++i)
1541  {
1542  m_fields[m_velocity[i]]->FwdTrans_IterPerExp(m_ForcingTerm[i], m_ForcingTerm_Coeffs[i]);
1543  }
1544  }
1545  else
1546  {
1547  cout << "'ForcingTerm' section has not been defined in the input file => forcing=0" << endl;
1548  }
1549  }
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 void EvaluateFunction(Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function as specified in the session file.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
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 convectioln terms
1358  EvaluateAdvectionTerms(inarray,outarray);
1359 
1360  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
1361  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
1362  {
1363  (*x)->Apply(m_fields, outarray, outarray, time);
1364  }
1365  }
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
void Nektar::CoupledLinearNS::EvaluateNewtonRHS ( Array< OneD, Array< OneD, NekDouble > > &  Velocity,
Array< OneD, Array< OneD, NekDouble > > &  outarray 
)

Definition at line 1721 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().

1723  {
1724  Array<OneD, Array<OneD, NekDouble> > Eval_Adv(m_velocity.num_elements());
1725  Array<OneD, Array<OneD, NekDouble> > tmp_DerVel(m_velocity.num_elements());
1726  Array<OneD, Array<OneD, NekDouble> > AdvTerm(m_velocity.num_elements());
1727  Array<OneD, Array<OneD, NekDouble> > ViscTerm(m_velocity.num_elements());
1728  Array<OneD, Array<OneD, NekDouble> > Forc(m_velocity.num_elements());
1729 
1730  for(int i = 0; i < m_velocity.num_elements(); ++i)
1731  {
1732  Eval_Adv[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1733  tmp_DerVel[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1734 
1735  AdvTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1736  ViscTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1737  Forc[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1738  outarray[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1739 
1740  m_fields[m_velocity[i]]->PhysDeriv(i, Velocity[i], tmp_DerVel[i]);
1741 
1742  Vmath::Smul(tmp_DerVel[i].num_elements(), m_kinvis, tmp_DerVel[i], 1, tmp_DerVel[i], 1);
1743  }
1744 
1745  EvaluateAdvectionTerms(Velocity, Eval_Adv);
1746 
1747  for(int i = 0; i < m_velocity.num_elements(); ++i)
1748  {
1749  m_fields[m_velocity[i]]->IProductWRTBase(Eval_Adv[i], AdvTerm[i]); //(w, (u.grad)u)
1750  m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_DerVel[i], ViscTerm[i]); //(grad w, grad u)
1751  m_fields[m_velocity[i]]->IProductWRTBase(m_ForcingTerm[i], Forc[i]); //(w, f)
1752 
1753  Vmath::Vsub(outarray[i].num_elements(), outarray[i], 1, AdvTerm[i], 1, outarray[i], 1);
1754  Vmath::Vsub(outarray[i].num_elements(), outarray[i], 1, ViscTerm[i], 1, outarray[i], 1);
1755 
1756  Vmath::Vadd(outarray[i].num_elements(), outarray[i], 1, Forc[i], 1, outarray[i], 1);
1757  }
1758  }
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 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:199
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:329
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
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:285
const SpatialDomains::ExpansionMap & Nektar::CoupledLinearNS::GenPressureExp ( const SpatialDomains::ExpansionMap VelExp)

Definition at line 1762 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().

1763  {
1764  int i;
1766 
1768 
1769  SpatialDomains::ExpansionMap::const_iterator expMapIter;
1770  int nummodes;
1771 
1772  for (expMapIter = VelExp.begin(); expMapIter != VelExp.end(); ++expMapIter)
1773  {
1775 
1776  for(i = 0; i < expMapIter->second->m_basisKeyVector.size(); ++i)
1777  {
1778  LibUtilities::BasisKey B = expMapIter->second->m_basisKeyVector[i];
1779  nummodes = B.GetNumModes();
1780  ASSERTL0(nummodes > 3,"Velocity polynomial space not sufficiently high (>= 4)");
1781  // Should probably set to be an orthogonal basis.
1782  LibUtilities::BasisKey newB(B.GetBasisType(),nummodes-2,B.GetPointsKey());
1783  BasisVec.push_back(newB);
1784  }
1785 
1786  // Put new expansion into list.
1787  SpatialDomains::ExpansionShPtr expansionElementShPtr =
1788  MemoryManager<SpatialDomains::Expansion>::AllocateSharedPtr(expMapIter->second->m_geomShPtr, BasisVec);
1789  (*returnval)[expMapIter->first] = expansionElementShPtr;
1790  }
1791 
1792  // Save expansion into graph.
1793  m_graph->SetExpansions("p",returnval);
1794 
1795  return *returnval;
1796  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:135
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:139
boost::shared_ptr< ExpansionMap > ExpansionMapShPtr
Definition: MeshGraph.h:175
std::vector< BasisKey > BasisKeyVector
Name for a vector of BasisKeys.
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:145
boost::shared_ptr< Expansion > ExpansionShPtr
Definition: MeshGraph.h:170
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
Describes the specification for a Basis.
Definition: Basis.h:50
void Nektar::CoupledLinearNS::InfNorm ( Array< OneD, Array< OneD, NekDouble > > &  inarray,
Array< OneD, NekDouble > &  outarray 
)

Definition at line 1688 of file CoupledLinearNS.cpp.

References Nektar::IncNavierStokes::m_velocity.

1690  {
1691  for(int i = 0; i < m_velocity.num_elements(); ++i)
1692  {
1693  outarray[i] = 0.0;
1694  for(int j = 0; j < inarray[i].num_elements(); ++j)
1695  {
1696  if(inarray[i][j] > outarray[i])
1697  {
1698  outarray[i] = inarray[i][j];
1699  }
1700  }
1701  cout << "InfNorm["<<i<<"] = "<< outarray[i] <<endl;
1702  }
1703  }
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
void Nektar::CoupledLinearNS::L2Norm ( Array< OneD, Array< OneD, NekDouble > > &  inarray,
Array< OneD, NekDouble > &  outarray 
)

Definition at line 1705 of file CoupledLinearNS.cpp.

References Nektar::IncNavierStokes::m_velocity.

Referenced by SolveSteadyNavierStokes().

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  outarray[i] += inarray[i][j]*inarray[i][j];
1714  }
1715  outarray[i]=sqrt(outarray[i]);
1716  cout << "L2Norm["<<i<<"] = "<< outarray[i] <<endl;
1717  }
1718  }
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
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 180 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, m_singleMode, and m_zeroMode.

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

181  {
182 
183  int nz;
184  if(m_singleMode)
185  {
186 
187  NekDouble lambda_imag;
188 
189  // load imaginary componont of any potential shift
190  // Probably should be called from DriverArpack but not yet
191  // clear how to do this
192  m_session->LoadParameter("imagShift",lambda_imag,NekConstants::kNekUnsetDouble);
193 
194  nz = 1;
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 
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  }
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)
bool m_singleMode
Identify if a single mode is required for stability analysis.
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:165
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(), 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, m_singleMode, Nektar::IncNavierStokes::m_velocity, Vmath::Neg(), Nektar::NullNekDoubleArrayofArray, Vmath::Smul(), Nektar::Timer::Start(), Nektar::Timer::Stop(), Nektar::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 = m_fields[m_velocity[0]]->GetOffset_Elmt_Id(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);
438  ::AllocateSharedPtr(nsize_bndry,nsize_int,blkmatStorage);
440  ::AllocateSharedPtr(nsize_bndry,nsize_int,blkmatStorage);
442  ::AllocateSharedPtr(nsize_int,nsize_int,blkmatStorage);
443 
445  ::AllocateSharedPtr(nsize_p,nsize_bndry,blkmatStorage);
446 
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  Timer timer;
460  timer.Start();
461  for(n = 0; n < nel; ++n)
462  {
463  eid = m_fields[m_velocity[0]]->GetOffset_Elmt_Id(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;
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  {
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
1198  mat.m_CoupledBndSys->Initialise(locToGloMap);
1199  timer.Stop();
1200  cout << "Multilevel condensation: " << timer.TimePerTest(1) << endl;
1201  }
DNekScalBlkMatSharedPtr m_Btilde
Interior-boundary Laplacian plus linearised convective terms .
DNekScalBlkMatSharedPtr m_D_bnd
Inner product of pressure system with divergence of the boundary velocity space . ...
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
NekDouble m_kinvis
Kinematic viscosity.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
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:248
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
bool m_singleMode
Identify if a single mode is required for stability analysis.
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
SOLVER_UTILS_EXPORT int GetTotPoints()
void Stop()
Definition: Timer.cpp:62
DNekScalBlkMatSharedPtr m_D_int
Inner product of pressure system with divergence of the interior velocity space . ...
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:199
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
static const NekDouble kNekUnsetDouble
Describe a linear system.
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
DNekScalBlkMatSharedPtr m_Cinv
Interior-Interior Laplaican plus linearised convective terms inverted, i.e. the inverse of ...
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
void Start()
Definition: Timer.cpp:51
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:359
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition: Timer.cpp:108
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:169
DNekScalBlkMatSharedPtr m_BCinv
Boundary-interior Laplacian plus linearised convective terms pre-multiplying Cinv: ...
MultiRegions::GlobalLinSysSharedPtr m_CoupledBndSys
void Nektar::CoupledLinearNS::Solve ( void  )

Definition at line 1495 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().

1496  {
1497  const unsigned int ncmpt = m_velocity.num_elements();
1498  Array<OneD, Array<OneD, NekDouble> > forcing_phys(ncmpt);
1499  Array<OneD, Array<OneD, NekDouble> > forcing (ncmpt);
1500 
1501  for(int i = 0; i < ncmpt; ++i)
1502  {
1503  forcing_phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[0]]->GetNpoints(), 0.0);
1504  forcing[i] = Array<OneD, NekDouble> (m_fields[m_velocity[0]]->GetNcoeffs(),0.0);
1505  }
1506 
1507  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
1508  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
1509  {
1510  const NekDouble time = 0;
1511  (*x)->Apply(m_fields, forcing_phys, forcing_phys, time);
1512  }
1513  for (unsigned int i = 0; i < ncmpt; ++i)
1514  {
1515  m_fields[i]->IProductWRTBase(forcing_phys[i], forcing[i]);
1516  }
1517 
1518  SolveLinearNS(forcing);
1519  }
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.
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 1862 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().

1863  {
1864  int i,n;
1865  Array<OneD, MultiRegions::ExpListSharedPtr> vel_fields(m_velocity.num_elements());
1866  Array<OneD, Array<OneD, NekDouble> > force(m_velocity.num_elements());
1867 
1869  {
1870  int ncoeffsplane = m_fields[m_velocity[0]]->GetPlane(0)->GetNcoeffs();
1871  for(n = 0; n < m_npointsZ/2; ++n)
1872  {
1873  // Get the a Fourier mode of velocity and forcing components.
1874  for(i = 0; i < m_velocity.num_elements(); ++i)
1875  {
1876  vel_fields[i] = m_fields[m_velocity[i]]->GetPlane(2*n);
1877  // Note this needs to correlate with how we pass forcing
1878  force[i] = forcing[i] + 2*n*ncoeffsplane;
1879  }
1880 
1881  SolveLinearNS(force,vel_fields,m_pressure->GetPlane(2*n),n);
1882  }
1883  for(i = 0; i < m_velocity.num_elements(); ++i)
1884  {
1885  m_fields[m_velocity[i]]->SetPhysState(false);
1886  }
1887  m_pressure->SetPhysState(false);
1888  }
1889  else
1890  {
1891  for(i = 0; i < m_velocity.num_elements(); ++i)
1892  {
1893  vel_fields[i] = m_fields[m_velocity[i]];
1894  // Note this needs to correlate with how we pass forcing
1895  force[i] = forcing[i];
1896  }
1897  SolveLinearNS(force,vel_fields,m_pressure);
1898  }
1899  }
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
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 1901 of file CoupledLinearNS.cpp.

References Nektar::SpatialDomains::eDirichlet, Nektar::SolverUtils::EquationSystem::eHomogeneous1D, 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, m_singleMode, Nektar::IncNavierStokes::m_velocity, Nektar::Transpose(), and Vmath::Zero().

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

Definition at line 1551 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::Timer::Start(), Nektar::Timer::Stop(), Nektar::Timer::TimePerTest(), and Vmath::Vadd().

Referenced by v_DoSolve().

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

References Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::IncNavierStokes::m_nConvectiveFields, Nektar::IncNavierStokes::m_velocity, Nektar::IncNavierStokes::SetBoundaryConditions(), SetUpCoupledMatrix(), Vmath::Smul(), and SolveLinearNS().

Referenced by v_DoInitialise().

1371  {
1372  int i;
1374  NekDouble lambda = 1.0/aii_Dt;
1375  static NekDouble lambda_store;
1376  Array <OneD, Array<OneD, NekDouble> > forcing(m_velocity.num_elements());
1377  // Matrix solution
1378  if(fabs(lambda_store - lambda) > 1e-10)
1379  {
1380  cout << "Setting up Stokes matrix problem [.";
1381  fflush(stdout);
1382  SetUpCoupledMatrix(lambda);
1383  cout << "]" << endl;
1384  lambda_store = lambda;
1385  }
1386 
1387  SetBoundaryConditions(time);
1388 
1389  // Forcing for advection solve
1390  for(i = 0; i < m_velocity.num_elements(); ++i)
1391  {
1392  m_fields[m_velocity[i]]->IProductWRTBase(inarray[i],m_fields[m_velocity[i]]->UpdateCoeffs());
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  }
void SetBoundaryConditions(NekDouble time)
time dependent boundary conditions updating
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:199
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()
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::EvaluateFunction(), 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  EvaluateFunction(fieldStr,AdvField,"AdvectionVelocity");
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  EvaluateFunction(fieldStr, Restart, "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  EvaluateFunction(fieldStr,AdvField,"AdvectionVelocity");
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:135
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 void EvaluateFunction(Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function as specified in the session file.
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)
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::Timer::Start(), Nektar::Timer::Stop(), and Nektar::Timer::TimePerTest().

1431  {
1432  switch(m_equationType)
1433  {
1434  case eUnsteadyStokes:
1435  case eUnsteadyNavierStokes:
1436  //AdvanceInTime(m_steps);
1437  UnsteadySystem::v_DoSolve();
1438  break;
1439  case eSteadyStokes:
1440  case eSteadyOseen:
1441  case eSteadyLinearisedNS:
1442  {
1443  Solve();
1444  break;
1445  }
1446  case eSteadyNavierStokes:
1447  {
1448  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  }
EquationType m_equationType
equation type;
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:135
NekDouble m_kinvis
Kinematic viscosity.
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
void Stop()
Definition: Timer.cpp:62
void Start()
Definition: Timer.cpp:51
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition: Timer.cpp:108
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:50
int Nektar::CoupledLinearNS::v_GetForceDimension ( void  )
privatevirtual

Implements Nektar::IncNavierStokes.

Definition at line 2216 of file CoupledLinearNS.cpp.

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

2217  {
2218  return m_session->GetVariables().size();
2219  }
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
void Nektar::CoupledLinearNS::v_InitObject ( )
protectedvirtual

Init object for UnsteadySystem class.

Initialization object for UnsteadySystem class.

Reimplemented from Nektar::IncNavierStokes.

Definition at line 65 of file CoupledLinearNS.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, ASSERTL1, Nektar::LibUtilities::NekFactory< tKey, tBase, >::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, m_singleMode, Nektar::SolverUtils::EquationSystem::m_useFFT, Nektar::IncNavierStokes::m_velocity, m_zeroMode, and Nektar::IncNavierStokes::v_InitObject().

66  {
68 
69  int i;
70  int expdim = m_graph->GetMeshDimension();
71 
72  // Get Expansion list for orthogonal expansion at p-2
73  const SpatialDomains::ExpansionMap &pressure_exp = GenPressureExp(m_graph->GetExpansions("u"));
74 
75  m_nConvectiveFields = m_fields.num_elements();
76  if(boost::iequals(m_boundaryConditions->GetVariable(m_nConvectiveFields-1), "p"))
77  {
78  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");
79  }
80  // Decide how to declare explist for pressure.
81  if(expdim == 2)
82  {
83  int nz;
84 
86  {
87  ASSERTL0(m_fields.num_elements() > 2,"Expect to have three at least three components of velocity variables");
88  LibUtilities::BasisKey Homo1DKey = m_fields[0]->GetHomogeneousBasis()->GetBasisKey();
89 
91 
92  ASSERTL1(m_npointsZ%2==0,"Non binary number of planes have been specified");
93  nz = m_npointsZ/2;
94 
95  }
96  else
97  {
98  //m_pressure2 = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(m_session, pressure_exp);
100  nz = 1;
101  }
102 
104  for(i =0 ; i < m_velocity.num_elements(); ++i)
105  {
106  velocity[i] = m_fields[m_velocity[i]];
107  }
108 
109  // Set up Array of mappings
111 
112  if(m_session->DefinesSolverInfo("SingleMode"))
113  {
114 
115  ASSERTL0(nz <=2 ,"For single mode calculation can only have nz <= 2");
116  m_singleMode = true;
117  if(m_session->DefinesSolverInfo("BetaZero"))
118  {
119  m_zeroMode = true;
120  }
121  int nz_loc = 2;
123  }
124  else
125  {
126  // base mode
127  int nz_loc = 1;
129 
130  if(nz > 1)
131  {
132  nz_loc = 2;
133  // Assume all higher modes have the same boundary conditions and re-use mapping
135  // Note high order modes cannot be singular
136  for(i = 2; i < nz; ++i)
137  {
139  }
140  }
141  }
142  }
143  else if (expdim == 3)
144  {
145  //m_pressure = MemoryManager<MultiRegions::ExpList3D>::AllocateSharedPtr(pressure_exp);
146  ASSERTL0(false,"Setup mapping aray");
147  }
148  else
149  {
150  ASSERTL0(false,"Exp dimension not recognised");
151  }
152 
153  // creation of the extrapolation object
155  {
156  std::string vExtrapolation = "Standard";
157 
158  if (m_session->DefinesSolverInfo("Extrapolation"))
159  {
160  vExtrapolation = m_session->GetSolverInfo("Extrapolation");
161  }
162 
164  vExtrapolation,
165  m_session,
166  m_fields,
167  m_pressure,
168  m_velocity,
169  m_advObject);
170  }
171 
172  }
EquationType m_equationType
equation type;
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:135
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:48
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;.
bool m_singleMode
Identify if a single mode is required for stability analysis.
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
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:165
virtual void v_InitObject()
Init object for UnsteadySystem class.
Describes the specification for a Basis.
Definition: Basis.h:50
enum HomogeneousType m_HomogeneousType
std::map< int, ExpansionShPtr > ExpansionMap
Definition: MeshGraph.h:171
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 2178 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, m_singleMode, and Nektar::SolverUtils::EquationSystem::WriteFld().

2179  {
2180  std::vector<Array<OneD, NekDouble> > fieldcoeffs(m_fields.num_elements()+1);
2181  std::vector<std::string> variables(m_fields.num_elements()+1);
2182  int i;
2183 
2184  for(i = 0; i < m_fields.num_elements(); ++i)
2185  {
2186  fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
2187  variables[i] = m_boundaryConditions->GetVariable(i);
2188  }
2189 
2190  fieldcoeffs[i] = Array<OneD, NekDouble>(m_fields[0]->GetNcoeffs());
2191  // project pressure field to velocity space
2192  if(m_singleMode==true)
2193  {
2194  Array<OneD, NekDouble > tmpfieldcoeffs (m_fields[0]->GetNcoeffs()/2);
2195  m_pressure->GetPlane(0)->BwdTrans_IterPerExp(m_pressure->GetPlane(0)->GetCoeffs(), m_pressure->GetPlane(0)->UpdatePhys());
2196  m_pressure->GetPlane(1)->BwdTrans_IterPerExp(m_pressure->GetPlane(1)->GetCoeffs(), m_pressure->GetPlane(1)->UpdatePhys());
2197  m_fields[0]->GetPlane(0)->FwdTrans_IterPerExp(m_pressure->GetPlane(0)->GetPhys(),fieldcoeffs[i]);
2198  m_fields[0]->GetPlane(1)->FwdTrans_IterPerExp(m_pressure->GetPlane(1)->GetPhys(),tmpfieldcoeffs);
2199  for(int e=0; e<m_fields[0]->GetNcoeffs()/2; e++)
2200  {
2201  fieldcoeffs[i][e+m_fields[0]->GetNcoeffs()/2] = tmpfieldcoeffs[e];
2202  }
2203  }
2204  else
2205  {
2206  m_pressure->BwdTrans_IterPerExp(m_pressure->GetCoeffs(),m_pressure->UpdatePhys());
2207  m_fields[0]->FwdTrans_IterPerExp(m_pressure->GetPhys(),fieldcoeffs[i]);
2208  }
2209  variables[i] = "p";
2210 
2211  std::string outname = m_sessionName + ".fld";
2212 
2213  WriteFld(outname,m_fields[0],fieldcoeffs,variables);
2214  }
std::string m_sessionName
Name of the session.
bool m_singleMode
Identify if a single mode is required for stability analysis.
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.
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.
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

friend class MemoryManager< CoupledLinearNS >
friend

Definition at line 92 of file CoupledLinearNS.h.

Member Data Documentation

string Nektar::CoupledLinearNS::className = SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction("CoupledLinearisedNS", CoupledLinearNS::create)
static

Name of class.

Definition at line 103 of file CoupledLinearNS.h.

int Nektar::CoupledLinearNS::m_counter
private

Definition at line 174 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

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

Definition at line 158 of file CoupledLinearNS.h.

Referenced by DefineForcingTerm(), and EvaluateNewtonRHS().

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

Definition at line 159 of file CoupledLinearNS.h.

Referenced by DefineForcingTerm().

bool Nektar::CoupledLinearNS::m_initialStep
private

Definition at line 175 of file CoupledLinearNS.h.

Referenced by SolveLinearNS(), SolveSteadyNavierStokes(), and v_DoInitialise().

NekDouble Nektar::CoupledLinearNS::m_kinvisMin
private

Definition at line 180 of file CoupledLinearNS.h.

Referenced by v_DoInitialise(), and v_DoSolve().

NekDouble Nektar::CoupledLinearNS::m_KinvisPercentage
private

Definition at line 182 of file CoupledLinearNS.h.

Referenced by Continuation(), and v_DoInitialise().

NekDouble Nektar::CoupledLinearNS::m_kinvisStep
private

Definition at line 181 of file CoupledLinearNS.h.

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

Definition at line 161 of file CoupledLinearNS.h.

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

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

Definition at line 187 of file CoupledLinearNS.h.

Referenced by SetUpCoupledMatrix(), and SolveLinearNS().

int Nektar::CoupledLinearNS::m_MatrixSetUpStep
private

Definition at line 179 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

int Nektar::CoupledLinearNS::m_maxIt
private

Definition at line 177 of file CoupledLinearNS.h.

Referenced by v_DoInitialise().

int Nektar::CoupledLinearNS::m_Restart
private

Definition at line 178 of file CoupledLinearNS.h.

Referenced by v_DoInitialise().

bool Nektar::CoupledLinearNS::m_singleMode
private

Identify if a single mode is required for stability analysis.

Definition at line 170 of file CoupledLinearNS.h.

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

NekDouble Nektar::CoupledLinearNS::m_tol
private

Definition at line 176 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

bool Nektar::CoupledLinearNS::m_zeroMode
private

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

Definition at line 172 of file CoupledLinearNS.h.

Referenced by SetUpCoupledMatrix(), and v_InitObject().