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.
NekDouble GetStabilityLimit (int n)
 Function to calculate the stability limit for DG/CG.
Array< OneD, NekDouble > GetStabilityLimitVector (const Array< OneD, int > &ExpOrder)
 Function to calculate the stability limit for DG/CG (a vector of them).
- Public Member Functions inherited from Nektar::SolverUtils::UnsteadySystem
virtual SOLVER_UTILS_EXPORT ~UnsteadySystem ()
 Destructor.
SOLVER_UTILS_EXPORT NekDouble GetTimeStep (const Array< OneD, const Array< OneD, NekDouble > > &inarray)
 Calculate the larger time-step mantaining the problem stable.
- Public Member Functions inherited from Nektar::SolverUtils::EquationSystem
virtual SOLVER_UTILS_EXPORT ~EquationSystem ()
 Destructor.
SOLVER_UTILS_EXPORT void SetUpTraceNormals (void)
SOLVER_UTILS_EXPORT void InitObject ()
 Initialises the members of this object.
SOLVER_UTILS_EXPORT void DoInitialise ()
 Perform any initialisation necessary before solving the problem.
SOLVER_UTILS_EXPORT void DoSolve ()
 Solve the problem.
SOLVER_UTILS_EXPORT void TransCoeffToPhys ()
 Transform from coefficient to physical space.
SOLVER_UTILS_EXPORT void TransPhysToCoeff ()
 Transform from physical to coefficient space.
SOLVER_UTILS_EXPORT void Output ()
 Perform output operations after solve.
