Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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::ExpansionMap
GenPressureExp (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 NekDouble &pTime=0.0, 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 NekDouble &pTime=0.0, 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, NekDouble
ErrorExtraPoints (unsigned int field)
 Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf]. More...
 
SOLVER_UTILS_EXPORT void 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::FieldMetaDataMap
UpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated. More...
 
SOLVER_UTILS_EXPORT NekDouble GetFinalTime ()
 Return final time. More...
 
SOLVER_UTILS_EXPORT int GetNcoeffs ()
 
SOLVER_UTILS_EXPORT int GetNcoeffs (const int eid)
 
SOLVER_UTILS_EXPORT int GetNumExpModes ()
 
SOLVER_UTILS_EXPORT const
Array< OneD, int > 
GetNumExpModesPerExp ()
 
SOLVER_UTILS_EXPORT int GetNvariables ()
 
SOLVER_UTILS_EXPORT const
std::string 
GetVariable (unsigned int i)
 
SOLVER_UTILS_EXPORT int GetTraceTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTraceNpoints ()
 
SOLVER_UTILS_EXPORT int GetExpSize ()
 
SOLVER_UTILS_EXPORT int GetPhys_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetCoeff_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTotPoints (int n)
 
SOLVER_UTILS_EXPORT int GetNpoints ()
 
SOLVER_UTILS_EXPORT int 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 SetSteps (const int steps)
 
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 std::string &s1, const std::string &s2)
 Perform a case-insensitive string comparison. More...
 
SOLVER_UTILS_EXPORT int GetCheckpointNumber ()
 
SOLVER_UTILS_EXPORT void SetCheckpointNumber (int num)
 
SOLVER_UTILS_EXPORT int GetCheckpointSteps ()
 
SOLVER_UTILS_EXPORT void SetCheckpointSteps (int num)
 
SOLVER_UTILS_EXPORT void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 

Static Public Member Functions

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

Public Attributes

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

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...
 
void SetZeroNormalVelocity ()
 Set Normal Velocity Component to Zero. 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)
 
SOLVER_UTILS_EXPORT void SVVVarDiffCoeff (const Array< OneD, Array< OneD, NekDouble > > vel, StdRegions::VarCoeffMap &varCoeffMap)
 Evaluate the SVV diffusion coefficient according to Moura's paper where it should proportional to h time velocity. More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession)
 Initialises EquationSystem class members. More...
 
int nocase_cmp (const std::string &s1, const std::string &s2)
 
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 bool v_NegatedOp (void)
 
virtual void v_TransCoeffToPhys (void)
 Virtual function for transformation to physical space. More...
 
virtual void v_TransPhysToCoeff (void)
 Virtual function for transformation to coefficient space. More...
 
virtual void v_Output (void)
 
virtual int v_GetForceDimension (void)
 

Private Attributes

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

Friends

class MemoryManager< CoupledLinearNS >
 

Additional Inherited Members

- Protected Types inherited from Nektar::SolverUtils::EquationSystem
enum  HomogeneousType { eHomogeneous1D, eHomogeneous2D, eHomogeneous3D, eNotHomogeneous }
 Parameter for homogeneous expansions. More...
 
- Protected Attributes inherited from Nektar::IncNavierStokes
ExtrapolateSharedPtr m_extrapolation
 
std::ofstream m_mdlFile
 modal energy file More...
 
bool m_SmoothAdvection
 bool to identify if advection term smoothing is requested More...
 
std::vector
< SolverUtils::ForcingSharedPtr
m_forcing
 Forcing terms. More...
 
int m_nConvectiveFields
 Number of fields to be convected;. More...
 
Array< OneD, int > m_velocity
 int which identifies which components of m_fields contains the velocity (u,v,w); More...
 
MultiRegions::ExpListSharedPtr m_pressure
 Pointer to field holding pressure field. More...
 
NekDouble m_kinvis
 Kinematic viscosity. More...
 
int m_energysteps
 dump energy to file at steps time More...
 
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...
 
std::map< std::string, Array
< OneD, Array< OneD, float > > > 
m_interpWeights
 Map of the interpolation weights for a specific filename. More...
 
