Nektar++
Static Public Member Functions | Public Attributes | Static Public Attributes | Protected Member Functions | Protected Attributes | List of all members
Nektar::SolverUtils::DiffusionLFR Class Reference

#include <DiffusionLFR.h>

Inheritance diagram for Nektar::SolverUtils::DiffusionLFR:
[legend]

Static Public Member Functions

static DiffusionSharedPtr create (std::string diffType)
 

Public Attributes

Array< OneD, NekDouble > m_jac
 
Array< OneD, Array< OneD, NekDouble > > m_gmat
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
 
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
 
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
 
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
 
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
 
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
 
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
 
DNekMatSharedPtr m_Ixm
 
DNekMatSharedPtr m_Ixp
 
- Public Attributes inherited from Nektar::SolverUtils::Diffusion
Array< OneD, NekDouble > m_divVel
 Params for Ducros sensor. More...
 
Array< OneD, NekDouble > m_divVelSquare
 
Array< OneD, NekDouble > m_curlVelSquare
 

Static Public Attributes

static std::string type []
 

Protected Member Functions

 DiffusionLFR (std::string diffType)
 DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion term. The implementation is only for segments, quadrilaterals and hexahedra at the moment. More...
 
virtual void v_InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Initiliase DiffusionLFR objects and store them before starting the time-stepping. More...
 
virtual void v_SetupMetrics (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform the fluxes at the interfaces of each element from the physical space to the standard space). More...
 