SOLVER_UTILS_EXPORT NekDouble LinfError (unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
 Linf error computation.
SOLVER_UTILS_EXPORT std::string GetSessionName ()
 Get Session name.
SOLVER_UTILS_EXPORT void ResetSessionName (std::string newname)
 Reset Session name.
SOLVER_UTILS_EXPORT
LibUtilities::SessionReaderSharedPtr 		GetSession ()
 Get Session name.
SOLVER_UTILS_EXPORT
MultiRegions::ExpListSharedPtr 		GetPressure ()
 Get pressure field if available.
SOLVER_UTILS_EXPORT void PrintSummary (std::ostream &out)
 Print a summary of parameters and solver characteristics.
SOLVER_UTILS_EXPORT void SetLambda (NekDouble lambda)
 Set parameter m_lambda.
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.
SOLVER_UTILS_EXPORT void EvaluateFunction (std::vector< std::string > pFieldNames, Array< OneD, Array< OneD, NekDouble > > &pFields, const std::string &pName, const int domain=0)
 Populate given fields with the function from session.
SOLVER_UTILS_EXPORT void EvaluateFunction (std::vector< std::string > pFieldNames, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const std::string &pName, const int domain=0)
 Populate given fields with the function from session.
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.
SOLVER_UTILS_EXPORT void InitialiseBaseFlow (Array< OneD, Array< OneD, NekDouble > > &base)
 Perform initialisation of the base flow.
SOLVER_UTILS_EXPORT void SetInitialConditions (NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
 Initialise the data in the dependent fields.
SOLVER_UTILS_EXPORT void EvaluateExactSolution (int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
 Evaluates an exact solution.
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.
SOLVER_UTILS_EXPORT NekDouble L2Error (unsigned int field, bool Normalised=false)
 Compute the L2 error of the fields.
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].
SOLVER_UTILS_EXPORT void WeakAdvectionGreensDivergenceForm (const Array< OneD, Array< OneD, NekDouble > > &F, Array< OneD, NekDouble > &outarray)
 Compute the inner product $ (\nabla \phi \cdot F) $.
SOLVER_UTILS_EXPORT void WeakAdvectionDivergenceForm (const Array< OneD, Array< OneD, NekDouble > > &F, Array< OneD, NekDouble > &outarray)
 Compute the inner product $ (\phi, \nabla \cdot F) $.
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) $.
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.
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.
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.
SOLVER_UTILS_EXPORT void Checkpoint_Output (const int n)
 Write checkpoint file of m_fields.
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.
SOLVER_UTILS_EXPORT void WriteFld (const std::string &outname)
 Write field data to the given filename.
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.
SOLVER_UTILS_EXPORT void ImportFld (const std::string &infile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
 Input field data from the given file.
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.
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.
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.
SOLVER_UTILS_EXPORT void ScanForHistoryPoints ()
 Builds map of which element holds each history point.
SOLVER_UTILS_EXPORT void WriteHistoryData (std::ostream &out)
 Probe each history point and write to file.
SOLVER_UTILS_EXPORT void SessionSummary (SummaryList &vSummary)
 Write out a session summary.
SOLVER_UTILS_EXPORT Array
< OneD,
MultiRegions::ExpListSharedPtr > & 		UpdateFields ()
SOLVER_UTILS_EXPORT
LibUtilities::FieldMetaDataMap & 		UpdateFieldMetaDataMap ()
 Get hold of FieldInfoMap so it can be updated.
SOLVER_UTILS_EXPORT NekDouble GetFinalTime ()
 Return final time.
SOLVER_UTILS_EXPORT int GetNcoeffs ()
SOLVER_UTILS_EXPORT int GetNcoeffs (const int eid)
SOLVER_UTILS_EXPORT int GetNumExpModes ()
SOLVER_UTILS_EXPORT const
Array< OneD, int > 		GetNumExpModesPerExp ()
SOLVER_UTILS_EXPORT int GetNvariables ()
SOLVER_UTILS_EXPORT const
std::string 		GetVariable (unsigned int i)
SOLVER_UTILS_EXPORT int GetTraceTotPoints ()
SOLVER_UTILS_EXPORT int GetTraceNpoints ()
SOLVER_UTILS_EXPORT int GetExpSize ()
SOLVER_UTILS_EXPORT int GetPhys_Offset (int n)
SOLVER_UTILS_EXPORT int GetCoeff_Offset (int n)
SOLVER_UTILS_EXPORT int GetTotPoints ()
SOLVER_UTILS_EXPORT int GetTotPoints (int n)
SOLVER_UTILS_EXPORT int GetNpoints ()
SOLVER_UTILS_EXPORT int GetNumElmVelocity ()
SOLVER_UTILS_EXPORT int GetSteps ()
SOLVER_UTILS_EXPORT NekDouble GetTimeStep ()
SOLVER_UTILS_EXPORT void CopyFromPhysField (const int i, Array< OneD, NekDouble > &output)
SOLVER_UTILS_EXPORT void CopyToPhysField (const int i, Array< OneD, NekDouble > &output)
SOLVER_UTILS_EXPORT void SetStepsToOne ()
SOLVER_UTILS_EXPORT void ZeroPhysFields ()
SOLVER_UTILS_EXPORT void FwdTransFields ()
SOLVER_UTILS_EXPORT void GetFluxVector (const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
SOLVER_UTILS_EXPORT void GetFluxVector (const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &fluxX, Array< OneD, Array< OneD, NekDouble > > &fluxY)
SOLVER_UTILS_EXPORT void GetFluxVector (const int i, const int j, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
SOLVER_UTILS_EXPORT void NumericalFlux (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
SOLVER_UTILS_EXPORT void NumericalFlux (Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxX, Array< OneD, Array< OneD, NekDouble > > &numfluxY)
SOLVER_UTILS_EXPORT void NumFluxforScalar (const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
SOLVER_UTILS_EXPORT void NumFluxforVector (const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
SOLVER_UTILS_EXPORT void SetModifiedBasis (const bool modbasis)
SOLVER_UTILS_EXPORT int NoCaseStringCompare (const string &s1, const string &s2)
 Perform a case-insensitive string comparison.

Static Public Member Functions

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

Static Public Attributes

static std::string className
 Name of class.

Protected Member Functions

 CompressibleFlowSystem (const LibUtilities::SessionReaderSharedPtr &pSession)
virtual void v_InitObject ()
 Initialization object for CompressibleFlowSystem class.
virtual void v_GenerateSummary (SolverUtils::SummaryList &s)
 Print a summary of time stepping parameters.
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.
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.
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.
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.
void WallBC (int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
 Wall boundary conditions for compressible flow problems.
void WallViscousBC (int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
 Wall boundary conditions for viscous compressible flow problems.
void SymmetryBC (int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
 Symmetry boundary conditions for compressible flow problems.
void RiemannInvariantBC (int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
 Outflow characteristic boundary conditions for compressible flow problems.
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.
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.
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)
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.
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.
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator ()
 Get the maximum timestep estimator for cfl control.
virtual SOLVER_UTILS_EXPORT void v_DoSolve ()
 Solves an unsteady problem.
virtual SOLVER_UTILS_EXPORT void v_DoInitialise ()
 Sets up initial conditions.
virtual SOLVER_UTILS_EXPORT void v_AppendOutput1D (Array< OneD, Array< OneD, NekDouble > > &solution1D)
 Print the solution at each solution point in a txt file.
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_PostIntegrate (int step)
SOLVER_UTILS_EXPORT void CheckForRestartTime (NekDouble &time)
- Protected Member Functions inherited from Nektar::SolverUtils::EquationSystem
SOLVER_UTILS_EXPORT EquationSystem (const LibUtilities::SessionReaderSharedPtr &pSession)
 Initialises EquationSystem class members.
int nocase_cmp (const string &s1, const string &s2)
SOLVER_UTILS_EXPORT void SetBoundaryConditions (NekDouble time)
 Evaluates the boundary conditions at the given time.
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.
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.
virtual SOLVER_UTILS_EXPORT void v_TransCoeffToPhys ()
 Virtual function for transformation to physical space.
virtual SOLVER_UTILS_EXPORT void v_TransPhysToCoeff ()
 Virtual function for transformation to coefficient space.
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_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
StdRegions::StdQuadExpSharedPtr m_OrthoQuadExp
StdRegions::StdHexExpSharedPtr m_OrthoHexExp
bool m_smoothDiffusion
- Protected Attributes inherited from Nektar::SolverUtils::UnsteadySystem
int m_infosteps
 Number of time steps between outputting status information.
LibUtilities::TimeIntegrationWrapperSharedPtr m_intScheme
 Wrapper to the time integration scheme.
LibUtilities::TimeIntegrationSchemeOperators m_ode
 The time integration scheme operators to use.
LibUtilities::TimeIntegrationSolutionSharedPtr m_intSoln
NekDouble m_epsilon
bool m_explicitDiffusion
 Indicates if explicit or implicit treatment of diffusion is used.
bool m_explicitAdvection
 Indicates if explicit or implicit treatment of advection is used.
bool m_explicitReaction
 Indicates if explicit or implicit treatment of reaction is used.
bool m_homoInitialFwd
 Flag to determine if simulation should start in homogeneous forward transformed state.
std::vector< int > m_intVariables
std::vector< FilterSharedPtr > m_filters
- Protected Attributes inherited from Nektar::SolverUtils::EquationSystem
LibUtilities::CommSharedPtr m_comm
 Communicator.
LibUtilities::SessionReaderSharedPtr m_session
 The session reader.
LibUtilities::FieldIOSharedPtr m_fld
 Field input/output.
Array< OneD,
MultiRegions::ExpListSharedPtr > 		m_fields
 Array holding all dependent variables.
Array< OneD,
MultiRegions::ExpListSharedPtr > 		m_base
 Base fields.
Array< OneD,
MultiRegions::ExpListSharedPtr > 		m_derivedfields
 Array holding all dependent variables.
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
 Pointer to boundary conditions object.
SpatialDomains::MeshGraphSharedPtr m_graph
 Pointer to graph defining mesh.
std::string m_filename
 Filename.
std::string m_sessionName
 Name of the session.
NekDouble m_time
 Current time of simulation.
NekDouble m_fintime
 Finish time of the simulation.
NekDouble m_timestep
 Time step size.
NekDouble m_lambda
 Lambda constant in real system if one required.
NekDouble m_checktime
 Time between checkpoints.
int m_steps
 Number of steps to take.
int m_checksteps
 Number of steps between checkpoints.
int m_spacedim
 Spatial dimension (>= expansion dim).
int m_expdim
 Expansion dimension.
bool m_SingleMode
 Flag to determine if single homogeneous mode is used.
bool m_HalfMode
 Flag to determine if half homogeneous mode is used.
bool m_MultipleModes
 Flag to determine if use multiple homogenenous modes are used.
bool m_useFFT
 Flag to determine if FFT is used for homogeneous transform.
bool m_homogen_dealiasing
 Flag to determine if dealiasing is used for homogeneous simulations.
bool m_specHP_dealiasing
 Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
enum MultiRegions::ProjectionType m_projectionType
 Type of projection; e.g continuous or discontinuous.
Array< OneD, Array< OneD,
NekDouble > > 		m_traceNormals
 Array holding trace normals for DG simulations in the forwards direction.
Array< OneD, Array< OneD,
Array< OneD, NekDouble > > > 		m_gradtan
 1 x nvariable x nq
Array< OneD, Array< OneD,
Array< OneD, NekDouble > > > 		m_tanbasis
 2 x m_spacedim x nq
Array< OneD, bool > m_checkIfSystemSingular
 Flag to indicate if the fields should be checked for singularity.
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
 Map to identify relevant solver info to dump in output fields.
int m_NumQuadPointsError
 Number of Quadrature points used to work out the error.
enum HomogeneousType m_HomogeneousType
NekDouble m_LhomX
 physical length in X direction (if homogeneous)
NekDouble m_LhomY
 physical length in Y direction (if homogeneous)
NekDouble m_LhomZ
 physical length in Z direction (if homogeneous)
int m_npointsX
 number of points in X direction (if homogeneous)
int m_npointsY
 number of points in Y direction (if homogeneous)
int m_npointsZ
 number of points in Z direction (if homogeneous)
int m_HomoDirec
 number of homogenous directions
int m_NumMode
 Mode to use in case of single mode analysis.

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).
- Protected Types inherited from Nektar::SolverUtils::EquationSystem
enum  HomogeneousType { eHomogeneous1D, eHomogeneous2D, eHomogeneous3D, eNotHomogeneous }
 Parameter for homogeneous expansions. More...

Detailed Description

Definition at line 66 of file CompressibleFlowSystem.h.

Constructor & Destructor Documentation

Nektar::CompressibleFlowSystem::~CompressibleFlowSystem ( )
virtual

Destructor for CompressibleFlowSystem class.

Definition at line 241 of file CompressibleFlowSystem.cpp.

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

Definition at line 53 of file CompressibleFlowSystem.cpp.

: UnsteadySystem(pSession)
{
}

Member Function Documentation

static SolverUtils::EquationSystemSharedPtr Nektar::CompressibleFlowSystem::create ( const LibUtilities::SessionReaderSharedPtr &  pSession)
inlinestatic

Creates an instance of this class.

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

Definition at line 73 of file CompressibleFlowSystem.h.

{
return MemoryManager<CompressibleFlowSystem>::AllocateSharedPtr(pSession);
}
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 781 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 Nektar::NavierStokesCFE::SetBoundaryConditions(), Nektar::EulerADCFE::SetBoundaryConditions(), and Nektar::EulerCFE::SetBoundaryConditions().

{
int i, j;
int e, pnt;
int id1, id2, nBCEdgePts;
int nTracePts = GetTraceTotPoints();
int nVariables = physarray.num_elements();
int nDimensions = m_spacedim;
const Array<OneD, const int> &traceBndMap
= m_fields[0]->GetTraceBndMap();
// Get physical values of the forward trace
Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
for (i = 0; i < nVariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
int eMax;
eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
// Loop on bcRegions
for (e = 0; e < eMax; ++e)
{
nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
GetExp(e)->GetNumPoints(0);
id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
GetPhys_Offset(e) ;
id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
// Loop on points of bcRegion 'e'
for (i = 0; i < nBCEdgePts; i++)
{
pnt = id2+i;
// Setting up bcs for density
(m_fields[0]->GetBndCondExpansions()[bcRegion]->
UpdatePhys())[id1+i] = Fwd[0][pnt];
// Setting up bcs for velocities
for (j = 1; j <=nDimensions; ++j)
{
(m_fields[j]->GetBndCondExpansions()[bcRegion]->
UpdatePhys())[id1+i] = Fwd[j][pnt];
}
// Setting up bcs for energy
(m_fields[nVariables-1]->GetBndCondExpansions()[bcRegion]->
UpdatePhys())[id1+i] = Fwd[nVariables-1][pnt];
}
}
}
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 847 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().

{
int i, j;
int nq = physfield[0].num_elements();
int nVariables = m_fields.num_elements();
Array<OneD, NekDouble> pressure(nq);
Array<OneD, Array<OneD, NekDouble> > velocity(m_spacedim);
// Flux vector for the rho equation
for (i = 0; i < m_spacedim; ++i)
{
velocity[i] = Array<OneD, NekDouble>(nq);
Vmath::Vcopy(nq, physfield[i+1], 1, flux[0][i], 1);
}
GetVelocityVector(physfield, velocity);
GetPressure (physfield, velocity, pressure);
// Flux vector for the velocity fields
for (i = 0; i < m_spacedim; ++i)
{
for (j = 0; j < m_spacedim; ++j)
{
Vmath::Vmul(nq, velocity[j], 1, physfield[i+1], 1,
flux[i+1][j], 1);
}
// Add pressure to appropriate field
Vmath::Vadd(nq, flux[i+1][i], 1, pressure, 1, flux[i+1][i], 1);
}
// Flux vector for energy.
Vmath::Vadd(nq, physfield[m_spacedim+1], 1, pressure, 1,
pressure, 1);
for (j = 0; j < m_spacedim; ++j)
{
Vmath::Vmul(nq, velocity[j], 1, pressure, 1,
flux[m_spacedim+1][j], 1);
}
// For the smooth viscosity model
if (nVariables == m_spacedim+3)
{
// Add a zero row for the advective fluxes
for (j = 0; j < m_spacedim; ++j)
{
Vmath::Zero(nq, flux[m_spacedim+2][j], 1);
}
}
}
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 910 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().

{
int i, j;
int nq = physfield[0].num_elements();
int nVariables = m_fields.num_elements();
// Factor to rescale 1d points in dealiasing
NekDouble OneDptscale = 2;
nq = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
Array<OneD, NekDouble> pressure(nq);
Array<OneD, Array<OneD, NekDouble> > velocity(m_spacedim);
Array<OneD, Array<OneD, NekDouble> > physfield_interp(nVariables);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > flux_interp(
nVariables);
for (i = 0; i < nVariables; ++ i)
{
physfield_interp[i] = Array<OneD, NekDouble>(nq);
flux_interp[i] = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
m_fields[0]->PhysInterp1DScaled(
OneDptscale, physfield[i], physfield_interp[i]);
for (j = 0; j < m_spacedim; ++j)
{
flux_interp[i][j] = Array<OneD, NekDouble>(nq);
}
}
// Flux vector for the rho equation
for (i = 0; i < m_spacedim; ++i)
{
velocity[i] = Array<OneD, NekDouble>(nq);
// Galerkin project solution back to original space
m_fields[0]->PhysGalerkinProjection1DScaled(
OneDptscale, physfield_interp[i+1], flux[0][i]);
}
GetVelocityVector(physfield_interp, velocity);
GetPressure (physfield_interp, velocity, pressure);
// Evaluation of flux vector for the velocity fields
for (i = 0; i < m_spacedim; ++i)
{
for (j = 0; j < m_spacedim; ++j)
{
Vmath::Vmul(nq, velocity[j], 1, physfield_interp[i+1], 1,
flux_interp[i+1][j], 1);
}
// Add pressure to appropriate field
Vmath::Vadd(nq, flux_interp[i+1][i], 1, pressure,1,
flux_interp[i+1][i], 1);
}
// Galerkin project solution back to origianl space
for (i = 0; i < m_spacedim; ++i)
{
for (j = 0; j < m_spacedim; ++j)
{
m_fields[0]->PhysGalerkinProjection1DScaled(
OneDptscale, flux_interp[i+1][j], flux[i+1][j]);
}
}
// Evaluation of flux vector for energy
Vmath::Vadd(nq, physfield_interp[m_spacedim+1], 1, pressure, 1,
pressure, 1);
for (j = 0; j < m_spacedim; ++j)
{
Vmath::Vmul(nq, velocity[j], 1, pressure, 1,
flux_interp[m_spacedim+1][j], 1);
// Galerkin project solution back to origianl space
m_fields[0]->PhysGalerkinProjection1DScaled(
OneDptscale,
flux_interp[m_spacedim+1][j],
flux[m_spacedim+1][j]);
}
}
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 243 of file CompressibleFlowSystem.h.

References m_gamma.

Referenced by v_InitObject().

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

Definition at line 238 of file CompressibleFlowSystem.h.

References m_gasConstant.

{
return m_gasConstant;
}
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 253 of file CompressibleFlowSystem.h.

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

Referenced by v_InitObject().

{
return m_traceNormals;
}
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 1759 of file CompressibleFlowSystem.cpp.

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

{
int nBCEdgePts = physfield[0].num_elements();
NekDouble alpha = -0.5;
// Calculate ||rho v||^2
Vmath::Vmul(nBCEdgePts, physfield[1], 1, physfield[1], 1, pressure, 1);
for (int i = 1; i < m_spacedim; ++i)
{
Vmath::Vvtvp(nBCEdgePts, physfield[1+i], 1, physfield[1+i], 1,
pressure, 1, pressure, 1);
}
// Divide by rho to get rho*||v||^2
Vmath::Vdiv (nBCEdgePts, pressure, 1, physfield[0], 1, pressure, 1);
// pressure <- E - 0.5*pressure
Vmath::Svtvp(nBCEdgePts, alpha,
pressure, 1, physfield[m_spacedim+1], 1, pressure, 1);
// Multiply by (gamma-1)
Vmath::Smul (nBCEdgePts, m_gamma-1, pressure, 1, pressure, 1);
}
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 248 of file CompressibleFlowSystem.h.

References m_vecLocs.

Referenced by v_InitObject().

{
return m_vecLocs;
}
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 1001 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().

{
int j, k;
int nVariables = m_fields.num_elements();
int nPts = physfield[0].num_elements();
// Stokes hypotesis
const NekDouble lambda = -2.0/3.0;
// Auxiliary variables
Array<OneD, NekDouble > mu (nPts, 0.0);
Array<OneD, NekDouble > thermalConductivity(nPts, 0.0);
Array<OneD, NekDouble > mu2 (nPts, 0.0);
Array<OneD, NekDouble > divVel (nPts, 0.0);
// Variable viscosity through the Sutherland's law
if (m_ViscosityType == "Variable")
{
GetDynamicViscosity(physfield[nVariables-2], mu);
NekDouble tRa = m_Cp / m_Prandtl;
Vmath::Smul(nPts, tRa, &mu[0], 1, &thermalConductivity[0], 1);
}
else
{
Vmath::Fill(nPts, m_mu, &mu[0], 1);
Vmath::Fill(nPts, m_thermalConductivity,
&thermalConductivity[0], 1);
}
// Computing diagonal terms of viscous stress tensor
Array<OneD, Array<OneD, NekDouble> > tmp(m_spacedim);
Array<OneD, Array<OneD, NekDouble> > Sgg(m_spacedim);
// mu2 = 2 * mu
Vmath::Smul(nPts, 2.0, &mu[0], 1, &mu2[0], 1);
// Velocity divergence
for (j = 0; j < m_spacedim; ++j)
{
Vmath::Vadd(nPts, &divVel[0], 1, &derivativesO1[j][j][0], 1,
&divVel[0], 1);
}
// Velocity divergence scaled by lambda * mu
Vmath::Smul(nPts, lambda, &divVel[0], 1, &divVel[0], 1);
Vmath::Vmul(nPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
// Diagonal terms of viscous stress tensor (Sxx, Syy, Szz)
// Sjj = 2 * mu * du_j/dx_j - (2 / 3) * mu * sum_j(du_j/dx_j)
for (j = 0; j < m_spacedim; ++j)
{
tmp[j] = Array<OneD, NekDouble>(nPts, 0.0);
Sgg[j] = Array<OneD, NekDouble>(nPts, 0.0);
Vmath::Vmul(nPts, &mu2[0], 1, &derivativesO1[j][j][0], 1,
&tmp[j][0], 1);
Vmath::Vadd(nPts, &tmp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
}
// Extra diagonal terms of viscous stress tensor (Sxy, Sxz, Syz)
// Note: they exist for 2D and 3D problems only
Array<OneD, NekDouble > Sxy(nPts, 0.0);
Array<OneD, NekDouble > Sxz(nPts, 0.0);
Array<OneD, NekDouble > Syz(nPts, 0.0);
if (m_spacedim == 2)
{
// Sxy = (du/dy + dv/dx)
Vmath::Vadd(nPts, &derivativesO1[0][1][0], 1,
&derivativesO1[1][0][0], 1, &Sxy[0], 1);
// Sxy = mu * (du/dy + dv/dx)
Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
}
else if (m_spacedim == 3)
{
// Sxy = (du/dy + dv/dx)
Vmath::Vadd(nPts, &derivativesO1[0][1][0], 1,
&derivativesO1[1][0][0], 1, &Sxy[0], 1);
// Sxz = (du/dz + dw/dx)
Vmath::Vadd(nPts, &derivativesO1[0][2][0], 1,
&derivativesO1[2][0][0], 1, &Sxz[0], 1);
// Syz = (dv/dz + dw/dy)
Vmath::Vadd(nPts, &derivativesO1[1][2][0], 1,
&derivativesO1[2][1][0], 1, &Syz[0], 1);
// Sxy = mu * (du/dy + dv/dx)
Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
// Sxz = mu * (du/dy + dv/dx)
Vmath::Vmul(nPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
// Syz = mu * (du/dy + dv/dx)
Vmath::Vmul(nPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
}
// Energy-related terms
Array<OneD, NekDouble > STx(nPts, 0.0);
Array<OneD, NekDouble > STy(nPts, 0.0);
Array<OneD, NekDouble > STz(nPts, 0.0);
// Building the viscous flux vector
// Viscous flux vector for the rho equation
for (k = 0; k < m_spacedim; ++k)
{
Vmath::Zero(nPts, viscousTensor[k][0], 1);
}
if (m_spacedim == 1)
{
Array<OneD, NekDouble > tmp1(nPts, 0.0);
// u * Sxx
Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
// k * dT/dx
Vmath::Vmul(nPts, &thermalConductivity[0], 1,
&derivativesO1[0][1][0], 1, &tmp1[0], 1);
// STx = u * Sxx + (K / mu) * dT/dx
Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
}
else if (m_spacedim == 2)
{
Array<OneD, NekDouble > tmp1(nPts, 0.0);
Array<OneD, NekDouble > tmp2(nPts, 0.0);
// Computation of STx
// u * Sxx
Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
// v * Sxy
Vmath::Vmul(nPts, &physfield[1][0], 1, &Sxy[0], 1, &tmp1[0], 1);
// k * dT/dx
Vmath::Vmul(nPts, &thermalConductivity[0], 1,
&derivativesO1[0][2][0], 1, &tmp2[0], 1);
// STx = u * Sxx + v * Sxy + K * dT/dx
Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
// Computation of STy
// Re-initialise temporary arrays
Vmath::Zero(nPts, &tmp1[0], 1);
Vmath::Zero(nPts, &tmp2[0], 1);
// v * Syy
Vmath::Vmul(nPts, &physfield[1][0], 1, &Sgg[1][0], 1, &STy[0], 1);
// u * Sxy
Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxy[0], 1, &tmp1[0], 1);
// k * dT/dy
Vmath::Vmul(nPts, &thermalConductivity[0], 1,
&derivativesO1[1][2][0], 1, &tmp2[0], 1);
// STy = v * Syy + u * Sxy + K * dT/dy
Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
}
else if (m_spacedim == 3)
{
Array<OneD, NekDouble > tmp1(nPts, 0.0);
Array<OneD, NekDouble > tmp2(nPts, 0.0);
Array<OneD, NekDouble > tmp3(nPts, 0.0);
// Computation of STx
// u * Sxx
Vmath::Vmul(nPts, &physfield[0][0], 1, &Sgg[0][0], 1, &STx[0], 1);
// v * Sxy
Vmath::Vmul(nPts, &physfield[1][0], 1, &Sxy[0], 1, &tmp1[0], 1);
// v * Sxz
Vmath::Vmul(nPts, &physfield[2][0], 1, &Sxz[0], 1, &tmp2[0], 1);
// k * dT/dx
Vmath::Vmul(nPts, &thermalConductivity[0], 1,
&derivativesO1[0][3][0], 1, &tmp3[0], 1);
// STx = u * Sxx + v * Sxy + w * Sxz + (K / mu) * dT/dx
Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
Vmath::Vadd(nPts, &STx[0], 1, &tmp3[0], 1, &STx[0], 1);
// Computation of STy
// Re-initialise temporary arrays
Vmath::Zero(nPts, &tmp1[0], 1);
Vmath::Zero(nPts, &tmp2[0], 1);
Vmath::Zero(nPts, &tmp3[0], 1);
// v * Syy
Vmath::Vmul(nPts, &physfield[1][0], 1, &Sgg[1][0], 1, &STy[0], 1);
// u * Sxy
Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxy[0], 1, &tmp1[0], 1);
// w * Syz
Vmath::Vmul(nPts, &physfield[2][0], 1, &Syz[0], 1, &tmp2[0], 1);
// k * dT/dy
Vmath::Vmul(nPts, &thermalConductivity[0], 1,
&derivativesO1[1][3][0], 1, &tmp3[0], 1);
// STy = v * Syy + u * Sxy + w * Syz + K * dT/dy
Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
Vmath::Vadd(nPts, &STy[0], 1, &tmp3[0], 1, &STy[0], 1);
// Computation of STz
// Re-initialise temporary arrays
Vmath::Zero(nPts, &tmp1[0], 1);
Vmath::Zero(nPts, &tmp2[0], 1);
Vmath::Zero(nPts, &tmp3[0], 1);
// w * Szz
Vmath::Vmul(nPts, &physfield[2][0], 1, &Sgg[2][0], 1, &STz[0], 1);
// u * Sxz
Vmath::Vmul(nPts, &physfield[0][0], 1, &Sxz[0], 1, &tmp1[0], 1);
// v * Syz
Vmath::Vmul(nPts, &physfield[1][0], 1, &Syz[0], 1, &tmp2[0], 1);
// k * dT/dz
Vmath::Vmul(nPts, &thermalConductivity[0], 1,
&derivativesO1[2][3][0], 1, &tmp3[0], 1);
// STz = w * Szz + u * Sxz + v * Syz + K * dT/dz
Vmath::Vadd(nPts, &STz[0], 1, &tmp1[0], 1, &STz[0], 1);
Vmath::Vadd(nPts, &STz[0], 1, &tmp2[0], 1, &STz[0], 1);
Vmath::Vadd(nPts, &STz[0], 1, &tmp3[0], 1, &STz[0], 1);
}
switch (m_spacedim)
{
case 1:
{
// f_11v = f_rho = 0
Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
// f_21v = f_rhou
Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
// f_31v = f_E
Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][2][0], 1);
break;
}
case 2:
{
// f_11v = f_rho1 = 0
Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
// f_12v = f_rho2 = 0
Vmath::Zero(nPts, &viscousTensor[1][0][0], 1);
// f_21v = f_rhou1
Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
// f_22v = f_rhou2
Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[1][1][0], 1);
// f_31v = f_rhov1
Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[0][2][0], 1);
// f_32v = f_rhov2
Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor[1][2][0], 1);
// f_41v = f_E1
Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][3][0], 1);
// f_42v = f_E2
Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor[1][3][0], 1);
break;
}
case 3:
{
// f_11v = f_rho1 = 0
Vmath::Zero(nPts, &viscousTensor[0][0][0], 1);
// f_12v = f_rho2 = 0
Vmath::Zero(nPts, &viscousTensor[1][0][0], 1);
// f_13v = f_rho3 = 0
Vmath::Zero(nPts, &viscousTensor[2][0][0], 1);
// f_21v = f_rhou1
Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor[0][1][0], 1);
// f_22v = f_rhou2
Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[1][1][0], 1);
// f_23v = f_rhou3
Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor[2][1][0], 1);
// f_31v = f_rhov1
Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor[0][2][0], 1);
// f_32v = f_rhov2
Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor[1][2][0], 1);
// f_33v = f_rhov3
Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor[2][2][0], 1);
// f_31v = f_rhow1
Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor[0][3][0], 1);
// f_32v = f_rhow2
Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor[1][3][0], 1);
// f_33v = f_rhow3
Vmath::Vcopy(nPts, &Sgg[2][0], 1, &viscousTensor[2][3][0], 1);
// f_41v = f_E1
Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor[0][4][0], 1);
// f_42v = f_E2
Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor[1][4][0], 1);
// f_43v = f_E3
Vmath::Vcopy(nPts, &STz[0], 1, &viscousTensor[2][4][0], 1);
break;
}
default:
{
ASSERTL0(false, "Illegal expansion dimension");
}
}
}
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 1334 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().