std::map< std::string, Array
< OneD, Array< OneD, unsigned
int > > > 
m_interpInds
 Map of the interpolation indices for a specific filename. More...
 
Array< OneD,
MultiRegions::ExpListSharedPtr
m_fields
 Array holding all dependent variables. More...
 
Array< OneD,
MultiRegions::ExpListSharedPtr
m_base
 Base fields. More...
 
Array< OneD,
MultiRegions::ExpListSharedPtr
m_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...
 
int m_initialStep
 Number of the step where the simulation should begin. More...
 
NekDouble m_fintime
 Finish time of the simulation. More...
 
NekDouble m_timestep
 Time step size. More...
 
NekDouble m_lambda
 Lambda constant in real system if one required. More...
 
std::set< std::string > m_loadedFields
 
NekDouble m_checktime
 Time between checkpoints. More...
 
int m_nchk
 Number of checkpoints written so far. More...
 
int m_steps
 Number of steps to take. More...
 
int m_checksteps
 Number of steps between checkpoints. More...
 
int m_spacedim
 Spatial dimension (>= expansion dim). More...
 
int m_expdim
 Expansion dimension. More...
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform. More...
 
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations. More...
 
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More...
 
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous. More...
 
Array< OneD, Array< OneD,
NekDouble > > 
m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction. More...
 
Array< OneD, 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...
 

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

59  :
60  UnsteadySystem(pSession),
61  IncNavierStokes(pSession),
62  m_zeroMode(false)
63  {
64  }
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 1657 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().

