Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | Protected Attributes | Friends | List of all members
Nektar::CompressibleFlowSystem Class Reference

#include <CompressibleFlowSystem.h>

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

Public Member Functions

virtual ~CompressibleFlowSystem ()
 Destructor for CompressibleFlowSystem class. More...
 
NekDouble GetStabilityLimit (int n)
 Function to calculate the stability limit for DG/CG. More...
 
Array< OneD, NekDouble > GetStabilityLimitVector (const Array< OneD, int > &ExpOrder)
 Function to calculate the stability limit for DG/CG (a vector of them). More...
 
- Public Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
virtual SOLVER_UTILS_EXPORT ~UnsteadySystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep (const Array< OneD, const Array< OneD, NekDouble > > &inarray)
 Calculate the larger time-step mantaining the problem stable. More...
 
- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor. More...
 
SOLVER_UTILS_EXPORT void SetUpTraceNormals (void)
 
SOLVER_UTILS_EXPORT void InitObject ()
 Initialises the members of this object. More...
 
SOLVER_UTILS_EXPORT void DoInitialise ()
 Perform any initialisation necessary before solving the problem. More...
 
SOLVER_UTILS_EXPORT void DoSolve ()
 Solve the problem. More...
 
SOLVER_UTILS_EXPORT void TransCoeffToPhys ()
 Transform from coefficient to physical space. More...
 
SOLVER_UTILS_EXPORT void TransPhysToCoeff ()
 Transform from physical to coefficient space. More...
 
SOLVER_UTILS_EXPORT void Output ()
 Perform output operations after solve. More...
 
SOLVER_UTILS_EXPORT NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation. More...
 
SOLVER_UTILS_EXPORT std::string GetSessionName ()
 Get Session name. More...
 
template<class T >
boost::shared_ptr< T > as ()
 
SOLVER_UTILS_EXPORT void ResetSessionName (std::string newname)
 Reset Session name. More...
 
SOLVER_UTILS_EXPORT
LibUtilities::SessionReaderSharedPtr 		GetSession ()
 Get Session name. More...
 
SOLVER_UTILS_EXPORT
MultiRegions::ExpListSharedPtr 		GetPressure ()
 Get pressure field if available. More...
 
SOLVER_UTILS_EXPORT void PrintSummary (std::ostream &out)
 Print a summary of parameters and solver characteristics. More...
 
SOLVER_UTILS_EXPORT void SetLambda (NekDouble lambda)
 Set parameter m_lambda. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
 Evaluates a function as specified in the session file. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (std::vector< std::string > pFieldNames, Array< OneD, Array< OneD, NekDouble > > &pFields, const std::string &pName, const NekDouble &pTime=0.0, const int domain=0)
 Populate given fields with the function from session. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (std::vector< std::string > pFieldNames, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const std::string &pName, const NekDouble &pTime=0.0, const int domain=0)
 Populate given fields with the function from session. More...
 
SOLVER_UTILS_EXPORT void EvaluateFunction (std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, const NekDouble &pTime=0.0, const int domain=0)
 
SOLVER_UTILS_EXPORT std::string DescribeFunction (std::string pFieldName, const std::string &pFunctionName, const int domain)
 Provide a description of a function for a given field name. More...
 
SOLVER_UTILS_EXPORT void InitialiseBaseFlow (Array< OneD, Array< OneD, NekDouble > > &base)
 Perform initialisation of the base flow. More...
 
SOLVER_UTILS_EXPORT void SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 Initialise the data in the dependent fields. More...
 
SOLVER_UTILS_EXPORT void EvaluateExactSolution (int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 Evaluates an exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised=false)
 Compute the L2 error between fields and a given exact solution. More...
 
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, bool Normalised=false)
 Compute the L2 error of the fields. More...
 
SOLVER_UTILS_EXPORT Array
< OneD, NekDouble > 		ErrorExtraPoints (unsigned int field)
 Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf]. More...
 
SOLVER_UTILS_EXPORT void WeakAdvectionGreensDivergenceForm (const Array< OneD, Array< OneD, NekDouble > > &F, Array< OneD, NekDouble > &outarray)
 Compute the inner product $ (\nabla \phi \cdot F) $. More...
 
SOLVER_UTILS_EXPORT void WeakAdvectionDivergenceForm (const Array< OneD, Array< OneD, NekDouble > > &F, Array< OneD, NekDouble > &outarray)
 Compute the inner product $ (\phi, \nabla \cdot F) $. More...
 
SOLVER_UTILS_EXPORT void WeakAdvectionNonConservativeForm (const Array< OneD, Array< OneD, NekDouble > > &V, const Array< OneD, const NekDouble > &u, Array< OneD, NekDouble > &outarray, bool UseContCoeffs=false)
 Compute the inner product $ (\phi, V\cdot \nabla u) $. More...
 
f SOLVER_UTILS_EXPORT void AdvectionNonConservativeForm (const Array< OneD, Array< OneD, NekDouble > > &V, const Array< OneD, const NekDouble > &u, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
 Compute the non-conservative advection. More...
 
SOLVER_UTILS_EXPORT void WeakDGAdvection (const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false, int nvariables=0)
 Calculate the weak discontinuous Galerkin advection. More...
 
SOLVER_UTILS_EXPORT void WeakDGDiffusion (const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false)
 Calculate weak DG Diffusion in the LDG form. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n)
 Write checkpoint file of m_fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write checkpoint file of custom data fields. More...
 
SOLVER_UTILS_EXPORT void Checkpoint_BaseFlow (const int n)
 Write base flow file of m_fields. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname)
 Write field data to the given filename. More...
 
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname, MultiRegions::ExpListSharedPtr &field, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 Write input fields to the given filename. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 Input field data from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFldToMultiDomains (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int ndomains)
 Input field data from the given file to multiple domains. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, std::vector< std::string > &fieldStr, Array< OneD, Array< OneD, NekDouble > > &coeffs)
 Output a field. Input field data into array from the given file. More...
 
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, MultiRegions::ExpListSharedPtr &pField, std::string &pFieldName)
 Output a field. Input field data into ExpList from the given file. More...
 
SOLVER_UTILS_EXPORT void ScanForHistoryPoints ()
 Builds map of which element holds each history point. More...
 
SOLVER_UTILS_EXPORT void WriteHistoryData (std::ostream &out)
 Probe each history point and write to file. More...
 
SOLVER_UTILS_EXPORT void SessionSummary (SummaryList &vSummary)
 Write out a session summary. More...
 
SOLVER_UTILS_EXPORT Array
< OneD,
MultiRegions::ExpListSharedPtr > & 		UpdateFields ()
 
SOLVER_UTILS_EXPORT
LibUtilities::FieldMetaDataMap & 		UpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated. More...
 
SOLVER_UTILS_EXPORT NekDouble GetFinalTime ()
 Return final time. More...
 
SOLVER_UTILS_EXPORT int GetNcoeffs ()
 
SOLVER_UTILS_EXPORT int GetNcoeffs (const int eid)
 
SOLVER_UTILS_EXPORT int GetNumExpModes ()
 
SOLVER_UTILS_EXPORT const
Array< OneD, int > 		GetNumExpModesPerExp ()
 
SOLVER_UTILS_EXPORT int GetNvariables ()
 
SOLVER_UTILS_EXPORT const
std::string 		GetVariable (unsigned int i)
 
SOLVER_UTILS_EXPORT int GetTraceTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTraceNpoints ()
 
SOLVER_UTILS_EXPORT int GetExpSize ()
 
SOLVER_UTILS_EXPORT int GetPhys_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetCoeff_Offset (int n)
 
SOLVER_UTILS_EXPORT int GetTotPoints ()
 
SOLVER_UTILS_EXPORT int GetTotPoints (int n)
 
SOLVER_UTILS_EXPORT int GetNpoints ()
 
SOLVER_UTILS_EXPORT int GetNumElmVelocity ()
 
SOLVER_UTILS_EXPORT int GetSteps ()
 
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
 
SOLVER_UTILS_EXPORT void CopyFromPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void CopyToPhysField (const int i, Array< OneD, NekDouble > &output)
 
SOLVER_UTILS_EXPORT void SetSteps (const int steps)
 
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
 
SOLVER_UTILS_EXPORT void FwdTransFields ()
 