{
#if 0
int i, j, k;
int nVariables = m_fields.num_elements();
int nPts = physfield[0].num_elements();
int variables_phys = physfield.num_elements();
// Factor to rescale 1d points in dealiasing.
NekDouble OneDptscale = 2;
// Get number of points to dealias a cubic non-linearity
nPts = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
int nVariables_aux = derivativesO1[0].num_elements();
Array<OneD, Array<OneD, NekDouble> > physfield_interp(variables_phys);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > derivativesO1_interp(
m_spacedim);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >viscousTensor_interp(
m_spacedim);
for (i = 0; i < m_spacedim; ++ i)
{
viscousTensor_interp[i] = Array<OneD, Array<OneD, NekDouble> >(
nVariables);
for (j = 0; j < nVariables; ++j)
{
viscousTensor_interp[i][j] = Array<OneD, NekDouble>(nPts);
}
}
// Stokes hypotesis
NekDouble lambda = -2.0/3.0;
// Auxiliary variables
Array<OneD, NekDouble > mu (nPts, 0.0);
Array<OneD, NekDouble > mu2 (nPts, 0.0);
Array<OneD, NekDouble > divVel (nPts, 0.0);
Array<OneD, NekDouble > pressure (nPts, 0.0);
Array<OneD, NekDouble > temperature(nPts, 0.0);
for (i = 0; i < nVariables; ++i)
{
m_fields[0]->PhysInterp1DScaled(
OneDptscale, physfield[i], fields_interp[i]);
}
for (i = 0; i < variables_phys; ++i)
{
physfield_interp[i] = Array<OneD, NekDouble> (nPts);
// Interpolation to higher space
m_fields[0]->PhysInterp1DScaled(OneDptscale, physfield[i],
physfield_interp[i]);
}
for (i = 0; i < m_spacedim; ++i)
{
derivativesO1_interp[i] = Array<OneD, Array<OneD, NekDouble> >(
nVariables_aux);
for (j = 0; j < nVariables_aux; ++j)
{
derivativesO1_interp[i][j] = Array<OneD, NekDouble>(nPts);
m_fields[0]->PhysInterp1DScaled(
OneDptscale, derivativesO1[i][j],
derivativesO1_interp[i][j]);
}
}
// Thermodynamic related quantities
GetPressure(fields_interp, pressure);
GetTemperature(fields_interp, pressure, temperature);
// Variable viscosity through the Sutherland's law
if (m_ViscosityType == "Variable")
{
GetDynamicViscosity(fields_interp[variables_phys-1], mu);
}
else
{
Vmath::Sadd(nPts, m_mu, &mu[0], 1, &mu[0], 1);
}
// Computing diagonal terms of viscous stress tensor
Array<OneD, Array<OneD, NekDouble> > tmp(m_spacedim);
Array<OneD, Array<OneD, NekDouble> > Sgg(m_spacedim);
// mu2 = 2 * mu
Vmath::Smul(nPts, 2.0, &mu[0], 1, &mu2[0], 1);
// Velocity divergence
for (j = 0; j < m_spacedim; ++j)
{
Vmath::Vadd(nPts, &divVel[0], 1, &derivativesO1_interp[j][j][0], 1,
&divVel[0], 1);
}
// Velocity divergence scaled by lambda * mu
Vmath::Smul(nPts, lambda, &divVel[0], 1, &divVel[0], 1);
Vmath::Vmul(nPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
// Digonal terms of viscous stress tensor (Sxx, Syy, Szz)
// Sjj = 2 * mu * du_j/dx_j - (2 / 3) * mu * sum_j(du_j/dx_j)
for (j = 0; j < m_spacedim; ++j)
{
tmp[j] = Array<OneD, NekDouble>(nPts, 0.0);
Sgg[j] = Array<OneD, NekDouble>(nPts, 0.0);
Vmath::Vmul(nPts, &mu2[0], 1, &derivativesO1_interp[j][j][0], 1,
&tmp[j][0], 1);
Vmath::Vadd(nPts, &tmp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
}
// Extra diagonal terms of viscous stress tensor (Sxy, Sxz, Syz)
// Note: they exist for 2D and 3D problems only
Array<OneD, NekDouble > Sxy(nPts, 0.0);
Array<OneD, NekDouble > Sxz(nPts, 0.0);
Array<OneD, NekDouble > Syz(nPts, 0.0);
if (m_spacedim == 2)
{
// Sxy = (du/dy + dv/dx)
Vmath::Vadd(nPts, &derivativesO1_interp[0][1][0], 1,
&derivativesO1_interp[1][0][0], 1, &Sxy[0], 1);
// Sxy = mu * (du/dy + dv/dx)
Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
}
else if (m_spacedim == 3)
{
// Sxy = (du/dy + dv/dx)
Vmath::Vadd(nPts, &derivativesO1_interp[0][1][0], 1,
&derivativesO1_interp[1][0][0], 1, &Sxy[0], 1);
// Sxz = (du/dz + dw/dx)
Vmath::Vadd(nPts, &derivativesO1_interp[0][2][0], 1,
&derivativesO1_interp[2][0][0], 1, &Sxz[0], 1);
// Syz = (dv/dz + dw/dy)
Vmath::Vadd(nPts, &derivativesO1_interp[1][2][0], 1,
&derivativesO1_interp[2][1][0], 1, &Syz[0], 1);
// Sxy = mu * (du/dy + dv/dx)
Vmath::Vmul(nPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
// Sxz = mu * (du/dy + dv/dx)
Vmath::Vmul(nPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
// Syz = mu * (du/dy + dv/dx)
Vmath::Vmul(nPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
}
// Energy-related terms
Array<OneD, NekDouble > STx(nPts, 0.0);
Array<OneD, NekDouble > STy(nPts, 0.0);
Array<OneD, NekDouble > STz(nPts, 0.0);
// Building the viscous flux vector
if (i == 0)
{
// Viscous flux vector for the rho equation
for (k = 0; k < m_spacedim; ++k)
{
Vmath::Zero(nPts, viscousTensor_interp[k][i], 1);
}
}
if (m_spacedim == 1)
{
Array<OneD, NekDouble > tmp1(nPts, 0.0);
// u * Sxx
Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
&Sgg[0][0], 1, &STx[0], 1);
// k * dT/dx
Vmath::Smul(nPts, m_thermalConductivity,
&derivativesO1_interp[0][1][0], 1, &tmp1[0], 1);
// STx = u * Sxx + (K / mu) * dT/dx
Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
}
else if (m_spacedim == 2)
{
Array<OneD, NekDouble > tmp1(nPts, 0.0);
Array<OneD, NekDouble > tmp2(nPts, 0.0);
// Computation of STx
// u * Sxx
Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
&Sgg[0][0], 1, &STx[0], 1);
// v * Sxy
Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
&Sxy[0], 1, &tmp1[0], 1);
// k * dT/dx
Vmath::Smul(nPts, m_thermalConductivity,
&derivativesO1_interp[0][2][0], 1, &tmp2[0], 1);
// STx = u * Sxx + v * Sxy + K * dT/dx
Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
// Computation of STy
// Re-initialise temporary arrays
Vmath::Zero(nPts, &tmp1[0], 1);
Vmath::Zero(nPts, &tmp2[0], 1);
// v * Syy
Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
&Sgg[1][0], 1, &STy[0], 1);
// u * Sxy
Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
&Sxy[0], 1, &tmp1[0], 1);
// k * dT/dy
Vmath::Smul(nPts, m_thermalConductivity,
&derivativesO1_interp[1][2][0], 1, &tmp2[0], 1);
// STy = v * Syy + u * Sxy + K * dT/dy
Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
}
else if (m_spacedim == 3)
{
Array<OneD, NekDouble > tmp1(nPts, 0.0);
Array<OneD, NekDouble > tmp2(nPts, 0.0);
Array<OneD, NekDouble > tmp3(nPts, 0.0);
// Computation of STx
// u * Sxx
Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
&Sgg[0][0], 1, &STx[0], 1);
// v * Sxy
Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
&Sxy[0], 1, &tmp1[0], 1);
// v * Sxy
Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
&Sxz[0], 1, &tmp2[0], 1);
// k * dT/dx
Vmath::Smul(nPts, m_thermalConductivity,
&derivativesO1_interp[0][3][0], 1, &tmp3[0], 1);
// STx = u * Sxx + v * Sxy + w * Sxz + (K / mu) * dT/dx
Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
Vmath::Vadd(nPts, &STx[0], 1, &tmp3[0], 1, &STx[0], 1);
// Computation of STy
// Re-initialise temporary arrays
Vmath::Zero(nPts, &tmp1[0], 1);
Vmath::Zero(nPts, &tmp2[0], 1);
Vmath::Zero(nPts, &tmp3[0], 1);
// v * Syy
Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
&Sgg[1][0], 1, &STy[0], 1);
// u * Sxy
Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
&Sxy[0], 1, &tmp1[0], 1);
// w * Syz
Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
&Syz[0], 1, &tmp2[0], 1);
// k * dT/dy
Vmath::Smul(nPts, m_thermalConductivity,
&derivativesO1_interp[1][3][0], 1, &tmp3[0], 1);
// STy = v * Syy + u * Sxy + w * Syz + K * dT/dy
Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
Vmath::Vadd(nPts, &STy[0], 1, &tmp3[0], 1, &STy[0], 1);
// Computation of STz
// Re-initialise temporary arrays
Vmath::Zero(nPts, &tmp1[0], 1);
Vmath::Zero(nPts, &tmp2[0], 1);
Vmath::Zero(nPts, &tmp3[0], 1);
// w * Szz
Vmath::Vmul(nPts, &physfield_interp[2][0], 1,
&Sgg[2][0], 1, &STz[0], 1);
// u * Sxz
Vmath::Vmul(nPts, &physfield_interp[0][0], 1,
&Sxz[0], 1, &tmp1[0], 1);
// v * Syz
Vmath::Vmul(nPts, &physfield_interp[1][0], 1,
&Syz[0], 1, &tmp2[0], 1);
// k * dT/dz
Vmath::Smul(nPts, m_thermalConductivity,
&derivativesO1_interp[2][3][0], 1, &tmp3[0], 1);
// STz = w * Szz + u * Sxz + v * Syz + K * dT/dz
Vmath::Vadd(nPts, &STz[0], 1, &tmp1[0], 1, &STz[0], 1);
Vmath::Vadd(nPts, &STz[0], 1, &tmp2[0], 1, &STz[0], 1);
Vmath::Vadd(nPts, &STz[0], 1, &tmp3[0], 1, &STz[0], 1);
}
switch (m_spacedim)
{
case 1:
{
int nq = physfield[0].num_elements();
// f_11v = f_rho = 0
Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
// f_21v = f_rhou
Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
// f_31v = f_E
Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][2][0], 1);
break;
}
case 2:
{
int nq = physfield[0].num_elements();
// f_11v = f_rho1 = 0
Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
// f_12v = f_rho2 = 0
Vmath::Zero(nq, &viscousTensor_interp[1][0][0], 1);
// f_21v = f_rhou1
Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
// f_22v = f_rhou2
Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[1][1][0], 1);
// f_31v = f_rhov1
Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[0][2][0], 1);
// f_32v = f_rhov2
Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor_interp[1][2][0], 1);
// f_41v = f_E1
Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][3][0], 1);
// f_42v = f_E2
Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor_interp[1][3][0], 1);
break;
}
case 3:
{
int nq = physfield[0].num_elements();
// f_11v = f_rho1 = 0
Vmath::Zero(nq, &viscousTensor_interp[0][0][0], 1);
// f_12v = f_rho2 = 0
Vmath::Zero(nq, &viscousTensor_interp[1][0][0], 1);
// f_13v = f_rho3 = 0
Vmath::Zero(nq, &viscousTensor_interp[2][0][0], 1);
// f_21v = f_rhou1
Vmath::Vcopy(nPts, &Sgg[0][0], 1, &viscousTensor_interp[0][1][0], 1);
// f_22v = f_rhou2
Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[1][1][0], 1);
// f_23v = f_rhou3
Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor_interp[2][1][0], 1);
// f_31v = f_rhov1
Vmath::Vcopy(nPts, &Sxy[0], 1, &viscousTensor_interp[0][2][0], 1);
// f_32v = f_rhov2
Vmath::Vcopy(nPts, &Sgg[1][0], 1, &viscousTensor_interp[1][2][0], 1);
// f_33v = f_rhov3
Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor_interp[2][2][0], 1);
// f_31v = f_rhow1
Vmath::Vcopy(nPts, &Sxz[0], 1, &viscousTensor_interp[0][3][0], 1);
// f_32v = f_rhow2
Vmath::Vcopy(nPts, &Syz[0], 1, &viscousTensor_interp[1][3][0], 1);
// f_33v = f_rhow3
Vmath::Vcopy(nPts, &Sgg[2][0], 1, &viscousTensor_interp[2][3][0], 1);
// f_41v = f_E1
Vmath::Vcopy(nPts, &STx[0], 1, &viscousTensor_interp[0][4][0], 1);
// f_42v = f_E2
Vmath::Vcopy(nPts, &STy[0], 1, &viscousTensor_interp[1][4][0], 1);
// f_43v = f_E3
Vmath::Vcopy(nPts, &STz[0], 1, &viscousTensor_interp[2][4][0], 1);
break;
}
default:
{
ASSERTL0(false, "Illegal expansion dimension");
}
}
for (i = 0; i < m_spacedim; ++i)
{
for (j = 1; j < nVariables; ++j)
{
// Galerkin project solution back to origianl space
m_fields[0]->PhysGalerkinProjection1DScaled(
OneDptscale,
viscousTensor_interp[i][j],
viscousTensor[i][j]);
}
}
#endif
}
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 525 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 Nektar::NavierStokesCFE::SetBoundaryConditions(), Nektar::EulerADCFE::SetBoundaryConditions(), and Nektar::EulerCFE::SetBoundaryConditions().

