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 > > &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 > > &physarray)
 Wall boundary conditions for compressible flow problems. More...
 
void WallViscousBC (int bcRegion, int cnt, 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 > > &physarray)
 Symmetry boundary conditions for compressible flow problems. More...
 
void RiemannInvariantBC (int bcRegion, int cnt, 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 > > &physarray)
 Pressure outflow non-reflective boundary conditions for compressible flow problems. More...
 
void PressureOutflowBC (int bcRegion, int cnt, 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 > > &physarray)
 Pressure outflow boundary conditions for compressible flow problems. More...
 
void PressureInflowFileBC (int bcRegion, int cnt, 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 > > &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 > > &  physarray 
)
protected

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

Definition at line 1574 of file CompressibleFlowSystem.cpp.

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

Referenced by SetCommonBC().

1578  {
1579  int i, j;
1580  int e, pnt;
1581  int id1, id2, nBCEdgePts;
1582  int nTracePts = GetTraceTotPoints();
1583  int nVariables = physarray.num_elements();
1584  int nDimensions = m_spacedim;
1585 
1586  const Array<OneD, const int> &traceBndMap
1587  = m_fields[0]->GetTraceBndMap();
1588 
1589  // Get physical values of the forward trace
1590  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
1591  for (i = 0; i < nVariables; ++i)
1592  {
1593  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
1594  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
1595  }
1596 
1597  int eMax;
1598 
1599  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
1600 
1601  // Loop on bcRegions
1602  for (e = 0; e < eMax; ++e)
1603  {
1604  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1605  GetExp(e)->GetTotPoints();
1606  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1607  GetPhys_Offset(e) ;
1608  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1609 
1610  // Loop on points of bcRegion 'e'
1611  for (i = 0; i < nBCEdgePts; i++)
1612  {
1613  pnt = id2+i;
1614 
1615  // Setting up bcs for density
1616  (m_fields[0]->GetBndCondExpansions()[bcRegion]->
1617  UpdatePhys())[id1+i] = Fwd[0][pnt];
1618 
1619  // Setting up bcs for velocities
1620  for (j = 1; j <=nDimensions; ++j)
1621  {
1622  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1623  UpdatePhys())[id1+i] = Fwd[j][pnt];
1624  }
1625 
1626  // Setting up bcs for energy
1627  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1628  UpdatePhys())[id1+i] = Fwd[nVariables-1][pnt];
1629  }
1630  }
1631  }
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:492
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
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 1640 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().

1643  {
1644  int i, j;
1645  int nq = physfield[0].num_elements();
1646  int nVariables = m_fields.num_elements();
1647 
1648  Array<OneD, NekDouble> pressure(nq);
1649  Array<OneD, Array<OneD, NekDouble> > velocity(m_spacedim);
1650 
1651  // Flux vector for the rho equation
1652  for (i = 0; i < m_spacedim; ++i)
1653  {
1654  velocity[i] = Array<OneD, NekDouble>(nq);
1655  Vmath::Vcopy(nq, physfield[i+1], 1, flux[0][i], 1);
1656  }
1657 
1658  GetVelocityVector(physfield, velocity);
1659  GetPressure (physfield, velocity, pressure);
1660 
1661  // Flux vector for the velocity fields
1662  for (i = 0; i < m_spacedim; ++i)
1663  {
1664  for (j = 0; j < m_spacedim; ++j)
1665  {
1666  Vmath::Vmul(nq, velocity[j], 1, physfield[i+1], 1,
1667  flux[i+1][j], 1);
1668  }
1669 
1670  // Add pressure to appropriate field
1671  Vmath::Vadd(nq, flux[i+1][i], 1, pressure, 1, flux[i+1][i], 1);
1672  }
1673 
1674  // Flux vector for energy.
1675  Vmath::Vadd(nq, physfield[m_spacedim+1], 1, pressure, 1,
1676  pressure, 1);
1677 
1678  for (j = 0; j < m_spacedim; ++j)
1679  {
1680  Vmath::Vmul(nq, velocity[j], 1, pressure, 1,
1681  flux[m_spacedim+1][j], 1);
1682  }
1683 
1684  // For the smooth viscosity model
1685  if (nVariables == m_spacedim+3)
1686  {
1687  // Add a zero row for the advective fluxes
1688  for (j = 0; j < m_spacedim; ++j)
1689  {
1690  Vmath::Zero(nq, flux[m_spacedim+2][j], 1);
1691  }
1692  }
1693  }
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 1703 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().

1706  {
1707  int i, j;
1708  int nq = physfield[0].num_elements();
1709  int nVariables = m_fields.num_elements();
1710 
1711  // Factor to rescale 1d points in dealiasing
1712  NekDouble OneDptscale = 2;
1713  nq = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
1714 
1715  Array<OneD, NekDouble> pressure(nq);
1716  Array<OneD, Array<OneD, NekDouble> > velocity(m_spacedim);
1717 
1718  Array<OneD, Array<OneD, NekDouble> > physfield_interp(nVariables);
1719  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > flux_interp(
1720  nVariables);
1721 
1722  for (i = 0; i < nVariables; ++ i)
1723  {
1724  physfield_interp[i] = Array<OneD, NekDouble>(nq);
1725  flux_interp[i] = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
1726  m_fields[0]->PhysInterp1DScaled(
1727  OneDptscale, physfield[i], physfield_interp[i]);
1728 
1729  for (j = 0; j < m_spacedim; ++j)
1730  {
1731  flux_interp[i][j] = Array<OneD, NekDouble>(nq);
1732  }
1733  }
1734 
1735  // Flux vector for the rho equation
1736  for (i = 0; i < m_spacedim; ++i)
1737  {
1738  velocity[i] = Array<OneD, NekDouble>(nq);
1739 
1740  // Galerkin project solution back to original space
1741  m_fields[0]->PhysGalerkinProjection1DScaled(
1742  OneDptscale, physfield_interp[i+1], flux[0][i]);
1743  }
1744 
1745  GetVelocityVector(physfield_interp, velocity);
1746  GetPressure (physfield_interp, velocity, pressure);
1747 
1748  // Evaluation of flux vector for the velocity fields
1749  for (i = 0; i < m_spacedim; ++i)
1750  {
1751  for (j = 0; j < m_spacedim; ++j)
1752  {
1753  Vmath::Vmul(nq, velocity[j], 1, physfield_interp[i+1], 1,
1754  flux_interp[i+1][j], 1);
1755  }
1756 
1757  // Add pressure to appropriate field
1758  Vmath::Vadd(nq, flux_interp[i+1][i], 1, pressure,1,
1759  flux_interp[i+1][i], 1);
1760  }
1761 
1762  // Galerkin project solution back to origianl space
1763  for (i = 0; i < m_spacedim; ++i)
1764  {
1765  for (j = 0; j < m_spacedim; ++j)
1766  {
1767  m_fields[0]->PhysGalerkinProjection1DScaled(
1768  OneDptscale, flux_interp[i+1][j], flux[i+1][j]);
1769  }
1770  }
1771 
1772  // Evaluation of flux vector for energy
1773  Vmath::Vadd(nq, physfield_interp[m_spacedim+1], 1, pressure, 1,
1774  pressure, 1);
1775 
1776  for (j = 0; j < m_spacedim; ++j)
1777  {
1778  Vmath::Vmul(nq, velocity[j], 1, pressure, 1,
1779  flux_interp[m_spacedim+1][j], 1);
1780 
1781  // Galerkin project solution back to origianl space
1782  m_fields[0]->PhysGalerkinProjection1DScaled(
1783  OneDptscale,
1784  flux_interp[m_spacedim+1][j],
1785  flux[m_spacedim+1][j]);
1786  }
1787  }
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 280 of file CompressibleFlowSystem.h.

References m_gamma.

Referenced by v_InitObject().

281  {
282  return m_gamma;
283  }
NekDouble Nektar::CompressibleFlowSystem::GetGasConstant ( )
inlineprotected

Definition at line 275 of file CompressibleFlowSystem.h.

References m_gasConstant.

276  {
277  return m_gasConstant;
278  }
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 290 of file CompressibleFlowSystem.h.

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

Referenced by v_InitObject().

291  {
292  return m_traceNormals;
293  }
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 2552 of file CompressibleFlowSystem.cpp.

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

2555  {
2556  int nBCEdgePts = physfield[0].num_elements();
2557  NekDouble alpha = -0.5;
2558 
2559  // Calculate ||rho v||^2
2560  Vmath::Vmul(nBCEdgePts, physfield[1], 1, physfield[1], 1, pressure, 1);
2561  for (int i = 1; i < m_spacedim; ++i)
2562  {
2563  Vmath::Vvtvp(nBCEdgePts, physfield[1+i], 1, physfield[1+i], 1,
2564  pressure, 1, pressure, 1);
2565  }
2566  // Divide by rho to get rho*||v||^2
2567  Vmath::Vdiv (nBCEdgePts, pressure, 1, physfield[0], 1, pressure, 1);
2568  // pressure <- E - 0.5*pressure
2569  Vmath::Svtvp(nBCEdgePts, alpha,
2570  pressure, 1, physfield[m_spacedim+1], 1, pressure, 1);
2571  // Multiply by (gamma-1)
2572  Vmath::Smul (nBCEdgePts, m_gamma-1, pressure, 1, pressure, 1);
2573  }
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 285 of file CompressibleFlowSystem.h.