SOLVER_UTILS_EXPORT void GetFluxVector (const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
 
SOLVER_UTILS_EXPORT void GetFluxVector (const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &fluxX, Array< OneD, Array< OneD, NekDouble > > &fluxY)
 
SOLVER_UTILS_EXPORT void GetFluxVector (const int i, const int j, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
 
SOLVER_UTILS_EXPORT void NumericalFlux (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
 
SOLVER_UTILS_EXPORT void NumericalFlux (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxX, Array< OneD, Array< OneD, NekDouble > > &numfluxY)
 
SOLVER_UTILS_EXPORT void NumFluxforScalar (const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
 
SOLVER_UTILS_EXPORT void NumFluxforVector (const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
 
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
 
SOLVER_UTILS_EXPORT int NoCaseStringCompare (const std::string &s1, const std::string &s2)
 Perform a case-insensitive string comparison. More...
 
SOLVER_UTILS_EXPORT int GetCheckpointNumber ()
 
SOLVER_UTILS_EXPORT void SetCheckpointNumber (int num)
 
SOLVER_UTILS_EXPORT int GetCheckpointSteps ()
 
SOLVER_UTILS_EXPORT void SetCheckpointSteps (int num)
 
SOLVER_UTILS_EXPORT void SetTime (const NekDouble time)
 
SOLVER_UTILS_EXPORT void SetInitialStep (const int step)
 
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time. More...
 
virtual SOLVER_UTILS_EXPORT bool v_NegatedOp ()
 Virtual function to identify if operator is negated in DoSolve. More...
 

Static Public Member Functions

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

Static Public Attributes

static std::string className
 Name of class. More...
 

Protected Member Functions

 CompressibleFlowSystem (const LibUtilities::SessionReaderSharedPtr &pSession)
 
virtual void v_InitObject ()
 Initialization object for CompressibleFlowSystem class. More...
 
virtual void v_GenerateSummary (SolverUtils::SummaryList &s)
 Print a summary of time stepping parameters. More...
 
void GetFluxVector (const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
 Return the flux vector for the compressible Euler equations. More...
 
void GetFluxVectorDeAlias (const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
 Return the flux vector for the compressible Euler equations by using the de-aliasing technique. More...
 
void GetViscousFluxVector (const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &derivatives, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &viscousTensor)
 Return the flux vector for the LDG diffusion problem. More...
 
void GetFluxVectorPDESC (const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
 
void GetViscousFluxVectorDeAlias (const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &derivatives, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &viscousTensor)
 Return the flux vector for the LDG diffusion problem. More...
 
void SetCommonBC (const std::string &userDefStr, const int n, const NekDouble time, int &cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &inarray)
 Set boundary conditions which can be: a) Wall and Symmerty BCs implemented at CompressibleFlowSystem level since they are compressible solver specific; b) Time dependent BCs. More...
 
void WallBC (int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
 Wall boundary conditions for compressible flow problems. More...
 
void WallViscousBC (int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
 Wall boundary conditions for viscous compressible flow problems. More...
 
void SymmetryBC (int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
 Symmetry boundary conditions for compressible flow problems. More...
 
void RiemannInvariantBC (int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
 Outflow characteristic boundary conditions for compressible flow problems. More...
 
void PressureOutflowNonReflectiveBC (int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
 Pressure outflow non-reflective boundary conditions for compressible flow problems. More...
 
void PressureOutflowBC (int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
 Pressure outflow boundary conditions for compressible flow problems. More...
 
void PressureOutflowFileBC (int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
 Pressure outflow boundary conditions for compressible flow problems. More...
 
void PressureInflowFileBC (int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
 Pressure inflow boundary conditions for compressible flow problems where either the density and the velocities are assigned from a file or the full state is assigned from a file (depending on the problem type, either subsonic or supersonic). More...
 
void ExtrapOrder0BC (int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
 Extrapolation of order 0 for all the variables such that, at the boundaries, a trivial Riemann problem is solved. More...
 
void GetVelocityVector (const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
 
void GetSoundSpeed (const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure, Array< OneD, NekDouble > &soundspeed)
 
void GetMach (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &soundspeed, Array< OneD, NekDouble > &mach)
 
void GetTemperature (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure, Array< OneD, NekDouble > &temperature)
 
void GetPressure (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure)
 Calculate the pressure field $ p = (\gamma-1)(E-\frac{1}{2}\rho\| \mathbf{v} \|^2) $ assuming an ideal gas law. More...
 
void GetPressure (const Array< OneD, const Array< OneD, NekDouble > > &physfield, const Array< OneD, const Array< OneD, NekDouble > > &velocity, Array< OneD, NekDouble > &pressure)
 
void GetEnthalpy (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure, Array< OneD, NekDouble > &enthalpy)
 
void GetEntropy (const Array< OneD, const Array< OneD, NekDouble > > &physfield, const Array< OneD, const NekDouble > &pressure, const Array< OneD, const NekDouble > &temperature, Array< OneD, NekDouble > &entropy)
 
void GetSmoothArtificialViscosity (const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &eps_bar)
 
void GetDynamicViscosity (const Array< OneD, const NekDouble > &temperature, Array< OneD, NekDouble > &mu)
 
void GetStdVelocity (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &stdV)
 
virtual bool v_PostIntegrate (int step)
 
bool CalcSteadyState (bool output)
 
void GetSensor (const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, NekDouble > &Sensor, Array< OneD, NekDouble > &SensorKappa)
 
void GetElementDimensions (Array< OneD, Array< OneD, NekDouble > > &outarray, Array< OneD, NekDouble > &hmin)
 
void GetAbsoluteVelocity (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &Vtot)
 
void GetArtificialDynamicViscosity (const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &mu_var)
 
void SetVarPOrderElmt (const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &PolyOrder)
 
void GetForcingTerm (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > outarrayForcing)
 
virtual 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 void v_SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 
NekDouble GetGasConstant ()
 
NekDouble GetGamma ()
 
const Array< OneD, const Array
< OneD, NekDouble > > & 		GetVecLocs ()
 
const Array< OneD, const Array
< OneD, NekDouble > > & 		GetNormals ()
 
virtual void v_ExtraFldOutput (std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
 
- 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_DoSolve ()
 Solves an unsteady problem. More...
 
virtual SOLVER_UTILS_EXPORT void v_DoInitialise ()
 Sets up initial conditions. 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 > > &numflux)
 
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 bool v_PreIntegrate (int step)
 
virtual SOLVER_UTILS_EXPORT bool v_SteadyStateCheck (int step)
 
SOLVER_UTILS_EXPORT void CheckForRestartTime (NekDouble &time)
 
SOLVER_UTILS_EXPORT void SVVVarDiffCoeff (const Array< OneD, Array< OneD, NekDouble > > vel, StdRegions::VarCoeffMap &varCoeffMap)
 Evaluate the SVV diffusion coefficient according to Moura's paper where it should proportional to h time velocity. More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession)
 Initialises EquationSystem class members. More...
 
int nocase_cmp (const std::string &s1, const std::string &s2)
 
virtual SOLVER_UTILS_EXPORT
NekDouble 		v_LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Virtual function for the L_inf error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT
NekDouble 		v_L2Error (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray, bool Normalised=false)
 Virtual function for the L_2 error computation between fields and a given exact solution. More...
 
virtual SOLVER_UTILS_EXPORT void v_TransCoeffToPhys ()
 Virtual function for transformation to physical space. More...
 
virtual SOLVER_UTILS_EXPORT void v_TransPhysToCoeff ()
 Virtual function for transformation to coefficient space. More...
 
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_Output (void)
 
virtual SOLVER_UTILS_EXPORT
MultiRegions::ExpListSharedPtr 		v_GetPressure (void)
 

Protected Attributes

SolverUtils::RiemannSolverSharedPtr m_riemannSolver
 
SolverUtils::RiemannSolverSharedPtr m_riemannSolverLDG
 
SolverUtils::AdvectionSharedPtr m_advection
 
SolverUtils::DiffusionSharedPtr m_diffusion
 
Array< OneD, Array< OneD,
NekDouble > > 		m_vecLocs
 
NekDouble m_gamma
 
NekDouble m_pInf
 
NekDouble m_rhoInf
 
NekDouble m_uInf
 
NekDouble m_vInf
 
NekDouble m_wInf
 
NekDouble m_UInf
 
NekDouble m_gasConstant
 
NekDouble m_Twall
 
std::string m_ViscosityType
 
std::string m_shockCaptureType
 
std::string m_EqTypeStr
 
NekDouble m_mu
 
NekDouble m_Skappa
 
NekDouble m_Kappa
 
NekDouble m_mu0
 
NekDouble m_FacL
 
NekDouble m_FacH
 
NekDouble m_eps_max
 
NekDouble m_thermalConductivity
 
NekDouble m_Cp
 
NekDouble m_C1
 
NekDouble m_C2
 
NekDouble m_hFactor
 
NekDouble m_Prandtl
 
NekDouble m_amplitude
 
NekDouble m_omega
 
std::ofstream m_errFile
 
int m_steadyStateSteps
 
NekDouble m_steadyStateTol
 
std::vector
< SolverUtils::ForcingSharedPtr > 		m_forcing
 
StdRegions::StdQuadExpSharedPtr m_OrthoQuadExp
 
StdRegions::StdHexExpSharedPtr m_OrthoHexExp
 
bool m_smoothDiffusion
 
Array< OneD, NekDouble > m_pressureStorage
 
Array< OneD, Array< OneD,
NekDouble > > 		m_fieldStorage
 
Array< OneD, Array< OneD,
NekDouble > > 		m_un
 
- 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< FilterSharedPtr > m_filters
 
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator. More...
 
LibUtilities::SessionReaderSharedPtr m_session
 The session reader. More...
 
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output. More...
 
std::map< std::string, Array
< OneD, Array< OneD, float > > > 		m_interpWeights
 Map of the interpolation weights for a specific filename. More...
 
std::map< std::string, Array
< OneD, Array< OneD, unsigned
int > > > 		m_interpInds
 Map of the interpolation indices for a specific filename. More...
 
Array< OneD,
MultiRegions::ExpListSharedPtr > 		m_fields
 Array holding all dependent variables. More...
 
Array< OneD,
MultiRegions::ExpListSharedPtr > 		m_base
 Base fields. More...
 
Array< OneD,
MultiRegions::ExpListSharedPtr > 		m_derivedfields
 Array holding all dependent variables. More...
 
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
 Pointer to boundary conditions object. More...
 
SpatialDomains::MeshGraphSharedPtr m_graph
 Pointer to graph defining mesh. More...
 
std::string m_sessionName
 Name of the session. More...
 
NekDouble m_time
 Current time of simulation. More...
 
int m_initialStep
 Number of the step where the simulation should begin. More...
 
NekDouble m_fintime
 Finish time of the simulation. More...
 
NekDouble m_timestep
 Time step size. More...
 
NekDouble m_lambda
 Lambda constant in real system if one required. More...
 
std::set< std::string > m_loadedFields
 
NekDouble m_checktime
 Time between checkpoints. More...
 
int m_nchk
 Number of checkpoints written so far. More...
 
int m_steps
 Number of steps to take. More...
 
int m_checksteps
 Number of steps between checkpoints. More...
 
int m_spacedim
 Spatial dimension (>= expansion dim). More...
 
int m_expdim
 Expansion dimension. More...
 
bool m_singleMode
 Flag to determine if single homogeneous mode is used. More...
 
bool m_halfMode
 Flag to determine if half homogeneous mode is used. More...
 
bool m_multipleModes
 Flag to determine if use multiple homogenenous modes are used. More...
 
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform. More...
 
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations. More...
 
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More...
 
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous. More...
 
Array< OneD, Array< OneD,
NekDouble > > 		m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction. More...
 
Array< OneD, Array< OneD,
Array< OneD, NekDouble > > > 		m_gradtan
 1 x nvariable x nq More...
 
Array< OneD, Array< OneD,
Array< OneD, NekDouble > > > 		m_tanbasis
 2 x m_spacedim x nq More...
 
Array< OneD, bool > m_checkIfSystemSingular
 Flag to indicate if the fields should be checked for singularity. More...
 
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
 Map to identify relevant solver info to dump in output fields. More...
 
int m_NumQuadPointsError
 Number of Quadrature points used to work out the error. More...
 
enum HomogeneousType m_HomogeneousType
 
NekDouble m_LhomX
 physical length in X direction (if homogeneous) More...
 
NekDouble m_LhomY
 physical length in Y direction (if homogeneous) More...
 
NekDouble m_LhomZ
 physical length in Z direction (if homogeneous) More...
 
int m_npointsX
 number of points in X direction (if homogeneous) More...
 
int m_npointsY
 number of points in Y direction (if homogeneous) More...
 
int m_npointsZ
 number of points in Z direction (if homogeneous) More...
 
int m_HomoDirec
 number of homogenous directions More...
 

Friends

class MemoryManager< CompressibleFlowSystem >
 

Additional Inherited Members

- Public Attributes inherited from Nektar::SolverUtils::UnsteadySystem
NekDouble m_cflSafetyFactor
 CFL safety factor (comprise between 0 to 1). More...
 
- Protected Types inherited from Nektar::SolverUtils::EquationSystem
enum  HomogeneousType { eHomogeneous1D, eHomogeneous2D, eHomogeneous3D, eNotHomogeneous }
 Parameter for homogeneous expansions. More...
 

Detailed Description

Definition at line 52 of file CompressibleFlowSystem.h.

Constructor & Destructor Documentation

Nektar::CompressibleFlowSystem::~CompressibleFlowSystem ( )
virtual

Destructor for CompressibleFlowSystem class.

Definition at line 302 of file CompressibleFlowSystem.cpp.

303  {
304 
305  }
Nektar::CompressibleFlowSystem::CompressibleFlowSystem ( const LibUtilities::SessionReaderSharedPtr &  pSession)
protected

Definition at line 56 of file CompressibleFlowSystem.cpp.

58  : UnsteadySystem(pSession)
59  {
60  }
Nektar::SolverUtils::UnsteadySystem::UnsteadySystem
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr &pSession)
Initialises UnsteadySystem class members.
Definition: UnsteadySystem.cpp:67

Member Function Documentation

bool Nektar::CompressibleFlowSystem::CalcSteadyState ( bool  output)
protected
static SolverUtils::EquationSystemSharedPtr Nektar::CompressibleFlowSystem::create ( const LibUtilities::SessionReaderSharedPtr &  pSession)
inlinestatic

Creates an instance of this class.

Definition at line 59 of file CompressibleFlowSystem.h.

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

61  {
62  return MemoryManager<CompressibleFlowSystem>::
63  AllocateSharedPtr(pSession);
64  }
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
void Nektar::CompressibleFlowSystem::ExtrapOrder0BC ( int  bcRegion,
int  cnt,
Array< OneD, Array< OneD, NekDouble > > &  Fwd,
Array< OneD, Array< OneD, NekDouble > > &  physarray 
)
protected

Extrapolation of order 0 for all the variables such that, at the boundaries, a trivial Riemann problem is solved.

Definition at line 1513 of file CompressibleFlowSystem.cpp.

References Nektar::SolverUtils::EquationSystem::GetPhys_Offset(), Nektar::SolverUtils::EquationSystem::m_fields, and Nektar::SolverUtils::EquationSystem::m_spacedim.

Referenced by SetCommonBC().

1518  {
1519  int i, j;
1520  int e, pnt;
1521  int id1, id2, nBCEdgePts;
1522  int nVariables = physarray.num_elements();
1523  int nDimensions = m_spacedim;
1524 
1525  const Array<OneD, const int> &traceBndMap
1526  = m_fields[0]->GetTraceBndMap();
1527 
1528  int eMax;
1529 
1530  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
1531 
1532  // Loop on bcRegions
1533  for (e = 0; e < eMax; ++e)
1534  {
1535  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1536  GetExp(e)->GetTotPoints();
1537  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1538  GetPhys_Offset(e) ;
1539  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1540 
1541  // Loop on points of bcRegion 'e'
1542  for (i = 0; i < nBCEdgePts; i++)
1543  {
1544  pnt = id2+i;
1545 
1546  // Setting up bcs for density
1547  (m_fields[0]->GetBndCondExpansions()[bcRegion]->
1548  UpdatePhys())[id1+i] = Fwd[0][pnt];
1549 
1550  // Setting up bcs for velocities
1551  for (j = 1; j <=nDimensions; ++j)
1552  {
1553  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1554  UpdatePhys())[id1+i] = Fwd[j][pnt];
1555  }
1556 
1557  // Setting up bcs for energy
1558  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1559  UpdatePhys())[id1+i] = Fwd[nVariables-1][pnt];
1560  }
1561  }
1562  }
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::SolverUtils::EquationSystem::GetPhys_Offset
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
Definition: EquationSystem.h:856
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
void Nektar::CompressibleFlowSystem::GetAbsoluteVelocity ( const Array< OneD, const Array< OneD, NekDouble > > &  inarray,
Array< OneD, NekDouble > &  Vtot 
)
protected

Referenced by Nektar::EulerADCFE::DoOdeRhs().

void Nektar::CompressibleFlowSystem::GetArtificialDynamicViscosity ( const Array< OneD, Array< OneD, NekDouble > > &  physfield,
Array< OneD, NekDouble > &  mu_var 
)
protected

Referenced by Nektar::EulerADCFE::v_InitObject(), and v_InitObject().

void Nektar::CompressibleFlowSystem::GetDynamicViscosity ( const Array< OneD, const NekDouble > &  temperature,
Array< OneD, NekDouble > &  mu 
)
protected

Referenced by GetViscousFluxVector(), and GetViscousFluxVectorDeAlias().

void Nektar::CompressibleFlowSystem::GetElementDimensions ( Array< OneD, Array< OneD, NekDouble > > &  outarray,
Array< OneD, NekDouble > &  hmin 
)
protected
void Nektar::CompressibleFlowSystem::GetEnthalpy ( const Array< OneD, const Array< OneD, NekDouble > > &  physfield,
Array< OneD, NekDouble > &  pressure,
Array< OneD, NekDouble > &  enthalpy 
)
protected
void Nektar::CompressibleFlowSystem::GetEntropy ( const Array< OneD, const Array< OneD, NekDouble > > &  physfield,
const Array< OneD, const NekDouble > &  pressure,
const Array< OneD, const NekDouble > &  temperature,
Array< OneD, NekDouble > &  entropy 
)
protected
void Nektar::CompressibleFlowSystem::GetFluxVector ( const Array< OneD, Array< OneD, NekDouble > > &  physfield,
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  flux 
)
protected

Return the flux vector for the compressible Euler equations.

Parameters
iComponent of the flux vector to calculate.
physfieldFields.
fluxResulting flux.

Definition at line 1571 of file CompressibleFlowSystem.cpp.

References Nektar::SolverUtils::EquationSystem::GetPressure(), GetVelocityVector(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vmul(), and Vmath::Zero().

Referenced by v_InitObject().

1574  {
1575  int i, j;
1576  int nq = physfield[0].num_elements();
1577  int nVariables = m_fields.num_elements();
1578 
1579  Array<OneD, NekDouble> pressure(nq);
1580  Array<OneD, Array<OneD, NekDouble> > velocity(m_spacedim);
1581 
1582  // Flux vector for the rho equation
1583  for (i = 0; i < m_spacedim; ++i)
1584  {
1585  velocity[i] = Array<OneD, NekDouble>(nq);
1586  Vmath::Vcopy(nq, physfield[i+1], 1, flux[0][i], 1);
1587  }
1588 
1589  GetVelocityVector(physfield, velocity);
1590  GetPressure (physfield, velocity, pressure);
1591 
1592  // Flux vector for the velocity fields
1593  for (i = 0; i < m_spacedim; ++i)
1594  {
1595  for (j = 0; j < m_spacedim; ++j)
1596  {
1597  Vmath::Vmul(nq, velocity[j], 1, physfield[i+1], 1,
1598  flux[i+1][j], 1);
1599  }
1600 
1601  // Add pressure to appropriate field
1602  Vmath::Vadd(nq, flux[i+1][i], 1, pressure, 1, flux[i+1][i], 1);
1603  }
1604 
1605  // Flux vector for energy.
1606  Vmath::Vadd(nq, physfield[m_spacedim+1], 1, pressure, 1,
1607  pressure, 1);
1608 
1609  for (j = 0; j < m_spacedim; ++j)
1610  {
1611  Vmath::Vmul(nq, velocity[j], 1, pressure, 1,
1612  flux[m_spacedim+1][j], 1);
1613  }
1614 
1615  // For the smooth viscosity model
1616  if (nVariables == m_spacedim+3)
1617  {
1618  // Add a zero row for the advective fluxes
1619  for (j = 0; j < m_spacedim; ++j)
1620  {
1621  Vmath::Zero(nq, flux[m_spacedim+2][j], 1);
1622  }
1623  }
1624  }
Nektar::CompressibleFlowSystem::GetVelocityVector
void GetVelocityVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Nektar::SolverUtils::EquationSystem::GetPressure
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure()
Get pressure field if available.
Definition: EquationSystem.h:748
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Vmath::Vadd
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
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
void Nektar::CompressibleFlowSystem::GetFluxVectorDeAlias ( const Array< OneD, Array< OneD, NekDouble > > &  physfield,
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  flux 
)
protected

Return the flux vector for the compressible Euler equations by using the de-aliasing technique.

Parameters
iComponent of the flux vector to calculate.
physfieldFields.
fluxResulting flux.

Definition at line 1634 of file CompressibleFlowSystem.cpp.

References Nektar::SolverUtils::EquationSystem::GetPressure(), GetVelocityVector(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Vadd(), and Vmath::Vmul().

Referenced by v_InitObject().

1637  {
1638  int i, j;
1639  int nq = physfield[0].num_elements();
1640  int nVariables = m_fields.num_elements();
1641 
1642  // Factor to rescale 1d points in dealiasing
1643  NekDouble OneDptscale = 2;
1644  nq = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
1645 
1646  Array<OneD, NekDouble> pressure(nq);
1647  Array<OneD, Array<OneD, NekDouble> > velocity(m_spacedim);
1648 
1649  Array<OneD, Array<OneD, NekDouble> > physfield_interp(nVariables);
1650  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > flux_interp(
1651  nVariables);
1652 
1653  for (i = 0; i < nVariables; ++ i)
1654  {
1655  physfield_interp[i] = Array<OneD, NekDouble>(nq);
1656  flux_interp[i] = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
1657  m_fields[0]->PhysInterp1DScaled(
1658  OneDptscale, physfield[i], physfield_interp[i]);
1659 
1660  for (j = 0; j < m_spacedim; ++j)
1661  {
1662  flux_interp[i][j] = Array<OneD, NekDouble>(nq);
1663  }
1664  }
1665 
1666  // Flux vector for the rho equation
1667  for (i = 0; i < m_spacedim; ++i)
1668  {
1669  velocity[i] = Array<OneD, NekDouble>(nq);
1670 
1671  // Galerkin project solution back to original space
1672  m_fields[0]->PhysGalerkinProjection1DScaled(
1673  OneDptscale, physfield_interp[i+1], flux[0][i]);
1674  }
1675 
1676  GetVelocityVector(physfield_interp, velocity);
1677  GetPressure (physfield_interp, velocity, pressure);
1678 
1679  // Evaluation of flux vector for the velocity fields
1680  for (i = 0; i < m_spacedim; ++i)
1681  {
1682  for (j = 0; j < m_spacedim; ++j)
1683  {
1684  Vmath::Vmul(nq, velocity[j], 1, physfield_interp[i+1], 1,
1685  flux_interp[i+1][j], 1);
1686  }
1687 
1688  // Add pressure to appropriate field
1689  Vmath::Vadd(nq, flux_interp[i+1][i], 1, pressure,1,
1690  flux_interp[i+1][i], 1);
1691  }
1692 
1693  // Galerkin project solution back to origianl space
1694  for (i = 0; i < m_spacedim; ++i)
1695  {
1696  for (j = 0; j < m_spacedim; ++j)
1697  {
1698  m_fields[0]->PhysGalerkinProjection1DScaled(
1699  OneDptscale, flux_interp[i+1][j], flux[i+1][j]);
1700  }
1701  }
1702 
1703  // Evaluation of flux vector for energy
1704  Vmath::Vadd(nq, physfield_interp[m_spacedim+1], 1, pressure, 1,
1705  pressure, 1);
1706 
1707  for (j = 0; j < m_spacedim; ++j)
1708  {
1709  Vmath::Vmul(nq, velocity[j], 1, pressure, 1,
1710  flux_interp[m_spacedim+1][j], 1);
1711 
1712  // Galerkin project solution back to origianl space
1713  m_fields[0]->PhysGalerkinProjection1DScaled(
1714  OneDptscale,
1715  flux_interp[m_spacedim+1][j],
1716  flux[m_spacedim+1][j]);
1717  }
1718  }
Nektar::CompressibleFlowSystem::GetVelocityVector
void GetVelocityVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Nektar::SolverUtils::EquationSystem::GetPressure
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure()
Get pressure field if available.
Definition: EquationSystem.h:748
Vmath::Vadd
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
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
void Nektar::CompressibleFlowSystem::GetFluxVectorPDESC ( const Array< OneD, Array< OneD, NekDouble > > &  physfield,
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  flux 
)
protected
void Nektar::CompressibleFlowSystem::GetForcingTerm ( const Array< OneD, const Array< OneD, NekDouble > > &  inarray,
Array< OneD, Array< OneD, NekDouble > >  outarrayForcing 
)
protected

Referenced by Nektar::EulerADCFE::DoOdeRhs().

NekDouble Nektar::CompressibleFlowSystem::GetGamma ( )
inlineprotected

Definition at line 289 of file CompressibleFlowSystem.h.

References m_gamma.

Referenced by v_InitObject().

290  {
291  return m_gamma;
292  }
NekDouble Nektar::CompressibleFlowSystem::GetGasConstant ( )
inlineprotected

Definition at line 284 of file CompressibleFlowSystem.h.

References m_gasConstant.

285  {
286  return m_gasConstant;
287  }
void Nektar::CompressibleFlowSystem::GetMach ( Array< OneD, Array< OneD, NekDouble > > &  physfield,
Array< OneD, NekDouble > &  soundspeed,
Array< OneD, NekDouble > &  mach 
)
protected
const Array<OneD, const Array<OneD, NekDouble> >& Nektar::CompressibleFlowSystem::GetNormals ( )
inlineprotected

Definition at line 299 of file CompressibleFlowSystem.h.

References Nektar::SolverUtils::EquationSystem::m_traceNormals.

Referenced by v_InitObject().

300  {
301  return m_traceNormals;
302  }
Nektar::SolverUtils::EquationSystem::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
Definition: EquationSystem.h:516
void Nektar::CompressibleFlowSystem::GetPressure ( const Array< OneD, const Array< OneD, NekDouble > > &  physfield,
Array< OneD, NekDouble > &  pressure 
)
protected

Calculate the pressure field $ p = (\gamma-1)(E-\frac{1}{2}\rho\| \mathbf{v} \|^2) $ assuming an ideal gas law.

Parameters
physfieldInput momentum.
pressureComputed pressure field.

Definition at line 2483 of file CompressibleFlowSystem.cpp.

References m_gamma, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Smul(), Vmath::Svtvp(), Vmath::Vdiv(), Vmath::Vmul(), and Vmath::Vvtvp().

2486  {
2487  int nBCEdgePts = physfield[0].num_elements();
2488  NekDouble alpha = -0.5;
2489 
2490  // Calculate ||rho v||^2
2491  Vmath::Vmul(nBCEdgePts, physfield[1], 1, physfield[1], 1, pressure, 1);
2492  for (int i = 1; i < m_spacedim; ++i)
2493  {
2494  Vmath::Vvtvp(nBCEdgePts, physfield[1+i], 1, physfield[1+i], 1,
2495  pressure, 1, pressure, 1);
2496  }
2497  // Divide by rho to get rho*||v||^2
2498  Vmath::Vdiv (nBCEdgePts, pressure, 1, physfield[0], 1, pressure, 1);
2499  // pressure <- E - 0.5*pressure
2500  Vmath::Svtvp(nBCEdgePts, alpha,
2501  pressure, 1, physfield[m_spacedim+1], 1, pressure, 1);
2502  // Multiply by (gamma-1)
2503  Vmath::Smul (nBCEdgePts, m_gamma-1, pressure, 1, pressure, 1);
2504  }
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:471
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Vmath::Vdiv
void Vdiv(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:227
Vmath::Smul
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
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
void Nektar::CompressibleFlowSystem::GetPressure ( const Array< OneD, const Array< OneD, NekDouble > > &  physfield,
const Array< OneD, const Array< OneD, NekDouble > > &  velocity,
Array< OneD, NekDouble > &  pressure 
)
protected
void Nektar::CompressibleFlowSystem::GetSensor ( const Array< OneD, const Array< OneD, NekDouble > > &  physarray,
Array< OneD, NekDouble > &  Sensor,
Array< OneD, NekDouble > &  SensorKappa 
)
protected
void Nektar::CompressibleFlowSystem::GetSmoothArtificialViscosity ( const Array< OneD, Array< OneD, NekDouble > > &  physfield,
Array< OneD, NekDouble > &  eps_bar 
)
protected

Referenced by v_InitObject().

void Nektar::CompressibleFlowSystem::GetSoundSpeed ( const Array< OneD, Array< OneD, NekDouble > > &  physfield,
Array< OneD, NekDouble > &  pressure,
Array< OneD, NekDouble > &  soundspeed 
)
protected

Referenced by Nektar::EulerADCFE::DoOdeRhs().

NekDouble Nektar::CompressibleFlowSystem::GetStabilityLimit ( int  n)

Function to calculate the stability limit for DG/CG.

Array<OneD, NekDouble> Nektar::CompressibleFlowSystem::GetStabilityLimitVector ( const Array< OneD, int > &  ExpOrder)

Function to calculate the stability limit for DG/CG (a vector of them).

void Nektar::CompressibleFlowSystem::GetStdVelocity ( const Array< OneD, const Array< OneD, NekDouble > > &  inarray,
Array< OneD, NekDouble > &  stdV 
)
protected
void Nektar::CompressibleFlowSystem::GetTemperature ( const Array< OneD, const Array< OneD, NekDouble > > &  physfield,
Array< OneD, NekDouble > &  pressure,
Array< OneD, NekDouble > &  temperature 
)
protected

Referenced by Nektar::NavierStokesCFE::DoOdeRhs(), and GetViscousFluxVectorDeAlias().

const Array<OneD, const Array<OneD, NekDouble> >& Nektar::CompressibleFlowSystem::GetVecLocs ( )
inlineprotected

Definition at line 294 of file CompressibleFlowSystem.h.

References m_vecLocs.

Referenced by v_InitObject().

295  {
296  return m_vecLocs;
297  }
Nektar::CompressibleFlowSystem::m_vecLocs
Array< OneD, Array< OneD, NekDouble > > m_vecLocs
Definition: CompressibleFlowSystem.h:83
void Nektar::CompressibleFlowSystem::GetVelocityVector ( const Array< OneD, Array< OneD, NekDouble > > &  physfield,
Array< OneD, Array< OneD, NekDouble > > &  velocity 
)
protected

Referenced by GetFluxVector(), and GetFluxVectorDeAlias().

void Nektar::CompressibleFlowSystem::GetViscousFluxVector ( const Array< OneD, Array< OneD, NekDouble > > &  physfield,
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  derivativesO1,
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  viscousTensor 
)
protected

Return the flux vector for the LDG diffusion problem.

Todo:
Complete the viscous flux vector

Definition at line 1725 of file CompressibleFlowSystem.cpp.

References ASSERTL0, Vmath::Fill(), GetDynamicViscosity(), m_Cp, Nektar::SolverUtils::EquationSystem::m_fields, m_mu, m_Prandtl, Nektar::SolverUtils::EquationSystem::m_spacedim, m_thermalConductivity, m_ViscosityType, Vmath::Smul(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vmul(), and Vmath::Zero().

Referenced by v_InitObject().

1729  {
1730  int j, k;
1731  int nVariables = m_fields.num_elements();
1732  int nPts = physfield[0].num_elements();
1733 
1734  // Stokes hypotesis
1735  const NekDouble lambda = -2.0/3.0;
1736 
1737  // Auxiliary variables
1738  Array<OneD, NekDouble > mu (nPts, 0.0);
1739  Array<OneD, NekDouble > thermalConductivity(nPts, 0.0);
1740  Array<OneD, NekDouble > mu2 (nPts, 0.0);
1741  Array<OneD, NekDouble > divVel (nPts, 0.0);
1742 
1743  // Variable viscosity through the Sutherland's law
1744  if (m_ViscosityType == "Variable")
1745  {
1746  GetDynamicViscosity(physfield[nVariables-2], mu);
1747  NekDouble tRa = m_Cp / m_Prandtl;
1748  Vmath::Smul(nPts, tRa, &mu[0], 1, &thermalConductivity[0], 1);
1749  }
1750  else
1751  {
1752  Vmath::Fill(nPts, m_mu, &mu[0], 1);
1753  Vmath::Fill(nPts, m_thermalConductivity,
1754  &thermalConductivity[0], 1);
1755  }
1756 
1757  // Computing diagonal terms of viscous stress tensor
1758  Array<OneD, Array<OneD, NekDouble> > tmp(m_spacedim);
1759  Array<OneD, Array<OneD, NekDouble> > Sgg(m_spacedim);
1760 
1761  // mu2 = 2 * mu
1762  Vmath::Smul(nPts, 2.0, &mu[0], 1, &mu2[0], 1);
1763 
1764  // Velocity divergence
1765  for (j = 0; j < m_spacedim; ++j)
1766  {
1767  Vmath::Vadd(nPts, &divVel[0], 1, &derivativesO1[j][j][0], 1,
1768  &divVel[0], 1);
1769  }
1770 
1771  // Velocity divergence scaled by lambda * mu
1772  Vmath::Smul(nPts, lambda, &divVel[0], 1, &divVel[0], 1);
1773  Vmath::Vmul(nPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
1774 
1775  // Diagonal terms of viscous stress tensor (Sxx, Syy, Szz)
1776  // Sjj = 2 * mu * du_j/dx_j - (2 / 3) * mu * sum_j(du_j/dx_j)
1777  for (j = 0; j < m_spacedim; ++j)
1778  {
1779  tmp[j] = Array<OneD, NekDouble>(nPts, 0.0);
1780  Sgg[j] = Array<OneD, NekDouble>(nPts, 0.0);
1781 
1782  Vmath::Vmul(nPts, &mu2[0], 1, &derivativesO1[j][j][0], 1,
1783  &tmp[j][0], 1);
1784 
1785  Vmath::Vadd(nPts, &tmp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
1786  }
1787 
1788  // Extra diagonal terms of viscous stress tensor (Sxy, Sxz, Syz)
1789  // Note: they exist for 2D and 3D problems only
1790  Array<OneD, NekDouble > Sxy(nPts, 0.0);
1791  Array<OneD, NekDouble > Sxz(nPts, 0.0);
1792  Array<OneD, NekDouble > Syz(nPts, 0.0);
1793 
1794  if (m_spacedim == 2)
1795  {
1796  // Sxy = (du/dy + dv/dx)
1797  Vmath::Vadd(nPts, &derivativesO1[0][1][0], 1,
1798  &derivativesO1[1][0][0], 1, &Sxy[0], 1);
1799 
1800  // Sxy = mu * (du/dy + dv/dx)
1801  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
1802  }
1803  else if (m_spacedim == 3)
1804  {
1805  // Sxy = (du/dy + dv/dx)
1806  Vmath::Vadd(nPts, &derivativesO1[0][1][0], 1,
1807  &derivativesO1[1][0][0], 1, &Sxy[0], 1);
1808 
1809  // Sxz = (du/dz + dw/dx)
1810  Vmath::Vadd(nPts, &derivativesO1[0][2][0], 1,
1811  &derivativesO1[2][0][0], 1, &Sxz[0], 1);
1812 
1813  // Syz = (dv/dz + dw/dy)
1814  Vmath::Vadd(nPts, &derivativesO1[1][2][0], 1,
1815  &derivativesO1[2][1][0], 1, &Syz[0], 1);
1816 
1817  // Sxy = mu * (du/dy + dv/dx)
1818  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
1819 
1820  // Sxz = mu * (du/dy + dv/dx)
1821  Vmath::Vmul(nPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
1822 
1823  // Syz = mu * (du/dy + dv/dx)
1824  Vmath::Vmul(nPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
1825  }
1826 
1827  // Energy-related terms
1828  Array<OneD, NekDouble > STx(nPts, 0.0);
1829  Array<OneD, NekDouble > STy(nPts, 0.0);
1830  Array<OneD, NekDouble > STz(nPts, 0.0);
1831 
1832  // Building the viscous flux vector
1833 
1834  // Viscous flux vector for the rho equation
1835  for (k = 0; k < m_spacedim; ++k)
1836  {
1837  Vmath::Zero(nPts, viscousTensor[k][0], 1);
1838  }
1839 
1840  if (m_spacedim == 1)
1841  {
1842  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1843 
1844  // u * Sxx
1845  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
1846 
1847  // k * dT/dx
1848  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1849  &derivativesO1[0][1][0], 1, &tmp1[0], 1);
1850 
1851  // STx = u * Sxx + (K / mu) * dT/dx
1852  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1853  }
1854  else if (m_spacedim == 2)
1855  {
1856  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1857  Array<OneD, NekDouble > tmp2(nPts, 0.0);
1858 
1859  // Computation of STx
1860 
1861  // u * Sxx
1862  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
1863 
1864  // v * Sxy
1865  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1866 
1867  // k * dT/dx
1868  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1869  &derivativesO1[0][2][0], 1, &tmp2[0], 1);
1870 
1871  // STx = u * Sxx + v * Sxy + K * dT/dx
1872  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1873  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
1874 
1875  // Computation of STy
1876 
1877  // Re-initialise temporary arrays
1878  Vmath::Zero(nPts, &tmp1[0], 1);
1879  Vmath::Zero(nPts, &tmp2[0], 1);
1880 
1881  // v * Syy
1882  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sgg[1][0], 1, &STy[0], 1);
1883 
1884  // u * Sxy
1885  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1886 
1887  // k * dT/dy
1888  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1889  &derivativesO1[1][2][0], 1, &tmp2[0], 1);
1890 
1891  // STy = v * Syy + u * Sxy + K * dT/dy
1892  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
1893  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
1894  }
1895  else if (m_spacedim == 3)
1896  {
1897  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1898  Array<OneD, NekDouble > tmp2(nPts, 0.0);
1899  Array<OneD, NekDouble > tmp3(nPts, 0.0);
1900 
1901  // Computation of STx
1902 
1903  // u * Sxx
1904  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
1905 
1906  // v * Sxy
1907  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1908 
1909  // v * Sxz
1910  Vmath::Vmul(nPts, &physfield[2][0], 1, &Sxz[0], 1, &tmp2[0], 1);
1911 
1912  // k * dT/dx
1913  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1914  &derivativesO1[0][3][0], 1, &tmp3[0], 1);
1915 
1916  // STx = u * Sxx + v * Sxy + w * Sxz + (K / mu) * dT/dx
1917  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1918  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
1919  Vmath::Vadd(nPts, &STx[0], 1, &tmp3[0], 1, &STx[0], 1);
1920 
1921  // Computation of STy
1922 
1923  // Re-initialise temporary arrays
1924  Vmath::Zero(nPts, &tmp1[0], 1);
1925  Vmath::Zero(nPts, &tmp2[0], 1);
1926  Vmath::Zero(nPts, &tmp3[0], 1);
1927 
1928  // v * Syy
1929  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sgg[1][0], 1, &STy[0], 1);
1930 
1931  // u * Sxy
1932  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1933 
1934  // w * Syz
1935  Vmath::Vmul(nPts, &physfield[2][0], 1, &Syz[0], 1, &tmp2[0], 1);
1936 
1937  // k * dT/dy
1938  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1939  &derivativesO1[1][3][0], 1, &tmp3[0], 1);
1940 
1941  // STy = v * Syy + u * Sxy + w * Syz + K * dT/dy
1942  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
1943  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
1944  Vmath::Vadd(nPts, &STy[0], 1, &tmp3[0], 1, &STy[0], 1);
1945 
1946  // Computation of STz
1947 
1948  // Re-initialise temporary arrays
1949  Vmath::Zero(nPts, &tmp1[0], 1);
1950  Vmath::Zero(nPts, &tmp2[0], 1);
1951  Vmath::Zero(nPts, &tmp3[0], 1);
1952 
1953  // w * Szz
1954  Vmath::Vmul(nPts, &physfield[2][0], 1, &Sgg[2][0], 1, &STz[0], 1);
1955 
1956  // u * Sxz
1957  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxz[0], 1, &tmp1[0], 1);
1958 
1959  // v * Syz
1960  Vmath::Vmul(nPts, &physfield[1][0], 1, &Syz[0], 1, &tmp2[0], 1);
1961 
1962  // k * dT/dz
1963  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1964  &derivativesO1[2][3][0], 1, &tmp3[0], 1);
1965 
1966  // STz = w * Szz + u * Sxz + v * Syz + K * dT/dz
1967  Vmath::Vadd(nPts, &STz[0], 1, &tmp1[0], 1, &STz[0], 1);
1968  Vmath::Vadd(nPts, &STz[0], 1, &tmp2[0], 1, &STz[0], 1);
1969  Vmath::Vadd(nPts, &STz[0], 1, &tmp3[0], 1, &STz[0], 1);
1970  }
1971 
1972  switch (m_spacedim)
1973  {
1974  case 1:
1975  {
1976  // f_11v = f_rho = 0
1977  Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
1978 
1979  // f_21v = f_rhou
1980  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
1981 
1982  // f_31v = f_E
1983  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][2][0], 1);
1984  break;
1985  }
1986  case 2:
1987  {
1988  // f_11v = f_rho1 = 0
1989  Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
1990  // f_12v = f_rho2 = 0
1991  Vmath::Zero(nPts, &viscousTensor[1][0][0], 1);
1992 
1993  // f_21v = f_rhou1
1994  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
1995  // f_22v = f_rhou2
1996  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[1][1][0], 1);
1997 
1998  // f_31v = f_rhov1
1999  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[0][2][0], 1);
2000  // f_32v = f_rhov2
2001  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor[1][2][0], 1);
2002 
2003  // f_41v = f_E1
2004  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][3][0], 1);
2005  // f_42v = f_E2
2006  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor[1][3][0], 1);
2007  break;
2008  }
2009  case 3:
2010  {
2011  // f_11v = f_rho1 = 0
2012  Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
2013  // f_12v = f_rho2 = 0
2014  Vmath::Zero(nPts, &viscousTensor[1][0][0], 1);
2015  // f_13v = f_rho3 = 0
2016  Vmath::Zero(nPts, &viscousTensor[2][0][0], 1);
2017 
2018  // f_21v = f_rhou1
2019  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
2020  // f_22v = f_rhou2
2021  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[1][1][0], 1);
2022  // f_23v = f_rhou3
2023  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor[2][1][0], 1);
2024 
2025  // f_31v = f_rhov1
2026  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[0][2][0], 1);
2027  // f_32v = f_rhov2
2028  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor[1][2][0], 1);
2029  // f_33v = f_rhov3
2030  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor[2][2][0], 1);
2031 
2032  // f_31v = f_rhow1
2033  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor[0][3][0], 1);
2034  // f_32v = f_rhow2
2035  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor[1][3][0], 1);
2036  // f_33v = f_rhow3
2037  Vmath::Vcopy(nPts, &Sgg[2][0], 1, &viscousTensor[2][3][0], 1);
2038 
2039  // f_41v = f_E1
2040  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][4][0], 1);
2041  // f_42v = f_E2
2042  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor[1][4][0], 1);
2043  // f_43v = f_E3
2044  Vmath::Vcopy(nPts, &STz[0], 1, &viscousTensor[2][4][0], 1);
2045  break;
2046  }
2047  default:
2048  {
2049  ASSERTL0(false, "Illegal expansion dimension");
2050  }
2051  }
2052  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
Vmath::Smul
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
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Nektar::CompressibleFlowSystem::GetDynamicViscosity
void GetDynamicViscosity(const Array< OneD, const NekDouble > &temperature, Array< OneD, NekDouble > &mu)
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Vmath::Vadd
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
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
void Nektar::CompressibleFlowSystem::GetViscousFluxVectorDeAlias ( const Array< OneD, Array< OneD, NekDouble > > &  physfield,
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  derivativesO1,
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  viscousTensor 
)
protected