{
int i, j;
int nTracePts = GetTraceTotPoints();
int nVariables = physarray.num_elements();
int nDimensions = m_spacedim;
const Array<OneD, const int> &traceBndMap
= m_fields[0]->GetTraceBndMap();
NekDouble gamma = m_gamma;
NekDouble gammaInv = 1.0 / gamma;
NekDouble gammaMinusOne = gamma - 1.0;
NekDouble gammaMinusOneInv = 1.0 / gammaMinusOne;
Array<OneD, NekDouble> tmp1 (nTracePts, 0.0);
Array<OneD, NekDouble> tmp2 (nTracePts, 0.0);
Array<OneD, NekDouble> VnInf(nTracePts, 0.0);
Array<OneD, NekDouble> velInf(nDimensions, 0.0);
// Computing the normal velocity for characteristics coming
// from outside the computational domain
velInf[0] = m_uInf;
Vmath::Smul(nTracePts, m_uInf, m_traceNormals[0], 1, VnInf, 1);
if (nDimensions == 2 || nDimensions == 3)
{
velInf[1] = m_vInf;
Vmath::Smul(nTracePts, m_vInf, m_traceNormals[1], 1, tmp1, 1);
Vmath::Vadd(nTracePts, VnInf, 1, tmp1, 1, VnInf, 1);
}
if (nDimensions == 3)
{
velInf[2] = m_wInf;
Vmath::Smul(nTracePts, m_wInf, m_traceNormals[2], 1, tmp2, 1);
Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
}
// Get physical values of the forward trace
Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
for (i = 0; i < nVariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
// Computing the normal velocity for characteristics coming
// from inside the computational domain
Array<OneD, NekDouble > Vn (nTracePts, 0.0);
Array<OneD, NekDouble > Vel(nTracePts, 0.0);
for (i = 0; i < nDimensions; ++i)
{
Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, Vel, 1);
Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Vel, 1, Vn, 1, Vn, 1);
}
// Computing the absolute value of the velocity in order to compute the
// Mach number to decide whether supersonic or subsonic
Array<OneD, NekDouble > absVel(nTracePts, 0.0);
for (i = 0; i < nDimensions; ++i)
{
Vmath::Vdiv(nTracePts, Fwd[i+1], 1, Fwd[0], 1, tmp1, 1);
Vmath::Vmul(nTracePts, tmp1, 1, tmp1, 1, tmp1, 1);
Vmath::Vadd(nTracePts, tmp1, 1, absVel, 1, absVel, 1);
}
Vmath::Vsqrt(nTracePts, absVel, 1, absVel, 1);
// Get speed of sound
Array<OneD, NekDouble > soundSpeed(nTracePts);
Array<OneD, NekDouble > pressure (nTracePts);
for (i = 0; i < nTracePts; i++)
{
if (m_spacedim == 1)
{
pressure[i] = (gammaMinusOne) * (Fwd[2][i] -
0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i]));
}
else if (m_spacedim == 2)
{
pressure[i] = (gammaMinusOne) * (Fwd[3][i] -
0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
Fwd[2][i] * Fwd[2][i] / Fwd[0][i]));
}
else
{
pressure[i] = (gammaMinusOne) * (Fwd[4][i] -
0.5 * (Fwd[1][i] * Fwd[1][i] / Fwd[0][i] +
Fwd[2][i] * Fwd[2][i] / Fwd[0][i] +
Fwd[3][i] * Fwd[3][i] / Fwd[0][i]));
}
soundSpeed[i] = sqrt(gamma * pressure[i] / Fwd[0][i]);
}
// Get Mach
Array<OneD, NekDouble > Mach(nTracePts, 0.0);
Vmath::Vdiv(nTracePts, Vn, 1, soundSpeed, 1, Mach, 1);
Vmath::Vabs(nTracePts, Mach, 1, Mach, 1);
// Auxiliary variables
int eMax;
int e, id1, id2, nBCEdgePts, pnt;
NekDouble cPlus, rPlus, cMinus, rMinus, VDBC, VNBC;
Array<OneD, NekDouble> velBC(nDimensions, 0.0);
Array<OneD, NekDouble> rhoVelBC(nDimensions, 0.0);
NekDouble rhoBC, EBC, cBC, sBC, pBC;
eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
// Loop on bcRegions
for (e = 0; e < eMax; ++e)
{
nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
GetExp(e)->GetNumPoints(0);
id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
GetPhys_Offset(e);
id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
// Loop on the points of the bcRegion
for (i = 0; i < nBCEdgePts; i++)
{
pnt = id2+i;
// Impose inflow Riemann invariant
if (Vn[pnt] <= 0.0)
{
// Subsonic flows
if (Mach[pnt] < 1.00)
{
// + Characteristic from inside
cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
// - Characteristic from boundary
cMinus = sqrt(gamma * m_pInf / m_rhoInf);
rMinus = VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
}
else
{
// + Characteristic from inside
cPlus = sqrt(gamma * m_pInf / m_rhoInf);
rPlus = VnInf[pnt] + 2.0 * cPlus * gammaMinusOneInv;
// + Characteristic from inside
cMinus = sqrt(gamma * m_pInf / m_rhoInf);
rMinus = VnInf[pnt] - 2.0 * cPlus * gammaMinusOneInv;
}
// Riemann boundary variables
VNBC = 0.5 * (rPlus + rMinus);
cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
VDBC = VNBC - VnInf[pnt];
// Thermodynamic boundary variables
sBC = m_pInf / (pow(m_rhoInf, gamma));
rhoBC = pow((cBC * cBC) / (gamma * sBC), gammaMinusOneInv);
pBC = rhoBC * cBC * cBC * gammaInv;
// Kinetic energy initialiasation
NekDouble EkBC = 0.0;
// Boundary velocities
for ( j = 0; j < nDimensions; ++j)
{
velBC[j] = velInf[j] + VDBC * m_traceNormals[j][pnt];
rhoVelBC[j] = rhoBC * velBC[j];
EkBC += 0.5 * rhoBC * velBC[j]*velBC[j];
}
// Boundary energy
EBC = pBC * gammaMinusOneInv + EkBC;
// Imposing Riemann Invariant boundary conditions
(m_fields[0]->GetBndCondExpansions()[bcRegion]->
UpdatePhys())[id1+i] = rhoBC;
for (j = 0; j < nDimensions; ++j)
{
(m_fields[j+1]->GetBndCondExpansions()[bcRegion]->
UpdatePhys())[id1+i] = rhoVelBC[j];
}
(m_fields[nDimensions+1]->GetBndCondExpansions()[bcRegion]->
UpdatePhys())[id1+i] = EBC;
}
else // Impose outflow Riemann invariant
{
// Subsonic flows
if (Mach[pnt] < 1.00)
{
// + Characteristic from inside
cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
// - Characteristic from boundary
cMinus = sqrt(gamma * m_pInf / m_rhoInf);
rMinus = VnInf[pnt] - 2.0 * cMinus * gammaMinusOneInv;
}
else
{
// + Characteristic from inside
cPlus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
rPlus = Vn[pnt] + 2.0 * cPlus * gammaMinusOneInv;
// + Characteristic from inside
cMinus = sqrt(gamma * pressure[pnt] / Fwd[0][pnt]);
rMinus = Vn[pnt] - 2.0 * cPlus * gammaMinusOneInv;
}
// Riemann boundary variables
VNBC = 0.5 * (rPlus + rMinus);
cBC = 0.25 * gammaMinusOne * (rPlus - rMinus);
VDBC = VNBC - Vn[pnt];
// Thermodynamic boundary variables
sBC = pressure[pnt] / (pow(Fwd[0][pnt], gamma));
rhoBC = pow((cBC * cBC) / (gamma * sBC), gammaMinusOneInv);
pBC = rhoBC * cBC * cBC * gammaInv;
// Kinetic energy initialiasation
NekDouble EkBC = 0.0;
// Boundary velocities
for ( j = 0; j < nDimensions; ++j)
{
velBC[j] = Fwd[j+1][pnt] / Fwd[0][pnt] +
VDBC * m_traceNormals[j][pnt];
rhoVelBC[j] = rhoBC * velBC[j];
EkBC += 0.5 * rhoBC * velBC[j]*velBC[j];
}
// Boundary energy
EBC = pBC * gammaMinusOneInv + EkBC;
// Imposing Riemann Invariant boundary conditions
(m_fields[0]->GetBndCondExpansions()[bcRegion]->
UpdatePhys())[id1+i] = rhoBC;
for (j = 0; j < nDimensions; ++j)
{
(m_fields[j+1]->GetBndCondExpansions()[bcRegion]->
UpdatePhys())[id1+i] = rhoVelBC[j];
}
(m_fields[nDimensions+1]->GetBndCondExpansions()[bcRegion]->
UpdatePhys())[id1+i] = EBC;
}
}
}
}
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 429 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 Nektar::NavierStokesCFE::SetBoundaryConditions(), Nektar::EulerADCFE::SetBoundaryConditions(), and Nektar::EulerCFE::SetBoundaryConditions().