References m_vecLocs.

Referenced by v_InitObject().

286  {
287  return m_vecLocs;
288  }
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 1794 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().

1798  {
1799  int j, k;
1800  int nVariables = m_fields.num_elements();
1801  int nPts = physfield[0].num_elements();
1802 
1803  // Stokes hypotesis
1804  const NekDouble lambda = -2.0/3.0;
1805 
1806  // Auxiliary variables
1807  Array<OneD, NekDouble > mu (nPts, 0.0);
1808  Array<OneD, NekDouble > thermalConductivity(nPts, 0.0);
1809  Array<OneD, NekDouble > mu2 (nPts, 0.0);
1810  Array<OneD, NekDouble > divVel (nPts, 0.0);
1811 
1812  // Variable viscosity through the Sutherland's law
1813  if (m_ViscosityType == "Variable")
1814  {
1815  GetDynamicViscosity(physfield[nVariables-2], mu);
1816  NekDouble tRa = m_Cp / m_Prandtl;
1817  Vmath::Smul(nPts, tRa, &mu[0], 1, &thermalConductivity[0], 1);
1818  }
1819  else
1820  {
1821  Vmath::Fill(nPts, m_mu, &mu[0], 1);
1822  Vmath::Fill(nPts, m_thermalConductivity,
1823  &thermalConductivity[0], 1);
1824  }
1825 
1826  // Computing diagonal terms of viscous stress tensor
1827  Array<OneD, Array<OneD, NekDouble> > tmp(m_spacedim);
1828  Array<OneD, Array<OneD, NekDouble> > Sgg(m_spacedim);
1829 
1830  // mu2 = 2 * mu
1831  Vmath::Smul(nPts, 2.0, &mu[0], 1, &mu2[0], 1);
1832 
1833  // Velocity divergence
1834  for (j = 0; j < m_spacedim; ++j)
1835  {
1836  Vmath::Vadd(nPts, &divVel[0], 1, &derivativesO1[j][j][0], 1,
1837  &divVel[0], 1);
1838  }
1839 
1840  // Velocity divergence scaled by lambda * mu
1841  Vmath::Smul(nPts, lambda, &divVel[0], 1, &divVel[0], 1);
1842  Vmath::Vmul(nPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
1843 
1844  // Diagonal terms of viscous stress tensor (Sxx, Syy, Szz)
1845  // Sjj = 2 * mu * du_j/dx_j - (2 / 3) * mu * sum_j(du_j/dx_j)
1846  for (j = 0; j < m_spacedim; ++j)
1847  {
1848  tmp[j] = Array<OneD, NekDouble>(nPts, 0.0);
1849  Sgg[j] = Array<OneD, NekDouble>(nPts, 0.0);
1850 
1851  Vmath::Vmul(nPts, &mu2[0], 1, &derivativesO1[j][j][0], 1,
1852  &tmp[j][0], 1);
1853 
1854  Vmath::Vadd(nPts, &tmp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
1855  }
1856 
1857  // Extra diagonal terms of viscous stress tensor (Sxy, Sxz, Syz)
1858  // Note: they exist for 2D and 3D problems only
1859  Array<OneD, NekDouble > Sxy(nPts, 0.0);
1860  Array<OneD, NekDouble > Sxz(nPts, 0.0);
1861  Array<OneD, NekDouble > Syz(nPts, 0.0);
1862 
1863  if (m_spacedim == 2)
1864  {
1865  // Sxy = (du/dy + dv/dx)
1866  Vmath::Vadd(nPts, &derivativesO1[0][1][0], 1,
1867  &derivativesO1[1][0][0], 1, &Sxy[0], 1);
1868 
1869  // Sxy = mu * (du/dy + dv/dx)
1870  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
1871  }
1872  else if (m_spacedim == 3)
1873  {
1874  // Sxy = (du/dy + dv/dx)
1875  Vmath::Vadd(nPts, &derivativesO1[0][1][0], 1,
1876  &derivativesO1[1][0][0], 1, &Sxy[0], 1);
1877 
1878  // Sxz = (du/dz + dw/dx)
1879  Vmath::Vadd(nPts, &derivativesO1[0][2][0], 1,
1880  &derivativesO1[2][0][0], 1, &Sxz[0], 1);
1881 
1882  // Syz = (dv/dz + dw/dy)
1883  Vmath::Vadd(nPts, &derivativesO1[1][2][0], 1,
1884  &derivativesO1[2][1][0], 1, &Syz[0], 1);
1885 
1886  // Sxy = mu * (du/dy + dv/dx)
1887  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
1888 
1889  // Sxz = mu * (du/dy + dv/dx)
1890  Vmath::Vmul(nPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
1891 
1892  // Syz = mu * (du/dy + dv/dx)
1893  Vmath::Vmul(nPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
1894  }
1895 
1896  // Energy-related terms
1897  Array<OneD, NekDouble > STx(nPts, 0.0);
1898  Array<OneD, NekDouble > STy(nPts, 0.0);
1899  Array<OneD, NekDouble > STz(nPts, 0.0);
1900 
1901  // Building the viscous flux vector
1902 
1903  // Viscous flux vector for the rho equation
1904  for (k = 0; k < m_spacedim; ++k)
1905  {
1906  Vmath::Zero(nPts, viscousTensor[k][0], 1);
1907  }
1908 
1909  if (m_spacedim == 1)
1910  {
1911  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1912 
1913  // u * Sxx
1914  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
1915 
1916  // k * dT/dx
1917  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1918  &derivativesO1[0][1][0], 1, &tmp1[0], 1);
1919 
1920  // STx = u * Sxx + (K / mu) * dT/dx
1921  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1922  }
1923  else if (m_spacedim == 2)
1924  {
1925  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1926  Array<OneD, NekDouble > tmp2(nPts, 0.0);
1927 
1928  // Computation of STx
1929 
1930  // u * Sxx
1931  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
1932 
1933  // v * Sxy
1934  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1935 
1936  // k * dT/dx
1937  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1938  &derivativesO1[0][2][0], 1, &tmp2[0], 1);
1939 
1940  // STx = u * Sxx + v * Sxy + K * dT/dx
1941  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1942  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
1943 
1944  // Computation of STy
1945 
1946  // Re-initialise temporary arrays
1947  Vmath::Zero(nPts, &tmp1[0], 1);
1948  Vmath::Zero(nPts, &tmp2[0], 1);
1949 
1950  // v * Syy
1951  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sgg[1][0], 1, &STy[0], 1);
1952 
1953  // u * Sxy
1954  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1955 
1956  // k * dT/dy
1957  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1958  &derivativesO1[1][2][0], 1, &tmp2[0], 1);
1959 
1960  // STy = v * Syy + u * Sxy + K * dT/dy
1961  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
1962  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
1963  }
1964  else if (m_spacedim == 3)
1965  {
1966  Array<OneD, NekDouble > tmp1(nPts, 0.0);
1967  Array<OneD, NekDouble > tmp2(nPts, 0.0);
1968  Array<OneD, NekDouble > tmp3(nPts, 0.0);
1969 
1970  // Computation of STx
1971 
1972  // u * Sxx
1973  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
1974 
1975  // v * Sxy
1976  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sxy[0], 1, &tmp1[0], 1);
1977 
1978  // v * Sxz
1979  Vmath::Vmul(nPts, &physfield[2][0], 1, &Sxz[0], 1, &tmp2[0], 1);
1980 
1981  // k * dT/dx
1982  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
1983  &derivativesO1[0][3][0], 1, &tmp3[0], 1);
1984 
1985  // STx = u * Sxx + v * Sxy + w * Sxz + (K / mu) * dT/dx
1986  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
1987  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
1988  Vmath::Vadd(nPts, &STx[0], 1, &tmp3[0], 1, &STx[0], 1);
1989 
1990  // Computation of STy
1991 
1992  // Re-initialise temporary arrays
1993  Vmath::Zero(nPts, &tmp1[0], 1);
1994  Vmath::Zero(nPts, &tmp2[0], 1);
1995  Vmath::Zero(nPts, &tmp3[0], 1);
1996 
1997  // v * Syy
1998  Vmath::Vmul(nPts, &physfield[1][0], 1, &Sgg[1][0], 1, &STy[0], 1);
1999 
2000  // u * Sxy
2001  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxy[0], 1, &tmp1[0], 1);
2002 
2003  // w * Syz
2004  Vmath::Vmul(nPts, &physfield[2][0], 1, &Syz[0], 1, &tmp2[0], 1);
2005 
2006  // k * dT/dy
2007  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
2008  &derivativesO1[1][3][0], 1, &tmp3[0], 1);
2009 
2010  // STy = v * Syy + u * Sxy + w * Syz + K * dT/dy
2011  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
2012  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
2013  Vmath::Vadd(nPts, &STy[0], 1, &tmp3[0], 1, &STy[0], 1);
2014 
2015  // Computation of STz
2016 
2017  // Re-initialise temporary arrays
2018  Vmath::Zero(nPts, &tmp1[0], 1);
2019  Vmath::Zero(nPts, &tmp2[0], 1);
2020  Vmath::Zero(nPts, &tmp3[0], 1);
2021 
2022  // w * Szz
2023  Vmath::Vmul(nPts, &physfield[2][0], 1, &Sgg[2][0], 1, &STz[0], 1);
2024 
2025  // u * Sxz
2026  Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxz[0], 1, &tmp1[0], 1);
2027 
2028  // v * Syz
2029  Vmath::Vmul(nPts, &physfield[1][0], 1, &Syz[0], 1, &tmp2[0], 1);
2030 
2031  // k * dT/dz
2032  Vmath::Vmul(nPts, &thermalConductivity[0], 1,
2033  &derivativesO1[2][3][0], 1, &tmp3[0], 1);
2034 
2035  // STz = w * Szz + u * Sxz + v * Syz + K * dT/dz
2036  Vmath::Vadd(nPts, &STz[0], 1, &tmp1[0], 1, &STz[0], 1);
2037  Vmath::Vadd(nPts, &STz[0], 1, &tmp2[0], 1, &STz[0], 1);
2038  Vmath::Vadd(nPts, &STz[0], 1, &tmp3[0], 1, &STz[0], 1);
2039  }
2040 
2041  switch (m_spacedim)
2042  {
2043  case 1:
2044  {
2045  // f_11v = f_rho = 0
2046  Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
2047 
2048  // f_21v = f_rhou
2049  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
2050 
2051  // f_31v = f_E
2052  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][2][0], 1);
2053  break;
2054  }
2055  case 2:
2056  {
2057  // f_11v = f_rho1 = 0
2058  Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
2059  // f_12v = f_rho2 = 0
2060  Vmath::Zero(nPts, &viscousTensor[1][0][0], 1);
2061 
2062  // f_21v = f_rhou1
2063  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
2064  // f_22v = f_rhou2
2065  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[1][1][0], 1);
2066 
2067  // f_31v = f_rhov1
2068  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[0][2][0], 1);
2069  // f_32v = f_rhov2
2070  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor[1][2][0], 1);
2071 
2072  // f_41v = f_E1
2073  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][3][0], 1);
2074  // f_42v = f_E2
2075  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor[1][3][0], 1);
2076  break;
2077  }
2078  case 3:
2079  {
2080  // f_11v = f_rho1 = 0
2081  Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
2082  // f_12v = f_rho2 = 0
2083  Vmath::Zero(nPts, &viscousTensor[1][0][0], 1);
2084  // f_13v = f_rho3 = 0
2085  Vmath::Zero(nPts, &viscousTensor[2][0][0], 1);
2086 
2087  // f_21v = f_rhou1
2088  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
2089  // f_22v = f_rhou2
2090  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[1][1][0], 1);
2091  // f_23v = f_rhou3
2092  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor[2][1][0], 1);
2093 
2094  // f_31v = f_rhov1
2095  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[0][2][0], 1);
2096  // f_32v = f_rhov2
2097  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor[1][2][0], 1);
2098  // f_33v = f_rhov3
2099  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor[2][2][0], 1);
2100 
2101  // f_31v = f_rhow1
2102  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor[0][3][0], 1);
2103  // f_32v = f_rhow2
2104  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor[1][3][0], 1);
2105  // f_33v = f_rhow3
2106  Vmath::Vcopy(nPts, &Sgg[2][0], 1, &viscousTensor[2][3][0], 1);
2107 
2108  // f_41v = f_E1
2109  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][4][0], 1);
2110  // f_42v = f_E2
2111  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor[1][4][0], 1);
2112  // f_43v = f_E3
2113  Vmath::Vcopy(nPts, &STz[0], 1, &viscousTensor[2][4][0], 1);
2114  break;
2115  }
2116  default:
2117  {
2118  ASSERTL0(false, "Illegal expansion dimension");
2119  }
2120  }
2121  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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 2127 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().