virtual void v_SetupCFunctions (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–72. More...
 
virtual void v_Diffuse (const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd)
 Calculate FR Diffusion for the linear problems using an LDG interface flux. More...
 
virtual void v_NumFluxforScalar (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
 Builds the numerical flux for the 1st order derivatives. More...
 
virtual void v_WeakPenaltyforScalar (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const Array< OneD, const NekDouble > &ufield, Array< OneD, NekDouble > &penaltyflux)
 Imposes appropriate bcs for the 1st order derivatives. More...
 
virtual void v_NumFluxforVector (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
 Build the numerical flux for the 2nd order derivatives. More...
 
virtual void v_WeakPenaltyforVector (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux, NekDouble C11)
 Imposes appropriate bcs for the 2nd order derivatives. More...
 
virtual void v_DerCFlux_1D (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, const NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
 Compute the derivative of the corrective flux for 1D problems. More...
 
virtual void v_DerCFlux_2D (const int nConvectiveFields, const int direction, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
 Compute the derivative of the corrective flux wrt a given coordinate for 2D problems. More...
 
virtual void v_DivCFlux_2D (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
 Compute the divergence of the corrective flux for 2D problems. More...
 
virtual void v_DivCFlux_2D_Gauss (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
 Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre". More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::Diffusion
virtual SOLVER_UTILS_EXPORT void v_DiffuseCoeffs (const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd)
 
virtual SOLVER_UTILS_EXPORT void v_DiffuseCoeffs (const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &vFwd, const Array< OneD, Array< OneD, NekDouble > > &vBwd, TensorOfArray3D< NekDouble > &qfield, Array< OneD, int > &nonZeroIndex)
 
virtual SOLVER_UTILS_EXPORT void v_ConsVarAveJump (const std::size_t nConvectiveFields, const size_t npnts, const Array< OneD, const Array< OneD, NekDouble > > &vFwd, const Array< OneD, const Array< OneD, NekDouble > > &vBwd, Array< OneD, Array< OneD, NekDouble > > &aver, Array< OneD, Array< OneD, NekDouble > > &jump)
 
virtual SOLVER_UTILS_EXPORT void v_DiffuseCalcDerivative (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &qfields, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd)
 Diffusion Flux, calculate the physical derivatives. More...
 
virtual SOLVER_UTILS_EXPORT void v_DiffuseVolumeFlux (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfields, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, int > &nonZeroIndex)
 Diffusion Volume Flux. More...
 
virtual SOLVER_UTILS_EXPORT void v_DiffuseTraceFlux (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &Qfields, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, Array< OneD, NekDouble > > &TraceFlux, const Array< OneD, Array< OneD, NekDouble >> &pFwd, const Array< OneD, Array< OneD, NekDouble >> &pBwd, Array< OneD, int > &nonZeroIndex)
 Diffusion term Trace Flux. More...
 
virtual SOLVER_UTILS_EXPORT void v_DiffuseTraceSymmFlux (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, const TensorOfArray3D< NekDouble > &qfield, const TensorOfArray3D< NekDouble > &VolumeFlux, TensorOfArray3D< NekDouble > &SymmFlux, const Array< OneD, Array< OneD, NekDouble >> &pFwd, const Array< OneD, Array< OneD, NekDouble >> &pBwd, Array< OneD, int > &nonZeroIndex, Array< OneD, Array< OneD, NekDouble > > &solution_Aver, Array< OneD, Array< OneD, NekDouble > > &solution_jump)
 
virtual SOLVER_UTILS_EXPORT void v_AddDiffusionSymmFluxToCoeff (const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &qfield, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd)
 
virtual SOLVER_UTILS_EXPORT void v_AddDiffusionSymmFluxToPhys (const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &qfield, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd)
 
virtual void v_SetHomoDerivs (Array< OneD, Array< OneD, NekDouble > > &deriv)
 
virtual TensorOfArray3D< NekDouble > & v_GetFluxTensor ()
 
virtual SOLVER_UTILS_EXPORT const Array< OneD, const Array< OneD, NekDouble > > & v_GetTraceNormal ()
 
virtual SOLVER_UTILS_EXPORT void v_GetPrimVar (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &primVar)
 Compute primary derivatives. More...
 
SOLVER_UTILS_EXPORT void GetDivCurl (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const TensorOfArray3D< NekDouble > &pVarDer)
 Compute divergence and curl squared. More...
 

Protected Attributes

Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 
LibUtilities::SessionReaderSharedPtr m_session
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
 
Array< OneD, Array< OneD, NekDouble > > m_IF2
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
 
Array< OneD, Array< OneD, NekDouble > > m_divFD
 
Array< OneD, Array< OneD, NekDouble > > m_divFC
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
 
std::string m_diffType
 
- Protected Attributes inherited from Nektar::SolverUtils::Diffusion
DiffusionFluxVecCB m_fluxVector
 
DiffusionFluxVecCBNS m_fluxVectorNS
 
DiffusionFluxPenaltyNS m_fluxPenaltyNS
 
DiffusionArtificialDiffusion m_ArtificialDiffusionVector
 
DiffusionFluxCons m_FunctorDiffusionfluxCons
 
DiffusionFluxCons m_FunctorDiffusionfluxConsTrace
 
SpecialBndTreat m_SpecialBndTreat
 
CalcViscosity m_CalcViscosity
 
DiffusionSymmFluxCons m_FunctorSymmetricfluxCons
 
NekDouble m_time =0.0
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::SolverUtils::Diffusion
virtual SOLVER_UTILS_EXPORT ~Diffusion ()
 
SOLVER_UTILS_EXPORT void InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 
SOLVER_UTILS_EXPORT void Diffuse (const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray)
 
SOLVER_UTILS_EXPORT void DiffuseCoeffs (const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray)
 Similar with Diffusion::Diffuse(): calculate diffusion flux The difference is in the outarray: it is the coefficients of basis for DiffuseCoeffs() it is the physics on quadrature points for Diffuse() More...
 
SOLVER_UTILS_EXPORT void Diffuse (const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, NekDouble time, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray)
 
SOLVER_UTILS_EXPORT void DiffuseCoeffs (const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd, TensorOfArray3D< NekDouble > &qfield, Array< OneD, int > &nonZeroIndex)
 
SOLVER_UTILS_EXPORT void DiffuseCoeffs (const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, NekDouble time, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray)
 
SOLVER_UTILS_EXPORT void DiffuseCalcDerivative (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfields, const Array< OneD, Array< OneD, NekDouble >> &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble >> &pBwd=NullNekDoubleArrayOfArray)
 
SOLVER_UTILS_EXPORT void DiffuseVolumeFlux (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfields, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, int > &nonZeroIndex=NullInt1DArray)
 Diffusion Volume FLux. More...
 
SOLVER_UTILS_EXPORT void DiffuseTraceFlux (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfields, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, Array< OneD, NekDouble > > &TraceFlux, const Array< OneD, Array< OneD, NekDouble >> &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble >> &pBwd=NullNekDoubleArrayOfArray, Array< OneD, int > &nonZeroIndex=NullInt1DArray)
 Diffusion term Trace Flux. More...
 
SOLVER_UTILS_EXPORT void DiffuseTraceSymmFlux (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, const TensorOfArray3D< NekDouble > &qfield, const TensorOfArray3D< NekDouble > &VolumeFlux, TensorOfArray3D< NekDouble > &SymmFlux, const Array< OneD, Array< OneD, NekDouble >> &pFwd, const Array< OneD, Array< OneD, NekDouble >> &pBwd, Array< OneD, int > &nonZeroIndex, Array< OneD, Array< OneD, NekDouble > > &solution_Aver=NullNekDoubleArrayOfArray, Array< OneD, Array< OneD, NekDouble > > &solution_jump=NullNekDoubleArrayOfArray)
 
SOLVER_UTILS_EXPORT void AddDiffusionSymmFluxToCoeff (const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &qfield, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd)
 Add symmetric flux to field in coeff space. More...
 
SOLVER_UTILS_EXPORT void AddDiffusionSymmFluxToPhys (const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &qfield, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd)
 Add symmetric flux to field in coeff physical space. More...
 
SOLVER_UTILS_EXPORT void AddSymmFluxIntegralToOffDiag (const int nConvectiveFields, const int nDim, const int nPts, const int nTracePts, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const int > &nonZeroIndex, TensorOfArray3D< NekDouble > &Fwdflux, TensorOfArray3D< NekDouble > &Bwdflux, Array< OneD, Array< OneD, NekDouble > > &outarray)
 
SOLVER_UTILS_EXPORT void FluxVec (TensorOfArray3D< NekDouble > &fluxvector)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetFluxVector (FuncPointerT func, ObjectPointerT obj)
 
SOLVER_UTILS_EXPORT void SetFluxVector (DiffusionFluxVecCB fluxVector)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetFluxVectorNS (FuncPointerT func, ObjectPointerT obj)
 
void SetFluxVectorNS (DiffusionFluxVecCBNS fluxVector)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetFluxPenaltyNS (FuncPointerT func, ObjectPointerT obj)
 
void SetFluxPenaltyNS (DiffusionFluxPenaltyNS flux)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetDiffusionFluxCons (FuncPointerT func, ObjectPointerT obj)
 
void SetDiffusionFluxCons (DiffusionFluxCons flux)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetDiffusionFluxConsTrace (FuncPointerT func, ObjectPointerT obj)
 
void SetDiffusionFluxConsTrace (DiffusionFluxCons flux)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetArtificialDiffusionVector (FuncPointerT func, ObjectPointerT obj)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetSpecialBndTreat (FuncPointerT func, ObjectPointerT obj)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetCalcViscosity (FuncPointerT func, ObjectPointerT obj)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetDiffusionSymmFluxCons (FuncPointerT func, ObjectPointerT obj)
 
void SetHomoDerivs (Array< OneD, Array< OneD, NekDouble > > &deriv)
 
virtual TensorOfArray3D< NekDouble > & GetFluxTensor ()
 
SOLVER_UTILS_EXPORT void GetAVmu (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &muvar, Array< OneD, NekDouble > &MuVarTrace)
 Get the mu of artifical viscosity(AV) More...
 
SOLVER_UTILS_EXPORT void ConsVarAveJump (const std::size_t nConvectiveFields, const size_t npnts, const Array< OneD, const Array< OneD, NekDouble > > &vFwd, const Array< OneD, const Array< OneD, NekDouble > > &vBwd, Array< OneD, Array< OneD, NekDouble > > &aver, Array< OneD, Array< OneD, NekDouble > > &jump)
 Get the average and jump value of conservative variables on trace. More...
 
SOLVER_UTILS_EXPORT const Array< OneD, const Array< OneD, NekDouble > > & GetTraceNormal ()
 Get trace normal. More...
 

Detailed Description

Definition at line 44 of file DiffusionLFR.h.

Constructor & Destructor Documentation

◆ DiffusionLFR()

Nektar::SolverUtils::DiffusionLFR::DiffusionLFR ( std::string  diffType)
protected

DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion term. The implementation is only for segments, quadrilaterals and hexahedra at the moment.

Todo:
Extension to triangles, tetrahedra and other shapes. (Long term objective)

Definition at line 71 of file DiffusionLFR.cpp.

71  :m_diffType(diffType)
72  {
73  }

Referenced by create().

Member Function Documentation

◆ create()

static DiffusionSharedPtr Nektar::SolverUtils::DiffusionLFR::create ( std::string  diffType)
inlinestatic

Definition at line 47 of file DiffusionLFR.h.

48  {
49  return DiffusionSharedPtr(new DiffusionLFR(diffType));
50  }
Nektar::SolverUtils::DiffusionLFR::DiffusionLFR
DiffusionLFR(std::string diffType)
DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion term....
Definition: DiffusionLFR.cpp:71
Nektar::SolverUtils::DiffusionSharedPtr
std::shared_ptr< SolverUtils::Diffusion > DiffusionSharedPtr
A shared pointer to an EquationSystem object.
Definition: Diffusion.h:627

References DiffusionLFR().

◆ v_DerCFlux_1D()

void Nektar::SolverUtils::DiffusionLFR::v_DerCFlux_1D ( const int  nConvectiveFields,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, const NekDouble > &  flux,
const Array< OneD, const NekDouble > &  iFlux,
Array< OneD, NekDouble > &  derCFlux 
)
protectedvirtual

Compute the derivative of the corrective flux for 1D problems.

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxVolumetric flux in the physical space.
iFluxNumerical interface flux in physical space.
derCFluxDerivative of the corrective flux.

Definition at line 1474 of file DiffusionLFR.cpp.

1480  {
1481  boost::ignore_unused(nConvectiveFields);
1482 
1483  int n;
1484  int nLocalSolutionPts, phys_offset, t_offset;
1485 
1486  Array<OneD, NekDouble> auxArray1, auxArray2;
1487  Array<TwoD, const NekDouble> gmat;
1488  Array<OneD, const NekDouble> jac;
1489 
1490  LibUtilities::PointsKeyVector ptsKeys;
1491  LibUtilities::BasisSharedPtr Basis;
1492  Basis = fields[0]->GetExp(0)->GetBasis(0);
1493 
1494  int nElements = fields[0]->GetExpSize();
1495  int nSolutionPts = fields[0]->GetTotPoints();
1496 
1497  vector<bool> negatedFluxNormal =
1498  std::static_pointer_cast<MultiRegions::DisContField>(
1499  fields[0])->GetNegatedFluxNormal();
1500 
1501  // Arrays to store the derivatives of the correction flux
1502  Array<OneD, NekDouble> DCL(nSolutionPts/nElements, 0.0);
1503  Array<OneD, NekDouble> DCR(nSolutionPts/nElements, 0.0);
1504 
1505  // Arrays to store the intercell numerical flux jumps
1506  Array<OneD, NekDouble> JumpL(nElements);
1507  Array<OneD, NekDouble> JumpR(nElements);
1508 
1509  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1510  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1511 
1512  for (n = 0; n < nElements; ++n)
1513  {
1514  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1515  phys_offset = fields[0]->GetPhys_Offset(n);
1516 
1517  Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1518  Array<OneD, NekDouble> tmpFluxVertex(1,0.0);
1519  Vmath::Vcopy(nLocalSolutionPts,
1520  &flux[phys_offset], 1,
1521  &tmparrayX1[0], 1);
1522 
1523  fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0],
1524  tmparrayX1,
1525  tmpFluxVertex);
1526 
1527  t_offset = fields[0]->GetTrace()
1528  ->GetPhys_Offset(elmtToTrace[n][0]->GetElmtId());
1529 
1530  if(negatedFluxNormal[2*n])
1531  {
1532  JumpL[n] = iFlux[t_offset] - tmpFluxVertex[0];
1533  }
1534  else
1535  {
1536  JumpL[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1537  }
1538 
1539  t_offset = fields[0]->GetTrace()
1540  ->GetPhys_Offset(elmtToTrace[n][1]->GetElmtId());
1541 
1542  fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1],
1543  tmparrayX1,
1544  tmpFluxVertex);
1545  if(negatedFluxNormal[2*n+1])
1546  {
1547  JumpR[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1548  }
1549  else
1550  {
1551  JumpR[n] = iFlux[t_offset] - tmpFluxVertex[0];
1552  }
1553  }
1554 
1555  for (n = 0; n < nElements; ++n)
1556  {
1557  ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1558  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1559  phys_offset = fields[0]->GetPhys_Offset(n);
1560  jac = fields[0]->GetExp(n)
1561  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
1562  ->GetMetricInfo()->GetJac(ptsKeys);
1563 
1564  JumpL[n] = JumpL[n] * jac[0];
1565  JumpR[n] = JumpR[n] * jac[0];
1566 
1567  // Left jump multiplied by left derivative of C function
1568  Vmath::Smul(nLocalSolutionPts, JumpL[n],
1569  &m_dGL_xi1[n][0], 1, &DCL[0], 1);
1570 
1571  // Right jump multiplied by right derivative of C function
1572  Vmath::Smul(nLocalSolutionPts, JumpR[n],
1573  &m_dGR_xi1[n][0], 1, &DCR[0], 1);
1574 
1575  // Assembling derivative of the correction flux
1576  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1577  &derCFlux[phys_offset], 1);
1578  }
1579  }
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:63
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:62
Nektar::LibUtilities::BasisSharedPtr
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:74
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
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:322
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199

References Nektar::LocalRegions::Expansion::GetMetricInfo(), m_dGL_xi1, m_dGR_xi1, Vmath::Smul(), Vmath::Vadd(), and Vmath::Vcopy().

Referenced by v_Diffuse().

◆ v_DerCFlux_2D()

void Nektar::SolverUtils::DiffusionLFR::v_DerCFlux_2D ( const int  nConvectiveFields,
const int  direction,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, const NekDouble > &  flux,
const Array< OneD, NekDouble > &  iFlux,
Array< OneD, NekDouble > &  derCFlux 
)
protectedvirtual

Compute the derivative of the corrective flux wrt a given coordinate for 2D problems.

There could be a bug for deformed elements since first derivatives are not exactly the same results of DiffusionLDG as expected

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxVolumetric flux in the physical space.
numericalFluxNumerical interface flux in physical space.
derCFluxDerivative of the corrective flux.
Todo:
: Switch on shapes eventually here.

Definition at line 1596 of file DiffusionLFR.cpp.

1603  {
1604  boost::ignore_unused(nConvectiveFields);
1605 
1606  int n, e, i, j, cnt;
1607 
1608  Array<TwoD, const NekDouble> gmat;
1609  Array<OneD, const NekDouble> jac;
1610 
1611  int nElements = fields[0]->GetExpSize();
1612 
1613  int trace_offset, phys_offset;
1614  int nLocalSolutionPts;
1615  int nquad0, nquad1;
1616  int nEdgePts;
1617 
1618  LibUtilities::PointsKeyVector ptsKeys;
1619  Array<OneD, NekDouble> auxArray1, auxArray2;
1620  Array<OneD, LibUtilities::BasisSharedPtr> base;
1621  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1622  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1623 
1624  // Loop on the elements
1625  for (n = 0; n < nElements; ++n)
1626  {
1627  // Offset of the element on the global vector
1628  phys_offset = fields[0]->GetPhys_Offset(n);
1629  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1630  ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1631 
1632  jac = fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1633  ->GetGeom2D()->GetMetricInfo()->GetJac(ptsKeys);
1634 
1635  base = fields[0]->GetExp(n)->GetBase();
1636  nquad0 = base[0]->GetNumPoints();
1637  nquad1 = base[1]->GetNumPoints();
1638 
1639  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1640  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1641  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1642  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1643 
1644  // Loop on the edges
1645  for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1646  {
1647  // Number of edge points of edge 'e'
1648  nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1649 
1650  // Array for storing volumetric fluxes on each edge
1651  Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1652 
1653  // Offset of the trace space correspondent to edge 'e'
1654  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1655  elmtToTrace[n][e]->GetElmtId());
1656 
1657  // Get the normals of edge 'e'
1658  //const Array<OneD, const Array<OneD, NekDouble> > &normals =
1659  //fields[0]->GetExp(n)->GetTraceNormal(e);
1660 
1661  // Extract the edge values of the volumetric fluxes
1662  // on edge 'e' and order them accordingly to the order
1663  // of the trace space
1664  fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1665  flux + phys_offset,
1666  auxArray1 = tmparray);
1667 
1668  // Compute the fluxJumps per each edge 'e' and each
1669  // flux point
1670  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1671  Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1,
1672  &tmparray[0], 1, &fluxJumps[0], 1);
1673 
1674  // Check the ordering of the fluxJumps and reverse
1675  // it in case of backward definition of edge 'e'
1676  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1677  StdRegions::eBackwards)
1678  {
1679  Vmath::Reverse(nEdgePts,
1680  &fluxJumps[0], 1,
1681  &fluxJumps[0], 1);
1682  }
1683 
1684  // Deformed elements
1685  if (fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1686  ->GetGeom2D()->GetMetricInfo()->GetGtype()
1687  == SpatialDomains::eDeformed)
1688  {
1689  // Extract the Jacobians along edge 'e'
1690  Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1691 
1692  fields[0]->GetExp(n)->GetTracePhysVals(
1693  e, elmtToTrace[n][e],
1694  jac, auxArray1 = jacEdge);
1695 
1696  // Check the ordering of the fluxJumps and reverse
1697  // it in case of backward definition of edge 'e'
1698  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1699  StdRegions::eBackwards)
1700  {
1701  Vmath::Reverse(nEdgePts,
1702  &jacEdge[0], 1,
1703  &jacEdge[0], 1);
1704  }
1705 
1706  // Multiply the fluxJumps by the edge 'e' Jacobians
1707  // to bring the fluxJumps into the standard space
1708  for (j = 0; j < nEdgePts; j++)
1709  {
1710  fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1711  }
1712  }
1713  // Non-deformed elements
1714  else
1715  {
1716  // Multiply the fluxJumps by the edge 'e' Jacobians
1717  // to bring the fluxJumps into the standard space
1718  Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1,
1719  fluxJumps, 1);
1720  }
1721 
1722  // Multiply jumps by derivatives of the correction functions
1723  // All the quntities at this level should be defined into
1724  // the standard space
1725  switch (e)
1726  {
1727  case 0:
1728  for (i = 0; i < nquad0; ++i)
1729  {
1730  // Multiply fluxJumps by Q factors
1731  //fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1732 
1733  for (j = 0; j < nquad1; ++j)
1734  {
1735  cnt = i + j*nquad0;
1736  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1737  }
1738  }
1739  break;
1740  case 1:
1741  for (i = 0; i < nquad1; ++i)
1742  {
1743  // Multiply fluxJumps by Q factors
1744  //fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1745 
1746  for (j = 0; j < nquad0; ++j)
1747  {
1748  cnt = (nquad0)*i + j;
1749  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1750  }
1751  }
1752  break;
1753  case 2:
1754  for (i = 0; i < nquad0; ++i)
1755  {
1756  // Multiply fluxJumps by Q factors
1757  //fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1758 
1759  for (j = 0; j < nquad1; ++j)
1760  {
1761  cnt = j*nquad0 + i;
1762  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1763  }
1764  }
1765  break;
1766  case 3:
1767  for (i = 0; i < nquad1; ++i)
1768  {
1769  // Multiply fluxJumps by Q factors
1770  //fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1771  for (j = 0; j < nquad0; ++j)
1772  {
1773  cnt = j + i*nquad0;
1774  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1775  }
1776  }
1777  break;
1778 
1779  default:
1780  ASSERTL0(false, "edge value (< 3) is out of range");
1781  break;
1782  }
1783  }
1784 
1785  // Summing the component of the flux for calculating the
1786  // derivatives per each direction
1787  if (direction == 0)
1788  {
1789  Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1,
1790  &divCFluxE3[0], 1, &derCFlux[phys_offset], 1);
1791  }
1792  else if (direction == 1)
1793  {
1794  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1795  &divCFluxE2[0], 1, &derCFlux[phys_offset], 1);
1796  }
1797  }
1798  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFR.h:64
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFR.h:65
Nektar::SpatialDomains::eDeformed
@ eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:56
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1226
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:372

References ASSERTL0, Nektar::StdRegions::eBackwards, Nektar::SpatialDomains::eDeformed, Nektar::LocalRegions::Expansion::GetMetricInfo(), m_dGL_xi1, m_dGL_xi2, m_dGR_xi1, m_dGR_xi2, Vmath::Reverse(), Vmath::Smul(), Vmath::Vadd(), and Vmath::Vsub().

Referenced by v_Diffuse().

◆ v_Diffuse()

void Nektar::SolverUtils::DiffusionLFR::v_Diffuse ( const std::size_t  nConvectiveFields,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, Array< OneD, NekDouble > > &  inarray,
Array< OneD, Array< OneD, NekDouble > > &  outarray,
const Array< OneD, Array< OneD, NekDouble > > &  pFwd,
const Array< OneD, Array< OneD, NekDouble > > &  pBwd 
)
protectedvirtual

Calculate FR Diffusion for the linear problems using an LDG interface flux.

Reimplemented from Nektar::SolverUtils::Diffusion.

Definition at line 853 of file DiffusionLFR.cpp.