Return the flux vector for the LDG diffusion problem.

Todo:
Complete the viscous flux vector

Definition at line 2058 of file CompressibleFlowSystem.cpp.

References ASSERTL0, GetDynamicViscosity(), Nektar::SolverUtils::EquationSystem::GetPressure(), GetTemperature(), Nektar::SolverUtils::EquationSystem::m_fields, m_mu, Nektar::SolverUtils::EquationSystem::m_spacedim, m_thermalConductivity, m_ViscosityType, Vmath::Sadd(), Vmath::Smul(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vmul(), and Vmath::Zero().

Referenced by v_InitObject().

2062  {
2063 #if 0
2064  int i, j, k;
2065  int nVariables = m_fields.num_elements();
2066  int nPts = physfield[0].num_elements();
2067 
2068  int variables_phys = physfield.num_elements();
2069 
2070  // Factor to rescale 1d points in dealiasing.
2071  NekDouble OneDptscale = 2;
2072 
2073  // Get number of points to dealias a cubic non-linearity
2074  nPts = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
2075 
2076  int nVariables_aux = derivativesO1[0].num_elements();
2077 
2078  Array<OneD, Array<OneD, NekDouble> > physfield_interp(variables_phys);
2079  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > derivativesO1_interp(
2080  m_spacedim);
2081  Array<OneD, Array<OneD, Array<OneD, NekDouble> > >viscousTensor_interp(
2082  m_spacedim);
2083 
2084  for (i = 0; i < m_spacedim; ++ i)
2085  {
2086  viscousTensor_interp[i] = Array<OneD, Array<OneD, NekDouble> >(
2087  nVariables);
2088  for (j = 0; j < nVariables; ++j)
2089  {
2090  viscousTensor_interp[i][j] = Array<OneD, NekDouble>(nPts);
2091  }
2092  }
2093 
2094  // Stokes hypotesis
2095  NekDouble lambda = -2.0/3.0;
2096 
2097  // Auxiliary variables
2098  Array<OneD, NekDouble > mu (nPts, 0.0);
2099  Array<OneD, NekDouble > mu2 (nPts, 0.0);
2100  Array<OneD, NekDouble > divVel (nPts, 0.0);
2101  Array<OneD, NekDouble > pressure (nPts, 0.0);
2102  Array<OneD, NekDouble > temperature(nPts, 0.0);
2103 
2104  for (i = 0; i < nVariables; ++i)
2105  {
2106  m_fields[0]->PhysInterp1DScaled(
2107  OneDptscale, physfield[i], fields_interp[i]);
2108  }
2109 
2110  for (i = 0; i < variables_phys; ++i)
2111  {
2112  physfield_interp[i] = Array<OneD, NekDouble> (nPts);
2113 
2114  // Interpolation to higher space
2115  m_fields[0]->PhysInterp1DScaled(OneDptscale, physfield[i],
2116  physfield_interp[i]);
2117  }
2118 
2119  for (i = 0; i < m_spacedim; ++i)
2120  {
2121  derivativesO1_interp[i] = Array<OneD, Array<OneD, NekDouble> >(
2122  nVariables_aux);
2123  for (j = 0; j < nVariables_aux; ++j)
2124  {
2125  derivativesO1_interp[i][j] = Array<OneD, NekDouble>(nPts);
2126  m_fields[0]->PhysInterp1DScaled(
2127  OneDptscale, derivativesO1[i][j],
2128  derivativesO1_interp[i][j]);
2129  }
2130  }
2131 
2132  // Thermodynamic related quantities
2133  GetPressure(fields_interp, pressure);
2134  GetTemperature(fields_interp, pressure, temperature);
2135 
2136  // Variable viscosity through the Sutherland's law
2137  if (m_ViscosityType == "Variable")
2138  {
2139  GetDynamicViscosity(fields_interp[variables_phys-1], mu);
2140  }
2141  else
2142  {
2143  Vmath::Sadd(nPts, m_mu, &mu[0], 1, &mu[0], 1);
2144  }
2145 
2146  // Computing diagonal terms of viscous stress tensor
2147  Array<OneD, Array<OneD, NekDouble> > tmp(m_spacedim);
2148  Array<OneD, Array<OneD, NekDouble> > Sgg(m_spacedim);
2149 
2150  // mu2 = 2 * mu
2151  Vmath::Smul(nPts, 2.0, &mu[0], 1, &mu2[0], 1);
2152 
2153  // Velocity divergence
2154  for (j = 0; j < m_spacedim; ++j)
2155  {
2156  Vmath::Vadd(nPts, &divVel[0], 1, &derivativesO1_interp[j][j][0], 1,
2157  &divVel[0], 1);
2158  }
2159 
2160  // Velocity divergence scaled by lambda * mu
2161  Vmath::Smul(nPts, lambda, &divVel[0], 1, &divVel[0], 1);
2162  Vmath::Vmul(nPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
2163 
2164  // Digonal terms of viscous stress tensor (Sxx, Syy, Szz)
2165  // Sjj = 2 * mu * du_j/dx_j - (2 / 3) * mu * sum_j(du_j/dx_j)
2166  for (j = 0; j < m_spacedim; ++j)
2167  {
2168  tmp[j] = Array<OneD, NekDouble>(nPts, 0.0);
2169  Sgg[j] = Array<OneD, NekDouble>(nPts, 0.0);
2170 
2171  Vmath::Vmul(nPts, &mu2[0], 1, &derivativesO1_interp[j][j][0], 1,
2172  &tmp[j][0], 1);
2173 
2174  Vmath::Vadd(nPts, &tmp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
2175  }
2176 
2177  // Extra diagonal terms of viscous stress tensor (Sxy, Sxz, Syz)
2178  // Note: they exist for 2D and 3D problems only
2179  Array<OneD, NekDouble > Sxy(nPts, 0.0);
2180  Array<OneD, NekDouble > Sxz(nPts, 0.0);
2181  Array<OneD, NekDouble > Syz(nPts, 0.0);
2182 
2183  if (m_spacedim == 2)
2184  {
2185  // Sxy = (du/dy + dv/dx)
2186  Vmath::Vadd(nPts, &derivativesO1_interp[0][1][0], 1,
2187  &derivativesO1_interp[1][0][0], 1, &Sxy[0], 1);
2188 
2189  // Sxy = mu * (du/dy + dv/dx)
2190  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
2191  }
2192  else if (m_spacedim == 3)
2193  {
2194  // Sxy = (du/dy + dv/dx)
2195  Vmath::Vadd(nPts, &derivativesO1_interp[0][1][0], 1,
2196  &derivativesO1_interp[1][0][0], 1, &Sxy[0], 1);
2197 
2198  // Sxz = (du/dz + dw/dx)
2199  Vmath::Vadd(nPts, &derivativesO1_interp[0][2][0], 1,
2200  &derivativesO1_interp[2][0][0], 1, &Sxz[0], 1);
2201 
2202  // Syz = (dv/dz + dw/dy)
2203  Vmath::Vadd(nPts, &derivativesO1_interp[1][2][0], 1,
2204  &derivativesO1_interp[2][1][0], 1, &Syz[0], 1);
2205 
2206  // Sxy = mu * (du/dy + dv/dx)
2207  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
2208 
2209  // Sxz = mu * (du/dy + dv/dx)
2210  Vmath::Vmul(nPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
2211 
2212  // Syz = mu * (du/dy + dv/dx)
2213  Vmath::Vmul(nPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
2214  }
2215 
2216  // Energy-related terms
2217  Array<OneD, NekDouble > STx(nPts, 0.0);
2218  Array<OneD, NekDouble > STy(nPts, 0.0);
2219  Array<OneD, NekDouble > STz(nPts, 0.0);
2220  // Building the viscous flux vector
2221  if (i == 0)
2222  {
2223  // Viscous flux vector for the rho equation
2224  for (k = 0; k < m_spacedim; ++k)
2225  {
2226  Vmath::Zero(nPts, viscousTensor_interp[k][i], 1);
2227  }
2228  }
2229 
2230  if (m_spacedim == 1)
2231  {
2232  Array<OneD, NekDouble > tmp1(nPts, 0.0);
2233 
2234  // u * Sxx
2235  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2236  &Sgg[0][0], 1, &STx[0], 1);
2237 
2238  // k * dT/dx
2239  Vmath::Smul(nPts, m_thermalConductivity,
2240  &derivativesO1_interp[0][1][0], 1, &tmp1[0], 1);
2241 
2242  // STx = u * Sxx + (K / mu) * dT/dx
2243  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
2244  }
2245  else if (m_spacedim == 2)
2246  {
2247  Array<OneD, NekDouble > tmp1(nPts, 0.0);
2248  Array<OneD, NekDouble > tmp2(nPts, 0.0);
2249 
2250  // Computation of STx
2251 
2252  // u * Sxx
2253  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2254  &Sgg[0][0], 1, &STx[0], 1);
2255 
2256  // v * Sxy
2257  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2258  &Sxy[0], 1, &tmp1[0], 1);
2259 
2260  // k * dT/dx
2261  Vmath::Smul(nPts, m_thermalConductivity,
2262  &derivativesO1_interp[0][2][0], 1, &tmp2[0], 1);
2263 
2264  // STx = u * Sxx + v * Sxy + K * dT/dx
2265  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
2266  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
2267 
2268  // Computation of STy
2269 
2270  // Re-initialise temporary arrays
2271  Vmath::Zero(nPts, &tmp1[0], 1);
2272  Vmath::Zero(nPts, &tmp2[0], 1);
2273 
2274  // v * Syy
2275  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2276  &Sgg[1][0], 1, &STy[0], 1);
2277 
2278  // u * Sxy
2279  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2280  &Sxy[0], 1, &tmp1[0], 1);
2281 
2282  // k * dT/dy
2283  Vmath::Smul(nPts, m_thermalConductivity,
2284  &derivativesO1_interp[1][2][0], 1, &tmp2[0], 1);
2285 
2286  // STy = v * Syy + u * Sxy + K * dT/dy
2287  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
2288  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
2289  }
2290  else if (m_spacedim == 3)
2291  {
2292  Array<OneD, NekDouble > tmp1(nPts, 0.0);
2293  Array<OneD, NekDouble > tmp2(nPts, 0.0);
2294  Array<OneD, NekDouble > tmp3(nPts, 0.0);
2295 
2296  // Computation of STx
2297 
2298  // u * Sxx
2299  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2300  &Sgg[0][0], 1, &STx[0], 1);
2301 
2302  // v * Sxy
2303  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2304  &Sxy[0], 1, &tmp1[0], 1);
2305 
2306  // v * Sxy
2307  Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
2308  &Sxz[0], 1, &tmp2[0], 1);
2309 
2310  // k * dT/dx
2311  Vmath::Smul(nPts, m_thermalConductivity,
2312  &derivativesO1_interp[0][3][0], 1, &tmp3[0], 1);
2313 
2314  // STx = u * Sxx + v * Sxy + w * Sxz + (K / mu) * dT/dx
2315  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
2316  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
2317  Vmath::Vadd(nPts, &STx[0], 1, &tmp3[0], 1, &STx[0], 1);
2318 
2319  // Computation of STy
2320 
2321  // Re-initialise temporary arrays
2322  Vmath::Zero(nPts, &tmp1[0], 1);
2323  Vmath::Zero(nPts, &tmp2[0], 1);
2324  Vmath::Zero(nPts, &tmp3[0], 1);
2325 
2326  // v * Syy
2327  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2328  &Sgg[1][0], 1, &STy[0], 1);
2329 
2330  // u * Sxy
2331  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2332  &Sxy[0], 1, &tmp1[0], 1);
2333 
2334  // w * Syz
2335  Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
2336  &Syz[0], 1, &tmp2[0], 1);
2337 
2338  // k * dT/dy
2339  Vmath::Smul(nPts, m_thermalConductivity,
2340  &derivativesO1_interp[1][3][0], 1, &tmp3[0], 1);
2341 
2342  // STy = v * Syy + u * Sxy + w * Syz + K * dT/dy
2343  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
2344  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
2345  Vmath::Vadd(nPts, &STy[0], 1, &tmp3[0], 1, &STy[0], 1);
2346 
2347  // Computation of STz
2348 
2349  // Re-initialise temporary arrays
2350  Vmath::Zero(nPts, &tmp1[0], 1);
2351  Vmath::Zero(nPts, &tmp2[0], 1);
2352  Vmath::Zero(nPts, &tmp3[0], 1);
2353 
2354  // w * Szz
2355  Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
2356  &Sgg[2][0], 1, &STz[0], 1);
2357 
2358  // u * Sxz
2359  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2360  &Sxz[0], 1, &tmp1[0], 1);
2361 
2362  // v * Syz
2363  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2364  &Syz[0], 1, &tmp2[0], 1);
2365 
2366  // k * dT/dz
2367  Vmath::Smul(nPts, m_thermalConductivity,
2368  &derivativesO1_interp[2][3][0], 1, &tmp3[0], 1);
2369 
2370  // STz = w * Szz + u * Sxz + v * Syz + K * dT/dz
2371  Vmath::Vadd(nPts, &STz[0], 1, &tmp1[0], 1, &STz[0], 1);
2372  Vmath::Vadd(nPts, &STz[0], 1, &tmp2[0], 1, &STz[0], 1);
2373  Vmath::Vadd(nPts, &STz[0], 1, &tmp3[0], 1, &STz[0], 1);
2374  }
2375 
2376  switch (m_spacedim)
2377  {
2378  case 1:
2379  {
2380 
2381  int nq = physfield[0].num_elements();
2382  // f_11v = f_rho = 0
2383  Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
2384 
2385  // f_21v = f_rhou
2386  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
2387 
2388  // f_31v = f_E
2389  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][2][0], 1);
2390  break;
2391  }
2392  case 2:
2393  {
2394  int nq = physfield[0].num_elements();
2395  // f_11v = f_rho1 = 0
2396  Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
2397  // f_12v = f_rho2 = 0
2398  Vmath::Zero(nq, &viscousTensor_interp[1][0][0], 1);
2399 
2400  // f_21v = f_rhou1
2401  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
2402  // f_22v = f_rhou2
2403  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[1][1][0], 1);
2404 
2405  // f_31v = f_rhov1
2406  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[0][2][0], 1);
2407  // f_32v = f_rhov2
2408  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor_interp[1][2][0], 1);
2409 
2410  // f_41v = f_E1
2411  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][3][0], 1);
2412  // f_42v = f_E2
2413  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor_interp[1][3][0], 1);
2414  break;
2415  }
2416  case 3:
2417  {
2418  int nq = physfield[0].num_elements();
2419  // f_11v = f_rho1 = 0
2420  Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
2421  // f_12v = f_rho2 = 0
2422  Vmath::Zero(nq, &viscousTensor_interp[1][0][0], 1);
2423  // f_13v = f_rho3 = 0
2424  Vmath::Zero(nq, &viscousTensor_interp[2][0][0], 1);
2425 
2426  // f_21v = f_rhou1
2427  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
2428  // f_22v = f_rhou2
2429  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[1][1][0], 1);
2430  // f_23v = f_rhou3
2431  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor_interp[2][1][0], 1);
2432 
2433  // f_31v = f_rhov1
2434  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[0][2][0], 1);
2435  // f_32v = f_rhov2
2436  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor_interp[1][2][0], 1);
2437  // f_33v = f_rhov3
2438  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor_interp[2][2][0], 1);
2439 
2440  // f_31v = f_rhow1
2441  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor_interp[0][3][0], 1);
2442  // f_32v = f_rhow2
2443  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor_interp[1][3][0], 1);
2444  // f_33v = f_rhow3
2445  Vmath::Vcopy(nPts, &Sgg[2][0], 1, &viscousTensor_interp[2][3][0], 1);
2446 
2447  // f_41v = f_E1
2448  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][4][0], 1);
2449  // f_42v = f_E2
2450  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor_interp[1][4][0], 1);
2451  // f_43v = f_E3
2452  Vmath::Vcopy(nPts, &STz[0], 1, &viscousTensor_interp[2][4][0], 1);
2453  break;
2454  }
2455  default:
2456  {
2457  ASSERTL0(false, "Illegal expansion dimension");
2458  }
2459  }
2460 
2461  for (i = 0; i < m_spacedim; ++i)
2462  {
2463  for (j = 1; j < nVariables; ++j)
2464  {
2465  // Galerkin project solution back to origianl space
2466  m_fields[0]->PhysGalerkinProjection1DScaled(
2467  OneDptscale,
2468  viscousTensor_interp[i][j],
2469  viscousTensor[i][j]);
2470  }
2471  }
2472 #endif
2473 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::CompressibleFlowSystem::GetTemperature
void GetTemperature(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure, Array< OneD, NekDouble > &temperature)
Vmath::Smul
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
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:301
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Nektar::CompressibleFlowSystem::GetDynamicViscosity
void GetDynamicViscosity(const Array< OneD, const NekDouble > &temperature, Array< OneD, NekDouble > &mu)
Nektar::SolverUtils::EquationSystem::GetPressure
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure()
Get pressure field if available.
Definition: EquationSystem.h:748
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Vmath::Vadd
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
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
void Nektar::CompressibleFlowSystem::PressureInflowFileBC ( int  bcRegion,
int  cnt,
Array< OneD, Array< OneD, NekDouble > > &  Fwd,
Array< OneD, Array< OneD, NekDouble > > &  physarray 
)
protected