2131  {
2132 #if 0
2133  int i, j, k;
2134  int nVariables = m_fields.num_elements();
2135  int nPts = physfield[0].num_elements();
2136 
2137  int variables_phys = physfield.num_elements();
2138 
2139  // Factor to rescale 1d points in dealiasing.
2140  NekDouble OneDptscale = 2;
2141 
2142  // Get number of points to dealias a cubic non-linearity
2143  nPts = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
2144 
2145  int nVariables_aux = derivativesO1[0].num_elements();
2146 
2147  Array<OneD, Array<OneD, NekDouble> > physfield_interp(variables_phys);
2148  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > derivativesO1_interp(
2149  m_spacedim);
2150  Array<OneD, Array<OneD, Array<OneD, NekDouble> > >viscousTensor_interp(
2151  m_spacedim);
2152 
2153  for (i = 0; i < m_spacedim; ++ i)
2154  {
2155  viscousTensor_interp[i] = Array<OneD, Array<OneD, NekDouble> >(
2156  nVariables);
2157  for (j = 0; j < nVariables; ++j)
2158  {
2159  viscousTensor_interp[i][j] = Array<OneD, NekDouble>(nPts);
2160  }
2161  }
2162 
2163  // Stokes hypotesis
2164  NekDouble lambda = -2.0/3.0;
2165 
2166  // Auxiliary variables
2167  Array<OneD, NekDouble > mu (nPts, 0.0);
2168  Array<OneD, NekDouble > mu2 (nPts, 0.0);
2169  Array<OneD, NekDouble > divVel (nPts, 0.0);
2170  Array<OneD, NekDouble > pressure (nPts, 0.0);
2171  Array<OneD, NekDouble > temperature(nPts, 0.0);
2172 
2173  for (i = 0; i < nVariables; ++i)
2174  {
2175  m_fields[0]->PhysInterp1DScaled(
2176  OneDptscale, physfield[i], fields_interp[i]);
2177  }
2178 
2179  for (i = 0; i < variables_phys; ++i)
2180  {
2181  physfield_interp[i] = Array<OneD, NekDouble> (nPts);
2182 
2183  // Interpolation to higher space
2184  m_fields[0]->PhysInterp1DScaled(OneDptscale, physfield[i],
2185  physfield_interp[i]);
2186  }
2187 
2188  for (i = 0; i < m_spacedim; ++i)
2189  {
2190  derivativesO1_interp[i] = Array<OneD, Array<OneD, NekDouble> >(
2191  nVariables_aux);
2192  for (j = 0; j < nVariables_aux; ++j)
2193  {
2194  derivativesO1_interp[i][j] = Array<OneD, NekDouble>(nPts);
2195  m_fields[0]->PhysInterp1DScaled(
2196  OneDptscale, derivativesO1[i][j],
2197  derivativesO1_interp[i][j]);
2198  }
2199  }
2200 
2201  // Thermodynamic related quantities
2202  GetPressure(fields_interp, pressure);
2203  GetTemperature(fields_interp, pressure, temperature);
2204 
2205  // Variable viscosity through the Sutherland's law
2206  if (m_ViscosityType == "Variable")
2207  {
2208  GetDynamicViscosity(fields_interp[variables_phys-1], mu);
2209  }
2210  else
2211  {
2212  Vmath::Sadd(nPts, m_mu, &mu[0], 1, &mu[0], 1);
2213  }
2214 
2215  // Computing diagonal terms of viscous stress tensor
2216  Array<OneD, Array<OneD, NekDouble> > tmp(m_spacedim);
2217  Array<OneD, Array<OneD, NekDouble> > Sgg(m_spacedim);
2218 
2219  // mu2 = 2 * mu
2220  Vmath::Smul(nPts, 2.0, &mu[0], 1, &mu2[0], 1);
2221 
2222  // Velocity divergence
2223  for (j = 0; j < m_spacedim; ++j)
2224  {
2225  Vmath::Vadd(nPts, &divVel[0], 1, &derivativesO1_interp[j][j][0], 1,
2226  &divVel[0], 1);
2227  }
2228 
2229  // Velocity divergence scaled by lambda * mu
2230  Vmath::Smul(nPts, lambda, &divVel[0], 1, &divVel[0], 1);
2231  Vmath::Vmul(nPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
2232 
2233  // Digonal terms of viscous stress tensor (Sxx, Syy, Szz)
2234  // Sjj = 2 * mu * du_j/dx_j - (2 / 3) * mu * sum_j(du_j/dx_j)
2235  for (j = 0; j < m_spacedim; ++j)
2236  {
2237  tmp[j] = Array<OneD, NekDouble>(nPts, 0.0);
2238  Sgg[j] = Array<OneD, NekDouble>(nPts, 0.0);
2239 
2240  Vmath::Vmul(nPts, &mu2[0], 1, &derivativesO1_interp[j][j][0], 1,
2241  &tmp[j][0], 1);
2242 
2243  Vmath::Vadd(nPts, &tmp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
2244  }
2245 
2246  // Extra diagonal terms of viscous stress tensor (Sxy, Sxz, Syz)
2247  // Note: they exist for 2D and 3D problems only
2248  Array<OneD, NekDouble > Sxy(nPts, 0.0);
2249  Array<OneD, NekDouble > Sxz(nPts, 0.0);
2250  Array<OneD, NekDouble > Syz(nPts, 0.0);
2251 
2252  if (m_spacedim == 2)
2253  {
2254  // Sxy = (du/dy + dv/dx)
2255  Vmath::Vadd(nPts, &derivativesO1_interp[0][1][0], 1,
2256  &derivativesO1_interp[1][0][0], 1, &Sxy[0], 1);
2257 
2258  // Sxy = mu * (du/dy + dv/dx)
2259  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
2260  }
2261  else if (m_spacedim == 3)
2262  {
2263  // Sxy = (du/dy + dv/dx)
2264  Vmath::Vadd(nPts, &derivativesO1_interp[0][1][0], 1,
2265  &derivativesO1_interp[1][0][0], 1, &Sxy[0], 1);
2266 
2267  // Sxz = (du/dz + dw/dx)
2268  Vmath::Vadd(nPts, &derivativesO1_interp[0][2][0], 1,
2269  &derivativesO1_interp[2][0][0], 1, &Sxz[0], 1);
2270 
2271  // Syz = (dv/dz + dw/dy)
2272  Vmath::Vadd(nPts, &derivativesO1_interp[1][2][0], 1,
2273  &derivativesO1_interp[2][1][0], 1, &Syz[0], 1);
2274 
2275  // Sxy = mu * (du/dy + dv/dx)
2276  Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
2277 
2278  // Sxz = mu * (du/dy + dv/dx)
2279  Vmath::Vmul(nPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
2280 
2281  // Syz = mu * (du/dy + dv/dx)
2282  Vmath::Vmul(nPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
2283  }
2284 
2285  // Energy-related terms
2286  Array<OneD, NekDouble > STx(nPts, 0.0);
2287  Array<OneD, NekDouble > STy(nPts, 0.0);
2288  Array<OneD, NekDouble > STz(nPts, 0.0);
2289  // Building the viscous flux vector
2290  if (i == 0)
2291  {
2292  // Viscous flux vector for the rho equation
2293  for (k = 0; k < m_spacedim; ++k)
2294  {
2295  Vmath::Zero(nPts, viscousTensor_interp[k][i], 1);
2296  }
2297  }
2298 
2299  if (m_spacedim == 1)
2300  {
2301  Array<OneD, NekDouble > tmp1(nPts, 0.0);
2302 
2303  // u * Sxx
2304  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2305  &Sgg[0][0], 1, &STx[0], 1);
2306 
2307  // k * dT/dx
2308  Vmath::Smul(nPts, m_thermalConductivity,
2309  &derivativesO1_interp[0][1][0], 1, &tmp1[0], 1);
2310 
2311  // STx = u * Sxx + (K / mu) * dT/dx
2312  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
2313  }
2314  else if (m_spacedim == 2)
2315  {
2316  Array<OneD, NekDouble > tmp1(nPts, 0.0);
2317  Array<OneD, NekDouble > tmp2(nPts, 0.0);
2318 
2319  // Computation of STx
2320 
2321  // u * Sxx
2322  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2323  &Sgg[0][0], 1, &STx[0], 1);
2324 
2325  // v * Sxy
2326  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2327  &Sxy[0], 1, &tmp1[0], 1);
2328 
2329  // k * dT/dx
2330  Vmath::Smul(nPts, m_thermalConductivity,
2331  &derivativesO1_interp[0][2][0], 1, &tmp2[0], 1);
2332 
2333  // STx = u * Sxx + v * Sxy + K * dT/dx
2334  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
2335  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
2336 
2337  // Computation of STy
2338 
2339  // Re-initialise temporary arrays
2340  Vmath::Zero(nPts, &tmp1[0], 1);
2341  Vmath::Zero(nPts, &tmp2[0], 1);
2342 
2343  // v * Syy
2344  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2345  &Sgg[1][0], 1, &STy[0], 1);
2346 
2347  // u * Sxy
2348  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2349  &Sxy[0], 1, &tmp1[0], 1);
2350 
2351  // k * dT/dy
2352  Vmath::Smul(nPts, m_thermalConductivity,
2353  &derivativesO1_interp[1][2][0], 1, &tmp2[0], 1);
2354 
2355  // STy = v * Syy + u * Sxy + K * dT/dy
2356  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
2357  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
2358  }
2359  else if (m_spacedim == 3)
2360  {
2361  Array<OneD, NekDouble > tmp1(nPts, 0.0);
2362  Array<OneD, NekDouble > tmp2(nPts, 0.0);
2363  Array<OneD, NekDouble > tmp3(nPts, 0.0);
2364 
2365  // Computation of STx
2366 
2367  // u * Sxx
2368  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2369  &Sgg[0][0], 1, &STx[0], 1);
2370 
2371  // v * Sxy
2372  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2373  &Sxy[0], 1, &tmp1[0], 1);
2374 
2375  // v * Sxy
2376  Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
2377  &Sxz[0], 1, &tmp2[0], 1);
2378 
2379  // k * dT/dx
2380  Vmath::Smul(nPts, m_thermalConductivity,
2381  &derivativesO1_interp[0][3][0], 1, &tmp3[0], 1);
2382 
2383  // STx = u * Sxx + v * Sxy + w * Sxz + (K / mu) * dT/dx
2384  Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
2385  Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
2386  Vmath::Vadd(nPts, &STx[0], 1, &tmp3[0], 1, &STx[0], 1);
2387 
2388  // Computation of STy
2389 
2390  // Re-initialise temporary arrays
2391  Vmath::Zero(nPts, &tmp1[0], 1);
2392  Vmath::Zero(nPts, &tmp2[0], 1);
2393  Vmath::Zero(nPts, &tmp3[0], 1);
2394 
2395  // v * Syy
2396  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2397  &Sgg[1][0], 1, &STy[0], 1);
2398 
2399  // u * Sxy
2400  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2401  &Sxy[0], 1, &tmp1[0], 1);
2402 
2403  // w * Syz
2404  Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
2405  &Syz[0], 1, &tmp2[0], 1);
2406 
2407  // k * dT/dy
2408  Vmath::Smul(nPts, m_thermalConductivity,
2409  &derivativesO1_interp[1][3][0], 1, &tmp3[0], 1);
2410 
2411  // STy = v * Syy + u * Sxy + w * Syz + K * dT/dy
2412  Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
2413  Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
2414  Vmath::Vadd(nPts, &STy[0], 1, &tmp3[0], 1, &STy[0], 1);
2415 
2416  // Computation of STz
2417 
2418  // Re-initialise temporary arrays
2419  Vmath::Zero(nPts, &tmp1[0], 1);
2420  Vmath::Zero(nPts, &tmp2[0], 1);
2421  Vmath::Zero(nPts, &tmp3[0], 1);
2422 
2423  // w * Szz
2424  Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
2425  &Sgg[2][0], 1, &STz[0], 1);
2426 
2427  // u * Sxz
2428  Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
2429  &Sxz[0], 1, &tmp1[0], 1);
2430 
2431  // v * Syz
2432  Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
2433  &Syz[0], 1, &tmp2[0], 1);
2434 
2435  // k * dT/dz
2436  Vmath::Smul(nPts, m_thermalConductivity,
2437  &derivativesO1_interp[2][3][0], 1, &tmp3[0], 1);
2438 
2439  // STz = w * Szz + u * Sxz + v * Syz + K * dT/dz
2440  Vmath::Vadd(nPts, &STz[0], 1, &tmp1[0], 1, &STz[0], 1);
2441  Vmath::Vadd(nPts, &STz[0], 1, &tmp2[0], 1, &STz[0], 1);
2442  Vmath::Vadd(nPts, &STz[0], 1, &tmp3[0], 1, &STz[0], 1);
2443  }
2444 
2445  switch (m_spacedim)
2446  {
2447  case 1:
2448  {
2449 
2450  int nq = physfield[0].num_elements();
2451  // f_11v = f_rho = 0
2452  Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
2453 
2454  // f_21v = f_rhou
2455  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
2456 
2457  // f_31v = f_E
2458  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][2][0], 1);
2459  break;
2460  }
2461  case 2:
2462  {
2463  int nq = physfield[0].num_elements();
2464  // f_11v = f_rho1 = 0
2465  Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
2466  // f_12v = f_rho2 = 0
2467  Vmath::Zero(nq, &viscousTensor_interp[1][0][0], 1);
2468 
2469  // f_21v = f_rhou1
2470  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
2471  // f_22v = f_rhou2
2472  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[1][1][0], 1);
2473 
2474  // f_31v = f_rhov1
2475  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[0][2][0], 1);
2476  // f_32v = f_rhov2
2477  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor_interp[1][2][0], 1);
2478 
2479  // f_41v = f_E1
2480  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][3][0], 1);
2481  // f_42v = f_E2
2482  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor_interp[1][3][0], 1);
2483  break;
2484  }
2485  case 3:
2486  {
2487  int nq = physfield[0].num_elements();
2488  // f_11v = f_rho1 = 0
2489  Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
2490  // f_12v = f_rho2 = 0
2491  Vmath::Zero(nq, &viscousTensor_interp[1][0][0], 1);
2492  // f_13v = f_rho3 = 0
2493  Vmath::Zero(nq, &viscousTensor_interp[2][0][0], 1);
2494 
2495  // f_21v = f_rhou1
2496  Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
2497  // f_22v = f_rhou2
2498  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[1][1][0], 1);
2499  // f_23v = f_rhou3
2500  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor_interp[2][1][0], 1);
2501 
2502  // f_31v = f_rhov1
2503  Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[0][2][0], 1);
2504  // f_32v = f_rhov2
2505  Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor_interp[1][2][0], 1);
2506  // f_33v = f_rhov3
2507  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor_interp[2][2][0], 1);
2508 
2509  // f_31v = f_rhow1
2510  Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor_interp[0][3][0], 1);
2511  // f_32v = f_rhow2
2512  Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor_interp[1][3][0], 1);
2513  // f_33v = f_rhow3
2514  Vmath::Vcopy(nPts, &Sgg[2][0], 1, &viscousTensor_interp[2][3][0], 1);
2515 
2516  // f_41v = f_E1
2517  Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][4][0], 1);
2518  // f_42v = f_E2
2519  Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor_interp[1][4][0], 1);
2520  // f_43v = f_E3
2521  Vmath::Vcopy(nPts, &STz[0], 1, &viscousTensor_interp[2][4][0], 1);
2522  break;
2523  }
2524  default:
2525  {
2526  ASSERTL0(false, "Illegal expansion dimension");
2527  }
2528  }
2529 
2530  for (i = 0; i < m_spacedim; ++i)
2531  {
2532  for (j = 1; j < nVariables; ++j)
2533  {
2534  // Galerkin project solution back to origianl space
2535  m_fields[0]->PhysGalerkinProjection1DScaled(
2536  OneDptscale,
2537  viscousTensor_interp[i][j],
2538  viscousTensor[i][j]);
2539  }
2540  }
2541 #endif
2542 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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 > > &  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 1411 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().