860  {
861  boost::ignore_unused(pFwd, pBwd);
862 
863  int i, j, n;
864  int phys_offset;
865  //Array<TwoD, const NekDouble> gmat;
866  //Array<OneD, const NekDouble> jac;
867  Array<OneD, NekDouble> auxArray1, auxArray2;
868 
869  Array<OneD, LibUtilities::BasisSharedPtr> Basis;
870  Basis = fields[0]->GetExp(0)->GetBase();
871 
872  int nElements = fields[0]->GetExpSize();
873  int nDim = fields[0]->GetCoordim(0);
874  int nSolutionPts = fields[0]->GetTotPoints();
875  int nCoeffs = fields[0]->GetNcoeffs();
876 
877  Array<OneD, Array<OneD, NekDouble> > outarrayCoeff(nConvectiveFields);
878  for (i = 0; i < nConvectiveFields; ++i)
879  {
880  outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
881  }
882 
883  // Compute interface numerical fluxes for inarray in physical space
884  v_NumFluxforScalar(fields, inarray, m_IF1);
885 
886  switch(nDim)
887  {
888  // 1D problems
889  case 1:
890  {
891  for (i = 0; i < nConvectiveFields; ++i)
892  {
893  // Computing the physical first-order discountinuous
894  // derivative
895  for (n = 0; n < nElements; n++)
896  {
897  phys_offset = fields[0]->GetPhys_Offset(n);
898 
899  fields[i]->GetExp(n)->PhysDeriv(0,
900  auxArray1 = inarray[i] + phys_offset,
901  auxArray2 = m_DU1[i][0] + phys_offset);
902  }
903 
904  // Computing the standard first-order correction
905  // derivative
906  v_DerCFlux_1D(nConvectiveFields, fields, inarray[i],
907  m_IF1[i][0], m_DFC1[i][0]);
908 
909  // Back to the physical space using global operations
910  Vmath::Vdiv(nSolutionPts, &m_DFC1[i][0][0], 1,
911  &m_jac[0], 1, &m_DFC1[i][0][0], 1);
912  Vmath::Vdiv(nSolutionPts, &m_DFC1[i][0][0], 1,
913  &m_jac[0], 1, &m_DFC1[i][0][0], 1);
914 
915  // Computing total first order derivatives
916  Vmath::Vadd(nSolutionPts, &m_DFC1[i][0][0], 1,
917  &m_DU1[i][0][0], 1, &m_D1[i][0][0], 1);
918 
919  Vmath::Vcopy(nSolutionPts, &m_D1[i][0][0], 1,
920  &m_tmp1[i][0][0], 1);
921  }
922 
923  // Computing interface numerical fluxes for m_D1
924  // in physical space
925  v_NumFluxforVector(fields, inarray, m_tmp1, m_IF2);
926 
927  for (i = 0; i < nConvectiveFields; ++i)
928  {
929  // Computing the physical second-order discountinuous
930  // derivative
931  for (n = 0; n < nElements; n++)
932  {
933  phys_offset = fields[0]->GetPhys_Offset(n);
934 
935  fields[i]->GetExp(n)->PhysDeriv(0,
936  auxArray1 = m_D1[i][0] + phys_offset,
937  auxArray2 = m_DD1[i][0] + phys_offset);
938  }
939 
940  // Computing the standard second-order correction
941  // derivative
942  v_DerCFlux_1D(nConvectiveFields, fields, m_D1[i][0],
943  m_IF2[i], m_DFC2[i][0]);
944 
945  // Back to the physical space using global operations
946  Vmath::Vdiv(nSolutionPts, &m_DFC2[i][0][0], 1,
947  &m_jac[0], 1, &m_DFC2[i][0][0], 1);
948  Vmath::Vdiv(nSolutionPts, &m_DFC2[i][0][0], 1,
949  &m_jac[0], 1, &m_DFC2[i][0][0], 1);
950 
951  // Adding the total divergence to outarray (RHS)
952  Vmath::Vadd(nSolutionPts, &m_DFC2[i][0][0], 1,
953  &m_DD1[i][0][0], 1, &outarray[i][0], 1);
954 
955  // Primitive Dealiasing 1D
956  if (!(Basis[0]->Collocation()))
957  {
958  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
959  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
960  }
961  }
962  break;
963  }
964  // 2D problems
965  case 2:
966  {
967  for(i = 0; i < nConvectiveFields; ++i)
968  {
969  for (j = 0; j < nDim; ++j)
970  {
971  // Temporary vectors
972  Array<OneD, NekDouble> u1_hat(nSolutionPts, 0.0);
973  Array<OneD, NekDouble> u2_hat(nSolutionPts, 0.0);
974 
975  if (j == 0)
976  {
977  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
978  &m_gmat[0][0], 1, &u1_hat[0], 1);
979 
980  Vmath::Vmul(nSolutionPts, &u1_hat[0], 1,
981  &m_jac[0], 1, &u1_hat[0], 1);
982 
983  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
984  &m_gmat[1][0], 1, &u2_hat[0], 1);
985 
986  Vmath::Vmul(nSolutionPts, &u2_hat[0], 1,
987  &m_jac[0], 1, &u2_hat[0], 1);
988  }
989  else if (j == 1)
990  {
991  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
992  &m_gmat[2][0], 1, &u1_hat[0], 1);
993 
994  Vmath::Vmul(nSolutionPts, &u1_hat[0], 1,
995  &m_jac[0], 1, &u1_hat[0], 1);
996 
997  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
998  &m_gmat[3][0], 1, &u2_hat[0], 1);
999 
1000  Vmath::Vmul(nSolutionPts, &u2_hat[0], 1,
1001  &m_jac[0], 1, &u2_hat[0], 1);
1002  }
1003 
1004  for (n = 0; n < nElements; n++)
1005  {
1006  phys_offset = fields[0]->GetPhys_Offset(n);
1007 
1008  fields[i]->GetExp(n)->StdPhysDeriv(0,
1009  auxArray1 = u1_hat + phys_offset,
1010  auxArray2 = m_tmp1[i][j] + phys_offset);
1011 
1012  fields[i]->GetExp(n)->StdPhysDeriv(1,
1013  auxArray1 = u2_hat + phys_offset,
1014  auxArray2 = m_tmp2[i][j] + phys_offset);
1015  }
1016 
1017  Vmath::Vadd(nSolutionPts, &m_tmp1[i][j][0], 1,
1018  &m_tmp2[i][j][0], 1,
1019  &m_DU1[i][j][0], 1);
1020 
1021  // Divide by the metric jacobian
1022  Vmath::Vdiv(nSolutionPts, &m_DU1[i][j][0], 1,
1023  &m_jac[0], 1, &m_DU1[i][j][0], 1);
1024 
1025  // Computing the standard first-order correction
1026  // derivatives
1027  v_DerCFlux_2D(nConvectiveFields, j, fields,
1028  inarray[i], m_IF1[i][j],
1029  m_DFC1[i][j]);
1030  }
1031 
1032  // Multiplying derivatives by B matrix to get auxiliary
1033  // variables
1034  for (j = 0; j < nSolutionPts; ++j)
1035  {
1036  // std(q1)
1037  m_BD1[i][0][j] =
1038  (m_gmat[0][j]*m_gmat[0][j] +
1039  m_gmat[2][j]*m_gmat[2][j]) *
1040  m_DFC1[i][0][j] +
1041  (m_gmat[1][j]*m_gmat[0][j] +
1042  m_gmat[3][j]*m_gmat[2][j]) *
1043  m_DFC1[i][1][j];
1044 
1045  // std(q2)
1046  m_BD1[i][1][j] =
1047  (m_gmat[1][j]*m_gmat[0][j] +
1048  m_gmat[3][j]*m_gmat[2][j]) *
1049  m_DFC1[i][0][j] +
1050  (m_gmat[1][j]*m_gmat[1][j] +
1051  m_gmat[3][j]*m_gmat[3][j]) *
1052  m_DFC1[i][1][j];
1053  }
1054 
1055  // Multiplying derivatives by A^(-1) to get back
1056  // into the physical space
1057  for (j = 0; j < nSolutionPts; j++)
1058  {
1059  // q1 = A11^(-1)*std(q1) + A12^(-1)*std(q2)
1060  m_DFC1[i][0][j] =
1061  m_gmat[3][j] * m_BD1[i][0][j] -
1062  m_gmat[2][j] * m_BD1[i][1][j];
1063 
1064  // q2 = A21^(-1)*std(q1) + A22^(-1)*std(q2)
1065  m_DFC1[i][1][j] =
1066  m_gmat[0][j] * m_BD1[i][1][j] -
1067  m_gmat[1][j] * m_BD1[i][0][j];
1068  }
1069 
1070  // Computing the physical first-order derivatives
1071  for (j = 0; j < nDim; ++j)
1072  {
1073  Vmath::Vadd(nSolutionPts, &m_DU1[i][j][0], 1,
1074  &m_DFC1[i][j][0], 1,
1075  &m_D1[i][j][0], 1);
1076  }
1077  }
1078 
1079  // Computing interface numerical fluxes for m_D1
1080  // in physical space
1081  v_NumFluxforVector(fields, inarray, m_D1, m_IF2);
1082 
1083  // Computing the standard second-order discontinuous
1084  // derivatives
1085  for (i = 0; i < nConvectiveFields; ++i)
1086  {
1087  Array<OneD, NekDouble> f_hat(nSolutionPts, 0.0);
1088  Array<OneD, NekDouble> g_hat(nSolutionPts, 0.0);
1089 
1090  for (j = 0; j < nSolutionPts; j++)
1091  {
1092  f_hat[j] = (m_D1[i][0][j] * m_gmat[0][j] +
1093  m_D1[i][1][j] * m_gmat[2][j])*m_jac[j];
1094 
1095  g_hat[j] = (m_D1[i][0][j] * m_gmat[1][j] +
1096  m_D1[i][1][j] * m_gmat[3][j])*m_jac[j];
1097  }
1098 
1099  for (n = 0; n < nElements; n++)
1100  {
1101  phys_offset = fields[0]->GetPhys_Offset(n);
1102 
1103  fields[0]->GetExp(n)->StdPhysDeriv(0,
1104  auxArray1 = f_hat + phys_offset,
1105  auxArray2 = m_DD1[i][0] + phys_offset);
1106 
1107  fields[0]->GetExp(n)->StdPhysDeriv(1,
1108  auxArray1 = g_hat + phys_offset,
1109  auxArray2 = m_DD1[i][1] + phys_offset);
1110  }
1111 
1112  // Divergence of the standard discontinuous flux
1113  Vmath::Vadd(nSolutionPts,
1114  &m_DD1[i][0][0], 1,
1115  &m_DD1[i][1][0], 1,
1116  &m_divFD[i][0], 1);
1117 
1118  // Divergence of the standard correction flux
1119  if (Basis[0]->GetPointsType() ==
1120  LibUtilities::eGaussGaussLegendre &&
1121  Basis[1]->GetPointsType() ==
1122  LibUtilities::eGaussGaussLegendre)
1123  {
1124  v_DivCFlux_2D_Gauss(nConvectiveFields, fields,
1125  f_hat, g_hat,
1126  m_IF2[i], m_divFC[i]);
1127  }
1128  else
1129  {
1130  v_DivCFlux_2D(nConvectiveFields, fields,
1131  m_D1[i][0], m_D1[i][1],
1132  m_IF2[i], m_divFC[i]);
1133  }
1134 
1135  // Divergence of the standard final flux
1136  Vmath::Vadd(nSolutionPts, &m_divFD[i][0], 1,
1137  &m_divFC[i][0], 1, &outarray[i][0], 1);
1138 
1139  // Dividing by the jacobian to get back into
1140  // physical space
1141  Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1,
1142  &m_jac[0], 1, &outarray[i][0], 1);
1143 
1144  // Primitive Dealiasing 2D
1145  if (!(Basis[0]->Collocation()))
1146  {
1147  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1148  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1149  }
1150  }//Close loop on nConvectiveFields
1151  break;
1152  }
1153  // 3D problems
1154  case 3:
1155  {
1156  ASSERTL0(false, "3D FRDG case not implemented yet");
1157  break;
1158  }
1159  }
1160  }
Nektar::SolverUtils::DiffusionLFR::v_NumFluxforVector
virtual void v_NumFluxforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
Build the numerical flux for the 2nd order derivatives.
Definition: DiffusionLFR.cpp:1307
Nektar::SolverUtils::DiffusionLFR::m_tmp2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Definition: DiffusionLFR.h:89
Nektar::SolverUtils::DiffusionLFR::v_DerCFlux_1D
virtual void v_DerCFlux_1D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, const NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux for 1D problems.
Definition: DiffusionLFR.cpp:1474
Nektar::SolverUtils::DiffusionLFR::m_IF1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
Definition: DiffusionLFR.h:77
Nektar::SolverUtils::DiffusionLFR::m_jac
Array< OneD, NekDouble > m_jac
Definition: DiffusionLFR.h:54
Nektar::SolverUtils::DiffusionLFR::v_NumFluxforScalar
virtual void v_NumFluxforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
Builds the numerical flux for the 1st order derivatives.
Definition: DiffusionLFR.cpp:1166
Nektar::SolverUtils::DiffusionLFR::m_DFC1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Definition: DiffusionLFR.h:79
Nektar::SolverUtils::DiffusionLFR::m_DD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
Definition: DiffusionLFR.h:82
Nektar::SolverUtils::DiffusionLFR::v_DivCFlux_2D_Gauss
virtual void v_DivCFlux_2D_Gauss(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre".
Definition: DiffusionLFR.cpp:2008
Nektar::SolverUtils::DiffusionLFR::m_divFC
Array< OneD, Array< OneD, NekDouble > > m_divFC
Definition: DiffusionLFR.h:86
Nektar::SolverUtils::DiffusionLFR::v_DivCFlux_2D
virtual void v_DivCFlux_2D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems.
Definition: DiffusionLFR.cpp:1813
Nektar::SolverUtils::DiffusionLFR::v_DerCFlux_2D
virtual void v_DerCFlux_2D(const int nConvectiveFields, const int direction, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux wrt a given coordinate for 2D problems.
Definition: DiffusionLFR.cpp:1596
Nektar::SolverUtils::DiffusionLFR::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: DiffusionLFR.h:55
Nektar::SolverUtils::DiffusionLFR::m_DU1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Definition: DiffusionLFR.h:78
Nektar::SolverUtils::DiffusionLFR::m_DFC2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Definition: DiffusionLFR.h:84
Nektar::SolverUtils::DiffusionLFR::m_tmp1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Definition: DiffusionLFR.h:88
Nektar::SolverUtils::DiffusionLFR::m_D1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
Definition: DiffusionLFR.h:81
Nektar::SolverUtils::DiffusionLFR::m_IF2
Array< OneD, Array< OneD, NekDouble > > m_IF2
Definition: DiffusionLFR.h:83
Nektar::SolverUtils::DiffusionLFR::m_divFD
Array< OneD, Array< OneD, NekDouble > > m_divFD
Definition: DiffusionLFR.h:85
Nektar::SolverUtils::DiffusionLFR::m_BD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
Definition: DiffusionLFR.h:80
Nektar::LibUtilities::eGaussGaussLegendre
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:48
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:192
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:257

References ASSERTL0, Nektar::LibUtilities::eGaussGaussLegendre, m_BD1, m_D1, m_DD1, m_DFC1, m_DFC2, m_divFC, m_divFD, m_DU1, m_gmat, m_IF1, m_IF2, m_jac, m_tmp1, m_tmp2, v_DerCFlux_1D(), v_DerCFlux_2D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), v_NumFluxforScalar(), v_NumFluxforVector(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), and Vmath::Vmul().

◆ v_DivCFlux_2D()

void Nektar::SolverUtils::DiffusionLFR::v_DivCFlux_2D ( const int  nConvectiveFields,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, const NekDouble > &  fluxX1,
const Array< OneD, const NekDouble > &  fluxX2,
const Array< OneD, const NekDouble > &  numericalFlux,
Array< OneD, NekDouble > &  divCFlux 
)
protectedvirtual

Compute the divergence of the corrective flux for 2D problems.

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxX1X1 - volumetric flux in physical space.
fluxX2X2 - volumetric flux in physical space.
numericalFluxInterface flux in physical space.
divCFluxDivergence of the corrective flux for 2D problems.
Todo:
: Switch on shapes eventually here.

Definition at line 1813 of file DiffusionLFR.cpp.

1820  {
1821  boost::ignore_unused(nConvectiveFields);
1822 
1823  int n, e, i, j, cnt;
1824 
1825  int nElements = fields[0]->GetExpSize();
1826 
1827  int nLocalSolutionPts;
1828  int nEdgePts;
1829  int trace_offset;
1830  int phys_offset;
1831  int nquad0;
1832  int nquad1;
1833 
1834  Array<OneD, NekDouble> auxArray1, auxArray2;
1835  Array<OneD, LibUtilities::BasisSharedPtr> base;
1836 
1837  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1838  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1839 
1840  // Loop on the elements
1841  for(n = 0; n < nElements; ++n)
1842  {
1843  // Offset of the element on the global vector
1844  phys_offset = fields[0]->GetPhys_Offset(n);
1845  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1846 
1847  base = fields[0]->GetExp(n)->GetBase();
1848  nquad0 = base[0]->GetNumPoints();
1849  nquad1 = base[1]->GetNumPoints();
1850 
1851  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1852  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1853  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1854  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1855 
1856  // Loop on the edges
1857  for(e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1858  {
1859  // Number of edge points of edge e
1860  nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1861 
1862  Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1863  Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1864  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
1865  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
1866  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1867 
1868  // Offset of the trace space correspondent to edge e
1869  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1870  elmtToTrace[n][e]->GetElmtId());
1871 
1872  // Get the normals of edge e
1873  const Array<OneD, const Array<OneD, NekDouble> > &normals =
1874  fields[0]->GetExp(n)->GetTraceNormal(e);
1875 
1876  // Extract the edge values of flux-x on edge e and order
1877  // them accordingly to the order of the trace space
1878  fields[0]->GetExp(n)->GetTracePhysVals(
1879  e, elmtToTrace[n][e],
1880  fluxX1 + phys_offset,
1881  auxArray1 = tmparrayX1);
1882 
1883  // Extract the edge values of flux-y on edge e and order
1884  // them accordingly to the order of the trace space
1885  fields[0]->GetExp(n)->GetTracePhysVals(
1886  e, elmtToTrace[n][e],
1887  fluxX2 + phys_offset,
1888  auxArray1 = tmparrayX2);
1889 
1890  // Multiply the edge components of the flux by the normal
1891  for (i = 0; i < nEdgePts; ++i)
1892  {
1893  fluxN[i] =
1894  tmparrayX1[i]*m_traceNormals[0][trace_offset+i] +
1895  tmparrayX2[i]*m_traceNormals[1][trace_offset+i];
1896  }
1897 
1898  // Subtract to the Riemann flux the discontinuous flux
1899  Vmath::Vsub(nEdgePts,
1900  &numericalFlux[trace_offset], 1,
1901  &fluxN[0], 1, &fluxJumps[0], 1);
1902 
1903  // Check the ordering of the jump vectors
1904  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1905  StdRegions::eBackwards)
1906  {
1907  Vmath::Reverse(nEdgePts,
1908  auxArray1 = fluxJumps, 1,
1909  auxArray2 = fluxJumps, 1);
1910  }
1911 
1912  for (i = 0; i < nEdgePts; ++i)
1913  {
1914  if (m_traceNormals[0][trace_offset+i] != normals[0][i]
1915  || m_traceNormals[1][trace_offset+i] != normals[1][i])
1916  {
1917  fluxJumps[i] = -fluxJumps[i];
1918  }
1919  }
1920 
1921  // Multiply jumps by derivatives of the correction functions
1922  switch (e)
1923  {
1924  case 0:
1925  for (i = 0; i < nquad0; ++i)
1926  {
1927  // Multiply fluxJumps by Q factors
1928  fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1929 
1930  for (j = 0; j < nquad1; ++j)
1931  {
1932  cnt = i + j*nquad0;
1933  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1934  }
1935  }
1936  break;
1937  case 1:
1938  for (i = 0; i < nquad1; ++i)
1939  {
1940  // Multiply fluxJumps by Q factors
1941  fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1942 
1943  for (j = 0; j < nquad0; ++j)
1944  {
1945  cnt = (nquad0)*i + j;
1946  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1947  }
1948  }
1949  break;
1950  case 2:
1951  for (i = 0; i < nquad0; ++i)
1952  {
1953  // Multiply fluxJumps by Q factors
1954  fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1955 
1956  for (j = 0; j < nquad1; ++j)
1957  {
1958  cnt = j*nquad0 + i;
1959  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1960  }
1961  }
1962  break;
1963  case 3:
1964  for (i = 0; i < nquad1; ++i)
1965  {
1966  // Multiply fluxJumps by Q factors
1967  fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1968  for (j = 0; j < nquad0; ++j)
1969  {
1970  cnt = j + i*nquad0;
1971  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1972  }
1973  }
1974  break;
1975 
1976  default:
1977  ASSERTL0(false,"edge value (< 3) is out of range");
1978  break;
1979  }
1980  }
1981 
1982  // Sum each edge contribution
1983  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1984  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1985 
1986  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1987  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1988 
1989  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1990  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1991  }
1992  }
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:74
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFR.h:57
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFR.h:60
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFR.h:59
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFR.h:58