Pressure inflow boundary conditions for compressible flow problems where either the density and the velocities are assigned from a file or the full state is assigned from a file (depending on the problem type, either subsonic or supersonic).

Definition at line 1358 of file CompressibleFlowSystem.cpp.

References Nektar::SolverUtils::EquationSystem::GetExpSize(), Nektar::SolverUtils::EquationSystem::GetPhys_Offset(), Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_fieldStorage, m_gamma, Nektar::SolverUtils::EquationSystem::m_spacedim, Nektar::SolverUtils::EquationSystem::m_traceNormals, m_uInf, m_vInf, m_wInf, npts, Vmath::Smul(), Vmath::Vabs(), Vmath::Vadd(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by SetCommonBC().

1363  {
1364  int i, j;
1365  int nTracePts = GetTraceTotPoints();
1366  int nVariables = physarray.num_elements();
1367  int nDimensions = m_spacedim;
1368 
1369  const Array<OneD, const int> &traceBndMap
1370  = m_fields[0]->GetTraceBndMap();
1371 
1372  NekDouble gamma = m_gamma;
1373  NekDouble gammaMinusOne = gamma - 1.0;
1374  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
1375 
1376  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
1377  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
1378  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
1379  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
1380 
1381  // Computing the normal velocity for characteristics coming
1382  // from outside the computational domain
1383  velInf[0] = m_uInf;
1384  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
1385  if (nDimensions == 2 || nDimensions == 3)
1386  {
1387  velInf[1] = m_vInf;
1388  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
1389  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
1390  }
1391  if (nDimensions == 3)
1392  {
1393  velInf[2] = m_wInf;
1394  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
1395  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
1396  }
1397 
1398  // Computing the normal velocity for characteristics coming
1399  // from inside the computational domain
1400  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1401  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
1402  for (i = 0; i < nDimensions; ++i)
1403  {
1404  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
1405  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
1406  }
1407 
1408  // Computing the absolute value of the velocity in order to compute the
1409  // Mach number to decide whether supersonic or subsonic
1410  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
1411  for (i = 0; i < nDimensions; ++i)
1412  {
1413  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
1414  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
1415  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
1416  }
1417  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
1418 
1419  // Get speed of sound
1420  Array<OneD, NekDouble > soundSpeed (nTracePts, 0.0);
1421  Array<OneD, NekDouble > pressure (nTracePts, 0.0);
1422 
1423  for (i = 0; i < nTracePts; i++)
1424  {
1425  if (m_spacedim == 1)
1426  {
1427  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
1428  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
1429  }
1430  else if (m_spacedim == 2)
1431  {
1432  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
1433  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1434  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
1435  }
1436  else
1437  {
1438  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
1439  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1440  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
1441  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
1442  }
1443 
1444  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
1445  }
1446 
1447  // Get Mach
1448  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
1449  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
1450  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
1451 
1452  // Auxiliary variables
1453  int e, id1, id2, npts, pnt;
1454  NekDouble rhoeb;
1455 
1456  // Loop on the bcRegions
1457  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->
1458  GetExpSize(); ++e)
1459  {
1460  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1461  GetExp(e)->GetTotPoints();
1462  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1463  GetPhys_Offset(e);
1464  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1465 
1466  // Loop on points of bcRegion 'e'
1467  for (i = 0; i < npts; i++)
1468  {
1469  pnt = id2+i;
1470 
1471  // Subsonic flows
1472  if (Mach[pnt] < 0.99)
1473  {
1474  // Partial extrapolation for subsonic cases
1475  for (j = 0; j < nVariables-1; ++j)
1476  {
1477  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1478  UpdatePhys())[id1+i] = m_fieldStorage[j][id1+i];
1479  }
1480 
1481  // Kinetic energy calculation
1482  NekDouble Ek = 0.0;
1483  for (j = 1; j < nVariables-1; ++j)
1484  {
1485  Ek += 0.5 * (m_fieldStorage[j][id1+i] *
1486  m_fieldStorage[j][id1+i]) /
1487  m_fieldStorage[0][id1+i];
1488  }
1489  rhoeb = gammaMinusOneInv * pressure[pnt] + Ek;
1490 
1491  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1492  UpdatePhys())[id1+i] = rhoeb;
1493  }
1494  // Supersonic flows
1495  else
1496  {
1497  for (j = 0; j < nVariables; ++j)
1498  {
1499  // Extrapolation for supersonic cases
1500  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1501  UpdatePhys())[id1+i] = Fwd[j][pnt];
1502  }
1503  }
1504  }
1505  }
1506  }
Nektar::CompressibleFlowSystem::m_fieldStorage
Array< OneD, Array< OneD, NekDouble > > m_fieldStorage
Definition: CompressibleFlowSystem.h:132
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Vmath::Vdiv
void Vdiv(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:227
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:410
Nektar::SolverUtils::EquationSystem::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
Definition: EquationSystem.h:516
Vmath::Smul
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
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:841
Nektar::SolverUtils::EquationSystem::GetPhys_Offset
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
Definition: EquationSystem.h:856
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Nektar::SolverUtils::EquationSystem::GetExpSize
SOLVER_UTILS_EXPORT int GetExpSize()
Definition: EquationSystem.h:851
Vmath::Vadd
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
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
void Nektar::CompressibleFlowSystem::PressureOutflowBC ( int  bcRegion,
int  cnt,
Array< OneD, Array< OneD, NekDouble > > &  Fwd,
Array< OneD, Array< OneD, NekDouble > > &  physarray 
)
protected

Pressure outflow boundary conditions for compressible flow problems.

Definition at line 1048 of file CompressibleFlowSystem.cpp.

References Nektar::SolverUtils::EquationSystem::GetExpSize(), Nektar::SolverUtils::EquationSystem::GetPhys_Offset(), Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_gamma, m_pInf, Nektar::SolverUtils::EquationSystem::m_spacedim, Nektar::SolverUtils::EquationSystem::m_traceNormals, m_uInf, m_vInf, m_wInf, npts, Vmath::Smul(), Vmath::Vabs(), Vmath::Vadd(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by SetCommonBC().

1053  {
1054  int i, j;
1055  int nTracePts = GetTraceTotPoints();
1056  int nVariables = physarray.num_elements();
1057  int nDimensions = m_spacedim;
1058 
1059  const Array<OneD, const int> &traceBndMap
1060  = m_fields[0]->GetTraceBndMap();
1061 
1062  NekDouble gamma = m_gamma;
1063  NekDouble gammaMinusOne = gamma - 1.0;
1064  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
1065 
1066  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
1067  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
1068  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
1069  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
1070 
1071  // Computing the normal velocity for characteristics coming
1072  // from outside the computational domain
1073  velInf[0] = m_uInf;
1074  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
1075  if (nDimensions == 2 || nDimensions == 3)
1076  {
1077  velInf[1] = m_vInf;
1078  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
1079  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
1080  }
1081  if (nDimensions == 3)
1082  {
1083  velInf[2] = m_wInf;
1084  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
1085  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
1086  }
1087 
1088  // Computing the normal velocity for characteristics coming
1089  // from inside the computational domain
1090  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1091  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
1092  for (i = 0; i < nDimensions; ++i)
1093  {
1094  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
1095  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
1096  }
1097 
1098  // Computing the absolute value of the velocity in order to compute the
1099  // Mach number to decide whether supersonic or subsonic
1100  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
1101  for (i = 0; i < nDimensions; ++i)
1102  {
1103  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
1104  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
1105  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
1106  }
1107  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
1108 
1109  // Get speed of sound
1110  Array<OneD, NekDouble > soundSpeed(nTracePts);
1111  Array<OneD, NekDouble > pressure (nTracePts);
1112 
1113  for (i = 0; i < nTracePts; i++)
1114  {
1115  if (m_spacedim == 1)
1116  {
1117  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
1118  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
1119  }
1120  else if (m_spacedim == 2)
1121  {
1122  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
1123  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1124  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
1125  }
1126  else
1127  {
1128  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
1129  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1130  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
1131  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
1132  }
1133 
1134  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
1135  }
1136 
1137  // Get Mach
1138  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
1139  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
1140  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
1141 
1142  // Auxiliary variables
1143  int e, id1, id2, npts, pnt;
1144  NekDouble rhoeb;
1145 
1146  // Loop on the bcRegions
1147  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->
1148  GetExpSize(); ++e)
1149  {
1150  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1151  GetExp(e)->GetTotPoints();
1152  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1153  GetPhys_Offset(e) ;
1154  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1155 
1156  // Loop on points of bcRegion 'e'
1157  for (i = 0; i < npts; i++)
1158  {
1159  pnt = id2+i;
1160 
1161  // Subsonic flows
1162  if (Mach[pnt] < 0.99)
1163  {
1164  // Kinetic energy calculation
1165  NekDouble Ek = 0.0;
1166  for (j = 1; j < nVariables-1; ++j)
1167  {
1168  Ek += 0.5 * (Fwd[j][pnt] * Fwd[j][pnt]) / Fwd[0][pnt];
1169  }
1170 
1171  rhoeb = m_pInf * gammaMinusOneInv + Ek;
1172 
1173  // Partial extrapolation for subsonic cases
1174  for (j = 0; j < nVariables-1; ++j)
1175  {
1176  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1177  UpdatePhys())[id1+i] = Fwd[j][pnt];
1178  }
1179 
1180  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1181  UpdatePhys())[id1+i] = rhoeb;
1182  }
1183  // Supersonic flows
1184  else
1185  {
1186  for (j = 0; j < nVariables; ++j)
1187  {
1188  // Extrapolation for supersonic cases
1189  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1190  UpdatePhys())[id1+i] = Fwd[j][pnt];
1191  }
1192  }
1193  }
1194  }
1195  }
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Vmath::Vdiv
void Vdiv(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:227
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:410
Nektar::SolverUtils::EquationSystem::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
Definition: EquationSystem.h:516
Vmath::Smul
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
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:841
Nektar::SolverUtils::EquationSystem::GetPhys_Offset
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
Definition: EquationSystem.h:856
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Nektar::SolverUtils::EquationSystem::GetExpSize
SOLVER_UTILS_EXPORT int GetExpSize()
Definition: EquationSystem.h:851
Vmath::Vadd
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
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
void Nektar::CompressibleFlowSystem::PressureOutflowFileBC ( int  bcRegion,
int  cnt,
Array< OneD, Array< OneD, NekDouble > > &  Fwd,
Array< OneD, Array< OneD, NekDouble > > &  physarray 
)
protected

Pressure outflow boundary conditions for compressible flow problems.

Definition at line 1202 of file CompressibleFlowSystem.cpp.

References Nektar::SolverUtils::EquationSystem::GetExpSize(), Nektar::SolverUtils::EquationSystem::GetPhys_Offset(), Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_gamma, m_pressureStorage, Nektar::SolverUtils::EquationSystem::m_spacedim, Nektar::SolverUtils::EquationSystem::m_traceNormals, m_uInf, m_vInf, m_wInf, npts, Vmath::Smul(), Vmath::Vabs(), Vmath::Vadd(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by SetCommonBC().

1207  {
1208  int i, j;
1209  int nTracePts = GetTraceTotPoints();
1210  int nVariables = physarray.num_elements();
1211  int nDimensions = m_spacedim;
1212 
1213  const Array<OneD, const int> &traceBndMap
1214  = m_fields[0]->GetTraceBndMap();
1215 
1216  NekDouble gamma = m_gamma;
1217  NekDouble gammaMinusOne = gamma - 1.0;
1218  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
1219 
1220  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
1221  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
1222  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
1223  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
1224 
1225  // Computing the normal velocity for characteristics coming
1226  // from outside the computational domain
1227  velInf[0] = m_uInf;
1228  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
1229  if (nDimensions == 2 || nDimensions == 3)
1230  {
1231  velInf[1] = m_vInf;
1232  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
1233  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
1234  }
1235  if (nDimensions == 3)
1236  {
1237  velInf[2] = m_wInf;
1238  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
1239  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
1240  }
1241 
1242  // Computing the normal velocity for characteristics coming
1243  // from inside the computational domain
1244  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1245  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
1246  for (i = 0; i < nDimensions; ++i)
1247  {
1248  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
1249  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
1250  }
1251 
1252  // Computing the absolute value of the velocity in order to compute the
1253  // Mach number to decide whether supersonic or subsonic
1254  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
1255  for (i = 0; i < nDimensions; ++i)
1256  {
1257  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
1258  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
1259  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
1260  }
1261  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
1262 
1263  // Get speed of sound
1264  Array<OneD, NekDouble > soundSpeed (nTracePts, 0.0);
1265  Array<OneD, NekDouble > pressure (nTracePts, 0.0);
1266 
1267  for (i = 0; i < nTracePts; i++)
1268  {
1269  if (m_spacedim == 1)
1270  {
1271  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
1272  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
1273  }
1274  else if (m_spacedim == 2)
1275  {
1276  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
1277  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1278  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
1279  }
1280  else
1281  {
1282  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
1283  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1284  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
1285  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
1286  }
1287 
1288  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
1289  }
1290 
1291  // Get Mach
1292  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
1293  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
1294  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
1295 
1296  // Auxiliary variables
1297  int e, id1, id2, npts, pnt;
1298  NekDouble rhoeb;
1299 
1300  // Loop on the bcRegions
1301  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->
1302  GetExpSize(); ++e)
1303  {
1304  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1305  GetExp(e)->GetTotPoints();
1306  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1307  GetPhys_Offset(e);
1308  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1309 
1310  // Loop on points of bcRegion 'e'
1311  for (i = 0; i < npts; i++)
1312  {
1313  pnt = id2+i;
1314 
1315  // Subsonic flows
1316  if (Mach[pnt] < 0.99)
1317  {
1318  // Kinetic energy calculation
1319  NekDouble Ek = 0.0;
1320  for (j = 1; j < nVariables-1; ++j)
1321  {
1322  Ek += 0.5 * (Fwd[j][pnt] * Fwd[j][pnt]) / Fwd[0][pnt];
1323  }
1324 
1325  rhoeb = m_pressureStorage[id1+i] * gammaMinusOneInv + Ek;
1326 
1327  // Partial extrapolation for subsonic cases
1328  for (j = 0; j < nVariables-1; ++j)
1329  {
1330  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1331  UpdatePhys())[id1+i] = Fwd[j][pnt];
1332  }
1333 
1334  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1335  UpdatePhys())[id1+i] = rhoeb;
1336  }
1337  // Supersonic flows
1338  else
1339  {
1340  for (j = 0; j < nVariables; ++j)
1341  {
1342  // Extrapolation for supersonic cases
1343  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1344  UpdatePhys())[id1+i] = Fwd[j][pnt];
1345  }
1346  }
1347  }
1348  }
1349  }
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Vmath::Vdiv
void Vdiv(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:227
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:410
Nektar::SolverUtils::EquationSystem::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
Definition: EquationSystem.h:516
Vmath::Smul
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
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::CompressibleFlowSystem::m_pressureStorage
Array< OneD, NekDouble > m_pressureStorage
Definition: CompressibleFlowSystem.h:129
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:841
Nektar::SolverUtils::EquationSystem::GetPhys_Offset
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
Definition: EquationSystem.h:856
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Nektar::SolverUtils::EquationSystem::GetExpSize
SOLVER_UTILS_EXPORT int GetExpSize()
Definition: EquationSystem.h:851
Vmath::Vadd
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
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
void Nektar::CompressibleFlowSystem::PressureOutflowNonReflectiveBC ( int  bcRegion,
int  cnt,
Array< OneD, Array< OneD, NekDouble > > &  Fwd,
Array< OneD, Array< OneD, NekDouble > > &  physarray 
)
protected

Pressure outflow non-reflective boundary conditions for compressible flow problems.

Definition at line 895 of file CompressibleFlowSystem.cpp.

References Nektar::SolverUtils::EquationSystem::GetExpSize(), Nektar::SolverUtils::EquationSystem::GetPhys_Offset(), Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_gamma, m_pInf, Nektar::SolverUtils::EquationSystem::m_spacedim, Nektar::SolverUtils::EquationSystem::m_traceNormals, m_uInf, m_vInf, m_wInf, npts, Vmath::Smul(), Vmath::Vabs(), Vmath::Vadd(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by SetCommonBC().

900  {
901  int i, j;
902  int nTracePts = GetTraceTotPoints();
903  int nVariables = physarray.num_elements();
904  int nDimensions = m_spacedim;
905 
906  const Array<OneD, const int> &traceBndMap
907  = m_fields[0]->GetTraceBndMap();
908 
909  NekDouble gamma = m_gamma;
910  NekDouble gammaMinusOne = gamma - 1.0;
911  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
912 
913  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
914  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
915  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
916  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
917 
918  // Computing the normal velocity for characteristics coming
919  // from outside the computational domain
920  velInf[0] = m_uInf;
921  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
922  if (nDimensions == 2 || nDimensions == 3)
923  {
924  velInf[1] = m_vInf;
925  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
926  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
927  }
928  if (nDimensions == 3)
929  {
930  velInf[2] = m_wInf;
931  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
932  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
933  }
934 
935  // Computing the normal velocity for characteristics coming
936  // from inside the computational domain
937  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
938  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
939  for (i = 0; i < nDimensions; ++i)
940  {
941  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
942  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
943  }
944 
945  // Computing the absolute value of the velocity in order to compute the
946  // Mach number to decide whether supersonic or subsonic
947  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
948  for (i = 0; i < nDimensions; ++i)
949  {
950  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
951  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
952  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
953  }
954  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
955 
956  // Get speed of sound
957  Array<OneD, NekDouble > soundSpeed(nTracePts);
958  Array<OneD, NekDouble > pressure (nTracePts);
959 
960  for (i = 0; i < nTracePts; i++)
961  {
962  if (m_spacedim == 1)
963  {
964  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
965  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
966  }
967  else if (m_spacedim == 2)
968  {
969  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
970  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
971  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
972  }
973  else
974  {
975  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
976  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
977  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
978  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
979  }
980 
981  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
982  }
983 
984  // Get Mach
985  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
986  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
987  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
988 
989  // Auxiliary variables
990  int e, id1, id2, npts, pnt;
991  NekDouble rhoeb;
992 
993  // Loop on the bcRegions
994  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->
995  GetExpSize(); ++e)
996  {
997  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
998  GetExp(e)->GetTotPoints();
999  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1000  GetPhys_Offset(e) ;
1001  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1002 
1003  // Loop on points of bcRegion 'e'
1004  for (i = 0; i < npts; i++)
1005  {
1006  pnt = id2+i;
1007 
1008  // Subsonic flows
1009  if (Mach[pnt] < 0.99)
1010  {
1011  // Kinetic energy calculation
1012  NekDouble Ek = 0.0;
1013  for (j = 1; j < nVariables-1; ++j)
1014  {
1015  Ek += 0.5 * (Fwd[j][pnt] * Fwd[j][pnt]) / Fwd[0][pnt];
1016  }
1017 
1018  rhoeb = m_pInf * gammaMinusOneInv + Ek;
1019 
1020  // Partial extrapolation for subsonic cases
1021  for (j = 0; j < nVariables-1; ++j)
1022  {
1023  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1024  UpdatePhys())[id1+i] = Fwd[j][pnt];
1025  }
1026 
1027  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1028  UpdatePhys())[id1+i] = 2.0 * rhoeb - Fwd[nVariables-1][pnt];
1029  }
1030  // Supersonic flows
1031  else
1032  {
1033  for (j = 0; j < nVariables; ++j)
1034  {
1035  // Extrapolation for supersonic cases
1036  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1037  UpdatePhys())[id1+i] = Fwd[j][pnt];
1038  }
1039  }
1040  }
1041  }
1042  }
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Vmath::Vdiv
void Vdiv(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:227
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:410
Nektar::SolverUtils::EquationSystem::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
Definition: EquationSystem.h:516
Vmath::Smul
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
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:841
Nektar::SolverUtils::EquationSystem::GetPhys_Offset
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
Definition: EquationSystem.h:856
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Nektar::SolverUtils::EquationSystem::GetExpSize
SOLVER_UTILS_EXPORT int GetExpSize()
Definition: EquationSystem.h:851
Vmath::Vadd
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
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
void Nektar::CompressibleFlowSystem::RiemannInvariantBC ( int  bcRegion,
int  cnt,
Array< OneD, Array< OneD, NekDouble > > &  Fwd,
Array< OneD, Array< OneD, NekDouble > > &  physarray 
)
protected

Outflow characteristic boundary conditions for compressible flow problems.

Definition at line 646 of file CompressibleFlowSystem.cpp.

References Nektar::SolverUtils::EquationSystem::GetPhys_Offset(), Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_gamma, m_pInf, m_rhoInf, Nektar::SolverUtils::EquationSystem::m_spacedim, Nektar::SolverUtils::EquationSystem::m_traceNormals, m_uInf, m_vInf, m_wInf, Vmath::Smul(), Vmath::Vabs(), Vmath::Vadd(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by SetCommonBC().

651  {
652  int i, j;
653  int nTracePts = GetTraceTotPoints();
654  int nDimensions = m_spacedim;
655 
656  const Array<OneD, const int> &traceBndMap
657  = m_fields[0]->GetTraceBndMap();
658 
659  NekDouble gamma = m_gamma;
660  NekDouble gammaInv = 1.0 / gamma;
661  NekDouble gammaMinusOne = gamma - 1.0;
662  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
663 
664  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
665  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
666  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
667  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
668 
669  // Computing the normal velocity for characteristics coming
670  // from outside the computational domain
671  velInf[0] = m_uInf;
672  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
673  if (nDimensions == 2 || nDimensions == 3)
674  {
675  velInf[1] = m_vInf;
676  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
677  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
678  }
679  if (nDimensions == 3)
680  {
681  velInf[2] = m_wInf;
682  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
683  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
684  }
685 
686  // Computing the normal velocity for characteristics coming
687  // from inside the computational domain
688  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
689  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
690  for (i = 0; i < nDimensions; ++i)
691  {
692  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
693  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
694  }
695 
696  // Computing the absolute value of the velocity in order to compute the
697  // Mach number to decide whether supersonic or subsonic
698  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
699  for (i = 0; i < nDimensions; ++i)
700  {
701  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
702  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
703  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
704  }
705  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
706 
707  // Get speed of sound
708  Array<OneD, NekDouble > soundSpeed(nTracePts);
709  Array<OneD, NekDouble > pressure (nTracePts);
710 
711  for (i = 0; i < nTracePts; i++)
712  {
713  if (m_spacedim == 1)
714  {
715  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
716  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
717  }
718  else if (m_spacedim == 2)
719  {
720  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
721  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
722  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
723  }
724  else
725  {
726  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
727  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
728  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
729  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
730  }
731 
732  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
733  }
734 
735  // Get Mach
736  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
737  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
738  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
739 
740  // Auxiliary variables
741  int eMax;
742  int e, id1, id2, nBCEdgePts, pnt;
743  NekDouble cPlus, rPlus, cMinus, rMinus, VDBC, VNBC;
744  Array<OneD, NekDouble> velBC(nDimensions, 0.0);
745  Array<OneD, NekDouble> rhoVelBC(nDimensions, 0.0);
746  NekDouble rhoBC, EBC, cBC, sBC, pBC;
747 
748  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
749 
750  // Loop on bcRegions
751  for (e = 0; e < eMax; ++e)
752  {
753  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
754  GetExp(e)->GetTotPoints();
755 
756  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
757  GetPhys_Offset(e);
758  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
759 
760  // Loop on the points of the bcRegion
761  for (i = 0; i < nBCEdgePts; i++)
762  {
763  pnt = id2+i;
764 
765  // Impose inflow Riemann invariant
766  if (Vn[pnt] <= 0.0)
767  {
768  // Subsonic flows
769  if (Mach[pnt] < 1.00)
770  {
771  // + Characteristic from inside
772  cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
773  rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
774 
775  // - Characteristic from boundary
776  cMinus = sqrt(gamma * m_pInf / m_rhoInf);
777  rMinus = VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
778  }
779  else
780  {
781  // + Characteristic from inside
782  cPlus = sqrt(gamma * m_pInf / m_rhoInf);
783  rPlus = VnInf[pnt] + 2.0 * cPlus * gammaMinusOneInv;
784 
785  // + Characteristic from inside
786  cMinus = sqrt(gamma * m_pInf / m_rhoInf);
787  rMinus = VnInf[pnt] - 2.0 * cPlus * gammaMinusOneInv;
788  }
789 
790  // Riemann boundary variables
791  VNBC = 0.5 * (rPlus + rMinus);
792  cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
793  VDBC = VNBC - VnInf[pnt];
794 
795  // Thermodynamic boundary variables
796  sBC = m_pInf / (pow(m_rhoInf, gamma));
797  rhoBC = pow((cBC * cBC) / (gamma * sBC), gammaMinusOneInv);
798  pBC = rhoBC * cBC * cBC * gammaInv;
799 
800  // Kinetic energy initialiasation
801  NekDouble EkBC = 0.0;
802 
803  // Boundary velocities
804  for ( j = 0; j < nDimensions; ++j)
805  {
806  velBC[j] = velInf[j] + VDBC * m_traceNormals[j][pnt];
807  rhoVelBC[j] = rhoBC * velBC[j];
808  EkBC += 0.5 * rhoBC * velBC[j]*velBC[j];
809  }
810 
811  // Boundary energy
812  EBC = pBC * gammaMinusOneInv + EkBC;
813 
814  // Imposing Riemann Invariant boundary conditions
815  (m_fields[0]->GetBndCondExpansions()[bcRegion]->
816  UpdatePhys())[id1+i] = rhoBC;
817  for (j = 0; j < nDimensions; ++j)
818  {
819  (m_fields[j+1]->GetBndCondExpansions()[bcRegion]->
820  UpdatePhys())[id1+i] = rhoVelBC[j];
821  }
822  (m_fields[nDimensions+1]->GetBndCondExpansions()[bcRegion]->
823  UpdatePhys())[id1+i] = EBC;
824 
825  }
826  else // Impose outflow Riemann invariant
827  {
828  // Subsonic flows
829  if (Mach[pnt] < 1.00)
830  {
831  // + Characteristic from inside
832  cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
833  rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
834 
835  // - Characteristic from boundary
836  cMinus = sqrt(gamma * m_pInf / m_rhoInf);
837  rMinus = VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
838  }
839  else
840  {
841  // + Characteristic from inside
842  cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
843  rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
844 
845  // + Characteristic from inside
846  cMinus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
847  rMinus = Vn[pnt] - 2.0 * cPlus * gammaMinusOneInv;
848  }
849 
850  // Riemann boundary variables
851  VNBC = 0.5 * (rPlus + rMinus);
852  cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
853  VDBC = VNBC - Vn[pnt];
854 
855  // Thermodynamic boundary variables
856  sBC = pressure[pnt] / (pow(Fwd[0][pnt], gamma));
857  rhoBC = pow((cBC * cBC) / (gamma * sBC), gammaMinusOneInv);
858  pBC = rhoBC * cBC * cBC * gammaInv;
859 
860  // Kinetic energy initialiasation
861  NekDouble EkBC = 0.0;
862 
863  // Boundary velocities
864  for ( j = 0; j < nDimensions; ++j)
865  {
866  velBC[j] = Fwd[j+1][pnt] / Fwd[0][pnt] +
867  VDBC * m_traceNormals[j][pnt];
868  rhoVelBC[j] = rhoBC * velBC[j];
869  EkBC += 0.5 * rhoBC * velBC[j]*velBC[j];
870  }
871 
872  // Boundary energy
873  EBC = pBC * gammaMinusOneInv + EkBC;
874 
875  // Imposing Riemann Invariant boundary conditions
876  (m_fields[0]->GetBndCondExpansions()[bcRegion]->
877  UpdatePhys())[id1+i] = rhoBC;
878  for (j = 0; j < nDimensions; ++j)
879  {
880  (m_fields[j+1]->GetBndCondExpansions()[bcRegion]->
881  UpdatePhys())[id1+i] = rhoVelBC[j];
882  }
883  (m_fields[nDimensions+1]->GetBndCondExpansions()[bcRegion]->
884  UpdatePhys())[id1+i] = EBC;
885  }
886  }
887  }
888  }
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Vmath::Vdiv
void Vdiv(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:227
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:410
Nektar::SolverUtils::EquationSystem::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
Definition: EquationSystem.h:516
Vmath::Smul
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
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:841
Nektar::SolverUtils::EquationSystem::GetPhys_Offset
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
Definition: EquationSystem.h:856
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Vmath::Vadd
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
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
void Nektar::CompressibleFlowSystem::SetCommonBC ( const std::string &  userDefStr,
const int  n,
const NekDouble  time,
int &  cnt,
Array< OneD, Array< OneD, NekDouble > > &  Fwd,
Array< OneD, Array< OneD, NekDouble > > &  inarray 
)
protected

Set boundary conditions which can be: a) Wall and Symmerty BCs implemented at CompressibleFlowSystem level since they are compressible solver specific; b) Time dependent BCs.

Parameters
inarrayfields array.
timetime.

Definition at line 324 of file CompressibleFlowSystem.cpp.

References ASSERTL0, ExtrapOrder0BC(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_session, PressureInflowFileBC(), PressureOutflowBC(), PressureOutflowFileBC(), PressureOutflowNonReflectiveBC(), RiemannInvariantBC(), SymmetryBC(), WallBC(), and WallViscousBC().

Referenced by Nektar::NavierStokesCFE::SetBoundaryConditions(), Nektar::EulerADCFE::SetBoundaryConditions(), and Nektar::EulerCFE::SetBoundaryConditions().

331  {
332  std::string varName;
333  int nvariables = m_fields.num_elements();
334 
335  if(!userDefStr.empty())
336  {
337  if(boost::iequals(userDefStr,"Wall"))
338  {
339  WallBC(n, cnt, Fwd, inarray);
340  }
341  else if(boost::iequals(userDefStr,"WallViscous") ||
342  boost::iequals(userDefStr,"WallAdiabatic"))
343  {
344  // Wall Boundary Condition
345  WallViscousBC(n, cnt, Fwd, inarray);
346  }
347  else if(boost::iequals(userDefStr,"Symmetry"))
348  {
349  // Symmetric Boundary Condition
350  SymmetryBC(n, cnt, Fwd, inarray);
351  }
352  else if(boost::iequals(userDefStr,"RiemannInvariant"))
353  {
354  // Riemann invariant characteristic Boundary Condition
355  RiemannInvariantBC(n, cnt, Fwd, inarray);
356  }
357  else if(boost::iequals(userDefStr,"PressureOutflowNonReflective"))
358  {
359  // Pressure outflow non-reflective Boundary Condition
360  PressureOutflowNonReflectiveBC(n, cnt, Fwd, inarray);
361  }
362  else if(boost::iequals(userDefStr,"PressureOutflow"))
363  {
364  // Pressure outflow Boundary Condition
365  PressureOutflowBC(n, cnt, Fwd, inarray);
366  }
367  else if(boost::iequals(userDefStr,"PressureOutflowFile"))
368  {
369  // Pressure outflow Boundary Condition from file
370  PressureOutflowFileBC(n, cnt, Fwd, inarray);
371  }
372  else if(boost::iequals(userDefStr,"PressureInflowFile"))
373  {
374  // Pressure inflow Boundary Condition from file
375  PressureInflowFileBC(n, cnt, Fwd, inarray);
376  }
377  else if(boost::iequals(userDefStr,"ExtrapOrder0"))
378  {
379  // Extrapolation of the data at the boundaries
380  ExtrapOrder0BC(n, cnt, Fwd, inarray);
381  }
382  else if(boost::iequals(userDefStr,"TimeDependent"))
383  {
384  for (int i = 0; i < nvariables; ++i)
385  {
386  varName = m_session->GetVariable(i);
387  m_fields[i]->EvaluateBoundaryConditions(time, varName);
388  }
389  }
390  else
391  {
392  string errmsg = "Unrecognised boundary condition: ";
393  errmsg += userDefStr;
394  ASSERTL0(false,errmsg.c_str());
395  }
396  }
397  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::CompressibleFlowSystem::WallBC
void WallBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:402
Nektar::CompressibleFlowSystem::SymmetryBC
void SymmetryBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Symmetry boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:558
Nektar::CompressibleFlowSystem::RiemannInvariantBC
void RiemannInvariantBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Outflow characteristic boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:646
Nektar::CompressibleFlowSystem::PressureOutflowBC
void PressureOutflowBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Pressure outflow boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:1048
Nektar::CompressibleFlowSystem::PressureOutflowFileBC
void PressureOutflowFileBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Pressure outflow boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:1202
Nektar::CompressibleFlowSystem::PressureOutflowNonReflectiveBC
void PressureOutflowNonReflectiveBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Pressure outflow non-reflective boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:895
Nektar::CompressibleFlowSystem::WallViscousBC
void WallViscousBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary conditions for viscous compressible flow problems.
Definition: CompressibleFlowSystem.cpp:491
Nektar::CompressibleFlowSystem::PressureInflowFileBC
void PressureInflowFileBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Pressure inflow boundary conditions for compressible flow problems where either the density and the v...
Definition: CompressibleFlowSystem.cpp:1358
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:452
Nektar::CompressibleFlowSystem::ExtrapOrder0BC
void ExtrapOrder0BC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Extrapolation of order 0 for all the variables such that, at the boundaries, a trivial Riemann proble...
Definition: CompressibleFlowSystem.cpp:1513
void Nektar::CompressibleFlowSystem::SetVarPOrderElmt ( const Array< OneD, const Array< OneD, NekDouble > > &  physfield,
Array< OneD, NekDouble > &  PolyOrder 
)
protected
void Nektar::CompressibleFlowSystem::SymmetryBC ( int  bcRegion,
int  cnt,
Array< OneD, Array< OneD, NekDouble > > &  Fwd,
Array< OneD, Array< OneD, NekDouble > > &  physarray 
)
protected

Symmetry boundary conditions for compressible flow problems.

Definition at line 558 of file CompressibleFlowSystem.cpp.

References Nektar::SolverUtils::EquationSystem::GetPhys_Offset(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_spacedim, Nektar::SolverUtils::EquationSystem::m_traceNormals, Vmath::Smul(), Vmath::Vcopy(), Vmath::Vsub(), and Vmath::Vvtvp().

Referenced by SetCommonBC().

563  {
564  int i;
565  int nVariables = physarray.num_elements();
566 
567  const Array<OneD, const int> &traceBndMap
568  = m_fields[0]->GetTraceBndMap();
569 
570  // Take into account that for PDE based shock capturing, eps = 0 at the
571  // wall.
572  int e, id1, id2, nBCEdgePts, eMax;
573 
574  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
575 
576  for (e = 0; e < eMax; ++e)
577  {
578  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
579  GetExp(e)->GetTotPoints();
580  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
581  GetPhys_Offset(e);
582  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
583 
584  if (nVariables == m_spacedim+3)
585  {
586  NekDouble factor = 0.0;
587  NekDouble factor2 = 1.0;
588 
589  Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
590  Vmath::Smul(nBCEdgePts,
591  factor,
592  &Fwd[nVariables-1][id2], 1,
593  &tmp2[0], 1);
594 
595  Vmath::Vsub(nBCEdgePts,
596  &Fwd[nVariables-1][id2], 1,
597  &tmp2[0], 1,
598  &Fwd[nVariables-1][id2], 1);
599 
600  Vmath::Smul(nBCEdgePts,
601  factor2,
602  &Fwd[nVariables-1][id2], 1,
603  &Fwd[nVariables-1][id2], 1);
604  }
605 
606  // For 2D/3D, define: v* = v - 2(v.n)n
607  Array<OneD, NekDouble> tmp(nBCEdgePts, 0.0);
608 
609  // Calculate (v.n)
610  for (i = 0; i < m_spacedim; ++i)
611  {
612  Vmath::Vvtvp(nBCEdgePts,
613  &Fwd[1+i][id2], 1,
614  &m_traceNormals[i][id2], 1,
615  &tmp[0], 1,
616  &tmp[0], 1);
617  }
618 
619  // Calculate 2.0(v.n)
620  Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
621 
622  // Calculate v* = v - 2.0(v.n)n
623  for (i = 0; i < m_spacedim; ++i)
624  {
625  Vmath::Vvtvp(nBCEdgePts,
626  &tmp[0], 1,
627  &m_traceNormals[i][id2], 1,
628  &Fwd[1+i][id2], 1,
629  &Fwd[1+i][id2], 1);
630  }
631 
632  // Copy boundary adjusted values into the boundary expansion
633  for (i = 0; i < nVariables; ++i)
634  {
635  Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
636  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
637  UpdatePhys())[id1], 1);
638  }
639  }
640  }
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Nektar::SolverUtils::EquationSystem::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
Definition: EquationSystem.h:516
Vmath::Smul
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
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Vmath::Vsub
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
Nektar::SolverUtils::EquationSystem::GetPhys_Offset
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
Definition: EquationSystem.h:856
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
virtual void Nektar::CompressibleFlowSystem::v_ExtraFldOutput ( std::vector< Array< OneD, NekDouble > > &  fieldcoeffs,
std::vector< std::string > &  variables 
)
protectedvirtual

Reimplemented from Nektar::SolverUtils::EquationSystem.

void Nektar::CompressibleFlowSystem::v_GenerateSummary ( SolverUtils::SummaryList &  s)
protectedvirtual

Print a summary of time stepping parameters.

Print out a summary with some relevant information.

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Reimplemented in Nektar::EulerADCFE, Nektar::EulerCFE, and Nektar::NavierStokesCFE.

Definition at line 310 of file CompressibleFlowSystem.cpp.

References Nektar::SolverUtils::UnsteadySystem::v_GenerateSummary().

Referenced by Nektar::NavierStokesCFE::v_GenerateSummary(), Nektar::EulerCFE::v_GenerateSummary(), and Nektar::EulerADCFE::v_GenerateSummary().

311  {
312  UnsteadySystem::v_GenerateSummary(s);
313  }
Nektar::SolverUtils::UnsteadySystem::v_GenerateSummary
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
Definition: UnsteadySystem.cpp:460
virtual NekDouble Nektar::CompressibleFlowSystem::v_GetTimeStep ( const Array< OneD, const Array< OneD, NekDouble > > &  inarray)
protectedvirtual

Return the timestep to be used for the next step in the time-marching loop.

See also
UnsteadySystem::GetTimeStep

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

void Nektar::CompressibleFlowSystem::v_InitObject ( )
protectedvirtual

Initialization object for CompressibleFlowSystem class.

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

Reimplemented in Nektar::EulerADCFE, Nektar::NavierStokesCFE, and Nektar::EulerCFE.

Definition at line 65 of file CompressibleFlowSystem.cpp.

References ASSERTL0, Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::MultiRegions::eDiscontinuous, Nektar::MultiRegions::eGalerkin, Nektar::SolverUtils::GetAdvectionFactory(), GetArtificialDynamicViscosity(), Nektar::SolverUtils::GetDiffusionFactory(), GetFluxVector(), GetFluxVectorDeAlias(), GetGamma(), GetNormals(), Nektar::SolverUtils::EquationSystem::GetPressure(), Nektar::SolverUtils::GetRiemannSolverFactory(), GetSmoothArtificialViscosity(), GetVecLocs(), GetViscousFluxVector(), GetViscousFluxVectorDeAlias(), Nektar::SolverUtils::Forcing::Load(), m_advection, m_amplitude, m_C1, m_C2, m_Cp, m_diffusion, m_eps_max, m_EqTypeStr, m_FacH, m_FacL, Nektar::SolverUtils::EquationSystem::m_fields, m_fieldStorage, m_forcing, m_gamma, m_gasConstant, m_hFactor, Nektar::SolverUtils::UnsteadySystem::m_homoInitialFwd, m_Kappa, m_mu, m_mu0, m_omega, m_pInf, m_Prandtl, m_pressureStorage, Nektar::SolverUtils::EquationSystem::m_projectionType, m_rhoInf, m_riemannSolver, m_riemannSolverLDG, Nektar::SolverUtils::EquationSystem::m_session, m_shockCaptureType, m_Skappa, Nektar::SolverUtils::EquationSystem::m_spacedim, Nektar::SolverUtils::EquationSystem::m_specHP_dealiasing, m_steadyStateTol, m_thermalConductivity, m_Twall, m_uInf, m_UInf, m_vecLocs, m_vInf, m_ViscosityType, m_wInf, Nektar::SolverUtils::UnsteadySystem::v_InitObject(), and Vmath::Vcopy().

Referenced by Nektar::EulerCFE::v_InitObject(), Nektar::NavierStokesCFE::v_InitObject(), and Nektar::EulerADCFE::v_InitObject().

66  {
67  UnsteadySystem::v_InitObject();
68 
69  ASSERTL0(m_session->DefinesSolverInfo("UPWINDTYPE"),
70  "No UPWINDTYPE defined in session.");
71 
72  // Do not forwards transform initial condition
73  m_homoInitialFwd = false;
74 
75  // Set up locations of velocity vector.
76  m_vecLocs = Array<OneD, Array<OneD, NekDouble> >(1);
77  m_vecLocs[0] = Array<OneD, NekDouble>(m_spacedim);
78  for (int i = 0; i < m_spacedim; ++i)
79  {
80  m_vecLocs[0][i] = 1 + i;
81  }
82 
83  // Get gamma parameter from session file.
84  ASSERTL0(m_session->DefinesParameter("Gamma"),
85  "Compressible flow sessions must define a Gamma parameter.");
86  m_session->LoadParameter("Gamma", m_gamma, 1.4);
87 
88  // Get E0 parameter from session file.
89  ASSERTL0(m_session->DefinesParameter("pInf"),
90  "Compressible flow sessions must define a pInf parameter.");
91  m_session->LoadParameter("pInf", m_pInf, 101325);
92 
93  // Get rhoInf parameter from session file.
94  ASSERTL0(m_session->DefinesParameter("rhoInf"),
95  "Compressible flow sessions must define a rhoInf parameter.");
96  m_session->LoadParameter("rhoInf", m_rhoInf, 1.225);
97 
98  // Get uInf parameter from session file.
99  ASSERTL0(m_session->DefinesParameter("uInf"),
100  "Compressible flow sessions must define a uInf parameter.");
101  m_session->LoadParameter("uInf", m_uInf, 0.1);
102 
103  m_UInf = m_uInf;
104 
105  // Get vInf parameter from session file.
106  if (m_spacedim == 2 || m_spacedim == 3)
107  {
108  ASSERTL0(m_session->DefinesParameter("vInf"),
109  "Compressible flow sessions must define a vInf parameter"
110  "for 2D/3D problems.");
111  m_session->LoadParameter("vInf", m_vInf, 0.0);
112  m_UInf = sqrt(m_uInf*m_uInf + m_vInf*m_vInf);
113  }
114 
115  // Get wInf parameter from session file.
116  if (m_spacedim == 3)
117  {
118  ASSERTL0(m_session->DefinesParameter("wInf"),
119  "Compressible flow sessions must define a wInf parameter"
120  "for 3D problems.");
121  m_session->LoadParameter("wInf", m_wInf, 0.0);
122  m_UInf = sqrt(m_uInf*m_uInf + m_vInf*m_vInf + m_wInf*m_wInf);
123  }
124 
125  m_session->LoadParameter ("GasConstant", m_gasConstant, 287.058);
126  m_session->LoadParameter ("Twall", m_Twall, 300.15);
127  m_session->LoadSolverInfo("ViscosityType", m_ViscosityType, "Constant");
128  m_session->LoadParameter ("mu", m_mu, 1.78e-05);
129  m_session->LoadParameter ("Skappa", m_Skappa, -2.048);
130  m_session->LoadParameter ("Kappa", m_Kappa, 0.0);
131  m_session->LoadParameter ("mu0", m_mu0, 1.0);
132  m_session->LoadParameter ("FL", m_FacL, 0.0);
133  m_session->LoadParameter ("FH", m_FacH, 0.0);
134  m_session->LoadParameter ("hFactor", m_hFactor, 1.0);
135  m_session->LoadParameter ("epsMax", m_eps_max, 1.0);
136  m_session->LoadParameter ("C1", m_C1, 3.0);
137  m_session->LoadParameter ("C2", m_C2, 5.0);
138  m_session->LoadSolverInfo("ShockCaptureType",
139  m_shockCaptureType, "Off");
140  m_session->LoadParameter ("thermalConductivity",
141  m_thermalConductivity, 0.0257);
142 
143  m_EqTypeStr = m_session->GetSolverInfo("EQTYPE");
144 
145  m_Cp = m_gamma / (m_gamma - 1.0) * m_gasConstant;
146  m_Prandtl = m_Cp * m_mu / m_thermalConductivity;
147 
148  m_session->LoadParameter("amplitude", m_amplitude, 0.001);
149  m_session->LoadParameter("omega", m_omega, 1.0);
150  m_session->LoadParameter("SteadyStateTol", m_steadyStateTol, 0.0);
151 
152  // Forcing terms for the sponge region
153  m_forcing = SolverUtils::Forcing::Load(m_session, m_fields,
154  m_fields.num_elements());
155 
156  // Loop over Boundary Regions for PressureOutflowFileBC
157  int nvariables = m_fields.num_elements();
158  Array<OneD, Array<OneD, NekDouble> > tmpStorage(nvariables);
159  for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
160  {
161  // PressureOutflowFile Boundary Condition
162  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
163  "PressureOutflowFile")
164  {
165  int numBCPts = m_fields[0]->
166  GetBndCondExpansions()[n]->GetNpoints();
167  m_pressureStorage = Array<OneD, NekDouble>(numBCPts, 0.0);
168  for (int i = 0; i < nvariables; ++i)
169  {
170  tmpStorage[i] = Array<OneD, NekDouble>(numBCPts, 0.0);
171 
172  Vmath::Vcopy(
173  numBCPts,
174  m_fields[i]->GetBndCondExpansions()[n]->GetPhys(), 1,
175  tmpStorage[i], 1);
176  }
177  GetPressure(tmpStorage, m_pressureStorage);
178  }
179  }
180 
181  // Loop over Boundary Regions for PressureInflowFileBC
182  m_fieldStorage = Array<OneD, Array<OneD, NekDouble> > (nvariables);
183  for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
184  {
185  // PressureInflowFile Boundary Condition
186  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
187  "PressureInflowFile")
188  {
189  int numBCPts = m_fields[0]->
190  GetBndCondExpansions()[n]->GetNpoints();
191  for (int i = 0; i < nvariables; ++i)
192  {
193  m_fieldStorage[i] = Array<OneD, NekDouble>(numBCPts, 0.0);
194  Vmath::Vcopy(
195  numBCPts,
196  m_fields[i]->GetBndCondExpansions()[n]->GetPhys(), 1,
197  m_fieldStorage[i], 1);
198  }
199  }
200  }
201 
202  // Type of advection class to be used
203  switch(m_projectionType)
204  {
205  // Continuous field
206  case MultiRegions::eGalerkin:
207  {
208  ASSERTL0(false, "Continuous field not supported.");
209  break;
210  }
211  // Discontinuous field
212  case MultiRegions::eDiscontinuous:
213  {
214  string advName, diffName, riemName;
215 
216  // Setting up advection and diffusion operators
217  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
218  m_session->LoadSolverInfo("DiffusionType", diffName, "LDGNS");
219 
220  m_advection = SolverUtils::GetAdvectionFactory()
221  .CreateInstance(advName, advName);
222  m_diffusion = SolverUtils::GetDiffusionFactory()
223  .CreateInstance(diffName, diffName);
224 
225  if (m_specHP_dealiasing)
226  {
227  m_advection->SetFluxVector(&CompressibleFlowSystem::
228  GetFluxVectorDeAlias, this);
229  m_diffusion->SetFluxVectorNS(
230  &CompressibleFlowSystem::GetViscousFluxVectorDeAlias,
231  this);
232  }
233  else
234  {
235  m_advection->SetFluxVector (&CompressibleFlowSystem::
236  GetFluxVector, this);
237  m_diffusion->SetFluxVectorNS(&CompressibleFlowSystem::
238  GetViscousFluxVector, this);
239  }
240 
241  if (m_shockCaptureType=="Smooth" && m_EqTypeStr=="EulerADCFE")
242  {
243  m_advection->SetFluxVector(&CompressibleFlowSystem::
244  GetFluxVector, this);
245 
246  m_diffusion->SetArtificialDiffusionVector(
247  &CompressibleFlowSystem::GetSmoothArtificialViscosity, this);
248  }
249  if (m_shockCaptureType=="NonSmooth" && m_EqTypeStr=="EulerADCFE")
250  {
251  m_advection->SetFluxVector(&CompressibleFlowSystem::
252  GetFluxVector, this);
253 
254  m_diffusion->SetArtificialDiffusionVector(
255  &CompressibleFlowSystem::GetArtificialDynamicViscosity, this);
256  }
257 
258  // Setting up Riemann solver for advection operator
259  m_session->LoadSolverInfo("UpwindType", riemName, "Average");
260 
261  m_riemannSolver = SolverUtils::GetRiemannSolverFactory()
262  .CreateInstance(riemName);
263 
264  // Setting up upwind solver for diffusion operator
265  m_riemannSolverLDG = SolverUtils::GetRiemannSolverFactory()
266  .CreateInstance("UpwindLDG");
267 
268  // Setting up parameters for advection operator Riemann solver
269  m_riemannSolver->SetParam (
270  "gamma", &CompressibleFlowSystem::GetGamma, this);
271  m_riemannSolver->SetAuxVec(
272  "vecLocs", &CompressibleFlowSystem::GetVecLocs, this);
273  m_riemannSolver->SetVector(
274  "N", &CompressibleFlowSystem::GetNormals, this);
275 
276  // Setting up parameters for diffusion operator Riemann solver
277  m_riemannSolverLDG->SetParam (
278  "gamma", &CompressibleFlowSystem::GetGamma, this);
279  m_riemannSolverLDG->SetVector(
280  "vecLocs", &CompressibleFlowSystem::GetVecLocs, this);
281  m_riemannSolverLDG->SetVector(
282  "N", &CompressibleFlowSystem::GetNormals, this);
283 
284  // Concluding initialisation of advection / diffusion operators
285  m_advection->SetRiemannSolver (m_riemannSolver);
286  m_diffusion->SetRiemannSolver (m_riemannSolverLDG);
287  m_advection->InitObject (m_session, m_fields);
288  m_diffusion->InitObject (m_session, m_fields);
289  break;
290  }
291  default:
292  {
293  ASSERTL0(false, "Unsupported projection type.");
294  break;
295  }
296  }
297  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::CompressibleFlowSystem::m_fieldStorage
Array< OneD, Array< OneD, NekDouble > > m_fieldStorage
Definition: CompressibleFlowSystem.h:132
Nektar::SolverUtils::UnsteadySystem::m_homoInitialFwd
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
Definition: UnsteadySystem.h:80
Nektar::LibUtilities::NekFactory::CreateInstance
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
Nektar::CompressibleFlowSystem::GetViscousFluxVector
void GetViscousFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &derivatives, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &viscousTensor)
Return the flux vector for the LDG diffusion problem.
Definition: CompressibleFlowSystem.cpp:1725
Nektar::CompressibleFlowSystem::m_riemannSolverLDG
SolverUtils::RiemannSolverSharedPtr m_riemannSolverLDG
Definition: CompressibleFlowSystem.h:80
Nektar::SolverUtils::GetDiffusionFactory
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:42
Nektar::SolverUtils::EquationSystem::m_projectionType
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Definition: EquationSystem.h:514
Nektar::SolverUtils::EquationSystem::m_specHP_dealiasing
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
Definition: EquationSystem.h:512
Nektar::CompressibleFlowSystem::GetArtificialDynamicViscosity
void GetArtificialDynamicViscosity(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &mu_var)
Nektar::SolverUtils::Forcing::Load
static SOLVER_UTILS_EXPORT std::vector< ForcingSharedPtr > Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields=0)
Definition: Forcing.cpp:86
Nektar::CompressibleFlowSystem::m_vecLocs
Array< OneD, Array< OneD, NekDouble > > m_vecLocs
Definition: CompressibleFlowSystem.h:83
Nektar::CompressibleFlowSystem::GetViscousFluxVectorDeAlias
void GetViscousFluxVectorDeAlias(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &derivatives, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &viscousTensor)
Return the flux vector for the LDG diffusion problem.
Definition: CompressibleFlowSystem.cpp:2058
Nektar::CompressibleFlowSystem::GetVecLocs
const Array< OneD, const Array< OneD, NekDouble > > & GetVecLocs()
Definition: CompressibleFlowSystem.h:294
Nektar::CompressibleFlowSystem::GetFluxVectorDeAlias
void GetFluxVectorDeAlias(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Return the flux vector for the compressible Euler equations by using the de-aliasing technique...
Definition: CompressibleFlowSystem.cpp:1634
Nektar::SolverUtils::GetRiemannSolverFactory
RiemannSolverFactory & GetRiemannSolverFactory()
Definition: RiemannSolver.cpp:63
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::SolverUtils::GetAdvectionFactory
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:46
Nektar::SolverUtils::UnsteadySystem::v_InitObject
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
Definition: UnsteadySystem.cpp:78
Nektar::CompressibleFlowSystem::m_pressureStorage
Array< OneD, NekDouble > m_pressureStorage
Definition: CompressibleFlowSystem.h:129
Nektar::CompressibleFlowSystem::m_advection
SolverUtils::AdvectionSharedPtr m_advection
Definition: CompressibleFlowSystem.h:81
Nektar::CompressibleFlowSystem::m_diffusion
SolverUtils::DiffusionSharedPtr m_diffusion
Definition: CompressibleFlowSystem.h:82
Nektar::CompressibleFlowSystem::GetNormals
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
Definition: CompressibleFlowSystem.h:299
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:452
Nektar::CompressibleFlowSystem::GetSmoothArtificialViscosity
void GetSmoothArtificialViscosity(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &eps_bar)
Nektar::SolverUtils::EquationSystem::GetPressure
SOLVER_UTILS_EXPORT MultiRegions::ExpListSharedPtr GetPressure()
Get pressure field if available.
Definition: EquationSystem.h:748
Nektar::CompressibleFlowSystem::CompressibleFlowSystem
CompressibleFlowSystem(const LibUtilities::SessionReaderSharedPtr &pSession)
Definition: CompressibleFlowSystem.cpp:56
Nektar::CompressibleFlowSystem::m_riemannSolver
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
Definition: CompressibleFlowSystem.h:79
Nektar::CompressibleFlowSystem::GetFluxVector
void GetFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Return the flux vector for the compressible Euler equations.
Definition: CompressibleFlowSystem.cpp:1571
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Nektar::CompressibleFlowSystem::m_forcing
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Definition: CompressibleFlowSystem.h:122
virtual bool Nektar::CompressibleFlowSystem::v_PostIntegrate ( int  step)
protectedvirtual

Reimplemented from Nektar::SolverUtils::UnsteadySystem.

virtual void Nektar::CompressibleFlowSystem::v_SetInitialConditions ( NekDouble  initialtime = 0.0,
bool  dumpInitialConditions = true,
const int  domain = 0 
)
protectedvirtual

Set the physical fields based on a restart file, or a function describing the initial condition given in the session.

Parameters
initialtimeTime at which to evaluate the function.
dumpInitialConditionsWrite the initial condition to file?

Reimplemented from Nektar::SolverUtils::EquationSystem.

Reimplemented in Nektar::EulerADCFE, Nektar::EulerCFE, and Nektar::NavierStokesCFE.

Referenced by Nektar::NavierStokesCFE::v_SetInitialConditions(), Nektar::EulerCFE::v_SetInitialConditions(), and Nektar::EulerADCFE::v_SetInitialConditions().

void Nektar::CompressibleFlowSystem::WallBC ( int  bcRegion,
int  cnt,
Array< OneD, Array< OneD, NekDouble > > &  Fwd,
Array< OneD, Array< OneD, NekDouble > > &  physarray 
)
protected

Wall boundary conditions for compressible flow problems.

Definition at line 402 of file CompressibleFlowSystem.cpp.

References Nektar::SolverUtils::EquationSystem::GetPhys_Offset(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_spacedim, Nektar::SolverUtils::EquationSystem::m_traceNormals, Vmath::Smul(), Vmath::Vcopy(), Vmath::Vsub(), and Vmath::Vvtvp().

Referenced by SetCommonBC().

407  {
408  int i;
409  int nVariables = physarray.num_elements();
410 
411  const Array<OneD, const int> &traceBndMap
412  = m_fields[0]->GetTraceBndMap();
413 
414  // Adjust the physical values of the trace to take
415  // user defined boundaries into account
416  int e, id1, id2, nBCEdgePts, eMax;
417 
418  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
419 
420  for (e = 0; e < eMax; ++e)
421  {
422  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
423  GetExp(e)->GetTotPoints();
424  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
425  GetPhys_Offset(e);
426  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
427 
428  // Boundary condition for epsilon term.
429  if (nVariables == m_spacedim+3)
430  {
431  NekDouble factor = 1.0;
432  NekDouble factor2 = 1.0;
433 
434  Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
435  Vmath::Smul(nBCEdgePts,
436  factor,
437  &Fwd[nVariables-1][id2], 1,
438  &tmp2[0], 1);
439 
440  Vmath::Vsub(nBCEdgePts,
441  &Fwd[nVariables-1][id2], 1,
442  &tmp2[0], 1,
443  &Fwd[nVariables-1][id2], 1);
444 
445  Vmath::Smul(nBCEdgePts,
446  factor2,
447  &Fwd[nVariables-1][id2], 1,
448  &Fwd[nVariables-1][id2], 1);
449 
450  }
451 
452  // For 2D/3D, define: v* = v - 2(v.n)n
453  Array<OneD, NekDouble> tmp(nBCEdgePts, 0.0);
454 
455  // Calculate (v.n)
456  for (i = 0; i < m_spacedim; ++i)
457  {
458  Vmath::Vvtvp(nBCEdgePts,
459  &Fwd[1+i][id2], 1,
460  &m_traceNormals[i][id2], 1,
461  &tmp[0], 1,
462  &tmp[0], 1);
463  }
464 
465  // Calculate 2.0(v.n)
466  Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
467 
468  // Calculate v* = v - 2.0(v.n)n
469  for (i = 0; i < m_spacedim; ++i)
470  {
471  Vmath::Vvtvp(nBCEdgePts,
472  &tmp[0], 1,
473  &m_traceNormals[i][id2], 1,
474  &Fwd[1+i][id2], 1,
475  &Fwd[1+i][id2], 1);
476  }
477 
478  // Copy boundary adjusted values into the boundary expansion
479  for (i = 0; i < nVariables; ++i)
480  {
481  Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
482  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
483  UpdatePhys())[id1], 1);
484  }
485  }
486  }
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Nektar::SolverUtils::EquationSystem::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
Definition: EquationSystem.h:516
Vmath::Smul
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
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Vmath::Vsub
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
Nektar::SolverUtils::EquationSystem::GetPhys_Offset
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
Definition: EquationSystem.h:856
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
void Nektar::CompressibleFlowSystem::WallViscousBC ( int  bcRegion,
int  cnt,
Array< OneD, Array< OneD, NekDouble > > &  Fwd,
Array< OneD, Array< OneD, NekDouble > > &  physarray 
)
protected