1415  {
1416  int i, j;
1417  int nTracePts = GetTraceTotPoints();
1418  int nVariables = physarray.num_elements();
1419  int nDimensions = m_spacedim;
1420 
1421  const Array<OneD, const int> &traceBndMap
1422  = m_fields[0]->GetTraceBndMap();
1423 
1424  NekDouble gamma = m_gamma;
1425  NekDouble gammaMinusOne = gamma - 1.0;
1426  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
1427 
1428  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
1429  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
1430  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
1431  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
1432 
1433  // Computing the normal velocity for characteristics coming
1434  // from outside the computational domain
1435  velInf[0] = m_uInf;
1436  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
1437  if (nDimensions == 2 || nDimensions == 3)
1438  {
1439  velInf[1] = m_vInf;
1440  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
1441  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
1442  }
1443  if (nDimensions == 3)
1444  {
1445  velInf[2] = m_wInf;
1446  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
1447  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
1448  }
1449 
1450  // Get physical values of the forward trace
1451  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
1452 
1453  for (i = 0; i < nVariables; ++i)
1454  {
1455  Fwd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1456  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
1457  }
1458 
1459  // Computing the normal velocity for characteristics coming
1460  // from inside the computational domain
1461  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1462  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
1463  for (i = 0; i < nDimensions; ++i)
1464  {
1465  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
1466  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
1467  }
1468 
1469  // Computing the absolute value of the velocity in order to compute the
1470  // Mach number to decide whether supersonic or subsonic
1471  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
1472  for (i = 0; i < nDimensions; ++i)
1473  {
1474  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
1475  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
1476  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
1477  }
1478  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
1479 
1480  // Get speed of sound
1481  Array<OneD, NekDouble > soundSpeed (nTracePts, 0.0);
1482  Array<OneD, NekDouble > pressure (nTracePts, 0.0);
1483 
1484  for (i = 0; i < nTracePts; i++)
1485  {
1486  if (m_spacedim == 1)
1487  {
1488  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
1489  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
1490  }
1491  else if (m_spacedim == 2)
1492  {
1493  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
1494  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1495  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
1496  }
1497  else
1498  {
1499  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
1500  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1501  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
1502  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
1503  }
1504 
1505  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
1506  }
1507 
1508  // Get Mach
1509  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
1510  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
1511  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
1512 
1513  // Auxiliary variables
1514  int e, id1, id2, npts, pnt;
1515  NekDouble rhoeb;
1516 
1517  // Loop on the bcRegions
1518  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->
1519  GetExpSize(); ++e)
1520  {
1521  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1522  GetExp(e)->GetTotPoints();
1523  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1524  GetPhys_Offset(e);
1525  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1526 
1527  // Loop on points of bcRegion 'e'
1528  for (i = 0; i < npts; i++)
1529  {
1530  pnt = id2+i;
1531 
1532  // Subsonic flows
1533  if (Mach[pnt] < 0.99)
1534  {
1535  // Partial extrapolation for subsonic cases
1536  for (j = 0; j < nVariables-1; ++j)
1537  {
1538  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1539  UpdatePhys())[id1+i] = m_fieldStorage[j][id1+i];
1540  }
1541 
1542  // Kinetic energy calculation
1543  NekDouble Ek = 0.0;
1544  for (j = 1; j < nVariables-1; ++j)
1545  {
1546  Ek += 0.5 * (m_fieldStorage[j][id1+i] *
1547  m_fieldStorage[j][id1+i]) /
1548  m_fieldStorage[0][id1+i];
1549  }
1550  rhoeb = gammaMinusOneInv * pressure[pnt] + Ek;
1551 
1552  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1553  UpdatePhys())[id1+i] = rhoeb;
1554  }
1555  // Supersonic flows
1556  else
1557  {
1558  for (j = 0; j < nVariables; ++j)
1559  {
1560  // Extrapolation for supersonic cases
1561  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1562  UpdatePhys())[id1+i] = Fwd[j][pnt];
1563  }
1564  }
1565  }
1566  }
1567  }
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 > > &  physarray 
)
protected

