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

#include <CoupledLinearNS.h>

Inheritance diagram for Nektar::CoupledLinearNS:
Inheritance graph
[legend]
Collaboration diagram for Nektar::CoupledLinearNS:
Collaboration graph
[legend]

Public Member Functions

void SetUpCoupledMatrix (const NekDouble lambda=0.0, const Array< OneD, Array< OneD, NekDouble > > &Advfield=NullNekDoubleArrayofArray, bool IsLinearNSEquation=true)
 
const SpatialDomains::ExpansionMapGenPressureExp (const SpatialDomains::ExpansionMap &VelExp)
 
void Solve (void)
 
void SolveLinearNS (const Array< OneD, Array< OneD, NekDouble > > &forcing)
 
void SolveLinearNS (const Array< OneD, Array< OneD, NekDouble > > &forcing, Array< OneD, MultiRegions::ExpListSharedPtr > &fields, MultiRegions::ExpListSharedPtr &pressure, const int HomogeneousMode=0)
 
void SolveUnsteadyStokesSystem (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time, const NekDouble a_iixDt)
 
void EvaluateAdvection (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
 
void SolveSteadyNavierStokes (void)
 
void Continuation (void)
 
void EvaluateNewtonRHS (Array< OneD, Array< OneD, NekDouble > > &Velocity, Array< OneD, Array< OneD, NekDouble > > &outarray)
 
void InfNorm (Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
 
void L2Norm (Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
 
void DefineForcingTerm (void)
 
- Public Member Functions inherited from Nektar::IncNavierStokes
virtual ~IncNavierStokes ()
 
virtual void v_GetFluxVector (const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
 
virtual void v_NumericalFlux (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
 
int GetNConvectiveFields (void)
 
Array< OneD, int > & GetVelocity (void)
 
Array< OneD, NekDoubleGetElmtCFLVals (void)
 
NekDouble GetCFLEstimate (int &elmtid)
 
void AddForcing (const SolverUtils::ForcingSharedPtr &pForce)
 
- Public Member Functions inherited from Nektar::SolverUtils::AdvectionSystem
SOLVER_UTILS_EXPORT AdvectionSystem (const LibUtilities::SessionReaderSharedPtr &pSession)
 
virtual SOLVER_UTILS_EXPORT ~AdvectionSystem ()
 
AdvectionSharedPtr GetAdvObject ()
 Returns the advection object held by this instance. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
virtual SOLVER_UTILS_EXPORT ~UnsteadySystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep (const Array< OneD, const Array< OneD, NekDouble > > &inarray)
 Calculate the larger time-step mantaining the problem stable. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT void SetUpTraceNormals (void)
 
SOLVER_UTILS_EXPORT void InitObject ()
 Initialises the members of this object. More...
 
SOLVER_UTILS_EXPORT void DoInitialise ()
 Perform any initialisation necessary before solving the problem. More...
 
SOLVER_UTILS_EXPORT void DoSolve ()
 Solve the problem. More...
 
SOLVER_UTILS_EXPORT void TransCoeffToPhys ()
 Transform from coefficient to physical space. More...
 
SOLVER_UTILS_EXPORT void TransPhysToCoeff ()
 Transform from physical to coefficient space. More...
 
SOLVER_UTILS_EXPORT void Output ()
 Perform output operations after solve. More...
 
SOLVER_UTILS_EXPORT NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation. More...
 
SOLVER_UTILS_EXPORT std::string GetSessionName ()
 Get Session name. More...
 
template<class T >
boost::shared_ptr< T > as ()
 
SOLVER_UTILS_EXPORT void ResetSessionName (std::string newname)
 Reset Session name. More...
 
SOLVER_UTILS_EXPORT LibUtilities::SessionReaderSharedPtr GetSession ()
 Get Session name. More...
 
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure ()
 Get pressure field if available. More...
 
SOLVER_UTILS_EXPORT void PrintSummary (std::ostream &out)
 Print a summary of parameters and solver characteristics. More...
 
SOLVER_UTILS_EXPORT void SetLambda (NekDouble lambda)
 Set parameter m_lambda. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
 Evaluates a function as specified in the session file. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (std::vector< std::string > pFieldNames, Array< OneD, Array< OneD, NekDouble > > &pFields, const std::string &pName, const int domain=0)
 Populate given fields with the function from session. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (std::vector< std::string > pFieldNames, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const std::string &pName, const int domain=0)
 Populate given fields with the function from session. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, const NekDouble &pTime=0.0, const int domain=0)
 
SOLVER_UTILS_EXPORT std::string DescribeFunction (std::string pFieldName, const std::string &pFunctionName, const int domain)
 Provide a description of a function for a given field name. More...
 
SOLVER_UTILS_EXPORT void InitialiseBaseFlow (Array< OneD, Array< OneD, NekDouble > > &base)
 Perform initialisation of the base flow. More...
 
SOLVER_UTILS_EXPORT void SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 Initialise the data in the dependent fields. More...
 
SOLVER_UTILS_EXPORT void EvaluateExactSolution (int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 Evaluates an exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised=false)
 Compute the L2 error between fields and a given exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, bool Normalised=false)
 Compute the L2 error of the fields. More...
 
SOLVER_UTILS_EXPORT Array< OneD, NekDoubleErrorExtraPoints (unsigned int field)
 Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf]. More...
 
SOLVER_UTILS_EXPORT void WeakAdvectionGreensDivergenceForm (const Array< OneD, Array< OneD, NekDouble > > &F, Array< OneD, NekDouble > &outarray)
 Compute the inner product $ (\nabla \phi \cdot F) $. More...
 
SOLVER_UTILS_EXPORT void WeakAdvectionDivergenceForm (const Array< OneD, Array< OneD, NekDouble > > &F, Array< OneD, NekDouble > &outarray)
 Compute the inner product $ (\phi, \nabla \cdot F) $. More...
 
SOLVER_UTILS_EXPORT void WeakAdvectionNonConservativeForm (const Array< OneD, Array< OneD, NekDouble > > &V, const Array< OneD, const NekDouble > &u, Array< OneD, NekDouble > &outarray, bool UseContCoeffs=false)
 Compute the inner product $ (\phi, V\cdot \nabla u) $. More...
 
f SOLVER_UTILS_EXPORT void AdvectionNonConservativeForm (const Array< OneD, Array< OneD, NekDouble > > &V, const Array< OneD, const NekDouble > &u, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
 Compute the non-conservative advection. More...
 
SOLVER_UTILS_EXPORT void WeakDGAdvection (const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false, int nvariables=0)
 Calculate the weak discontinuous Galerkin advection. More...
 
SOLVER_UTILS_EXPORT void WeakDGDiffusion (const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false)
 Calculate weak DG Diffusion in the LDG form. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n)
 Write checkpoint file of m_fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write checkpoint file of custom data fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow (const int n)
 Write base flow file of m_fields. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname)
 Write field data to the given filename. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write input fields to the given filename. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 Input field data from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFldToMultiDomains (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int ndomains)
 Input field data from the given file to multiple domains. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, std::vector< std::string > &fieldStr, Array< OneD, Array< OneD, NekDouble > > &coeffs)
 Output a field. Input field data into array from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, MultiRegions::ExpListSharedPtr &pField, std::string &pFieldName)
 Output a field. Input field data into ExpList from the given file. More...
 