{
int i;
int nTracePts = GetTraceTotPoints();
int nVariables = physarray.num_elements();
const Array<OneD, const int> &traceBndMap
= m_fields[0]->GetTraceBndMap();
// Get physical values of the forward trace (from exp to phys)
Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
for (i = 0; i < nVariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
// Take into account that for PDE based shock capturing, eps = 0 at the
// wall.
int e, id1, id2, nBCEdgePts, eMax;
eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
for (e = 0; e < eMax; ++e)
{
nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
GetExp(e)->GetTotPoints();
id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
GetPhys_Offset(e);
id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
if (nVariables == m_spacedim+3)
{
NekDouble factor = 0.0;
NekDouble factor2 = 1.0;
Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
Vmath::Smul(nBCEdgePts,
factor,
&Fwd[nVariables-1][id2], 1,
&tmp2[0], 1);
Vmath::Vsub(nBCEdgePts,
&Fwd[nVariables-1][id2], 1,
&tmp2[0], 1,
&Fwd[nVariables-1][id2], 1);
Vmath::Smul(nBCEdgePts,
factor2,
&Fwd[nVariables-1][id2], 1,
&Fwd[nVariables-1][id2], 1);
}
// For 2D/3D, define: v* = v - 2(v.n)n
Array<OneD, NekDouble> tmp(nBCEdgePts, 0.0);
// Calculate (v.n)
for (i = 0; i < m_spacedim; ++i)
{
Vmath::Vvtvp(nBCEdgePts,
&Fwd[1+i][id2], 1,
&m_traceNormals[i][id2], 1,
&tmp[0], 1,
&tmp[0], 1);
}
// Calculate 2.0(v.n)
Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
// Calculate v* = v - 2.0(v.n)n
for (i = 0; i < m_spacedim; ++i)
{
Vmath::Vvtvp(nBCEdgePts,
&tmp[0], 1,
&m_traceNormals[i][id2], 1,
&Fwd[1+i][id2], 1,
&Fwd[1+i][id2], 1);
}
// Copy boundary adjusted values into the boundary expansion
for (i = 0; i < nVariables; ++i)
{
Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
&(m_fields[i]->GetBndCondExpansions()[bcRegion]->
UpdatePhys())[id1], 1);
}
}
}
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 249 of file CompressibleFlowSystem.cpp.