Pressure outflow boundary conditions for compressible flow problems.

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

1090  {
1091  int i, j;
1092  int nTracePts = GetTraceTotPoints();
1093  int nVariables = physarray.num_elements();
1094  int nDimensions = m_spacedim;
1095 
1096  const Array<OneD, const int> &traceBndMap
1097  = m_fields[0]->GetTraceBndMap();
1098 
1099  NekDouble gamma = m_gamma;
1100  NekDouble gammaMinusOne = gamma - 1.0;
1101  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
1102 
1103  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
1104  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
1105  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
1106  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
1107 
1108  // Computing the normal velocity for characteristics coming
1109  // from outside the computational domain
1110  velInf[0] = m_uInf;
1111  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
1112  if (nDimensions == 2 || nDimensions == 3)
1113  {
1114  velInf[1] = m_vInf;
1115  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
1116  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
1117  }
1118  if (nDimensions == 3)
1119  {
1120  velInf[2] = m_wInf;
1121  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
1122  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
1123  }
1124 
1125  // Get physical values of the forward trace
1126  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
1127  for (i = 0; i < nVariables; ++i)
1128  {
1129  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
1130  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
1131  }
1132 
1133  // Computing the normal velocity for characteristics coming
1134  // from inside the computational domain
1135  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1136  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
1137  for (i = 0; i < nDimensions; ++i)
1138  {
1139  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
1140  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
1141  }
1142 
1143  // Computing the absolute value of the velocity in order to compute the
1144  // Mach number to decide whether supersonic or subsonic
1145  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
1146  for (i = 0; i < nDimensions; ++i)
1147  {
1148  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
1149  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
1150  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
1151  }
1152  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
1153 
1154  // Get speed of sound
1155  Array<OneD, NekDouble > soundSpeed(nTracePts);
1156  Array<OneD, NekDouble > pressure (nTracePts);
1157 
1158  for (i = 0; i < nTracePts; i++)
1159  {
1160  if (m_spacedim == 1)
1161  {
1162  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
1163  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
1164  }
1165  else if (m_spacedim == 2)
1166  {
1167  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
1168  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1169  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
1170  }
1171  else
1172  {
1173  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
1174  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1175  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
1176  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
1177  }
1178 
1179  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
1180  }
1181 
1182  // Get Mach
1183  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
1184  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
1185  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
1186 
1187  // Auxiliary variables
1188  int e, id1, id2, npts, pnt;
1189  NekDouble rhoeb;
1190 
1191  // Loop on the bcRegions
1192  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->
1193  GetExpSize(); ++e)
1194  {
1195  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1196  GetExp(e)->GetTotPoints();
1197  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1198  GetPhys_Offset(e) ;
1199  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1200 
1201  // Loop on points of bcRegion 'e'
1202  for (i = 0; i < npts; i++)
1203  {
1204  pnt = id2+i;
1205 
1206  // Subsonic flows
1207  if (Mach[pnt] < 0.99)
1208  {
1209  // Kinetic energy calculation
1210  NekDouble Ek = 0.0;
1211  for (j = 1; j < nVariables-1; ++j)
1212  {
1213  Ek += 0.5 * (Fwd[j][pnt] * Fwd[j][pnt]) / Fwd[0][pnt];
1214  }
1215 
1216  rhoeb = m_pInf * gammaMinusOneInv + Ek;
1217 
1218  // Partial extrapolation for subsonic cases
1219  for (j = 0; j < nVariables-1; ++j)
1220  {
1221  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1222  UpdatePhys())[id1+i] = Fwd[j][pnt];
1223  }
1224 
1225  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1226  UpdatePhys())[id1+i] = rhoeb;
1227  }
1228  // Supersonic flows
1229  else
1230  {
1231  for (j = 0; j < nVariables; ++j)
1232  {
1233  // Extrapolation for supersonic cases
1234  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1235  UpdatePhys())[id1+i] = Fwd[j][pnt];
1236  }
1237  }
1238  }
1239  }
1240  }
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 > > &  physarray 
)
protected