References ASSERTL0, Nektar::StdRegions::eBackwards, m_dGL_xi1, m_dGL_xi2, m_dGR_xi1, m_dGR_xi2, m_Q2D_e0, m_Q2D_e1, m_Q2D_e2, m_Q2D_e3, m_traceNormals, Vmath::Reverse(), Vmath::Vadd(), and Vmath::Vsub().

Referenced by v_Diffuse().

◆ v_DivCFlux_2D_Gauss()

void Nektar::SolverUtils::DiffusionLFR::v_DivCFlux_2D_Gauss ( const int  nConvectiveFields,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, const NekDouble > &  fluxX1,
const Array< OneD, const NekDouble > &  fluxX2,
const Array< OneD, const NekDouble > &  numericalFlux,
Array< OneD, NekDouble > &  divCFlux 
)
protectedvirtual

Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre".

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxX1X1-volumetric flux in physical space.
fluxX2X2-volumetric flux in physical space.
numericalFluxInterface flux in physical space.
divCFluxDivergence of the corrective flux.
Todo:
: Switch on shapes eventually here.

Definition at line 2008 of file DiffusionLFR.cpp.

2015  {
2016  boost::ignore_unused(nConvectiveFields);
2017 
2018  int n, e, i, j, cnt;
2019 
2020  int nElements = fields[0]->GetExpSize();
2021  int nLocalSolutionPts;
2022  int nEdgePts;
2023  int trace_offset;
2024  int phys_offset;
2025  int nquad0;
2026  int nquad1;
2027 
2028  Array<OneD, NekDouble> auxArray1, auxArray2;
2029  Array<OneD, LibUtilities::BasisSharedPtr> base;
2030 
2031  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
2032  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
2033 
2034  // Loop on the elements
2035  for(n = 0; n < nElements; ++n)
2036  {
2037  // Offset of the element on the global vector
2038  phys_offset = fields[0]->GetPhys_Offset(n);
2039  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
2040 
2041  base = fields[0]->GetExp(n)->GetBase();
2042  nquad0 = base[0]->GetNumPoints();
2043  nquad1 = base[1]->GetNumPoints();
2044 
2045  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
2046  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
2047  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
2048  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
2049 
2050  // Loop on the edges
2051  for(e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
2052  {
2053  // Number of edge points of edge e
2054  nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
2055 
2056  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
2057  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
2058  Array<OneD, NekDouble> fluxN_R (nEdgePts, 0.0);
2059  Array<OneD, NekDouble> fluxN_D (nEdgePts, 0.0);
2060  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
2061 
2062  // Offset of the trace space correspondent to edge e
2063  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
2064  elmtToTrace[n][e]->GetElmtId());
2065 
2066  // Get the normals of edge e
2067  const Array<OneD, const Array<OneD, NekDouble> > &normals =
2068  fields[0]->GetExp(n)->GetTraceNormal(e);
2069 
2070  // Extract the trasformed normal flux at each edge
2071  switch (e)
2072  {
2073  case 0:
2074  // Extract the edge values of transformed flux-y on
2075  // edge e and order them accordingly to the order of
2076  // the trace space
2077  fields[0]->GetExp(n)->GetTracePhysVals(
2078  e, elmtToTrace[n][e],
2079  fluxX2 + phys_offset,
2080  auxArray1 = fluxN_D);
2081 
2082  Vmath::Neg (nEdgePts, fluxN_D, 1);
2083 
2084  // Extract the physical Rieamann flux at each edge
2085  Vmath::Vcopy(nEdgePts,
2086  &numericalFlux[trace_offset], 1,
2087  &fluxN[0], 1);
2088 
2089  // Check the ordering of vectors
2090  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2091  StdRegions::eBackwards)
2092  {
2093  Vmath::Reverse(nEdgePts,
2094  auxArray1 = fluxN, 1,
2095  auxArray2 = fluxN, 1);
2096 
2097  Vmath::Reverse(nEdgePts,
2098  auxArray1 = fluxN_D, 1,
2099  auxArray2 = fluxN_D, 1);
2100  }
2101 
2102  // Transform Riemann Fluxes in the standard element
2103  for (i = 0; i < nquad0; ++i)
2104  {
2105  // Multiply Riemann Flux by Q factors
2106  fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
2107  }
2108 
2109  for (i = 0; i < nEdgePts; ++i)
2110  {
2111  if (m_traceNormals[0][trace_offset+i] !=
2112  normals[0][i] ||
2113  m_traceNormals[1][trace_offset+i] !=
2114  normals[1][i])
2115  {
2116  fluxN_R[i] = -fluxN_R[i];
2117  }
2118  }
2119 
2120  // Subtract to the Riemann flux the discontinuous
2121  // flux
2122  Vmath::Vsub(nEdgePts,
2123  &fluxN_R[0], 1,
2124  &fluxN_D[0], 1, &fluxJumps[0], 1);
2125 
2126  // Multiplicate the flux jumps for the correction
2127  // function
2128  for (i = 0; i < nquad0; ++i)
2129  {
2130  for (j = 0; j < nquad1; ++j)
2131  {
2132  cnt = i + j*nquad0;
2133  divCFluxE0[cnt] =
2134  -fluxJumps[i] * m_dGL_xi2[n][j];
2135  }
2136  }
2137  break;
2138  case 1:
2139  // Extract the edge values of transformed flux-x on
2140  // edge e and order them accordingly to the order of
2141  // the trace space
2142  fields[0]->GetExp(n)->GetTracePhysVals(
2143  e, elmtToTrace[n][e],
2144  fluxX1 + phys_offset,
2145  auxArray1 = fluxN_D);
2146 
2147  // Extract the physical Rieamann flux at each edge
2148  Vmath::Vcopy(nEdgePts,
2149  &numericalFlux[trace_offset], 1,
2150  &fluxN[0], 1);
2151 
2152  // Check the ordering of vectors
2153  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2154  StdRegions::eBackwards)
2155  {
2156  Vmath::Reverse(nEdgePts,
2157  auxArray1 = fluxN, 1,
2158  auxArray2 = fluxN, 1);
2159 
2160  Vmath::Reverse(nEdgePts,
2161  auxArray1 = fluxN_D, 1,
2162  auxArray2 = fluxN_D, 1);
2163  }
2164 
2165  // Transform Riemann Fluxes in the standard element
2166  for (i = 0; i < nquad1; ++i)
2167  {
2168  // Multiply Riemann Flux by Q factors
2169  fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
2170  }
2171 
2172  for (i = 0; i < nEdgePts; ++i)
2173  {
2174  if (m_traceNormals[0][trace_offset+i] !=
2175  normals[0][i] ||
2176  m_traceNormals[1][trace_offset+i] !=
2177  normals[1][i])
2178  {
2179  fluxN_R[i] = -fluxN_R[i];
2180  }
2181  }
2182 
2183  // Subtract to the Riemann flux the discontinuous
2184  // flux
2185  Vmath::Vsub(nEdgePts,
2186  &fluxN_R[0], 1,
2187  &fluxN_D[0], 1, &fluxJumps[0], 1);
2188 
2189  // Multiplicate the flux jumps for the correction
2190  // function
2191  for (i = 0; i < nquad1; ++i)
2192  {
2193  for (j = 0; j < nquad0; ++j)
2194  {
2195  cnt = (nquad0)*i + j;
2196  divCFluxE1[cnt] =
2197  fluxJumps[i] * m_dGR_xi1[n][j];
2198  }
2199  }
2200  break;
2201  case 2:
2202 
2203  // Extract the edge values of transformed flux-y on
2204  // edge e and order them accordingly to the order of
2205  // the trace space
2206 
2207  fields[0]->GetExp(n)->GetTracePhysVals(
2208  e, elmtToTrace[n][e],
2209  fluxX2 + phys_offset,
2210  auxArray1 = fluxN_D);
2211 
2212  // Extract the physical Rieamann flux at each edge
2213  Vmath::Vcopy(nEdgePts,
2214  &numericalFlux[trace_offset], 1,
2215  &fluxN[0], 1);
2216 
2217  // Check the ordering of vectors
2218  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2219  StdRegions::eBackwards)
2220  {
2221  Vmath::Reverse(nEdgePts,
2222  auxArray1 = fluxN, 1,
2223  auxArray2 = fluxN, 1);
2224 
2225  Vmath::Reverse(nEdgePts,
2226  auxArray1 = fluxN_D, 1,
2227  auxArray2 = fluxN_D, 1);
2228  }
2229 
2230  // Transform Riemann Fluxes in the standard element
2231  for (i = 0; i < nquad0; ++i)
2232  {
2233  // Multiply Riemann Flux by Q factors
2234  fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
2235  }
2236 
2237  for (i = 0; i < nEdgePts; ++i)
2238  {
2239  if (m_traceNormals[0][trace_offset+i] !=
2240  normals[0][i] ||
2241  m_traceNormals[1][trace_offset+i] !=
2242  normals[1][i])
2243  {
2244  fluxN_R[i] = -fluxN_R[i];
2245  }
2246  }
2247 
2248  // Subtract to the Riemann flux the discontinuous
2249  // flux
2250 
2251  Vmath::Vsub(nEdgePts,
2252  &fluxN_R[0], 1,
2253  &fluxN_D[0], 1, &fluxJumps[0], 1);
2254 
2255  // Multiplicate the flux jumps for the correction
2256  // function
2257  for (i = 0; i < nquad0; ++i)
2258  {
2259  for (j = 0; j < nquad1; ++j)
2260  {
2261  cnt = j*nquad0 + i;
2262  divCFluxE2[cnt] =
2263  fluxJumps[i] * m_dGR_xi2[n][j];
2264  }
2265  }
2266  break;
2267  case 3:
2268  // Extract the edge values of transformed flux-x on
2269  // edge e and order them accordingly to the order of
2270  // the trace space
2271 
2272  fields[0]->GetExp(n)->GetTracePhysVals(
2273  e, elmtToTrace[n][e],
2274  fluxX1 + phys_offset,
2275  auxArray1 = fluxN_D);
2276  Vmath::Neg (nEdgePts, fluxN_D, 1);
2277 
2278  // Extract the physical Rieamann flux at each edge
2279  Vmath::Vcopy(nEdgePts,
2280  &numericalFlux[trace_offset], 1,
2281  &fluxN[0], 1);
2282 
2283  // Check the ordering of vectors
2284  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2285  StdRegions::eBackwards)
2286  {
2287  Vmath::Reverse(nEdgePts,
2288  auxArray1 = fluxN, 1,
2289  auxArray2 = fluxN, 1);
2290 
2291  Vmath::Reverse(nEdgePts,
2292  auxArray1 = fluxN_D, 1,
2293  auxArray2 = fluxN_D, 1);
2294  }
2295 
2296  // Transform Riemann Fluxes in the standard element
2297  for (i = 0; i < nquad1; ++i)
2298  {
2299  // Multiply Riemann Flux by Q factors
2300  fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
2301  }
2302 
2303  for (i = 0; i < nEdgePts; ++i)
2304  {
2305  if (m_traceNormals[0][trace_offset+i] !=
2306  normals[0][i] ||
2307  m_traceNormals[1][trace_offset+i] !=
2308  normals[1][i])
2309  {
2310  fluxN_R[i] = -fluxN_R[i];
2311  }
2312  }
2313 
2314  // Subtract to the Riemann flux the discontinuous
2315  // flux
2316 
2317  Vmath::Vsub(nEdgePts,
2318  &fluxN_R[0], 1,
2319  &fluxN_D[0], 1, &fluxJumps[0], 1);
2320 
2321  // Multiplicate the flux jumps for the correction
2322  // function
2323  for (i = 0; i < nquad1; ++i)
2324  {
2325  for (j = 0; j < nquad0; ++j)
2326  {
2327  cnt = j + i*nquad0;
2328  divCFluxE3[cnt] =
2329  -fluxJumps[i] * m_dGL_xi1[n][j];
2330  }
2331  }
2332  break;
2333  default:
2334  ASSERTL0(false,"edge value (< 3) is out of range");
2335  break;
2336  }
2337  }
2338 
2339 
2340  // Sum each edge contribution
2341  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
2342  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
2343 
2344  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2345  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2346 
2347  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2348  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2349  }
2350  }
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:461