{
UnsteadySystem::v_GenerateSummary(s);
}
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::EulerCFE, and Nektar::NavierStokesCFE.

Definition at line 62 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::GetRiemannSolverFactory(), GetSmoothArtificialViscosity(), GetVecLocs(), GetViscousFluxVector(), GetViscousFluxVectorDeAlias(), m_advection, m_C1, m_C2, m_Cp, m_diffusion, m_eps_max, m_EqTypeStr, m_FacH, m_FacL, Nektar::SolverUtils::EquationSystem::m_fields, m_gamma, m_gasConstant, m_hFactor, Nektar::SolverUtils::UnsteadySystem::m_homoInitialFwd, m_Kappa, m_mu, m_mu0, m_pInf, m_Prandtl, 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_thermalConductivity, m_Twall, m_uInf, m_vecLocs, m_vInf, m_ViscosityType, and m_wInf.

{
UnsteadySystem::v_InitObject();
ASSERTL0(m_session->DefinesSolverInfo("UPWINDTYPE"),
"No UPWINDTYPE defined in session.");
// Do not forwards transform initial condition
m_homoInitialFwd = false;
// Set up locations of velocity vector.
m_vecLocs = Array<OneD, Array<OneD, NekDouble> >(1);
m_vecLocs[0] = Array<OneD, NekDouble>(m_spacedim);
for (int i = 0; i < m_spacedim; ++i)
{
m_vecLocs[0][i] = 1 + i;
}
// Get gamma parameter from session file.
ASSERTL0(m_session->DefinesParameter("Gamma"),
"Compressible flow sessions must define a Gamma parameter.");
m_session->LoadParameter("Gamma", m_gamma, 1.4);
// Get E0 parameter from session file.
ASSERTL0(m_session->DefinesParameter("pInf"),
"Compressible flow sessions must define a pInf parameter.");
m_session->LoadParameter("pInf", m_pInf, 101325);
// Get rhoInf parameter from session file.
ASSERTL0(m_session->DefinesParameter("rhoInf"),
"Compressible flow sessions must define a rhoInf parameter.");
m_session->LoadParameter("rhoInf", m_rhoInf, 1.225);
// Get uInf parameter from session file.
ASSERTL0(m_session->DefinesParameter("uInf"),
"Compressible flow sessions must define a uInf parameter.");
m_session->LoadParameter("uInf", m_uInf, 0.1);
// Get vInf parameter from session file.
if (m_spacedim == 2 || m_spacedim == 3)
{
ASSERTL0(m_session->DefinesParameter("vInf"),
"Compressible flow sessions must define a vInf parameter"
"for 2D/3D problems.");
m_session->LoadParameter("vInf", m_vInf, 0.0);
}
// Get wInf parameter from session file.
if (m_spacedim == 3)
{
ASSERTL0(m_session->DefinesParameter("wInf"),
"Compressible flow sessions must define a wInf parameter"
"for 3D problems.");
m_session->LoadParameter("wInf", m_wInf, 0.0);
}
m_session->LoadParameter ("GasConstant", m_gasConstant, 287.058);
m_session->LoadParameter ("Twall", m_Twall, 300.15);
m_session->LoadSolverInfo("ViscosityType", m_ViscosityType, "Constant");
m_session->LoadParameter ("mu", m_mu, 1.78e-05);
m_session->LoadParameter ("Skappa", m_Skappa, -2.048);
m_session->LoadParameter ("Kappa", m_Kappa, 0.0);
m_session->LoadParameter ("mu0", m_mu0, 1.0);
m_session->LoadParameter ("FL", m_FacL, 0.0);
m_session->LoadParameter ("FH", m_FacH, 0.0);
m_session->LoadParameter ("hFactor", m_hFactor, 1.0);
m_session->LoadParameter ("epsMax", m_eps_max, 1.0);
m_session->LoadParameter ("C1", m_C1, 3.0);
m_session->LoadParameter ("C2", m_C2, 5.0);
m_session->LoadSolverInfo("ShockCaptureType",
m_shockCaptureType, "Off");
m_session->LoadParameter ("thermalConductivity",
m_thermalConductivity, 0.0257);
m_EqTypeStr = m_session->GetSolverInfo("EQTYPE");
m_Cp = m_gamma / (m_gamma - 1.0) * m_gasConstant;
m_Prandtl = m_Cp * m_mu / m_thermalConductivity;
// Type of advection class to be used
switch(m_projectionType)
{
// Continuous field
case MultiRegions::eGalerkin:
{
ASSERTL0(false, "Continuous field not supported.");
break;
}
// Discontinuous field
case MultiRegions::eDiscontinuous:
{
string advName, diffName, riemName;
// Setting up advection and diffusion operators
m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
m_session->LoadSolverInfo("DiffusionType", diffName, "LDGNS");
m_advection = SolverUtils::GetAdvectionFactory()
.CreateInstance(advName, advName);
m_diffusion = SolverUtils::GetDiffusionFactory()
.CreateInstance(diffName, diffName);
if (m_specHP_dealiasing)
{
m_advection->SetFluxVector(&CompressibleFlowSystem::
GetFluxVectorDeAlias, this);
m_diffusion->SetFluxVectorNS(
&CompressibleFlowSystem::GetViscousFluxVectorDeAlias,
this);
}
else
{
m_advection->SetFluxVector (&CompressibleFlowSystem::
GetFluxVector, this);
m_diffusion->SetFluxVectorNS(&CompressibleFlowSystem::
GetViscousFluxVector, this);
}
if (m_shockCaptureType=="Smooth" && m_EqTypeStr=="EulerADCFE")
{
m_advection->SetFluxVector(&CompressibleFlowSystem::
GetFluxVector, this);
m_diffusion->SetArtificialDiffusionVector(
&CompressibleFlowSystem::GetSmoothArtificialViscosity, this);
}
if (m_shockCaptureType=="NonSmooth" && m_EqTypeStr=="EulerADCFE")
{
m_advection->SetFluxVector(&CompressibleFlowSystem::
GetFluxVector, this);
m_diffusion->SetArtificialDiffusionVector(
&CompressibleFlowSystem::GetArtificialDynamicViscosity, this);
}
// Setting up Riemann solver for advection operator
m_session->LoadSolverInfo("UpwindType", riemName, "Average");
m_riemannSolver = SolverUtils::GetRiemannSolverFactory()
.CreateInstance(riemName);
// Setting up upwind solver for diffusion operator
m_riemannSolverLDG = SolverUtils::GetRiemannSolverFactory()
.CreateInstance("UpwindLDG");
// Setting up parameters for advection operator Riemann solver
m_riemannSolver->SetParam (
"gamma", &CompressibleFlowSystem::GetGamma, this);
m_riemannSolver->SetAuxVec(
"vecLocs", &CompressibleFlowSystem::GetVecLocs, this);
m_riemannSolver->SetVector(
"N", &CompressibleFlowSystem::GetNormals, this);
// Setting up parameters for diffusion operator Riemann solver
m_riemannSolverLDG->SetParam (
"gamma", &CompressibleFlowSystem::GetGamma, this);
m_riemannSolverLDG->SetVector(
"vecLocs", &CompressibleFlowSystem::GetVecLocs, this);
m_riemannSolverLDG->SetVector(
"N", &CompressibleFlowSystem::GetNormals, this);
// Concluding initialisation of advection / diffusion operators
m_advection->SetRiemannSolver (m_riemannSolver);
m_diffusion->SetRiemannSolver (m_riemannSolverLDG);
m_advection->InitObject (m_session, m_fields);
m_diffusion->InitObject (m_session, m_fields);
break;
}
default:
{
ASSERTL0(false, "Unsupported projection type.");
break;
}
}
}
virtual void Nektar::CompressibleFlowSystem::v_SetInitialConditions ( NekDouble  initialtime = 0.0,
bool  dumpInitialConditions = true,
const int  domain = 0 
)
inlineprotectedvirtual

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.