Pressure outflow boundary conditions for compressible flow problems.

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

1251  {
1252  int i, j;
1253  int nTracePts = GetTraceTotPoints();
1254  int nVariables = physarray.num_elements();
1255  int nDimensions = m_spacedim;
1256 
1257  const Array<OneD, const int> &traceBndMap
1258  = m_fields[0]->GetTraceBndMap();
1259 
1260  NekDouble gamma = m_gamma;
1261  NekDouble gammaMinusOne = gamma - 1.0;
1262  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
1263 
1264  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
1265  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
1266  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
1267  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
1268 
1269  // Computing the normal velocity for characteristics coming
1270  // from outside the computational domain
1271  velInf[0] = m_uInf;
1272  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
1273  if (nDimensions == 2 || nDimensions == 3)
1274  {
1275  velInf[1] = m_vInf;
1276  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
1277  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
1278  }
1279  if (nDimensions == 3)
1280  {
1281  velInf[2] = m_wInf;
1282  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
1283  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
1284  }
1285 
1286  // Get physical values of the forward trace
1287  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
1288 
1289  for (i = 0; i < nVariables; ++i)
1290  {
1291  Fwd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1292  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
1293  }
1294 
1295  // Computing the normal velocity for characteristics coming
1296  // from inside the computational domain
1297  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1298  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
1299  for (i = 0; i < nDimensions; ++i)
1300  {
1301  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
1302  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
1303  }
1304 
1305  // Computing the absolute value of the velocity in order to compute the
1306  // Mach number to decide whether supersonic or subsonic
1307  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
1308  for (i = 0; i < nDimensions; ++i)
1309  {
1310  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
1311  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
1312  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
1313  }
1314  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
1315 
1316  // Get speed of sound
1317  Array<OneD, NekDouble > soundSpeed (nTracePts, 0.0);
1318  Array<OneD, NekDouble > pressure (nTracePts, 0.0);
1319 
1320  for (i = 0; i < nTracePts; i++)
1321  {
1322  if (m_spacedim == 1)
1323  {
1324  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
1325  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
1326  }
1327  else if (m_spacedim == 2)
1328  {
1329  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
1330  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1331  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
1332  }
1333  else
1334  {
1335  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
1336  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1337  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
1338  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
1339  }
1340 
1341  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
1342  }
1343 
1344  // Get Mach
1345  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
1346  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
1347  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
1348 
1349  // Auxiliary variables
1350  int e, id1, id2, npts, pnt;
1351  NekDouble rhoeb;
1352 
1353  // Loop on the bcRegions
1354  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->
1355  GetExpSize(); ++e)
1356  {
1357  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1358  GetExp(e)->GetTotPoints();
1359  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1360  GetPhys_Offset(e);
1361  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1362 
1363  // Loop on points of bcRegion 'e'
1364  for (i = 0; i < npts; i++)
1365  {
1366  pnt = id2+i;
1367 
1368  // Subsonic flows
1369  if (Mach[pnt] < 0.99)
1370  {
1371  // Kinetic energy calculation
1372  NekDouble Ek = 0.0;
1373  for (j = 1; j < nVariables-1; ++j)
1374  {
1375  Ek += 0.5 * (Fwd[j][pnt] * Fwd[j][pnt]) / Fwd[0][pnt];
1376  }
1377 
1378  rhoeb = m_pressureStorage[id1+i] * gammaMinusOneInv + Ek;
1379 
1380  // Partial extrapolation for subsonic cases
1381  for (j = 0; j < nVariables-1; ++j)
1382  {
1383  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1384  UpdatePhys())[id1+i] = Fwd[j][pnt];
1385  }
1386 
1387  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1388  UpdatePhys())[id1+i] = rhoeb;
1389  }
1390  // Supersonic flows
1391  else
1392  {
1393  for (j = 0; j < nVariables; ++j)
1394  {
1395  // Extrapolation for supersonic cases
1396  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1397  UpdatePhys())[id1+i] = Fwd[j][pnt];
1398  }
1399  }
1400  }
1401  }
1402  }
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 > > &  physarray 
)
protected

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

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

930  {
931  int i, j;
932  int nTracePts = GetTraceTotPoints();
933  int nVariables = physarray.num_elements();
934  int nDimensions = m_spacedim;
935 
936  const Array<OneD, const int> &traceBndMap
937  = m_fields[0]->GetTraceBndMap();
938 
939  NekDouble gamma = m_gamma;
940  NekDouble gammaMinusOne = gamma - 1.0;
941  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
942 
943  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
944  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
945  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
946  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
947 
948  // Computing the normal velocity for characteristics coming
949  // from outside the computational domain
950  velInf[0] = m_uInf;
951  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
952  if (nDimensions == 2 || nDimensions == 3)
953  {
954  velInf[1] = m_vInf;
955  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
956  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
957  }
958  if (nDimensions == 3)
959  {
960  velInf[2] = m_wInf;
961  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
962  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
963  }
964 
965  // Get physical values of the forward trace
966  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
967  for (i = 0; i < nVariables; ++i)
968  {
969  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
970  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
971  }
972 
973  // Computing the normal velocity for characteristics coming
974  // from inside the computational domain
975  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
976  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
977  for (i = 0; i < nDimensions; ++i)
978  {
979  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
980  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
981  }
982 
983  // Computing the absolute value of the velocity in order to compute the
984  // Mach number to decide whether supersonic or subsonic
985  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
986  for (i = 0; i < nDimensions; ++i)
987  {
988  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
989  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
990  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
991  }
992  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
993 
994  // Get speed of sound
995  Array<OneD, NekDouble > soundSpeed(nTracePts);
996  Array<OneD, NekDouble > pressure (nTracePts);
997 
998  for (i = 0; i < nTracePts; i++)
999  {
1000  if (m_spacedim == 1)
1001  {
1002  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
1003  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
1004  }
1005  else if (m_spacedim == 2)
1006  {
1007  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
1008  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1009  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
1010  }
1011  else
1012  {
1013  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
1014  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
1015  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
1016  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
1017  }
1018 
1019  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
1020  }
1021 
1022  // Get Mach
1023  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
1024  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
1025  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
1026 
1027  // Auxiliary variables
1028  int e, id1, id2, npts, pnt;
1029  NekDouble rhoeb;
1030 
1031  // Loop on the bcRegions
1032  for (e = 0; e < m_fields[0]->GetBndCondExpansions()[bcRegion]->
1033  GetExpSize(); ++e)
1034  {
1035  npts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1036  GetExp(e)->GetTotPoints();
1037  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
1038  GetPhys_Offset(e) ;
1039  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
1040 
1041  // Loop on points of bcRegion 'e'
1042  for (i = 0; i < npts; i++)
1043  {
1044  pnt = id2+i;
1045 
1046  // Subsonic flows
1047  if (Mach[pnt] < 0.99)
1048  {
1049  // Kinetic energy calculation
1050  NekDouble Ek = 0.0;
1051  for (j = 1; j < nVariables-1; ++j)
1052  {
1053  Ek += 0.5 * (Fwd[j][pnt] * Fwd[j][pnt]) / Fwd[0][pnt];
1054  }
1055 
1056  rhoeb = m_pInf * gammaMinusOneInv + Ek;
1057 
1058  // Partial extrapolation for subsonic cases
1059  for (j = 0; j < nVariables-1; ++j)
1060  {
1061  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1062  UpdatePhys())[id1+i] = Fwd[j][pnt];
1063  }
1064 
1065  (m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
1066  UpdatePhys())[id1+i] = 2.0 * rhoeb - Fwd[nVariables-1][pnt];
1067  }
1068  // Supersonic flows
1069  else
1070  {
1071  for (j = 0; j < nVariables; ++j)
1072  {
1073  // Extrapolation for supersonic cases
1074  (m_fields[j]->GetBndCondExpansions()[bcRegion]->
1075  UpdatePhys())[id1+i] = Fwd[j][pnt];
1076  }
1077  }
1078  }
1079  }
1080  }
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 > > &  physarray 
)
protected

Outflow characteristic boundary conditions for compressible flow problems.

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