1658  {
1659  Array<OneD, Array<OneD, NekDouble> > u_N(m_velocity.num_elements());
1660  Array<OneD, Array<OneD, NekDouble> > tmp_RHS(m_velocity.num_elements());
1661  Array<OneD, Array<OneD, NekDouble> > RHS(m_velocity.num_elements());
1662  Array<OneD, Array<OneD, NekDouble> > u_star(m_velocity.num_elements());
1663 
1664  cout << "We apply the continuation method: " <<endl;
1665 
1666  for(int i = 0; i < m_velocity.num_elements(); ++i)
1667  {
1668  u_N[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1669  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), u_N[i]);
1670 
1671  RHS[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1672  tmp_RHS[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1673 
1674  m_fields[m_velocity[i]]->PhysDeriv(i, u_N[i], tmp_RHS[i]);
1675  Vmath::Smul(tmp_RHS[i].num_elements(), m_kinvis, tmp_RHS[i], 1, tmp_RHS[i], 1);
1676 
1677  m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_RHS[i], RHS[i]);
1678  }
1679 
1680  SetUpCoupledMatrix(0.0, u_N, true);
1681  SolveLinearNS(RHS);
1682 
1683  for(int i = 0; i < m_velocity.num_elements(); ++i)
1684  {
1685  u_star[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1686  m_fields[m_velocity[i]]->BwdTrans_IterPerExp(m_fields[m_velocity[i]]->GetCoeffs(), u_star[i]);
1687 
1688  //u_star(k+1) = u_N(k) + DeltaKinvis * u_star(k)
1689  Vmath::Smul(u_star[i].num_elements(), m_kinvis, u_star[i], 1, u_star[i], 1);
1690  Vmath::Vadd(u_star[i].num_elements(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1691 
1692  m_fields[m_velocity[i]]->FwdTrans(u_star[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1693  }
1694 
1696  }
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 1529 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().

1530  {
1531  m_ForcingTerm = Array<OneD, Array<OneD, NekDouble> > (m_velocity.num_elements());
1532  m_ForcingTerm_Coeffs = Array<OneD, Array<OneD, NekDouble> > (m_velocity.num_elements());
1533 
1534  for(int i = 0; i < m_velocity.num_elements(); ++i)
1535  {
1536  m_ForcingTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1537  m_ForcingTerm_Coeffs[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetNcoeffs(),0.0);
1538  }
1539 
1540  if(m_session->DefinesFunction("ForcingTerm"))
1541  {
1542  std::vector<std::string> fieldStr;
1543  for(int i = 0; i < m_velocity.num_elements(); ++i)
1544  {
1545  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1546  }
1547  EvaluateFunction(fieldStr, m_ForcingTerm, "ForcingTerm");
1548  for(int i = 0; i < m_velocity.num_elements(); ++i)
1549  {
1550  m_fields[m_velocity[i]]->FwdTrans_IterPerExp(m_ForcingTerm[i], m_ForcingTerm_Coeffs[i]);
1551  }
1552  }
1553  else
1554  {
1555  cout << "'ForcingTerm' section has not been defined in the input file => forcing=0" << endl;
1556  }
1557  }
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 1352 of file CoupledLinearNS.cpp.

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

Referenced by v_DoInitialise().

1355  {
1356  // evaluate convection terms
1357  EvaluateAdvectionTerms(inarray,outarray);
1358 
1359  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
1360  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
1361  {
1362  (*x)->Apply(m_fields, outarray, outarray, time);
1363  }
1364  }
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 1732 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().

1734  {
1735  Array<OneD, Array<OneD, NekDouble> > Eval_Adv(m_velocity.num_elements());
1736  Array<OneD, Array<OneD, NekDouble> > tmp_DerVel(m_velocity.num_elements());
1737  Array<OneD, Array<OneD, NekDouble> > AdvTerm(m_velocity.num_elements());
1738  Array<OneD, Array<OneD, NekDouble> > ViscTerm(m_velocity.num_elements());
1739  Array<OneD, Array<OneD, NekDouble> > Forc(m_velocity.num_elements());
1740 
1741  for(int i = 0; i < m_velocity.num_elements(); ++i)
1742  {
1743  Eval_Adv[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1744  tmp_DerVel[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1745 
1746  AdvTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1747  ViscTerm[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1748  Forc[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1749  outarray[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1750 
1751  m_fields[m_velocity[i]]->PhysDeriv(i, Velocity[i], tmp_DerVel[i]);
1752 
1753  Vmath::Smul(tmp_DerVel[i].num_elements(), m_kinvis, tmp_DerVel[i], 1, tmp_DerVel[i], 1);
1754  }
1755 
1756  EvaluateAdvectionTerms(Velocity, Eval_Adv);
1757 
1758  for(int i = 0; i < m_velocity.num_elements(); ++i)
1759  {
1760  m_fields[m_velocity[i]]->IProductWRTBase(Eval_Adv[i], AdvTerm[i]); //(w, (u.grad)u)
1761  m_fields[m_velocity[i]]->IProductWRTDerivBase(i, tmp_DerVel[i], ViscTerm[i]); //(grad w, grad u)
1762  m_fields[m_velocity[i]]->IProductWRTBase(m_ForcingTerm[i], Forc[i]); //(w, f)
1763 
1764  Vmath::Vsub(outarray[i].num_elements(), outarray[i], 1, AdvTerm[i], 1, outarray[i], 1);
1765  Vmath::Vsub(outarray[i].num_elements(), outarray[i], 1, ViscTerm[i], 1, outarray[i], 1);
1766 
1767  Vmath::Vadd(outarray[i].num_elements(), outarray[i], 1, Forc[i], 1, outarray[i], 1);
1768  }
1769  }
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 1773 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().

1774  {
1775  int i;
1777 
1779 
1780  SpatialDomains::ExpansionMap::const_iterator expMapIter;
1781  int nummodes;
1782 
1783  for (expMapIter = VelExp.begin(); expMapIter != VelExp.end(); ++expMapIter)
1784  {
1786 
1787  for(i = 0; i < expMapIter->second->m_basisKeyVector.size(); ++i)
1788  {
1789  LibUtilities::BasisKey B = expMapIter->second->m_basisKeyVector[i];
1790  nummodes = B.GetNumModes();
1791  ASSERTL0(nummodes > 3,"Velocity polynomial space not sufficiently high (>= 4)");
1792  // Should probably set to be an orthogonal basis.
1793  LibUtilities::BasisKey newB(B.GetBasisType(),nummodes-2,B.GetPointsKey());
1794  BasisVec.push_back(newB);
1795  }
1796 
1797  // Put new expansion into list.
1798  SpatialDomains::ExpansionShPtr expansionElementShPtr =
1799  MemoryManager<SpatialDomains::Expansion>::AllocateSharedPtr(expMapIter->second->m_geomShPtr, BasisVec);
1800  (*returnval)[expMapIter->first] = expansionElementShPtr;
1801  }
1802 
1803  // Save expansion into graph.
1804  m_graph->SetExpansions("p",returnval);
1805 
1806  return *returnval;
1807  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
boost::shared_ptr< ExpansionMap > ExpansionMapShPtr
Definition: MeshGraph.h:178
std::vector< BasisKey > BasisKeyVector
Name for a vector of BasisKeys.
boost::shared_ptr< Expansion > ExpansionShPtr
Definition: MeshGraph.h:173
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
void Nektar::CoupledLinearNS::InfNorm ( Array< OneD, Array< OneD, NekDouble > > &  inarray,
Array< OneD, NekDouble > &  outarray 
)

Definition at line 1699 of file CoupledLinearNS.cpp.

References Nektar::IncNavierStokes::m_velocity.

1701  {
1702  for(int i = 0; i < m_velocity.num_elements(); ++i)
1703  {
1704  outarray[i] = 0.0;
1705  for(int j = 0; j < inarray[i].num_elements(); ++j)
1706  {
1707  if(inarray[i][j] > outarray[i])
1708  {
1709  outarray[i] = inarray[i][j];
1710  }
1711  }
1712  cout << "InfNorm["<<i<<"] = "<< outarray[i] <<endl;
1713  }
1714  }
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 1716 of file CoupledLinearNS.cpp.

References Nektar::IncNavierStokes::m_velocity.

Referenced by SolveSteadyNavierStokes().

1718  {
1719  for(int i = 0; i < m_velocity.num_elements(); ++i)
1720  {
1721  outarray[i] = 0.0;
1722  for(int j = 0; j < inarray[i].num_elements(); ++j)
1723  {
1724  outarray[i] += inarray[i][j]*inarray[i][j];
1725  }
1726  outarray[i]=sqrt(outarray[i]);
1727  cout << "L2Norm["<<i<<"] = "<< outarray[i] <<endl;
1728  }
1729  }
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, Nektar::SolverUtils::EquationSystem::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 component 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  nz = 1;
194  m_mat = Array<OneD, CoupledSolverMatrices> (nz);
195 
196  ASSERTL1(m_npointsZ <=2,"Only expected a maxmimum of two planes in single mode linear NS solver");
197 
198  if(m_zeroMode)
199  {
200  SetUpCoupledMatrix(lambda,Advfield,IsLinearNSEquation,0,m_mat[0],m_locToGloMap[0],lambda_imag);
201  }
202  else
203  {
204  NekDouble beta = 2*M_PI/m_LhomZ;
205  NekDouble lam = lambda + m_kinvis*beta*beta;
206 
207  SetUpCoupledMatrix(lam,Advfield,IsLinearNSEquation,1,m_mat[0],m_locToGloMap[0],lambda_imag);
208  }
209  }
210  else
211  {
212  int n;
213  if(m_npointsZ > 1)
214  {
215  nz = m_npointsZ/2;
216  }
217  else
218  {
219  nz = 1;
220  }
221 
222  m_mat = Array<OneD, CoupledSolverMatrices> (nz);
223 
224  // mean mode or 2D mode.
225  SetUpCoupledMatrix(lambda,Advfield,IsLinearNSEquation,0,m_mat[0],m_locToGloMap[0]);
226 
227  for(n = 1; n < nz; ++n)
228  {
229  NekDouble beta = 2*M_PI*n/m_LhomZ;
230 
231  NekDouble lam = lambda + m_kinvis*beta*beta;
232 
233  SetUpCoupledMatrix(lam,Advfield,IsLinearNSEquation,n,m_mat[n],m_locToGloMap[n]);
234  }
235  }
236 
237  }
bool m_singleMode
Flag to determine if single homogeneous mode is used.
NekDouble m_kinvis
Kinematic viscosity.
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
Array< OneD, CoupledSolverMatrices > m_mat
int m_npointsZ
number of points in Z direction (if homogeneous)
double NekDouble
void SetUpCoupledMatrix(const NekDouble lambda=0.0, const Array< OneD, Array< OneD, NekDouble > > &Advfield=NullNekDoubleArrayofArray, bool IsLinearNSEquation=true)
static const NekDouble kNekUnsetDouble
Array< OneD, CoupledLocalToGlobalC0ContMapSharedPtr > m_locToGloMap
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
bool m_zeroMode
Id to identify when single mode is mean mode (i.e. beta=0);.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:218
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 376 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, Nektar::SolverUtils::EquationSystem::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().

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

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

1504  {
1505  const unsigned int ncmpt = m_velocity.num_elements();
1506  Array<OneD, Array<OneD, NekDouble> > forcing_phys(ncmpt);
1507  Array<OneD, Array<OneD, NekDouble> > forcing (ncmpt);
1508 
1509  for(int i = 0; i < ncmpt; ++i)
1510  {
1511  forcing_phys[i] = Array<OneD, NekDouble> (m_fields[m_velocity[0]]->GetNpoints(), 0.0);
1512  forcing[i] = Array<OneD, NekDouble> (m_fields[m_velocity[0]]->GetNcoeffs(),0.0);
1513  }
1514 
1515  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
1516  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
1517  {
1518  const NekDouble time = 0;
1519  (*x)->Apply(m_fields, forcing_phys, forcing_phys, time);
1520  }
1521  for (unsigned int i = 0; i < ncmpt; ++i)
1522  {
1523  m_fields[i]->IProductWRTBase(forcing_phys[i], forcing[i]);
1524  }
1525 
1526  SolveLinearNS(forcing);
1527  }
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 1873 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().

1874  {
1875  int i,n;
1876  Array<OneD, MultiRegions::ExpListSharedPtr> vel_fields(m_velocity.num_elements());
1877  Array<OneD, Array<OneD, NekDouble> > force(m_velocity.num_elements());
1878 
1880  {
1881  int ncoeffsplane = m_fields[m_velocity[0]]->GetPlane(0)->GetNcoeffs();
1882  for(n = 0; n < m_npointsZ/2; ++n)
1883  {
1884  // Get the a Fourier mode of velocity and forcing components.
1885  for(i = 0; i < m_velocity.num_elements(); ++i)
1886  {
1887  vel_fields[i] = m_fields[m_velocity[i]]->GetPlane(2*n);
1888  // Note this needs to correlate with how we pass forcing
1889  force[i] = forcing[i] + 2*n*ncoeffsplane;
1890  }
1891 
1892  SolveLinearNS(force,vel_fields,m_pressure->GetPlane(2*n),n);
1893  }
1894  for(i = 0; i < m_velocity.num_elements(); ++i)
1895  {
1896  m_fields[m_velocity[i]]->SetPhysState(false);
1897  }
1898  m_pressure->SetPhysState(false);
1899  }
1900  else
1901  {
1902  for(i = 0; i < m_velocity.num_elements(); ++i)
1903  {
1904  vel_fields[i] = m_fields[m_velocity[i]];
1905  // Note this needs to correlate with how we pass forcing
1906  force[i] = forcing[i];
1907  }
1908  SolveLinearNS(force,vel_fields,m_pressure);
1909  }
1910  }
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 1912 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, Nektar::SolverUtils::EquationSystem::m_singleMode, Nektar::IncNavierStokes::m_velocity, Nektar::Transpose(), and Vmath::Zero().

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

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

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

1370  {
1371  int i;
1372  Array<OneD, Array< OneD, NekDouble> > F(m_nConvectiveFields);
1373  NekDouble lambda = 1.0/aii_Dt;
1374  static NekDouble lambda_store;
1375  Array <OneD, Array<OneD, NekDouble> > forcing(m_velocity.num_elements());
1376  // Matrix solution
1377  if(fabs(lambda_store - lambda) > 1e-10)
1378  {
1379  cout << "Setting up Stokes matrix problem [.";
1380  fflush(stdout);
1381  SetUpCoupledMatrix(lambda);
1382  cout << "]" << endl;
1383  lambda_store = lambda;
1384  }
1385 
1386  SetBoundaryConditions(time);
1387 
1388  // Forcing for advection solve
1389  for(i = 0; i < m_velocity.num_elements(); ++i)
1390  {
1391  m_fields[m_velocity[i]]->IProductWRTBase(inarray[i],m_fields[m_velocity[i]]->UpdateCoeffs());
1392  Vmath::Smul(m_fields[m_velocity[i]]->GetNcoeffs(),lambda,m_fields[m_velocity[i]]->GetCoeffs(), 1,m_fields[m_velocity[i]]->UpdateCoeffs(),1);
1393  forcing[i] = m_fields[m_velocity[i]]->GetCoeffs();
1394  }
1395 
1396  SolveLinearNS(forcing);
1397 
1398  for(i = 0; i < m_velocity.num_elements(); ++i)
1399  {
1400  m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),outarray[i]);
1401  }
1402  }
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 1207 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().

1208  {
1209  switch(m_equationType)
1210  {
1211  case eUnsteadyStokes:
1212  case eUnsteadyNavierStokes:
1213  {
1214 
1215 // LibUtilities::TimeIntegrationMethod intMethod;
1216 // std::string TimeIntStr = m_session->GetSolverInfo("TIMEINTEGRATIONMETHOD");
1217 // int i;
1218 // for(i = 0; i < (int) LibUtilities::SIZE_TimeIntegrationMethod; ++i)
1219 // {
1220 // if(boost::iequals(LibUtilities::TimeIntegrationMethodMap[i],TimeIntStr))
1221 // {
1222 // intMethod = (LibUtilities::TimeIntegrationMethod)i;
1223 // break;
1224 // }
1225 // }
1226 //
1227 // ASSERTL0(i != (int) LibUtilities::SIZE_TimeIntegrationMethod, "Invalid time integration type.");
1228 //
1229 // m_integrationScheme = LibUtilities::GetTimeIntegrationWrapperFactory().CreateInstance(LibUtilities::TimeIntegrationMethodMap[intMethod]);
1230 
1231  // Could defind this from IncNavierStokes class?
1233 
1235 
1236  // Set initial condition using time t=0
1237 
1238  SetInitialConditions(0.0);
1239 
1240  }
1241  case eSteadyStokes:
1242  SetUpCoupledMatrix(0.0);
1243  break;
1244  case eSteadyOseen:
1245  {
1246  Array<OneD, Array<OneD, NekDouble> > AdvField(m_velocity.num_elements());
1247  for(int i = 0; i < m_velocity.num_elements(); ++i)
1248  {
1249  AdvField[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1250  }
1251 
1252  ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1253  "Advection Velocity section must be defined in "
1254  "session file.");
1255 
1256  std::vector<std::string> fieldStr;
1257  for(int i = 0; i < m_velocity.num_elements(); ++i)
1258  {
1259  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1260  }
1261  EvaluateFunction(fieldStr,AdvField,"AdvectionVelocity");
1262 
1263  SetUpCoupledMatrix(0.0,AdvField,false);
1264  }
1265  break;
1266  case eSteadyNavierStokes:
1267  {
1268  m_session->LoadParameter("KinvisMin", m_kinvisMin);
1269  m_session->LoadParameter("KinvisPercentage", m_KinvisPercentage);
1270  m_session->LoadParameter("Tolerence", m_tol);
1271  m_session->LoadParameter("MaxIteration", m_maxIt);
1272  m_session->LoadParameter("MatrixSetUpStep", m_MatrixSetUpStep);
1273  m_session->LoadParameter("Restart", m_Restart);
1274 
1275 
1277 
1278  if (m_Restart == 1)
1279  {
1280  ASSERTL0(m_session->DefinesFunction("Restart"),
1281  "Restart section must be defined in session file.");
1282 
1283  Array<OneD, Array<OneD, NekDouble> > Restart(m_velocity.num_elements());
1284  for(int i = 0; i < m_velocity.num_elements(); ++i)
1285  {
1286  Restart[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1287  }
1288  std::vector<std::string> fieldStr;
1289  for(int i = 0; i < m_velocity.num_elements(); ++i)
1290  {
1291  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1292  }
1293  EvaluateFunction(fieldStr, Restart, "Restart");
1294 
1295  for(int i = 0; i < m_velocity.num_elements(); ++i)
1296  {
1297  m_fields[m_velocity[i]]->FwdTrans_IterPerExp(Restart[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1298  }
1299  cout << "Saving the RESTART file for m_kinvis = "<< m_kinvis << " (<=> Re = " << 1/m_kinvis << ")" <<endl;
1300  }
1301  else //We solve the Stokes Problem
1302  {
1303 
1304  /*Array<OneD, Array<OneD, NekDouble> >ZERO(m_velocity.num_elements());
1305  *
1306  * for(int i = 0; i < m_velocity.num_elements(); ++i)
1307  * {
1308  * ZERO[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1309  * m_fields[m_velocity[i]]->FwdTrans(ZERO[i], m_fields[m_velocity[i]]->UpdateCoeffs());
1310  }*/
1311 
1312  SetUpCoupledMatrix(0.0);
1313  m_initialStep = true;
1314  m_counter=1;
1315  //SolveLinearNS(m_ForcingTerm_Coeffs);
1316  Solve();
1317  m_initialStep = false;
1318  cout << "Saving the Stokes Flow for m_kinvis = "<< m_kinvis << " (<=> Re = " << 1/m_kinvis << ")" <<endl;
1319  }
1320  }
1321  break;
1322  case eSteadyLinearisedNS:
1323  {
1324  SetInitialConditions(0.0);
1325 
1326  Array<OneD, Array<OneD, NekDouble> > AdvField(m_velocity.num_elements());
1327  for(int i = 0; i < m_velocity.num_elements(); ++i)
1328  {
1329  AdvField[i] = Array<OneD, NekDouble> (m_fields[m_velocity[i]]->GetTotPoints(),0.0);
1330  }
1331 
1332  ASSERTL0(m_session->DefinesFunction("AdvectionVelocity"),
1333  "Advection Velocity section must be defined in "
1334  "session file.");
1335 
1336  std::vector<std::string> fieldStr;
1337  for(int i = 0; i < m_velocity.num_elements(); ++i)
1338  {
1339  fieldStr.push_back(m_boundaryConditions->GetVariable(m_velocity[i]));
1340  }
1341  EvaluateFunction(fieldStr,AdvField,"AdvectionVelocity");
1342 
1343  SetUpCoupledMatrix(m_lambda,AdvField,true);
1344  }
1345  break;
1346  case eNoEquationType:
1347  default:
1348  ASSERTL0(false,"Unknown or undefined equation type for CoupledLinearNS");
1349  }
1350  }
EquationType m_equationType
equation type;
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
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 1429 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(), Nektar::Timer::TimePerTest(), and Nektar::SolverUtils::UnsteadySystem::v_DoSolve().

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

References Nektar::SolverUtils::AddSummaryItem().

1203  {
1204  SolverUtils::AddSummaryItem(s, "Solver Type", "Coupled Linearised NS");
1205  }
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 2227 of file CoupledLinearNS.cpp.

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

2228  {
2229  return m_session->GetVariables().size();
2230  }
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 66 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, Nektar::SolverUtils::EquationSystem::m_singleMode, Nektar::SolverUtils::EquationSystem::m_useFFT, Nektar::IncNavierStokes::m_velocity, m_zeroMode, and Nektar::IncNavierStokes::v_InitObject().

67  {
69 
70  int i;
71  int expdim = m_graph->GetMeshDimension();
72 
73  // Get Expansion list for orthogonal expansion at p-2
74  const SpatialDomains::ExpansionMap &pressure_exp = GenPressureExp(m_graph->GetExpansions("u"));
75 
76  m_nConvectiveFields = m_fields.num_elements();
77  if(boost::iequals(m_boundaryConditions->GetVariable(m_nConvectiveFields-1), "p"))
78  {
79  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");
80  }
81  // Decide how to declare explist for pressure.
82  if(expdim == 2)
83  {
84  int nz;
85 
87  {
88  ASSERTL0(m_fields.num_elements() > 2,"Expect to have three at least three components of velocity variables");
89  LibUtilities::BasisKey Homo1DKey = m_fields[0]->GetHomogeneousBasis()->GetBasisKey();
90 
92 
93  ASSERTL1(m_npointsZ%2==0,"Non binary number of planes have been specified");
94  nz = m_npointsZ/2;
95 
96  }
97  else
98  {
99  //m_pressure2 = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(m_session, pressure_exp);
101  nz = 1;
102  }
103 
104  Array<OneD, MultiRegions::ExpListSharedPtr> velocity(m_velocity.num_elements());
105  for(i =0 ; i < m_velocity.num_elements(); ++i)
106  {
107  velocity[i] = m_fields[m_velocity[i]];
108  }
109 
110  // Set up Array of mappings
111  m_locToGloMap = Array<OneD, CoupledLocalToGlobalC0ContMapSharedPtr> (nz);
112 
113  if(m_singleMode)
114  {
115 
116  ASSERTL0(nz <=2 ,"For single mode calculation can only have nz <= 2");
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;
bool m_singleMode
Flag to determine if single homogeneous mode is used.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
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:50
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;.
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:218
virtual void v_InitObject()
Init object for UnsteadySystem class.
enum HomogeneousType m_HomogeneousType
std::map< int, ExpansionShPtr > ExpansionMap
Definition: MeshGraph.h:174
bool Nektar::CoupledLinearNS::v_NegatedOp ( void  )
privatevirtual

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

Reimplemented from Nektar::SolverUtils::EquationSystem.

Definition at line 1498 of file CoupledLinearNS.cpp.

1499  {
1500  return true;
1501  }
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 2189 of file CoupledLinearNS.cpp.

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

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

Virtual function for transformation to physical space.

Reimplemented from Nektar::IncNavierStokes.

Definition at line 1405 of file CoupledLinearNS.cpp.

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

1406  {
1407  int nfields = m_fields.num_elements();
1408  for (int k=0 ; k < nfields; ++k)
1409  {
1410  //Backward Transformation in physical space for time evolution
1411  m_fields[k]->BwdTrans_IterPerExp(m_fields[k]->GetCoeffs(),
1412  m_fields[k]->UpdatePhys());
1413  }
1414 
1415  }
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 1417 of file CoupledLinearNS.cpp.

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

1418  {
1419  int nfields = m_fields.num_elements();
1420  for (int k=0 ; k < nfields; ++k)
1421  {
1422  //Forward Transformation in physical space for time evolution
1423  m_fields[k]->FwdTrans_IterPerExp(m_fields[k]->GetPhys(),
1424  m_fields[k]->UpdateCoeffs());
1425 
1426  }
1427  }
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 172 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 173 of file CoupledLinearNS.h.

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

NekDouble Nektar::CoupledLinearNS::m_kinvisMin
private

Definition at line 178 of file CoupledLinearNS.h.

Referenced by v_DoInitialise(), and v_DoSolve().

NekDouble Nektar::CoupledLinearNS::m_KinvisPercentage
private

Definition at line 180 of file CoupledLinearNS.h.

Referenced by Continuation(), and v_DoInitialise().

NekDouble Nektar::CoupledLinearNS::m_kinvisStep
private

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

Referenced by SetUpCoupledMatrix(), and SolveLinearNS().

int Nektar::CoupledLinearNS::m_MatrixSetUpStep
private

Definition at line 177 of file CoupledLinearNS.h.

Referenced by SolveSteadyNavierStokes(), and v_DoInitialise().

int Nektar::CoupledLinearNS::m_maxIt
private

Definition at line 175 of file CoupledLinearNS.h.

Referenced by v_DoInitialise().

int Nektar::CoupledLinearNS::m_Restart
private

Definition at line 176 of file CoupledLinearNS.h.

Referenced by v_DoInitialise().

NekDouble Nektar::CoupledLinearNS::m_tol
private

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

Referenced by SetUpCoupledMatrix(), and v_InitObject().