Wall boundary conditions for viscous compressible flow problems.

Definition at line 491 of file CompressibleFlowSystem.cpp.

References Nektar::SolverUtils::EquationSystem::GetPhys_Offset(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Neg(), Vmath::Smul(), Vmath::Vcopy(), and Vmath::Vsub().

Referenced by SetCommonBC().

496  {
497  int i;
498  int nVariables = physarray.num_elements();
499 
500  const Array<OneD, const int> &traceBndMap
501  = m_fields[0]->GetTraceBndMap();
502 
503  // Take into account that for PDE based shock capturing, eps = 0 at the
504  // wall. Adjust the physical values of the trace to take user defined
505  // boundaries into account
506  int e, id1, id2, nBCEdgePts, eMax;
507 
508  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
509 
510  for (e = 0; e < eMax; ++e)
511  {
512  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
513  GetExp(e)->GetTotPoints();
514  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
515  GetPhys_Offset(e);
516  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
517 
518  if (nVariables == m_spacedim+3)
519  {
520  NekDouble factor = 0.0;
521  NekDouble factor2 = 1.0;
522 
523  Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
524  Vmath::Smul(nBCEdgePts,
525  factor,
526  &Fwd[nVariables-1][id2], 1,
527  &tmp2[0], 1);
528 
529  Vmath::Vsub(nBCEdgePts,
530  &Fwd[nVariables-1][id2], 1,
531  &tmp2[0], 1,
532  &Fwd[nVariables-1][id2], 1);
533 
534  Vmath::Smul(nBCEdgePts,
535  factor2,
536  &Fwd[nVariables-1][id2], 1,
537  &Fwd[nVariables-1][id2], 1);
538  }
539 
540  for (i = 0; i < m_spacedim; i++)
541  {
542  Vmath::Neg(nBCEdgePts, &Fwd[i+1][id2], 1);
543  }
544 
545  // Copy boundary adjusted values into the boundary expansion
546  for (i = 0; i < nVariables; ++i)
547  {
548  Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
549  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
550  UpdatePhys())[id1], 1);
551  }
552  }
553  }
Vmath::Smul
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
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Vmath::Vsub
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
Nektar::SolverUtils::EquationSystem::GetPhys_Offset
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
Definition: EquationSystem.h:856
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:460
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047