673  {
674  int i, j;
675  int nTracePts = GetTraceTotPoints();
676  int nVariables = physarray.num_elements();
677  int nDimensions = m_spacedim;
678 
679  const Array<OneD, const int> &traceBndMap
680  = m_fields[0]->GetTraceBndMap();
681 
682  NekDouble gamma = m_gamma;
683  NekDouble gammaInv = 1.0 / gamma;
684  NekDouble gammaMinusOne = gamma - 1.0;
685  NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
686 
687  Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
688  Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
689  Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
690  Array<OneD, NekDouble> velInf(nDimensions, 0.0);
691 
692  // Computing the normal velocity for characteristics coming
693  // from outside the computational domain
694  velInf[0] = m_uInf;
695  Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
696  if (nDimensions == 2 || nDimensions == 3)
697  {
698  velInf[1] = m_vInf;
699  Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
700  Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
701  }
702  if (nDimensions == 3)
703  {
704  velInf[2] = m_wInf;
705  Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
706  Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
707  }
708 
709  // Get physical values of the forward trace
710  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
711  for (i = 0; i < nVariables; ++i)
712  {
713  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
714  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
715  }
716 
717  // Computing the normal velocity for characteristics coming
718  // from inside the computational domain
719  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
720  Array<OneD, NekDouble > Vel(nTracePts, 0.0);
721  for (i = 0; i < nDimensions; ++i)
722  {
723  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
724  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
725  }
726 
727  // Computing the absolute value of the velocity in order to compute the
728  // Mach number to decide whether supersonic or subsonic
729  Array<OneD, NekDouble > absVel(nTracePts, 0.0);
730  for (i = 0; i < nDimensions; ++i)
731  {
732  Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
733  Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
734  Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
735  }
736  Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
737 
738  // Get speed of sound
739  Array<OneD, NekDouble > soundSpeed(nTracePts);
740  Array<OneD, NekDouble > pressure (nTracePts);
741 
742  for (i = 0; i < nTracePts; i++)
743  {
744  if (m_spacedim == 1)
745  {
746  pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
747  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
748  }
749  else if (m_spacedim == 2)
750  {
751  pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
752  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
753  Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
754  }
755  else
756  {
757  pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
758  0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
759  Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
760  Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
761  }
762 
763  soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
764  }
765 
766  // Get Mach
767  Array<OneD, NekDouble > Mach(nTracePts, 0.0);
768  Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
769  Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
770 
771  // Auxiliary variables
772  int eMax;
773  int e, id1, id2, nBCEdgePts, pnt;
774  NekDouble cPlus, rPlus, cMinus, rMinus, VDBC, VNBC;
775  Array<OneD, NekDouble> velBC(nDimensions, 0.0);
776  Array<OneD, NekDouble> rhoVelBC(nDimensions, 0.0);
777  NekDouble rhoBC, EBC, cBC, sBC, pBC;
778 
779  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
780 
781  // Loop on bcRegions
782  for (e = 0; e < eMax; ++e)
783  {
784  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
785  GetExp(e)->GetTotPoints();
786 
787  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
788  GetPhys_Offset(e);
789  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
790 
791  // Loop on the points of the bcRegion
792  for (i = 0; i < nBCEdgePts; i++)
793  {
794  pnt = id2+i;
795 
796  // Impose inflow Riemann invariant
797  if (Vn[pnt] <= 0.0)
798  {
799  // Subsonic flows
800  if (Mach[pnt] < 1.00)
801  {
802  // + Characteristic from inside
803  cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
804  rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
805 
806  // - Characteristic from boundary
807  cMinus = sqrt(gamma * m_pInf / m_rhoInf);
808  rMinus = VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
809  }
810  else
811  {
812  // + Characteristic from inside
813  cPlus = sqrt(gamma * m_pInf / m_rhoInf);
814  rPlus = VnInf[pnt] + 2.0 * cPlus * gammaMinusOneInv;
815 
816  // + Characteristic from inside
817  cMinus = sqrt(gamma * m_pInf / m_rhoInf);
818  rMinus = VnInf[pnt] - 2.0 * cPlus * gammaMinusOneInv;
819  }
820 
821  // Riemann boundary variables
822  VNBC = 0.5 * (rPlus + rMinus);
823  cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
824  VDBC = VNBC - VnInf[pnt];
825 
826  // Thermodynamic boundary variables
827  sBC = m_pInf / (pow(m_rhoInf, gamma));
828  rhoBC = pow((cBC * cBC) / (gamma * sBC), gammaMinusOneInv);
829  pBC = rhoBC * cBC * cBC * gammaInv;
830 
831  // Kinetic energy initialiasation
832  NekDouble EkBC = 0.0;
833 
834  // Boundary velocities
835  for ( j = 0; j < nDimensions; ++j)
836  {
837  velBC[j] = velInf[j] + VDBC * m_traceNormals[j][pnt];
838  rhoVelBC[j] = rhoBC * velBC[j];
839  EkBC += 0.5 * rhoBC * velBC[j]*velBC[j];
840  }
841 
842  // Boundary energy
843  EBC = pBC * gammaMinusOneInv + EkBC;
844 
845  // Imposing Riemann Invariant boundary conditions
846  (m_fields[0]->GetBndCondExpansions()[bcRegion]->
847  UpdatePhys())[id1+i] = rhoBC;
848  for (j = 0; j < nDimensions; ++j)
849  {
850  (m_fields[j+1]->GetBndCondExpansions()[bcRegion]->
851  UpdatePhys())[id1+i] = rhoVelBC[j];
852  }
853  (m_fields[nDimensions+1]->GetBndCondExpansions()[bcRegion]->
854  UpdatePhys())[id1+i] = EBC;
855 
856  }
857  else // Impose outflow Riemann invariant
858  {
859  // Subsonic flows
860  if (Mach[pnt] < 1.00)
861  {
862  // + Characteristic from inside
863  cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
864  rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
865 
866  // - Characteristic from boundary
867  cMinus = sqrt(gamma * m_pInf / m_rhoInf);
868  rMinus = VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
869  }
870  else
871  {
872  // + Characteristic from inside
873  cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
874  rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
875 
876  // + Characteristic from inside
877  cMinus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
878  rMinus = Vn[pnt] - 2.0 * cPlus * gammaMinusOneInv;
879  }
880 
881  // Riemann boundary variables
882  VNBC = 0.5 * (rPlus + rMinus);
883  cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
884  VDBC = VNBC - Vn[pnt];
885 
886  // Thermodynamic boundary variables
887  sBC = pressure[pnt] / (pow(Fwd[0][pnt], gamma));
888  rhoBC = pow((cBC * cBC) / (gamma * sBC), gammaMinusOneInv);
889  pBC = rhoBC * cBC * cBC * gammaInv;
890 
891  // Kinetic energy initialiasation
892  NekDouble EkBC = 0.0;
893 
894  // Boundary velocities
895  for ( j = 0; j < nDimensions; ++j)
896  {
897  velBC[j] = Fwd[j+1][pnt] / Fwd[0][pnt] +
898  VDBC * m_traceNormals[j][pnt];
899  rhoVelBC[j] = rhoBC * velBC[j];
900  EkBC += 0.5 * rhoBC * velBC[j]*velBC[j];
901  }
902 
903  // Boundary energy
904  EBC = pBC * gammaMinusOneInv + EkBC;
905 
906  // Imposing Riemann Invariant boundary conditions
907  (m_fields[0]->GetBndCondExpansions()[bcRegion]->
908  UpdatePhys())[id1+i] = rhoBC;
909  for (j = 0; j < nDimensions; ++j)
910  {
911  (m_fields[j+1]->GetBndCondExpansions()[bcRegion]->
912  UpdatePhys())[id1+i] = rhoVelBC[j];
913  }
914  (m_fields[nDimensions+1]->GetBndCondExpansions()[bcRegion]->
915  UpdatePhys())[id1+i] = EBC;
916  }
917  }
918  }
919  }
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 > > &  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().

330  {
331  std::string varName;
332  int nvariables = m_fields.num_elements();
333 
334  if(!userDefStr.empty())
335  {
336  if(boost::iequals(userDefStr,"Wall"))
337  {
338  WallBC(n, cnt, inarray);
339  }
340  else if(boost::iequals(userDefStr,"WallViscous") ||
341  boost::iequals(userDefStr,"WallAdiabatic"))
342  {
343  // Wall Boundary Condition
344  WallViscousBC(n, cnt, inarray);
345  }
346  else if(boost::iequals(userDefStr,"Symmetry"))
347  {
348  // Symmetric Boundary Condition
349  SymmetryBC(n, cnt, inarray);
350  }
351  else if(boost::iequals(userDefStr,"RiemannInvariant"))
352  {
353  // Riemann invariant characteristic Boundary Condition
354  RiemannInvariantBC(n, cnt, inarray);
355  }
356  else if(boost::iequals(userDefStr,"PressureOutflowNonReflective"))
357  {
358  // Pressure outflow non-reflective Boundary Condition
359  PressureOutflowNonReflectiveBC(n, cnt, inarray);
360  }
361  else if(boost::iequals(userDefStr,"PressureOutflow"))
362  {
363  // Pressure outflow Boundary Condition
364  PressureOutflowBC(n, cnt, inarray);
365  }
366  else if(boost::iequals(userDefStr,"PressureOutflowFile"))
367  {
368  // Pressure outflow Boundary Condition from file
369  PressureOutflowFileBC(n, cnt, inarray);
370  }
371  else if(boost::iequals(userDefStr,"PressureInflowFile"))
372  {
373  // Pressure inflow Boundary Condition from file
374  PressureInflowFileBC(n, cnt, inarray);
375  }
376  else if(boost::iequals(userDefStr,"ExtrapOrder0"))
377  {
378  // Extrapolation of the data at the boundaries
379  ExtrapOrder0BC(n, cnt, inarray);
380  }
381  else if(boost::iequals(userDefStr,"TimeDependent"))
382  {
383  for (int i = 0; i < nvariables; ++i)
384  {
385  varName = m_session->GetVariable(i);
386  m_fields[i]->EvaluateBoundaryConditions(time, varName);
387  }
388  }
389  else
390  {
391  string errmsg = "Unrecognised boundary condition: ";
392  errmsg += userDefStr;
393  ASSERTL0(false,errmsg.c_str());
394  }
395  }
396  }
Nektar::CompressibleFlowSystem::ExtrapOrder0BC
void ExtrapOrder0BC(int bcRegion, int cnt, 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:1574
Nektar::CompressibleFlowSystem::RiemannInvariantBC
void RiemannInvariantBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Outflow characteristic boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:669
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::CompressibleFlowSystem::SymmetryBC
void SymmetryBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Symmetry boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:573
Nektar::CompressibleFlowSystem::PressureOutflowNonReflectiveBC
void PressureOutflowNonReflectiveBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Pressure outflow non-reflective boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:926
Nektar::CompressibleFlowSystem::PressureOutflowFileBC
void PressureOutflowFileBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Pressure outflow boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:1247
Nektar::CompressibleFlowSystem::PressureOutflowBC
void PressureOutflowBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Pressure outflow boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:1086
Nektar::CompressibleFlowSystem::WallViscousBC
void WallViscousBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary conditions for viscous compressible flow problems.
Definition: CompressibleFlowSystem.cpp:498
Nektar::CompressibleFlowSystem::WallBC
void WallBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary conditions for compressible flow problems.
Definition: CompressibleFlowSystem.cpp:401
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::PressureInflowFileBC
void PressureInflowFileBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Pressure inflow boundary conditions for compressible flow problems where either the density and the v...
Definition: CompressibleFlowSystem.cpp:1411
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 > > &  physarray 
)
protected