Definition at line 231 of file CompressibleFlowSystem.h.

{
}
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 257 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 Nektar::EulerADCFE::SetBoundaryConditions(), and Nektar::EulerCFE::SetBoundaryConditions().

{
int i;
int nTracePts = GetTraceTotPoints();
int nVariables = physarray.num_elements();
const Array<OneD, const int> &traceBndMap
= m_fields[0]->GetTraceBndMap();
// Get physical values of the forward trace
Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
for (i = 0; i < nVariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
// Adjust the physical values of the trace to take
// user defined boundaries into account
int e, id1, id2, nBCEdgePts, eMax;
eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
for (e = 0; e < eMax; ++e)
{
nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
GetExp(e)->GetTotPoints();
id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
GetPhys_Offset(e);
id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
// Boundary condition for epsilon term.
if (nVariables == m_spacedim+3)
{
NekDouble factor = 1.0;
NekDouble factor2 = 1.0;
Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
Vmath::Smul(nBCEdgePts,
factor,
&Fwd[nVariables-1][id2], 1,
&tmp2[0], 1);
Vmath::Vsub(nBCEdgePts,
&Fwd[nVariables-1][id2], 1,
&tmp2[0], 1,
&Fwd[nVariables-1][id2], 1);
Vmath::Smul(nBCEdgePts,
factor2,
&Fwd[nVariables-1][id2], 1,
&Fwd[nVariables-1][id2], 1);
}
// For 2D/3D, define: v* = v - 2(v.n)n
Array<OneD, NekDouble> tmp(nBCEdgePts, 0.0);
// Calculate (v.n)
for (i = 0; i < m_spacedim; ++i)
{
Vmath::Vvtvp(nBCEdgePts,
&Fwd[1+i][id2], 1,
&m_traceNormals[i][id2], 1,
&tmp[0], 1,
&tmp[0], 1);
}
// Calculate 2.0(v.n)
Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
// Calculate v* = v - 2.0(v.n)n
for (i = 0; i < m_spacedim; ++i)
{
Vmath::Vvtvp(nBCEdgePts,
&tmp[0], 1,
&m_traceNormals[i][id2], 1,
&Fwd[1+i][id2], 1,
&Fwd[1+i][id2], 1);
}
// Copy boundary adjusted values into the boundary expansion
for (i = 0; i < nVariables; ++i)
{
Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
&(m_fields[i]->GetBndCondExpansions()[bcRegion]->
UpdatePhys())[id1], 1);
}
}
}
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 354 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 Nektar::NavierStokesCFE::SetBoundaryConditions().