Friends And Related Function Documentation

friend class MemoryManager< CompressibleFlowSystem >
friend

Definition at line 56 of file CompressibleFlowSystem.h.

Member Data Documentation

string Nektar::CompressibleFlowSystem::className
static
Initial value:
=
SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
"CompressibleFlowSystem",
CompressibleFlowSystem::create,
"Auxiliary functions for the compressible flow system.")

Name of class.

Definition at line 66 of file CompressibleFlowSystem.h.

SolverUtils::AdvectionSharedPtr Nektar::CompressibleFlowSystem::m_advection
protected

Definition at line 81 of file CompressibleFlowSystem.h.

Referenced by Nektar::NavierStokesCFE::DoOdeRhs(), Nektar::EulerCFE::DoOdeRhs(), Nektar::EulerADCFE::DoOdeRhs(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_amplitude
protected

Definition at line 109 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_C1
protected

Definition at line 105 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_C2
protected

Definition at line 106 of file CompressibleFlowSystem.h.

Referenced by Nektar::EulerADCFE::DoOdeRhs(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_Cp
protected

Definition at line 104 of file CompressibleFlowSystem.h.

Referenced by GetViscousFluxVector(), and v_InitObject().

SolverUtils::DiffusionSharedPtr Nektar::CompressibleFlowSystem::m_diffusion
protected

Definition at line 82 of file CompressibleFlowSystem.h.

Referenced by Nektar::NavierStokesCFE::DoOdeRhs(), Nektar::EulerADCFE::DoOdeRhs(), Nektar::EulerADCFE::v_InitObject(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_eps_max
protected

Definition at line 102 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

std::string Nektar::CompressibleFlowSystem::m_EqTypeStr
protected

Definition at line 95 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

std::ofstream Nektar::CompressibleFlowSystem::m_errFile
protected

Definition at line 113 of file CompressibleFlowSystem.h.

NekDouble Nektar::CompressibleFlowSystem::m_FacH
protected

Definition at line 101 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_FacL
protected

Definition at line 100 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

Array<OneD, Array<OneD, NekDouble> > Nektar::CompressibleFlowSystem::m_fieldStorage
protected

Definition at line 132 of file CompressibleFlowSystem.h.

Referenced by PressureInflowFileBC(), and v_InitObject().

std::vector<SolverUtils::ForcingSharedPtr> Nektar::CompressibleFlowSystem::m_forcing
protected

Definition at line 122 of file CompressibleFlowSystem.h.

Referenced by Nektar::NavierStokesCFE::DoOdeRhs(), Nektar::EulerCFE::DoOdeRhs(), Nektar::EulerADCFE::DoOdeRhs(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_gamma
protected

Definition at line 84 of file CompressibleFlowSystem.h.

Referenced by Nektar::EulerCFE::EvaluateIsentropicVortex(), Nektar::EulerCFE::GetExactRinglebFlow(), GetGamma(), GetPressure(), PressureInflowFileBC(), PressureOutflowBC(), PressureOutflowFileBC(), PressureOutflowNonReflectiveBC(), RiemannInvariantBC(), Nektar::EulerCFE::SetBoundaryRinglebFlow(), Nektar::EulerCFE::SetInitialRinglebFlow(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_gasConstant
protected

Definition at line 91 of file CompressibleFlowSystem.h.

Referenced by GetGasConstant(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_hFactor
protected

Definition at line 107 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_Kappa
protected

Definition at line 98 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_mu
protected

Definition at line 96 of file CompressibleFlowSystem.h.

Referenced by GetViscousFluxVector(), GetViscousFluxVectorDeAlias(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_mu0
protected

Definition at line 99 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_omega
protected

Definition at line 110 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

StdRegions::StdHexExpSharedPtr Nektar::CompressibleFlowSystem::m_OrthoHexExp
protected

Definition at line 124 of file CompressibleFlowSystem.h.

StdRegions::StdQuadExpSharedPtr Nektar::CompressibleFlowSystem::m_OrthoQuadExp
protected

Definition at line 123 of file CompressibleFlowSystem.h.

NekDouble Nektar::CompressibleFlowSystem::m_pInf
protected

Definition at line 85 of file CompressibleFlowSystem.h.

Referenced by PressureOutflowBC(), PressureOutflowNonReflectiveBC(), RiemannInvariantBC(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_Prandtl
protected

Definition at line 108 of file CompressibleFlowSystem.h.

Referenced by GetViscousFluxVector(), and v_InitObject().

Array<OneD, NekDouble> Nektar::CompressibleFlowSystem::m_pressureStorage
protected

Definition at line 129 of file CompressibleFlowSystem.h.

Referenced by PressureOutflowFileBC(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_rhoInf
protected

Definition at line 86 of file CompressibleFlowSystem.h.

Referenced by RiemannInvariantBC(), and v_InitObject().

SolverUtils::RiemannSolverSharedPtr Nektar::CompressibleFlowSystem::m_riemannSolver
protected

Definition at line 79 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

SolverUtils::RiemannSolverSharedPtr Nektar::CompressibleFlowSystem::m_riemannSolverLDG
protected

Definition at line 80 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

std::string Nektar::CompressibleFlowSystem::m_shockCaptureType
protected

Definition at line 94 of file CompressibleFlowSystem.h.

Referenced by Nektar::EulerADCFE::DoOdeRhs(), Nektar::EulerADCFE::v_InitObject(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_Skappa
protected

Definition at line 97 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

bool Nektar::CompressibleFlowSystem::m_smoothDiffusion
protected

Definition at line 125 of file CompressibleFlowSystem.h.

Referenced by Nektar::EulerADCFE::v_InitObject().

int Nektar::CompressibleFlowSystem::m_steadyStateSteps
protected

Definition at line 116 of file CompressibleFlowSystem.h.

NekDouble Nektar::CompressibleFlowSystem::m_steadyStateTol
protected

Definition at line 119 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_thermalConductivity
protected

Definition at line 103 of file CompressibleFlowSystem.h.

Referenced by GetViscousFluxVector(), GetViscousFluxVectorDeAlias(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_Twall
protected

Definition at line 92 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_uInf
protected

Definition at line 87 of file CompressibleFlowSystem.h.

Referenced by PressureInflowFileBC(), PressureOutflowBC(), PressureOutflowFileBC(), PressureOutflowNonReflectiveBC(), RiemannInvariantBC(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_UInf
protected

Definition at line 90 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

Array<OneD, Array<OneD, NekDouble> > Nektar::CompressibleFlowSystem::m_un
protected

Definition at line 135 of file CompressibleFlowSystem.h.

Array<OneD, Array<OneD, NekDouble> > Nektar::CompressibleFlowSystem::m_vecLocs
protected

Definition at line 83 of file CompressibleFlowSystem.h.

Referenced by GetVecLocs(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_vInf
protected

Definition at line 88 of file CompressibleFlowSystem.h.

Referenced by PressureInflowFileBC(), PressureOutflowBC(), PressureOutflowFileBC(), PressureOutflowNonReflectiveBC(), RiemannInvariantBC(), and v_InitObject().

std::string Nektar::CompressibleFlowSystem::m_ViscosityType
protected

Definition at line 93 of file CompressibleFlowSystem.h.

Referenced by GetViscousFluxVector(), GetViscousFluxVectorDeAlias(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_wInf
protected

Definition at line 89 of file CompressibleFlowSystem.h.

Referenced by PressureInflowFileBC(), PressureOutflowBC(), PressureOutflowFileBC(), PressureOutflowNonReflectiveBC(), RiemannInvariantBC(), and v_InitObject().

Generated on Fri Dec 2 2016 13:10:58 for Nektar++ by   doxygen 1.8.8