References ASSERTL0, Nektar::StdRegions::eBackwards, m_dGL_xi1, m_dGL_xi2, m_dGR_xi1, m_dGR_xi2, m_Q2D_e0, m_Q2D_e1, m_Q2D_e2, m_Q2D_e3, m_traceNormals, Vmath::Neg(), Vmath::Reverse(), Vmath::Vadd(), Vmath::Vcopy(), and Vmath::Vsub().

Referenced by v_Diffuse().

◆ v_InitObject()

void Nektar::SolverUtils::DiffusionLFR::v_InitObject ( LibUtilities::SessionReaderSharedPtr  pSession,
Array< OneD, MultiRegions::ExpListSharedPtr >  pFields 
)
protectedvirtual

Initiliase DiffusionLFR objects and store them before starting the time-stepping.

This routine calls the virtual functions v_SetupMetrics and v_SetupCFunctions to initialise the objects needed by DiffusionLFR.

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.

Reimplemented from Nektar::SolverUtils::Diffusion.

Definition at line 85 of file DiffusionLFR.cpp.

88  {
89  m_session = pSession;
90  v_SetupMetrics(pSession, pFields);
91  v_SetupCFunctions(pSession, pFields);
92 
93  // Initialising arrays
94  int i, j;
95  int nConvectiveFields = pFields.size();
96  int nDim = pFields[0]->GetCoordim(0);
97  int nSolutionPts = pFields[0]->GetTotPoints();
98  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
99 
100  m_IF1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
101  nConvectiveFields);
102  m_DU1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
103  nConvectiveFields);
104  m_DFC1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
105  nConvectiveFields);
106  m_BD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
107  nConvectiveFields);
108  m_D1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
109  nConvectiveFields);
110  m_DD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
111  nConvectiveFields);
112  m_IF2 = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
113  m_DFC2 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
114  nConvectiveFields);
115  m_divFD = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
116  m_divFC = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
117 
118  m_tmp1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
119  nConvectiveFields);
120  m_tmp2 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
121  nConvectiveFields);
122 
123 
124 
125  for (i = 0; i < nConvectiveFields; ++i)
126  {
127  m_IF1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
128  m_DU1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
129  m_DFC1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
130  m_BD1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
131  m_D1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
132  m_DD1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
133  m_IF2[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
134  m_DFC2[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
135  m_divFD[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
136  m_divFC[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
137 
138  m_tmp1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
139  m_tmp2[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
140 
141  for (j = 0; j < nDim; ++j)
142  {
143  m_IF1[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
144  m_DU1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
145  m_DFC1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
146  m_BD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
147  m_D1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
148  m_DD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
149  m_DFC2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
150 
151  m_tmp1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
152  m_tmp2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
153  }
154  }
155 
156  }
Nektar::SolverUtils::DiffusionLFR::v_SetupMetrics
virtual void v_SetupMetrics(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform...
Definition: DiffusionLFR.cpp:174
Nektar::SolverUtils::DiffusionLFR::m_session
LibUtilities::SessionReaderSharedPtr m_session
Definition: DiffusionLFR.h:75
Nektar::SolverUtils::DiffusionLFR::v_SetupCFunctions
virtual void v_SetupCFunctions(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–72...
Definition: DiffusionLFR.cpp:326

References m_BD1, m_D1, m_DD1, m_DFC1, m_DFC2, m_divFC, m_divFD, m_DU1, m_IF1, m_IF2, m_session, m_tmp1, m_tmp2, v_SetupCFunctions(), and v_SetupMetrics().

◆ v_NumFluxforScalar()

void Nektar::SolverUtils::DiffusionLFR::v_NumFluxforScalar ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, Array< OneD, NekDouble > > &  ufield,
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  uflux 
)
protectedvirtual

Builds the numerical flux for the 1st order derivatives.

Definition at line 1166 of file DiffusionLFR.cpp.

1170  {
1171  int i, j;
1172  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1173  int nvariables = fields.size();
1174  int nDim = fields[0]->GetCoordim(0);
1175 
1176  Array<OneD, NekDouble > Fwd (nTracePts);
1177  Array<OneD, NekDouble > Bwd (nTracePts);
1178  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1179  Array<OneD, NekDouble > fluxtemp(nTracePts, 0.0);
1180 
1181  // Get the normal velocity Vn
1182  for (i = 0; i < nDim; ++i)
1183  {
1184  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
1185  Vn, 1, Vn, 1);
1186  }
1187 
1188  // Get the sign of (v \cdot n), v = an arbitrary vector
1189  // Evaluate upwind flux:
1190  // uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
1191  for (j = 0; j < nDim; ++j)
1192  {
1193  for (i = 0; i < nvariables ; ++i)
1194  {
1195  // Compute Fwd and Bwd value of ufield of i direction
1196  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1197 
1198  // if Vn >= 0, flux = uFwd, i.e.,
1199  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
1200  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
1201 
1202  // else if Vn < 0, flux = uBwd, i.e.,
1203  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
1204  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
1205 
1206  fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, fluxtemp);
1207 
1208  // Imposing weak boundary condition with flux
1209  // if Vn >= 0, uflux = uBwd at Neumann, i.e.,
1210  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
1211  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
1212 
1213  // if Vn >= 0, uflux = uFwd at Neumann, i.e.,
1214  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
1215  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
1216 
1217  if(fields[0]->GetBndCondExpansions().size())
1218  {
1219  v_WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
1220  }
1221 
1222  // if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
1223  // i.e,
1224  // edge::eForward, uFwd \‍(\tan_{\xi}^Fwd \cdot \vec{n})
1225  // edge::eBackward, uFwd \‍(\tan_{\xi}^Bwd \cdot \vec{n})
1226 
1227  // else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}),
1228  // i.e,
1229  // edge::eForward, uBwd \‍(\tan_{\xi}^Fwd \cdot \vec{n})
1230  // edge::eBackward, uBwd \‍(\tan_{\xi}^Bwd \cdot \vec{n})
1231 
1232  Vmath::Vcopy(nTracePts, fluxtemp, 1, uflux[i][j], 1);
1233  }
1234  }
1235  }
Nektar::SolverUtils::DiffusionLFR::v_WeakPenaltyforScalar
virtual void v_WeakPenaltyforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const Array< OneD, const NekDouble > &ufield, Array< OneD, NekDouble > &penaltyflux)
Imposes appropriate bcs for the 1st order derivatives.
Definition: DiffusionLFR.cpp:1241
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:565

References m_traceNormals, Vmath::Svtvp(), v_WeakPenaltyforScalar(), and Vmath::Vcopy().

Referenced by v_Diffuse().

◆ v_NumFluxforVector()

void Nektar::SolverUtils::DiffusionLFR::v_NumFluxforVector ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, Array< OneD, NekDouble > > &  ufield,
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  qfield,
Array< OneD, Array< OneD, NekDouble > > &  qflux 
)
protectedvirtual

Build the numerical flux for the 2nd order derivatives.

Definition at line 1307 of file DiffusionLFR.cpp.

1312  {
1313  boost::ignore_unused(ufield);
1314 
1315  int i, j;
1316  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1317  int nvariables = fields.size();
1318  int nDim = fields[0]->GetCoordim(0);
1319 
1320  NekDouble C11 = 0.0;
1321  Array<OneD, NekDouble > Fwd(nTracePts);
1322  Array<OneD, NekDouble > Bwd(nTracePts);
1323  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1324 
1325  Array<OneD, NekDouble > qFwd (nTracePts);
1326  Array<OneD, NekDouble > qBwd (nTracePts);
1327  Array<OneD, NekDouble > qfluxtemp(nTracePts, 0.0);
1328  Array<OneD, NekDouble > uterm(nTracePts);
1329 
1330  // Get the normal velocity Vn
1331  for(i = 0; i < nDim; ++i)
1332  {
1333  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
1334  Vn, 1, Vn, 1);
1335  }
1336 
1337  // Evaulate upwind flux:
1338  // qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
1339  for (i = 0; i < nvariables; ++i)
1340  {
1341  qflux[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
1342  for (j = 0; j < nDim; ++j)
1343  {
1344  // Compute Fwd and Bwd value of ufield of jth direction
1345  fields[i]->GetFwdBwdTracePhys(qfield[i][j], qFwd, qBwd);
1346 
1347  // if Vn >= 0, flux = uFwd, i.e.,
1348  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick
1349  // qflux = qBwd = q+
1350  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick
1351  // qflux = qBwd = q-
1352 
1353  // else if Vn < 0, flux = uBwd, i.e.,
1354  // edge::eForward, if V*n<0 <=> V*n_F<0, pick
1355  // qflux = qFwd = q-
1356  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
1357  // qflux = qFwd = q+
1358 
1359  fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1360 
1361  Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
1362  qfluxtemp, 1);
1363 
1364  /*
1365  // Generate Stability term = - C11 ( u- - u+ )
1366  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1367 
1368  Vmath::Vsub(nTracePts,
1369  Fwd, 1, Bwd, 1,
1370  uterm, 1);
1371 
1372  Vmath::Smul(nTracePts,
1373  -1.0 * C11, uterm, 1,
1374  uterm, 1);
1375 
1376  // Flux = {Fwd, Bwd} * (nx, ny, nz) + uterm * (nx, ny)
1377  Vmath::Vadd(nTracePts,
1378  uterm, 1,
1379  qfluxtemp, 1,
1380  qfluxtemp, 1);
1381  */
1382 
1383  // Imposing weak boundary condition with flux
1384  if (fields[0]->GetBndCondExpansions().size())
1385  {
1386  v_WeakPenaltyforVector(fields, i, j, qfield[i][j],
1387  qfluxtemp, C11);
1388  }
1389 
1390  // q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
1391  // n_xi = n_x * tan_xi_x + n_y * tan_xi_y + n_z * tan_xi_z
1392  // n_xi = n_x * tan_eta_x + n_y * tan_eta_y + n_z*tan_eta_z
1393  Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1,
1394  qflux[i], 1);
1395  }
1396  }
1397  }
Nektar::SolverUtils::DiffusionLFR::v_WeakPenaltyforVector
virtual void v_WeakPenaltyforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux, NekDouble C11)
Imposes appropriate bcs for the 2nd order derivatives.
Definition: DiffusionLFR.cpp:1405
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43