{
int i;
int nTracePts = GetTraceTotPoints();
int nVariables = physarray.num_elements();
const Array<OneD, const int> &traceBndMap
= m_fields[0]->GetTraceBndMap();
// Get physical values of the forward trace
Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
for (i = 0; i < nVariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
// Take into account that for PDE based shock capturing, eps = 0 at the
// wall. Adjust the physical values of the trace to take user defined
// boundaries into account
int e, id1, id2, nBCEdgePts, eMax;
eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
for (e = 0; e < eMax; ++e)
{
nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
GetExp(e)->GetTotPoints();
id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->
GetPhys_Offset(e);
id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
if (nVariables == m_spacedim+3)
{
NekDouble factor = 0.0;
NekDouble factor2 = 1.0;
Array<OneD, NekDouble > tmp2(nBCEdgePts, 0.0);
Vmath::Smul(nBCEdgePts,
factor,
&Fwd[nVariables-1][id2], 1,
&tmp2[0], 1);
Vmath::Vsub(nBCEdgePts,
&Fwd[nVariables-1][id2], 1,
&tmp2[0], 1,
&Fwd[nVariables-1][id2], 1);
Vmath::Smul(nBCEdgePts,
factor2,
&Fwd[nVariables-1][id2], 1,
&Fwd[nVariables-1][id2], 1);
}
for (i = 0; i < m_spacedim; i++)
{
Vmath::Neg(nBCEdgePts, &Fwd[i+1][id2], 1);
}
// Copy boundary adjusted values into the boundary expansion
for (i = 0; i < nVariables; ++i)
{
Vmath::Vcopy(nBCEdgePts, &Fwd[i][id2], 1,
&(m_fields[i]->GetBndCondExpansions()[bcRegion]->
UpdatePhys())[id1], 1);
}
}
}

Friends And Related Function Documentation

friend class MemoryManager< CompressibleFlowSystem >
friend

Definition at line 70 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 79 of file CompressibleFlowSystem.h.

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

Definition at line 94 of file CompressibleFlowSystem.h.

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

NekDouble Nektar::CompressibleFlowSystem::m_C1
protected

Definition at line 117 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_C2
protected

Definition at line 118 of file CompressibleFlowSystem.h.

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

NekDouble Nektar::CompressibleFlowSystem::m_Cp
protected

Definition at line 116 of file CompressibleFlowSystem.h.

Referenced by GetViscousFluxVector(), and v_InitObject().

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

Definition at line 95 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 114 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

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

Definition at line 107 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_FacH
protected

Definition at line 113 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_FacL
protected

Definition at line 112 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_gamma
protected

Definition at line 97 of file CompressibleFlowSystem.h.

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

NekDouble Nektar::CompressibleFlowSystem::m_gasConstant
protected

Definition at line 103 of file CompressibleFlowSystem.h.

Referenced by GetGasConstant(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_hFactor
protected

Definition at line 119 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_Kappa
protected

Definition at line 110 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_mu
protected

Definition at line 108 of file CompressibleFlowSystem.h.

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

NekDouble Nektar::CompressibleFlowSystem::m_mu0
protected

Definition at line 111 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

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

Definition at line 122 of file CompressibleFlowSystem.h.

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

Definition at line 121 of file CompressibleFlowSystem.h.

NekDouble Nektar::CompressibleFlowSystem::m_pInf
protected

Definition at line 98 of file CompressibleFlowSystem.h.

Referenced by RiemannInvariantBC(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_Prandtl
protected

Definition at line 120 of file CompressibleFlowSystem.h.

Referenced by GetViscousFluxVector(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_rhoInf
protected

Definition at line 99 of file CompressibleFlowSystem.h.

Referenced by RiemannInvariantBC(), and v_InitObject().

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

Definition at line 92 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

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

Definition at line 93 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

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

Definition at line 106 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 109 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

bool Nektar::CompressibleFlowSystem::m_smoothDiffusion
protected

Definition at line 123 of file CompressibleFlowSystem.h.

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

NekDouble Nektar::CompressibleFlowSystem::m_thermalConductivity
protected

Definition at line 115 of file CompressibleFlowSystem.h.

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

NekDouble Nektar::CompressibleFlowSystem::m_Twall
protected

Definition at line 104 of file CompressibleFlowSystem.h.

Referenced by v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_uInf
protected

Definition at line 100 of file CompressibleFlowSystem.h.

Referenced by RiemannInvariantBC(), and v_InitObject().

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

Definition at line 96 of file CompressibleFlowSystem.h.

Referenced by GetVecLocs(), and v_InitObject().

NekDouble Nektar::CompressibleFlowSystem::m_vInf
protected

Definition at line 101 of file CompressibleFlowSystem.h.

Referenced by RiemannInvariantBC(), and v_InitObject().

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

Definition at line 105 of file CompressibleFlowSystem.h.

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

NekDouble Nektar::CompressibleFlowSystem::m_wInf
protected

Definition at line 102 of file CompressibleFlowSystem.h.

Referenced by RiemannInvariantBC(), and v_InitObject().

Generated on Sun Mar 15 2015 00:11:50 for Nektar++ by   doxygen 1.8.1.2