SOLVER_UTILS_EXPORT void ScanForHistoryPoints ()
 Builds map of which element holds each history point. More...
 
SOLVER_UTILS_EXPORT void WriteHistoryData (std::ostream &out)
 Probe each history point and write to file. More...
 
SOLVER_UTILS_EXPORT void SessionSummary (SummaryList &vSummary)
 Write out a session summary. More...
 
SOLVER_UTILS_EXPORT Array< OneD, MultiRegions::ExpListSharedPtr > & UpdateFields ()
 
SOLVER_UTILS_EXPORT LibUtilities::FieldMetaDataMapUpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated. More...
 
SOLVER_UTILS_EXPORT NekDouble GetFinalTime ()
 Return final time. More...
 
SOLVER_UTILS_EXPORT int GetNcoeffs ()
 
SOLVER_UTILS_EXPORT int GetNcoeffs (const int eid)
 
SOLVER_UTILS_EXPORT int GetNumExpModes ()
 
SOLVER_UTILS_EXPORT const Array< OneD, int > GetNumExpModesPerExp ()
 
SOLVER_UTILS_EXPORT int GetNvariables ()
 
SOLVER_UTILS_EXPORT const std::string GetVariable (unsigned int i)
 
SOLVER_UTILS_EXPORT int GetTraceTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTraceNpoints ()
 