References m_traceNormals, Vmath::Svtvp(), v_WeakPenaltyforVector(), Vmath::Vadd(), and Vmath::Vmul().

Referenced by v_Diffuse().

◆ v_SetupCFunctions()

void Nektar::SolverUtils::DiffusionLFR::v_SetupCFunctions ( LibUtilities::SessionReaderSharedPtr  pSession,
Array< OneD, MultiRegions::ExpListSharedPtr >  pFields 
)
protectedvirtual

Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–72.

This routine calls 3 different bases: #eDG_DG_Left - #eDG_DG_Left which recovers a nodal DG scheme, #eDG_SD_Left - #eDG_SD_Left which recovers the SD scheme, #eDG_HU_Left - #eDG_HU_Left which recovers the Huynh scheme. The values of the derivatives of the correction function are then stored into global variables and reused into the virtual functions #v_DivCFlux_1D, v_DivCFlux_2D, #v_DivCFlux_3D to compute the the divergence of the correction flux for 1D, 2D or 3D problems.

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.

Definition at line 326 of file DiffusionLFR.cpp.

329  {
330  boost::ignore_unused(pSession);
331 
332  int i, n;
333  NekDouble c0 = 0.0;
334  NekDouble c1 = 0.0;
335  NekDouble c2 = 0.0;
336  int nquad0, nquad1, nquad2;
337  int nmodes0, nmodes1, nmodes2;
338  Array<OneD, LibUtilities::BasisSharedPtr> base;
339 
340  int nElements = pFields[0]->GetExpSize();
341  int nDim = pFields[0]->GetCoordim(0);
342 
343  switch (nDim)
344  {
345  case 1:
346  {
347  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
348  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
349 
350  for (n = 0; n < nElements; ++n)
351  {
352  base = pFields[0]->GetExp(n)->GetBase();
353  nquad0 = base[0]->GetNumPoints();
354  nmodes0 = base[0]->GetNumModes();
355  Array<OneD, const NekDouble> z0;
356  Array<OneD, const NekDouble> w0;
357 
358  base[0]->GetZW(z0, w0);
359 
360  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
361  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
362 
363  // Auxiliary vectors to build up the auxiliary
364  // derivatives of the Legendre polynomials
365  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
366  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
367  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
368 
369  // Degree of the correction functions
370  int p0 = nmodes0 - 1;
371 
372  // Function sign to compute left correction function
373  NekDouble sign0 = pow(-1.0, p0);
374 
375  // Factorial factor to build the scheme
376  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
377  / (pow(2.0, p0)
378  * boost::math::tgamma(p0 + 1)
379  * boost::math::tgamma(p0 + 1));
380 
381  // Scalar parameter which recovers the FR schemes
382  if (m_diffType == "LFRDG")
383  {
384  c0 = 0.0;
385  }
386  else if (m_diffType == "LFRSD")
387  {
388  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
389  * (ap0 * boost::math::tgamma(p0 + 1))
390  * (ap0 * boost::math::tgamma(p0 + 1)));
391  }
392  else if (m_diffType == "LFRHU")
393  {
394  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
395  * (ap0 * boost::math::tgamma(p0 + 1))
396  * (ap0 * boost::math::tgamma(p0 + 1)));
397  }
398  else if (m_diffType == "LFRcmin")
399  {
400  c0 = -2.0 / ((2.0 * p0 + 1.0)
401  * (ap0 * boost::math::tgamma(p0 + 1))
402  * (ap0 * boost::math::tgamma(p0 + 1)));
403  }
404  else if (m_diffType == "LFRinf")
405  {
406  c0 = 10000000000000000.0;
407  }
408 
409  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
410  * (ap0 * boost::math::tgamma(p0 + 1))
411  * (ap0 * boost::math::tgamma(p0 + 1));
412 
413  NekDouble overeta0 = 1.0 / (1.0 + etap0);
414 
415  // Derivative of the Legendre polynomials
416  // dLp = derivative of the Legendre polynomial of order p
417  // dLpp = derivative of the Legendre polynomial of order p+1
418  // dLpm = derivative of the Legendre polynomial of order p-1
419  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
420  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
421  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
422 
423  // Building the DG_c_Left
424  for(i = 0; i < nquad0; ++i)
425  {
426  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
427  m_dGL_xi1[n][i] += dLpp0[i];
428  m_dGL_xi1[n][i] *= overeta0;
429  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
430  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
431  }
432 
433  // Building the DG_c_Right
434  for(i = 0; i < nquad0; ++i)
435  {
436  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
437  m_dGR_xi1[n][i] += dLpp0[i];
438  m_dGR_xi1[n][i] *= overeta0;
439  m_dGR_xi1[n][i] += dLp0[i];
440  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
441  }
442  }
443  break;
444  }
445  case 2:
446  {
447  LibUtilities::BasisSharedPtr BasisFR_Left0;
448  LibUtilities::BasisSharedPtr BasisFR_Right0;
449  LibUtilities::BasisSharedPtr BasisFR_Left1;
450  LibUtilities::BasisSharedPtr BasisFR_Right1;
451 
452  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
453  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
454  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
455  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
456 
457  for (n = 0; n < nElements; ++n)
458  {
459  base = pFields[0]->GetExp(n)->GetBase();
460  nquad0 = base[0]->GetNumPoints();
461  nquad1 = base[1]->GetNumPoints();
462  nmodes0 = base[0]->GetNumModes();
463  nmodes1 = base[1]->GetNumModes();
464 
465  Array<OneD, const NekDouble> z0;
466  Array<OneD, const NekDouble> w0;
467  Array<OneD, const NekDouble> z1;
468  Array<OneD, const NekDouble> w1;
469 
470  base[0]->GetZW(z0, w0);
471  base[1]->GetZW(z1, w1);
472 
473  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
474  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
475  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
476  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
477 
478  // Auxiliary vectors to build up the auxiliary
479  // derivatives of the Legendre polynomials
480  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
481  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
482  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
483  Array<OneD,NekDouble> dLp1 (nquad1, 0.0);
484  Array<OneD,NekDouble> dLpp1 (nquad1, 0.0);
485  Array<OneD,NekDouble> dLpm1 (nquad1, 0.0);
486 
487  // Degree of the correction functions
488  int p0 = nmodes0 - 1;
489  int p1 = nmodes1 - 1;
490 
491  // Function sign to compute left correction function
492  NekDouble sign0 = pow(-1.0, p0);
493  NekDouble sign1 = pow(-1.0, p1);
494 
495  // Factorial factor to build the scheme
496  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
497  / (pow(2.0, p0)
498  * boost::math::tgamma(p0 + 1)
499  * boost::math::tgamma(p0 + 1));
500 
501  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
502  / (pow(2.0, p1)
503  * boost::math::tgamma(p1 + 1)
504  * boost::math::tgamma(p1 + 1));
505 
506  // Scalar parameter which recovers the FR schemes
507  if (m_diffType == "LFRDG")
508  {
509  c0 = 0.0;
510  c1 = 0.0;
511  }
512  else if (m_diffType == "LFRSD")
513  {
514  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
515  * (ap0 * boost::math::tgamma(p0 + 1))
516  * (ap0 * boost::math::tgamma(p0 + 1)));
517 
518  c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
519  * (ap1 * boost::math::tgamma(p1 + 1))
520  * (ap1 * boost::math::tgamma(p1 + 1)));
521  }
522  else if (m_diffType == "LFRHU")
523  {
524  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
525  * (ap0 * boost::math::tgamma(p0 + 1))
526  * (ap0 * boost::math::tgamma(p0 + 1)));
527 
528  c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
529  * (ap1 * boost::math::tgamma(p1 + 1))
530  * (ap1 * boost::math::tgamma(p1 + 1)));
531  }
532  else if (m_diffType == "LFRcmin")
533  {
534  c0 = -2.0 / ((2.0 * p0 + 1.0)
535  * (ap0 * boost::math::tgamma(p0 + 1))
536  * (ap0 * boost::math::tgamma(p0 + 1)));
537 
538  c1 = -2.0 / ((2.0 * p1 + 1.0)
539  * (ap1 * boost::math::tgamma(p1 + 1))
540  * (ap1 * boost::math::tgamma(p1 + 1)));
541  }
542  else if (m_diffType == "LFRinf")
543  {
544  c0 = 10000000000000000.0;
545  c1 = 10000000000000000.0;
546  }
547 
548  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
549  * (ap0 * boost::math::tgamma(p0 + 1))
550  * (ap0 * boost::math::tgamma(p0 + 1));
551 
552  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
553  * (ap1 * boost::math::tgamma(p1 + 1))
554  * (ap1 * boost::math::tgamma(p1 + 1));
555 
556  NekDouble overeta0 = 1.0 / (1.0 + etap0);
557  NekDouble overeta1 = 1.0 / (1.0 + etap1);
558 
559  // Derivative of the Legendre polynomials
560  // dLp = derivative of the Legendre polynomial of order p
561  // dLpp = derivative of the Legendre polynomial of order p+1
562  // dLpm = derivative of the Legendre polynomial of order p-1
563  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
564  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
565  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
566 
567  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
568  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1+1, 0.0, 0.0);
569  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1-1, 0.0, 0.0);
570 
571  // Building the DG_c_Left
572  for(i = 0; i < nquad0; ++i)
573  {
574  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
575  m_dGL_xi1[n][i] += dLpp0[i];
576  m_dGL_xi1[n][i] *= overeta0;
577  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
578  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
579  }
580 
581  // Building the DG_c_Left
582  for(i = 0; i < nquad1; ++i)
583  {
584  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
585  m_dGL_xi2[n][i] += dLpp1[i];
586  m_dGL_xi2[n][i] *= overeta1;
587  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
588  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
589  }
590 
591  // Building the DG_c_Right
592  for(i = 0; i < nquad0; ++i)
593  {
594  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
595  m_dGR_xi1[n][i] += dLpp0[i];
596  m_dGR_xi1[n][i] *= overeta0;
597  m_dGR_xi1[n][i] += dLp0[i];
598  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
599  }
600 
601  // Building the DG_c_Right
602  for(i = 0; i < nquad1; ++i)
603  {
604  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
605  m_dGR_xi2[n][i] += dLpp1[i];
606  m_dGR_xi2[n][i] *= overeta1;
607  m_dGR_xi2[n][i] += dLp1[i];
608  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
609  }
610  }
611  break;
612  }
613  case 3:
614  {
615  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
616  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
617  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
618  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
619  m_dGL_xi3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
620  m_dGR_xi3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
621 
622  for (n = 0; n < nElements; ++n)
623  {
624  base = pFields[0]->GetExp(n)->GetBase();
625  nquad0 = base[0]->GetNumPoints();
626  nquad1 = base[1]->GetNumPoints();
627  nquad2 = base[2]->GetNumPoints();
628  nmodes0 = base[0]->GetNumModes();
629  nmodes1 = base[1]->GetNumModes();
630  nmodes2 = base[2]->GetNumModes();
631 
632  Array<OneD, const NekDouble> z0;
633  Array<OneD, const NekDouble> w0;
634  Array<OneD, const NekDouble> z1;
635  Array<OneD, const NekDouble> w1;
636  Array<OneD, const NekDouble> z2;
637  Array<OneD, const NekDouble> w2;
638 
639  base[0]->GetZW(z0, w0);
640  base[1]->GetZW(z1, w1);
641  base[1]->GetZW(z2, w2);
642 
643  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
644  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
645  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
646  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
647  m_dGL_xi3[n] = Array<OneD, NekDouble>(nquad2);
648  m_dGR_xi3[n] = Array<OneD, NekDouble>(nquad2);
649 
650  // Auxiliary vectors to build up the auxiliary
651  // derivatives of the Legendre polynomials
652  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
653  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
654  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
655  Array<OneD,NekDouble> dLp1 (nquad1, 0.0);
656  Array<OneD,NekDouble> dLpp1 (nquad1, 0.0);
657  Array<OneD,NekDouble> dLpm1 (nquad1, 0.0);
658  Array<OneD,NekDouble> dLp2 (nquad2, 0.0);
659  Array<OneD,NekDouble> dLpp2 (nquad2, 0.0);
660  Array<OneD,NekDouble> dLpm2 (nquad2, 0.0);
661 
662  // Degree of the correction functions
663  int p0 = nmodes0 - 1;
664  int p1 = nmodes1 - 1;
665  int p2 = nmodes2 - 1;
666 
667  // Function sign to compute left correction function
668  NekDouble sign0 = pow(-1.0, p0);
669  NekDouble sign1 = pow(-1.0, p1);
670  NekDouble sign2 = pow(-1.0, p2);
671 
672  // Factorial factor to build the scheme
673  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
674  / (pow(2.0, p0)
675  * boost::math::tgamma(p0 + 1)
676  * boost::math::tgamma(p0 + 1));
677 
678  // Factorial factor to build the scheme
679  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
680  / (pow(2.0, p1)
681  * boost::math::tgamma(p1 + 1)
682  * boost::math::tgamma(p1 + 1));
683 
684  // Factorial factor to build the scheme
685  NekDouble ap2 = boost::math::tgamma(2 * p2 + 1)
686  / (pow(2.0, p2)
687  * boost::math::tgamma(p2 + 1)
688  * boost::math::tgamma(p2 + 1));
689 
690  // Scalar parameter which recovers the FR schemes
691  if (m_diffType == "LFRDG")
692  {
693  c0 = 0.0;
694  c1 = 0.0;
695  c2 = 0.0;
696  }
697  else if (m_diffType == "LFRSD")
698  {
699  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
700  * (ap0 * boost::math::tgamma(p0 + 1))
701  * (ap0 * boost::math::tgamma(p0 + 1)));
702 
703  c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
704  * (ap1 * boost::math::tgamma(p1 + 1))
705  * (ap1 * boost::math::tgamma(p1 + 1)));
706 
707  c2 = 2.0 * p2 / ((2.0 * p2 + 1.0) * (p2 + 1.0)
708  * (ap2 * boost::math::tgamma(p2 + 1))
709  * (ap2 * boost::math::tgamma(p2 + 1)));
710  }
711  else if (m_diffType == "LFRHU")
712  {
713  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
714  * (ap0 * boost::math::tgamma(p0 + 1))
715  * (ap0 * boost::math::tgamma(p0 + 1)));
716 
717  c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
718  * (ap1 * boost::math::tgamma(p1 + 1))
719  * (ap1 * boost::math::tgamma(p1 + 1)));
720 
721  c2 = 2.0 * (p2 + 1.0) / ((2.0 * p2 + 1.0) * p2
722  * (ap2 * boost::math::tgamma(p2 + 1))
723  * (ap2 * boost::math::tgamma(p2 + 1)));
724  }
725  else if (m_diffType == "LFRcmin")
726  {
727  c0 = -2.0 / ((2.0 * p0 + 1.0)
728  * (ap0 * boost::math::tgamma(p0 + 1))
729  * (ap0 * boost::math::tgamma(p0 + 1)));
730 
731  c1 = -2.0 / ((2.0 * p1 + 1.0)
732  * (ap1 * boost::math::tgamma(p1 + 1))
733  * (ap1 * boost::math::tgamma(p1 + 1)));
734 
735  c2 = -2.0 / ((2.0 * p2 + 1.0)
736  * (ap2 * boost::math::tgamma(p2 + 1))
737  * (ap2 * boost::math::tgamma(p2 + 1)));
738  }
739  else if (m_diffType == "LFRinf")
740  {
741  c0 = 10000000000000000.0;
742  c1 = 10000000000000000.0;
743  c2 = 10000000000000000.0;
744  }
745 
746  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
747  * (ap0 * boost::math::tgamma(p0 + 1))
748  * (ap0 * boost::math::tgamma(p0 + 1));
749 
750  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
751  * (ap1 * boost::math::tgamma(p1 + 1))
752  * (ap1 * boost::math::tgamma(p1 + 1));
753 
754  NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0)
755  * (ap2 * boost::math::tgamma(p2 + 1))
756  * (ap2 * boost::math::tgamma(p2 + 1));
757 
758  NekDouble overeta0 = 1.0 / (1.0 + etap0);
759  NekDouble overeta1 = 1.0 / (1.0 + etap1);
760  NekDouble overeta2 = 1.0 / (1.0 + etap2);
761 
762  // Derivative of the Legendre polynomials
763  // dLp = derivative of the Legendre polynomial of order p
764  // dLpp = derivative of the Legendre polynomial of order p+1
765  // dLpm = derivative of the Legendre polynomial of order p-1
766  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
767  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
768  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
769 
770  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
771  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1+1, 0.0, 0.0);
772  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1-1, 0.0, 0.0);
773 
774  Polylib::jacobd(nquad2, z2.data(), &(dLp2[0]), p2, 0.0, 0.0);
775  Polylib::jacobd(nquad2, z2.data(), &(dLpp2[0]), p2+1, 0.0, 0.0);
776  Polylib::jacobd(nquad2, z2.data(), &(dLpm2[0]), p2-1, 0.0, 0.0);
777 
778  // Building the DG_c_Left
779  for(i = 0; i < nquad0; ++i)
780  {
781  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
782  m_dGL_xi1[n][i] += dLpp0[i];
783  m_dGL_xi1[n][i] *= overeta0;
784  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
785  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
786  }
787 
788  // Building the DG_c_Left
789  for(i = 0; i < nquad1; ++i)
790  {
791  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
792  m_dGL_xi2[n][i] += dLpp1[i];
793  m_dGL_xi2[n][i] *= overeta1;
794  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
795  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
796  }
797 
798  // Building the DG_c_Left
799  for(i = 0; i < nquad2; ++i)
800  {
801  m_dGL_xi3[n][i] = etap2 * dLpm2[i];
802  m_dGL_xi3[n][i] += dLpp2[i];
803  m_dGL_xi3[n][i] *= overeta2;
804  m_dGL_xi3[n][i] = dLp2[i] - m_dGL_xi3[n][i];
805  m_dGL_xi3[n][i] = 0.5 * sign2 * m_dGL_xi3[n][i];
806  }
807 
808  // Building the DG_c_Right
809  for(i = 0; i < nquad0; ++i)
810  {
811  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
812  m_dGR_xi1[n][i] += dLpp0[i];
813  m_dGR_xi1[n][i] *= overeta0;
814  m_dGR_xi1[n][i] += dLp0[i];
815  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
816  }
817 
818  // Building the DG_c_Right
819  for(i = 0; i < nquad1; ++i)
820  {
821  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
822  m_dGR_xi2[n][i] += dLpp1[i];
823  m_dGR_xi2[n][i] *= overeta1;
824  m_dGR_xi2[n][i] += dLp1[i];
825  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
826  }
827 
828  // Building the DG_c_Right
829  for(i = 0; i < nquad2; ++i)
830  {
831  m_dGR_xi3[n][i] = etap2 * dLpm2[i];
832  m_dGR_xi3[n][i] += dLpp2[i];
833  m_dGR_xi3[n][i] *= overeta2;
834  m_dGR_xi3[n][i] += dLp2[i];
835  m_dGR_xi3[n][i] = 0.5 * m_dGR_xi3[n][i];
836  }
837  }
838  break;
839  }
840  default:
841  {
842  ASSERTL0(false,"Expansion dimension not recognised");
843  break;
844  }
845  }
846  }
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition: DiffusionLFR.h:67
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition: DiffusionLFR.h:66
Polylib::jacobd
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition: Polylib.cpp:1134