Symmetry boundary conditions for compressible flow problems.

Definition at line 573 of file CompressibleFlowSystem.cpp.

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

577  {
578  int i;
579  int nTracePts = GetTraceTotPoints();
580  int nVariables = physarray.num_elements();
581 
582  const Array<OneD, const int> &traceBndMap
583  = m_fields[0]->GetTraceBndMap();
584 
585  // Get physical values of the forward trace (from exp to phys)
586  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
587  for (i = 0; i < nVariables; ++i)
588  {
589  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
590  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
591  }
592 
593  // Take into account that for PDE based shock capturing, eps = 0 at the
594  // wall.
595  int e, id1, id2, nBCEdgePts, eMax;
596 
597  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
598 
599  for (e = 0; e < eMax; ++e)
600  {
601  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
602  GetExp(e)->GetTotPoints();
603  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
604  GetPhys_Offset(e);
605  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
606 
607  if (nVariables == m_spacedim+3)
608  {
609  NekDouble factor = 0.0;
610  NekDouble factor2 = 1.0;
611 
612  Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
613  Vmath::Smul(nBCEdgePts,
614  factor,
615  &Fwd[nVariables-1][id2], 1,
616  &tmp2[0], 1);
617 
618  Vmath::Vsub(nBCEdgePts,
619  &Fwd[nVariables-1][id2], 1,
620  &tmp2[0], 1,
621  &Fwd[nVariables-1][id2], 1);
622 
623  Vmath::Smul(nBCEdgePts,
624  factor2,
625  &Fwd[nVariables-1][id2], 1,
626  &Fwd[nVariables-1][id2], 1);
627  }
628 
629  // For 2D/3D, define: v* = v - 2(v.n)n
630  Array<OneD, NekDouble> tmp(nBCEdgePts, 0.0);
631 
632  // Calculate (v.n)
633  for (i = 0; i < m_spacedim; ++i)
634  {
635  Vmath::Vvtvp(nBCEdgePts,
636  &Fwd[1+i][id2], 1,
637  &m_traceNormals[i][id2], 1,
638  &tmp[0], 1,
639  &tmp[0], 1);
640  }
641 
642  // Calculate 2.0(v.n)
643  Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
644 
645  // Calculate v* = v - 2.0(v.n)n
646  for (i = 0; i < m_spacedim; ++i)
647  {
648  Vmath::Vvtvp(nBCEdgePts,
649  &tmp[0], 1,
650  &m_traceNormals[i][id2], 1,
651  &Fwd[1+i][id2], 1,
652  &Fwd[1+i][id2], 1);
653  }
654 
655  // Copy boundary adjusted values into the boundary expansion
656  for (i = 0; i < nVariables; ++i)
657  {
658  Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
659  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
660  UpdatePhys())[id1], 1);
661  }
662  }
663  }
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::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::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:161
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:1794
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:2127
Nektar::CompressibleFlowSystem::GetVecLocs
const Array< OneD, const Array< OneD, NekDouble > > & GetVecLocs()
Definition: CompressibleFlowSystem.h:285
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:1703
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:290
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:1640
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 > > &  physarray 
)
protected

Wall boundary conditions for compressible flow problems.

Definition at line 401 of file CompressibleFlowSystem.cpp.

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

405  {
406  int i;
407  int nTracePts = GetTraceTotPoints();
408  int nVariables = physarray.num_elements();
409 
410  const Array<OneD, const int> &traceBndMap
411  = m_fields[0]->GetTraceBndMap();
412 
413  // Get physical values of the forward trace
414  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
415  for (i = 0; i < nVariables; ++i)
416  {
417  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
418  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
419  }
420 
421  // Adjust the physical values of the trace to take
422  // user defined boundaries into account
423  int e, id1, id2, nBCEdgePts, eMax;
424 
425  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
426 
427  for (e = 0; e < eMax; ++e)
428  {
429  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
430  GetExp(e)->GetTotPoints();
431  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
432  GetPhys_Offset(e);
433  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
434 
435  // Boundary condition for epsilon term.
436  if (nVariables == m_spacedim+3)
437  {
438  NekDouble factor = 1.0;
439  NekDouble factor2 = 1.0;
440 
441  Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
442  Vmath::Smul(nBCEdgePts,
443  factor,
444  &Fwd[nVariables-1][id2], 1,
445  &tmp2[0], 1);
446 
447  Vmath::Vsub(nBCEdgePts,
448  &Fwd[nVariables-1][id2], 1,
449  &tmp2[0], 1,
450  &Fwd[nVariables-1][id2], 1);
451 
452  Vmath::Smul(nBCEdgePts,
453  factor2,
454  &Fwd[nVariables-1][id2], 1,
455  &Fwd[nVariables-1][id2], 1);
456 
457  }
458 
459  // For 2D/3D, define: v* = v - 2(v.n)n
460  Array<OneD, NekDouble> tmp(nBCEdgePts, 0.0);
461 
462  // Calculate (v.n)
463  for (i = 0; i < m_spacedim; ++i)
464  {
465  Vmath::Vvtvp(nBCEdgePts,
466  &Fwd[1+i][id2], 1,
467  &m_traceNormals[i][id2], 1,
468  &tmp[0], 1,
469  &tmp[0], 1);
470  }
471 
472  // Calculate 2.0(v.n)
473  Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
474 
475  // Calculate v* = v - 2.0(v.n)n
476  for (i = 0; i < m_spacedim; ++i)
477  {
478  Vmath::Vvtvp(nBCEdgePts,
479  &tmp[0], 1,
480  &m_traceNormals[i][id2], 1,
481  &Fwd[1+i][id2], 1,
482  &Fwd[1+i][id2], 1);
483  }
484 
485  // Copy boundary adjusted values into the boundary expansion
486  for (i = 0; i < nVariables; ++i)
487  {
488  Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
489  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
490  UpdatePhys())[id1], 1);
491  }
492  }
493  }
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::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::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 > > &  physarray 
)
protected

Wall boundary conditions for viscous compressible flow problems.

Definition at line 498 of file CompressibleFlowSystem.cpp.

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

Referenced by SetCommonBC().

502  {
503  int i;
504  int nTracePts = GetTraceTotPoints();
505  int nVariables = physarray.num_elements();
506 
507  const Array<OneD, const int> &traceBndMap
508  = m_fields[0]->GetTraceBndMap();
509 
510  // Get physical values of the forward trace
511  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
512  for (i = 0; i < nVariables; ++i)
513  {
514  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
515  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
516  }
517 
518  // Take into account that for PDE based shock capturing, eps = 0 at the
519  // wall. Adjust the physical values of the trace to take user defined
520  // boundaries into account
521  int e, id1, id2, nBCEdgePts, eMax;
522 
523  eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
524 
525  for (e = 0; e < eMax; ++e)
526  {
527  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
528  GetExp(e)->GetTotPoints();
529  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
530  GetPhys_Offset(e);
531  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
532 
533  if (nVariables == m_spacedim+3)
534  {
535  NekDouble factor = 0.0;
536  NekDouble factor2 = 1.0;
537 
538  Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
539  Vmath::Smul(nBCEdgePts,
540  factor,
541  &Fwd[nVariables-1][id2], 1,
542  &tmp2[0], 1);
543 
544  Vmath::Vsub(nBCEdgePts,
545  &Fwd[nVariables-1][id2], 1,
546  &tmp2[0], 1,
547  &Fwd[nVariables-1][id2], 1);
548 
549  Vmath::Smul(nBCEdgePts,
550  factor2,
551  &Fwd[nVariables-1][id2], 1,
552  &Fwd[nVariables-1][id2], 1);
553  }
554 
555  for (i = 0; i < m_spacedim; i++)
556  {
557  Vmath::Neg(nBCEdgePts, &Fwd[i+1][id2], 1);
558  }
559 
560  // Copy boundary adjusted values into the boundary expansion
561  for (i = 0; i < nVariables; ++i)
562  {
563  Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
564  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
565  UpdatePhys())[id1], 1);
566  }
567  }
568  }
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::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::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 Wed May 4 2016 17:54:23 for Nektar++ by   doxygen 1.8.8