SOLVER_UTILS_EXPORT int GetExpSize ()
 
SOLVER_UTILS_EXPORT int GetPhys_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetCoeff_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTotPoints (int n)
 
SOLVER_UTILS_EXPORT int GetNpoints ()
 
SOLVER_UTILS_EXPORT int GetNumElmVelocity ()
 
SOLVER_UTILS_EXPORT int GetSteps ()
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void CopyFromPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void CopyToPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void SetStepsToOne ()
 
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
 
SOLVER_UTILS_EXPORT void FwdTransFields ()
 
SOLVER_UTILS_EXPORT void GetFluxVector (const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
 
SOLVER_UTILS_EXPORT void GetFluxVector (const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &fluxX, Array< OneD, Array< OneD, NekDouble > > &fluxY)
 
SOLVER_UTILS_EXPORT void GetFluxVector (const int i, const int j, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
 
SOLVER_UTILS_EXPORT void NumericalFlux (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
 
SOLVER_UTILS_EXPORT void NumericalFlux (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxX, Array< OneD, Array< OneD, NekDouble > > &numfluxY)
 
SOLVER_UTILS_EXPORT void NumFluxforScalar (const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
 
SOLVER_UTILS_EXPORT void NumFluxforVector (const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int NoCaseStringCompare (const string &s1, const string &s2)
 Perform a case-insensitive string comparison. More...
 

Static Public Member Functions

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

Public Attributes

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

Static Public Attributes

static std::string className = SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction("CoupledLinearisedNS", CoupledLinearNS::create)
 Name of class. More...
 

Protected Member Functions

 CoupledLinearNS (const LibUtilities::SessionReaderSharedPtr &pSesssion)
 
virtual void v_InitObject ()
 Init object for UnsteadySystem class. More...
 
- Protected Member Functions inherited from Nektar::IncNavierStokes
 IncNavierStokes (const LibUtilities::SessionReaderSharedPtr &pSession)
 Constructor. More...
 
EquationType GetEquationType (void)
 
void EvaluateAdvectionTerms (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
 
void WriteModalEnergy (void)
 
void SetBoundaryConditions (NekDouble time)
 time dependent boundary conditions updating More...
 
void SetRadiationBoundaryForcing (int fieldid)
 Set Radiation forcing term. More...
 
bool CalcSteadyState (void)
 evaluate steady state More...
 
virtual MultiRegions::ExpListSharedPtr v_GetPressure ()
 
virtual bool v_PreIntegrate (int step)
 
virtual bool v_PostIntegrate (int step)
 
- Protected Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
SOLVER_UTILS_EXPORT UnsteadySystem (const LibUtilities::SessionReaderSharedPtr &pSession)
 Initialises UnsteadySystem class members. More...
 
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator ()
 Get the maximum timestep estimator for cfl control. More...
 
virtual SOLVER_UTILS_EXPORT void v_AppendOutput1D (Array< OneD, Array< OneD, NekDouble > > &solution1D)
 Print the solution at each solution point in a txt file. More...
 
virtual SOLVER_UTILS_EXPORT void v_NumericalFlux (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxX, Array< OneD, Array< OneD, NekDouble > > &numfluxY)
 
virtual SOLVER_UTILS_EXPORT void v_NumFluxforScalar (const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
 
virtual SOLVER_UTILS_EXPORT void v_NumFluxforVector (const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
 
virtual SOLVER_UTILS_EXPORT NekDouble v_GetTimeStep (const Array< OneD, const Array< OneD, NekDouble > > &inarray)
 Return the timestep to be used for the next step in the time-marching loop. More...
 
virtual SOLVER_UTILS_EXPORT bool v_SteadyStateCheck (int step)
 
SOLVER_UTILS_EXPORT void CheckForRestartTime (NekDouble &time)
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession)
 Initialises EquationSystem class members. More...
 
int nocase_cmp (const string &s1, const string &s2)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Virtual function for the L_inf error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT NekDouble v_L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
 Virtual function for the L_2 error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 
virtual SOLVER_UTILS_EXPORT void v_EvaluateExactSolution (unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetUpBaseFields (SpatialDomains::MeshGraphSharedPtr &mesh)
 
SOLVER_UTILS_EXPORT void ImportFldBase (std::string pInfile, SpatialDomains::MeshGraphSharedPtr pGraph)
 
virtual SOLVER_UTILS_EXPORT void v_ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 

Private Member Functions

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

Private Attributes

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

Friends

class MemoryManager< CoupledLinearNS >
 

Additional Inherited Members

- Protected Types inherited from Nektar::SolverUtils::EquationSystem
enum  HomogeneousType { eHomogeneous1D, eHomogeneous2D, eHomogeneous3D, eNotHomogeneous }
 Parameter for homogeneous expansions. More...
 
- Protected Attributes inherited from Nektar::IncNavierStokes
ExtrapolateSharedPtr m_extrapolation
 
std::ofstream m_mdlFile
 modal energy file More...
 
bool m_subSteppingScheme
 bool to identify if using a substepping scheme More...
 
bool m_SmoothAdvection
 bool to identify if advection term smoothing is requested More...
 
LibUtilities::TimeIntegrationWrapperSharedPtr m_subStepIntegrationScheme
 
std::vector< SolverUtils::ForcingSharedPtrm_forcing
 Forcing terms. More...
 
int m_nConvectiveFields
 Number of fields to be convected;. More...
 
Array< OneD, int > m_velocity
 int which identifies which components of m_fields contains the velocity (u,v,w); More...
 
MultiRegions::ExpListSharedPtr m_pressure
 Pointer to field holding pressure field. More...
 
NekDouble m_kinvis
 Kinematic viscosity. More...
 
int m_energysteps
 dump energy to file at steps time More...
 
int m_cflsteps
 dump cfl estimate More...
 
int m_steadyStateSteps
 Check for steady state at step interval. More...
 
NekDouble m_steadyStateTol
 Tolerance to which steady state should be evaluated at. More...
 
EquationType m_equationType
 equation type; More...
 
Array< OneD, Array< OneD, int > > m_fieldsBCToElmtID
 Mapping from BCs to Elmt IDs. More...
 
Array< OneD, Array< OneD, int > > m_fieldsBCToTraceID
 Mapping from BCs to Elmt Edge IDs. More...
 
Array< OneD, Array< OneD, NekDouble > > m_fieldsRadiationFactor
 RHS Factor for Radiation Condition. More...
 
int m_intSteps
 Number of time integration steps AND Order of extrapolation for pressure boundary conditions. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::AdvectionSystem
SolverUtils::AdvectionSharedPtr m_advObject
 Advection term. More...
 
- Protected Attributes inherited from Nektar::SolverUtils::UnsteadySystem
int m_infosteps
 Number of time steps between outputting status information. More...
 
LibUtilities::TimeIntegrationWrapperSharedPtr m_intScheme
 Wrapper to the time integration scheme. More...
 
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use. More...
 
LibUtilities::TimeIntegrationSolutionSharedPtr m_intSoln
 
NekDouble m_epsilon
 
bool m_explicitDiffusion
 Indicates if explicit or implicit treatment of diffusion is used. More...
 
bool m_explicitAdvection
 Indicates if explicit or implicit treatment of advection is used. More...
 
bool m_explicitReaction
 Indicates if explicit or implicit treatment of reaction is used. More...
 
bool m_homoInitialFwd
 Flag to determine if simulation should start in homogeneous forward transformed state. More...
 
std::vector< int > m_intVariables
 
std::vector< FilterSharedPtrm_filters
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
map< std::string, Array< OneD, Array< OneD, float > > > m_interpWeights
 Map of the interpolation weights for a specific filename. More...
 
map< std::string, Array< OneD, Array< OneD, unsigned int > > > m_interpInds
 Map of the interpolation indices for a specific filename. More...
 
Array< OneD, MultiRegions::ExpListSharedPtrm_fields
 Array holding all dependent variables. More...
 
Array< OneD, MultiRegions::ExpListSharedPtrm_base
 Base fields. More...
 
Array< OneD, MultiRegions::ExpListSharedPtrm_derivedfields
 Array holding all dependent variables. More...
 
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
 Pointer to boundary conditions object. More...
 
SpatialDomains::MeshGraphSharedPtr m_graph
 Pointer to graph defining mesh. More...
 
std::string m_sessionName
 Name of the session. More...
 
NekDouble m_time
 Current time of simulation. More...
 
NekDouble m_fintime
 Finish time of the simulation. More...
 
NekDouble m_timestep
 Time step size. More...
 
NekDouble m_lambda
 Lambda constant in real system if one required. More...
 
NekDouble m_checktime
 Time between checkpoints. More...
 
int m_steps
 Number of steps to take. More...
 
int m_checksteps
 Number of steps between checkpoints. More...
 
int m_spacedim
 Spatial dimension (>= expansion dim). More...
 
int m_expdim
 Expansion dimension. More...
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform. More...
 
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations. More...
 
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More...
 
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous. More...
 
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction. More...
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_gradtan
 1 x nvariable x nq More...
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tanbasis
 2 x m_spacedim x nq More...
 
Array< OneD, bool > m_checkIfSystemSingular
 Flag to indicate if the fields should be checked for singularity. More...
 
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
 Map to identify relevant solver info to dump in output fields. More...
 
int m_NumQuadPointsError
 Number of Quadrature points used to work out the error. More...
 
enum HomogeneousType m_HomogeneousType
 
NekDouble m_LhomX
 physical length in X direction (if homogeneous) More...
 
NekDouble m_LhomY
 physical length in Y direction (if homogeneous) More...
 
NekDouble m_LhomZ
 physical length in Z direction (if homogeneous) More...
 
int m_npointsX
 number of points in X direction (if homogeneous) More...
 
int m_npointsY
 number of points in Y direction (if homogeneous) More...
 
int m_npointsZ
 number of points in Z direction (if homogeneous) More...
 
int m_HomoDirec
 number of homogenous directions More...
 
int m_NumMode
 Mode to use in case of single mode analysis. More...
 

Detailed Description

Set up expansion field for velocity and pressure, the local to global mapping arrays and the basic memory definitions for coupled matrix system

Definition at line 89 of file CoupledLinearNS.h.

Constructor & Destructor Documentation

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

Definition at line 57 of file CoupledLinearNS.cpp.

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

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

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

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

Referenced by v_DoInitialise().

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

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

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

Definition at line 1697 of file CoupledLinearNS.cpp.

References Nektar::IncNavierStokes::m_velocity.

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

References Nektar::IncNavierStokes::m_velocity.

Referenced by SolveSteadyNavierStokes().

1716  {
1717  for(int i = 0; i < m_velocity.num_elements(); ++i)
1718  {
1719  outarray[i] = 0.0;
1720  for(int j = 0; j < inarray[i].num_elements(); ++j)
1721  {
1722  outarray[i] += inarray[i][j]*inarray[i][j];
1723  }
1724  outarray[i]=sqrt(outarray[i]);
1725  cout << "L2Norm["<<i<<"] = "<< outarray[i] <<endl;
1726  }
1727  }
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 178 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().

179  {
180 
181  int nz;
182  if(m_singleMode)
183  {
184 
185  NekDouble lambda_imag;
186 
187  // load imaginary component of any potential shift
188  // Probably should be called from DriverArpack but not yet
189  // clear how to do this
190  m_session->LoadParameter("imagShift",lambda_imag,NekConstants::kNekUnsetDouble);
191  nz = 1;
193 
194  ASSERTL1(m_npointsZ <=2,"Only expected a maxmimum of two planes in single mode linear NS solver");
195 
196  if(m_zeroMode)
197  {
198  SetUpCoupledMatrix(lambda,Advfield,IsLinearNSEquation,0,m_mat[0],m_locToGloMap[0],lambda_imag);
199  }
200  else
201  {
202  NekDouble beta = 2*M_PI/m_LhomZ;
203  NekDouble lam = lambda + m_kinvis*beta*beta;
204 
205  SetUpCoupledMatrix(lam,Advfield,IsLinearNSEquation,1,m_mat[0],m_locToGloMap[0],lambda_imag);
206  }
207  }
208  else
209  {
210  int n;
211  if(m_npointsZ > 1)
212  {
213  nz = m_npointsZ/2;
214  }
215  else
216  {
217  nz = 1;
218  }
219 
221 
222  // mean mode or 2D mode.
223  SetUpCoupledMatrix(lambda,Advfield,IsLinearNSEquation,0,m_mat[0],m_locToGloMap[0]);
224 
225  for(n = 1; n < nz; ++n)
226  {
227  NekDouble beta = 2*M_PI*n/m_LhomZ;
228 
229  NekDouble lam = lambda + m_kinvis*beta*beta;
230 
231  SetUpCoupledMatrix(lam,Advfield,IsLinearNSEquation,n,m_mat[n],m_locToGloMap[n]);
232  }
233  }
234 
235  }
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:191
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 374 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().

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

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

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

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

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

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

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

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

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

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

1428  {
1429  switch(m_equationType)
1430  {
1431  case eUnsteadyStokes:
1432  case eUnsteadyNavierStokes:
1433  //AdvanceInTime(m_steps);
1434  UnsteadySystem::v_DoSolve();
1435  break;
1436  case eSteadyStokes:
1437  case eSteadyOseen:
1438  case eSteadyLinearisedNS:
1439  {
1440  Solve();
1441  break;
1442  }
1443  case eSteadyNavierStokes:
1444  {
1445  Timer Generaltimer;
1446  Generaltimer.Start();
1447 
1448  int Check(0);
1449 
1450  //Saving the init datas
1451  Checkpoint_Output(Check);
1452  Check++;
1453 
1454  cout<<"We execute INITIALLY SolveSteadyNavierStokes for m_kinvis = "<<m_kinvis<<" (<=> Re = "<<1/m_kinvis<<")"<<endl;
1456 
1457  while(m_kinvis > m_kinvisMin)
1458  {
1459  if (Check == 1)
1460  {
1461  cout<<"We execute SolveSteadyNavierStokes for m_kinvis = "<<m_kinvis<<" (<=> Re = "<<1/m_kinvis<<")"<<endl;
1463  Checkpoint_Output(Check);
1464  Check++;
1465  }
1466 
1467  Continuation();
1468 
1469  if (m_kinvis > m_kinvisMin)
1470  {
1471  cout<<"We execute SolveSteadyNavierStokes for m_kinvis = "<<m_kinvis<<" (<=> Re = "<<1/m_kinvis<<")"<<endl;
1473  Checkpoint_Output(Check);
1474  Check++;
1475  }
1476  }
1477 
1478 
1479  Generaltimer.Stop();
1480  cout<<"\nThe total calculation time is : " << Generaltimer.TimePerTest(1)/60 << " minute(s). \n\n";
1481 
1482  break;
1483  }
1484  case eNoEquationType:
1485  default:
1486  ASSERTL0(false,"Unknown or undefined equation type for CoupledLinearNS");
1487  }
1488  }
EquationType m_equationType
equation type;
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
NekDouble m_kinvis
Kinematic viscosity.
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
void Stop()
Definition: Timer.cpp:62
void Start()
Definition: Timer.cpp:51
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition: Timer.cpp:108
void Nektar::CoupledLinearNS::v_GenerateSummary ( SolverUtils::SummaryList s)
privatevirtual

Print a summary of time stepping parameters.

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

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Definition at line 1200 of file CoupledLinearNS.cpp.

References Nektar::SolverUtils::AddSummaryItem().

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

Implements Nektar::IncNavierStokes.

Definition at line 2225 of file CoupledLinearNS.cpp.

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

2226  {
2227  return m_session->GetVariables().size();
2228  }
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 64 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().

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

1497  {
1498  return true;
1499  }
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 2187 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().

2188  {
2189  std::vector<Array<OneD, NekDouble> > fieldcoeffs(m_fields.num_elements()+1);
2190  std::vector<std::string> variables(m_fields.num_elements()+1);
2191  int i;
2192 
2193  for(i = 0; i < m_fields.num_elements(); ++i)
2194  {
2195  fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
2196  variables[i] = m_boundaryConditions->GetVariable(i);
2197  }
2198 
2199  fieldcoeffs[i] = Array<OneD, NekDouble>(m_fields[0]->GetNcoeffs());
2200  // project pressure field to velocity space
2201  if(m_singleMode==true)
2202  {
2203  Array<OneD, NekDouble > tmpfieldcoeffs (m_fields[0]->GetNcoeffs()/2);
2204  m_pressure->GetPlane(0)->BwdTrans_IterPerExp(m_pressure->GetPlane(0)->GetCoeffs(), m_pressure->GetPlane(0)->UpdatePhys());
2205  m_pressure->GetPlane(1)->BwdTrans_IterPerExp(m_pressure->GetPlane(1)->GetCoeffs(), m_pressure->GetPlane(1)->UpdatePhys());
2206  m_fields[0]->GetPlane(0)->FwdTrans_IterPerExp(m_pressure->GetPlane(0)->GetPhys(),fieldcoeffs[i]);
2207  m_fields[0]->GetPlane(1)->FwdTrans_IterPerExp(m_pressure->GetPlane(1)->GetPhys(),tmpfieldcoeffs);
2208  for(int e=0; e<m_fields[0]->GetNcoeffs()/2; e++)
2209  {
2210  fieldcoeffs[i][e+m_fields[0]->GetNcoeffs()/2] = tmpfieldcoeffs[e];
2211  }
2212  }
2213  else
2214  {
2215  m_pressure->BwdTrans_IterPerExp(m_pressure->GetCoeffs(),m_pressure->UpdatePhys());
2216  m_fields[0]->FwdTrans_IterPerExp(m_pressure->GetPhys(),fieldcoeffs[i]);
2217  }
2218  variables[i] = "p";
2219 
2220  std::string outname = m_sessionName + ".fld";
2221 
2222  WriteFld(outname,m_fields[0],fieldcoeffs,variables);
2223  }
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 1403 of file CoupledLinearNS.cpp.

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

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

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

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