References ASSERTL0, Polylib::jacobd(), m_dGL_xi1, m_dGL_xi2, m_dGL_xi3, m_dGR_xi1, m_dGR_xi2, m_dGR_xi3, and m_diffType.

Referenced by v_InitObject().

◆ v_SetupMetrics()

void Nektar::SolverUtils::DiffusionLFR::v_SetupMetrics ( LibUtilities::SessionReaderSharedPtr  pSession,
Array< OneD, MultiRegions::ExpListSharedPtr >  pFields 
)
protectedvirtual

Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform the fluxes at the interfaces of each element from the physical space to the standard space).

This routine calls the function #GetTraceQFactors to compute and store the metric factors following an anticlockwise conventions along the edges/faces of each element. Note: for 1D problem the transformation is not needed.

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.
Todo:
Add the metric terms for 3D Hexahedra.

Definition at line 174 of file DiffusionLFR.cpp.

177  {
178  boost::ignore_unused(pSession);
179 
180  int i, n;
181  int nquad0, nquad1;
182  int phys_offset;
183  int nLocalSolutionPts;
184  int nElements = pFields[0]->GetExpSize();
185  int nDimensions = pFields[0]->GetCoordim(0);
186  int nSolutionPts = pFields[0]->GetTotPoints();
187  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
188 
189  m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(nDimensions);
190  for (i = 0; i < nDimensions; ++i)
191  {
192  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
193  }
194  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
195 
196  Array<TwoD, const NekDouble> gmat;
197  Array<OneD, const NekDouble> jac;
198 
199  m_jac = Array<OneD, NekDouble>(nSolutionPts);
200 
201  Array<OneD, NekDouble> auxArray1;
202  Array<OneD, LibUtilities::BasisSharedPtr> base;
203  LibUtilities::PointsKeyVector ptsKeys;
204 
205  switch (nDimensions)
206  {
207  case 1:
208  {
209  for (n = 0; n < nElements; ++n)
210  {
211  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
212  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
213  phys_offset = pFields[0]->GetPhys_Offset(n);
214  jac = pFields[0]->GetExp(n)
215  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
216  ->GetMetricInfo()->GetJac(ptsKeys);
217  for (i = 0; i < nLocalSolutionPts; ++i)
218  {
219  m_jac[i+phys_offset] = jac[0];
220  }
221  }
222  break;
223  }
224  case 2:
225  {
226  m_gmat = Array<OneD, Array<OneD, NekDouble> >(4);
227  m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
228  m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
229  m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
230  m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
231 
232  m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble> >(nElements);
233  m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
234  m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
235  m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
236 
237  for (n = 0; n < nElements; ++n)
238  {
239  base = pFields[0]->GetExp(n)->GetBase();
240  nquad0 = base[0]->GetNumPoints();
241  nquad1 = base[1]->GetNumPoints();
242 
243  m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
244  m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
245  m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
246  m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
247 
248  // Extract the Q factors at each edge point
249  pFields[0]->GetExp(n)->GetTraceQFactors(
250  0, auxArray1 = m_Q2D_e0[n]);
251  pFields[0]->GetExp(n)->GetTraceQFactors(
252  1, auxArray1 = m_Q2D_e1[n]);
253  pFields[0]->GetExp(n)->GetTraceQFactors(
254  2, auxArray1 = m_Q2D_e2[n]);
255  pFields[0]->GetExp(n)->GetTraceQFactors(
256  3, auxArray1 = m_Q2D_e3[n]);
257 
258  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
259  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
260  phys_offset = pFields[0]->GetPhys_Offset(n);
261 
262  jac = pFields[0]->GetExp(n)
263  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
264  ->GetMetricInfo()->GetJac(ptsKeys);
265  gmat = pFields[0]->GetExp(n)
266  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
267  ->GetMetricInfo()->GetDerivFactors(ptsKeys);
268 
269  if (pFields[0]->GetExp(n)
270  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
271  ->GetMetricInfo()->GetGtype()
272  == SpatialDomains::eDeformed)
273  {
274  for (i = 0; i < nLocalSolutionPts; ++i)
275  {
276  m_jac[i+phys_offset] = jac[i];
277  m_gmat[0][i+phys_offset] = gmat[0][i];
278  m_gmat[1][i+phys_offset] = gmat[1][i];
279  m_gmat[2][i+phys_offset] = gmat[2][i];
280  m_gmat[3][i+phys_offset] = gmat[3][i];
281  }
282  }
283  else
284  {
285  for (i = 0; i < nLocalSolutionPts; ++i)
286  {
287  m_jac[i+phys_offset] = jac[0];
288  m_gmat[0][i+phys_offset] = gmat[0][0];
289  m_gmat[1][i+phys_offset] = gmat[1][0];
290  m_gmat[2][i+phys_offset] = gmat[2][0];
291  m_gmat[3][i+phys_offset] = gmat[3][0];
292  }
293  }
294  }
295  break;
296  }
297  case 3:
298  {
299  ASSERTL0(false,"3DFR Metric terms not implemented yet");
300  break;
301  }
302  default:
303  {
304  ASSERTL0(false, "Expansion dimension not recognised");
305  break;
306  }
307  }
308  }

References ASSERTL0, Nektar::SpatialDomains::eDeformed, Nektar::LocalRegions::Expansion::GetMetricInfo(), m_gmat, m_jac, m_Q2D_e0, m_Q2D_e1, m_Q2D_e2, m_Q2D_e3, and m_traceNormals.

Referenced by v_InitObject().

◆ v_WeakPenaltyforScalar()

void Nektar::SolverUtils::DiffusionLFR::v_WeakPenaltyforScalar ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const int  var,
const Array< OneD, const NekDouble > &  ufield,
Array< OneD, NekDouble > &  penaltyflux 
)
protectedvirtual

Imposes appropriate bcs for the 1st order derivatives.

Definition at line 1241 of file DiffusionLFR.cpp.

1246  {
1247  // Variable initialisation
1248  int i, e, id1, id2;
1249  int nBndEdgePts, nBndEdges;
1250  int cnt = 0;
1251  int nBndRegions = fields[var]->GetBndCondExpansions().size();
1252  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1253 
1254  // Extract physical values of the fields at the trace space
1255  Array<OneD, NekDouble > uplus(nTracePts);
1256  fields[var]->ExtractTracePhys(ufield, uplus);
1257 
1258  // Impose boundary conditions
1259  for (i = 0; i < nBndRegions; ++i)
1260  {
1261  // Number of boundary edges related to bcRegion 'i'
1262  nBndEdges = fields[var]->
1263  GetBndCondExpansions()[i]->GetExpSize();
1264 
1265  // Weakly impose boundary conditions by modifying flux values
1266  for (e = 0; e < nBndEdges ; ++e)
1267  {
1268  // Number of points on the expansion
1269  nBndEdgePts = fields[var]->
1270  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
1271 
1272  // Offset of the boundary expansion
1273  id1 = fields[var]->
1274  GetBndCondExpansions()[i]->GetPhys_Offset(e);
1275 
1276  // Offset of the trace space related to boundary expansion
1277  id2 = fields[0]->GetTrace()->
1278  GetPhys_Offset(fields[0]->GetTraceMap()->
1279  GetBndCondIDToGlobalTraceID(cnt++));
1280 
1281  // Dirichlet bcs ==> uflux = gD
1282  if (fields[var]->GetBndConditions()[i]->
1283  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1284  {
1285  Vmath::Vcopy(nBndEdgePts,
1286  &(fields[var]->
1287  GetBndCondExpansions()[i]->
1288  GetPhys())[id1], 1,
1289  &penaltyflux[id2], 1);
1290  }
1291  // Neumann bcs ==> uflux = u+
1292  else if ((fields[var]->GetBndConditions()[i])->
1293  GetBoundaryConditionType() == SpatialDomains::eNeumann)
1294  {
1295  Vmath::Vcopy(nBndEdgePts,
1296  &uplus[id2], 1,
1297  &penaltyflux[id2], 1);
1298  }
1299  }
1300  }
1301  }

References Nektar::SpatialDomains::eDirichlet, Nektar::SpatialDomains::eNeumann, and Vmath::Vcopy().

Referenced by v_NumFluxforScalar().

◆ v_WeakPenaltyforVector()

void Nektar::SolverUtils::DiffusionLFR::v_WeakPenaltyforVector ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const int  var,
const int  dir,
const Array< OneD, const NekDouble > &  qfield,
Array< OneD, NekDouble > &  penaltyflux,
NekDouble  C11 
)
protectedvirtual

Imposes appropriate bcs for the 2nd order derivatives.

Definition at line 1405 of file DiffusionLFR.cpp.

1412  {
1413  boost::ignore_unused(C11);
1414 
1415  int i, e, id1, id2;
1416  int nBndEdges, nBndEdgePts;
1417  int nBndRegions = fields[var]->GetBndCondExpansions().size();
1418  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1419 
1420  Array<OneD, NekDouble > uterm(nTracePts);
1421  Array<OneD, NekDouble > qtemp(nTracePts);
1422  int cnt = 0;
1423 
1424  fields[var]->ExtractTracePhys(qfield, qtemp);
1425 
1426  for (i = 0; i < nBndRegions; ++i)
1427  {
1428  nBndEdges = fields[var]->
1429  GetBndCondExpansions()[i]->GetExpSize();
1430 
1431  // Weakly impose boundary conditions by modifying flux values
1432  for (e = 0; e < nBndEdges ; ++e)
1433  {
1434  nBndEdgePts = fields[var]->
1435  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
1436 
1437  id1 = fields[var]->
1438  GetBndCondExpansions()[i]->GetPhys_Offset(e);
1439 
1440  id2 = fields[0]->GetTrace()->
1441  GetPhys_Offset(fields[0]->GetTraceMap()->
1442  GetBndCondIDToGlobalTraceID(cnt++));
1443 
1444  // For Dirichlet boundary condition:
1445  //qflux = q+ - C_11 (u+ - g_D) (nx, ny)
1446  if (fields[var]->GetBndConditions()[i]->
1447  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1448  {
1449  Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1450  &qtemp[id2], 1, &penaltyflux[id2], 1);
1451  }
1452  // For Neumann boundary condition: qflux = g_N
1453  else if ((fields[var]->GetBndConditions()[i])->
1454  GetBoundaryConditionType() == SpatialDomains::eNeumann)
1455  {
1456  Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1457  &(fields[var]->GetBndCondExpansions()[i]->
1458  GetPhys())[id1], 1, &penaltyflux[id2], 1);
1459  }
1460  }
1461  }
1462  }

References Nektar::SpatialDomains::eDirichlet, Nektar::SpatialDomains::eNeumann, m_traceNormals, and Vmath::Vmul().

Referenced by v_NumFluxforVector().

Member Data Documentation

◆ m_BD1

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_BD1
protected

Definition at line 80 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_D1

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_D1
protected

Definition at line 81 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_DD1

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_DD1
protected

Definition at line 82 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_DFC1

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_DFC1
protected

Definition at line 79 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_DFC2

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_DFC2
protected

Definition at line 84 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_dGL_xi1

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1

Definition at line 62 of file DiffusionLFR.h.

Referenced by v_DerCFlux_1D(), v_DerCFlux_2D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupCFunctions().

◆ m_dGL_xi2

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_dGL_xi2

Definition at line 64 of file DiffusionLFR.h.

Referenced by v_DerCFlux_2D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupCFunctions().

◆ m_dGL_xi3

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_dGL_xi3

Definition at line 66 of file DiffusionLFR.h.

Referenced by v_SetupCFunctions().

◆ m_dGR_xi1

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1

Definition at line 63 of file DiffusionLFR.h.

Referenced by v_DerCFlux_1D(), v_DerCFlux_2D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupCFunctions().

◆ m_dGR_xi2

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_dGR_xi2

Definition at line 65 of file DiffusionLFR.h.

Referenced by v_DerCFlux_2D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupCFunctions().

◆ m_dGR_xi3

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_dGR_xi3

Definition at line 67 of file DiffusionLFR.h.

Referenced by v_SetupCFunctions().

◆ m_diffType

std::string Nektar::SolverUtils::DiffusionLFR::m_diffType
protected

Definition at line 91 of file DiffusionLFR.h.

Referenced by v_SetupCFunctions().

◆ m_divFC

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_divFC
protected

Definition at line 86 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_divFD

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_divFD
protected

Definition at line 85 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_DU1

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_DU1
protected

Definition at line 78 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_gmat

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_gmat

Definition at line 55 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_SetupMetrics().

◆ m_IF1

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_IF1
protected

Definition at line 77 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_IF2

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_IF2
protected

Definition at line 83 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_Ixm

DNekMatSharedPtr Nektar::SolverUtils::DiffusionLFR::m_Ixm

Definition at line 68 of file DiffusionLFR.h.

◆ m_Ixp

DNekMatSharedPtr Nektar::SolverUtils::DiffusionLFR::m_Ixp

Definition at line 69 of file DiffusionLFR.h.

◆ m_jac

Array<OneD, NekDouble> Nektar::SolverUtils::DiffusionLFR::m_jac

Definition at line 54 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_SetupMetrics().

◆ m_Q2D_e0

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_Q2D_e0

Definition at line 57 of file DiffusionLFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupMetrics().

◆ m_Q2D_e1

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_Q2D_e1

Definition at line 58 of file DiffusionLFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupMetrics().

◆ m_Q2D_e2

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_Q2D_e2

Definition at line 59 of file DiffusionLFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupMetrics().

◆ m_Q2D_e3

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_Q2D_e3

Definition at line 60 of file DiffusionLFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupMetrics().

◆ m_session

LibUtilities::SessionReaderSharedPtr Nektar::SolverUtils::DiffusionLFR::m_session
protected

Definition at line 75 of file DiffusionLFR.h.

Referenced by v_InitObject().

◆ m_tmp1

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_tmp1
protected

Definition at line 88 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_tmp2

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_tmp2
protected

Definition at line 89 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_traceNormals

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_traceNormals
protected

Definition at line 74 of file DiffusionLFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), v_NumFluxforScalar(), v_NumFluxforVector(), v_SetupMetrics(), and v_WeakPenaltyforVector().

◆ type

std::string Nektar::SolverUtils::DiffusionLFR::type
static
Initial value:
= {
GetDiffusionFactory().RegisterCreatorFunction(
"LFRDG", DiffusionLFR::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRSD", DiffusionLFR::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRHU", DiffusionLFR::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRcmin", DiffusionLFR::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRcinf", DiffusionLFR::create)}
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:200
Nektar::SolverUtils::DiffusionLFR::create
static DiffusionSharedPtr create(std::string diffType)
Definition: DiffusionLFR.h:47
Nektar::SolverUtils::GetDiffusionFactory
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:41

Definition at line 52 of file DiffusionLFR.h.

Generated on Wed Nov 24 2021 10:33:36 for Nektar++ by doxygen 1.9.1