Nektar++
Static Public Member Functions | Public Attributes | Static Public Attributes | Protected Member Functions | Protected Attributes | Private Member Functions | 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, NekDoublem_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, NekDoublem_divVel
 Params for Ducros sensor. More...
 
Array< OneD, NekDoublem_divVelSquare
 
Array< OneD, NekDoublem_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) override
 Initiliase DiffusionLFR objects and store them before starting the time-stepping. 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) override
 Calculate FR Diffusion for the linear problems using an LDG interface flux. 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
 
DiffusionFluxCons m_FunctorDiffusionfluxCons
 
DiffusionFluxCons m_FunctorDiffusionfluxConsTrace
 
SpecialBndTreat m_SpecialBndTreat
 
DiffusionSymmFluxCons m_FunctorSymmetricfluxCons
 
NekDouble m_time = 0.0
 

Private Member Functions

void 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...
 
void 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...
 
void 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...
 
void 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...
 
void 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...
 
void 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...
 
void 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...
 
void 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...
 
void 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...
 
void 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...
 

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)
 
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 SetSpecialBndTreat (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 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  }
DiffusionLFR(std::string diffType)
DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion term....
std::shared_ptr< SolverUtils::Diffusion > DiffusionSharedPtr
A shared pointer to an EquationSystem object.
Definition: Diffusion.h:550

References DiffusionLFR().

◆ DerCFlux_1D()

void Nektar::SolverUtils::DiffusionLFR::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 
)
private

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 1477 of file DiffusionLFR.cpp.

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

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

Referenced by v_Diffuse().

◆ DerCFlux_2D()

void Nektar::SolverUtils::DiffusionLFR::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 
)
private

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 1598 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 
1619  Array<OneD, NekDouble> auxArray1, auxArray2;
1620  Array<OneD, LibUtilities::BasisSharedPtr> base;
1621  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1622  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]
1633  ->GetExp(n)
1634  ->as<LocalRegions::Expansion2D>()
1635  ->GetGeom2D()
1636  ->GetMetricInfo()
1637  ->GetJac(ptsKeys);
1638 
1639  base = fields[0]->GetExp(n)->GetBase();
1640  nquad0 = base[0]->GetNumPoints();
1641  nquad1 = base[1]->GetNumPoints();
1642 
1643  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1644  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1645  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1646  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1647 
1648  // Loop on the edges
1649  for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1650  {
1651  // Number of edge points of edge 'e'
1652  nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1653 
1654  // Array for storing volumetric fluxes on each edge
1655  Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1656 
1657  // Offset of the trace space correspondent to edge 'e'
1658  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1659  elmtToTrace[n][e]->GetElmtId());
1660 
1661  // Get the normals of edge 'e'
1662  // const Array<OneD, const Array<OneD, NekDouble> > &normals =
1663  // fields[0]->GetExp(n)->GetTraceNormal(e);
1664 
1665  // Extract the edge values of the volumetric fluxes
1666  // on edge 'e' and order them accordingly to the order
1667  // of the trace space
1668  fields[0]->GetExp(n)->GetTracePhysVals(
1669  e, elmtToTrace[n][e], flux + phys_offset, auxArray1 = tmparray);
1670 
1671  // Compute the fluxJumps per each edge 'e' and each
1672  // flux point
1673  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1674  Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1, &tmparray[0], 1,
1675  &fluxJumps[0], 1);
1676 
1677  // Check the ordering of the fluxJumps and reverse
1678  // it in case of backward definition of edge 'e'
1679  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1681  {
1682  Vmath::Reverse(nEdgePts, &fluxJumps[0], 1, &fluxJumps[0], 1);
1683  }
1684 
1685  // Deformed elements
1686  if (fields[0]
1687  ->GetExp(n)
1688  ->as<LocalRegions::Expansion2D>()
1689  ->GetGeom2D()
1690  ->GetMetricInfo()
1691  ->GetGtype() == SpatialDomains::eDeformed)
1692  {
1693  // Extract the Jacobians along edge 'e'
1694  Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1695 
1696  fields[0]->GetExp(n)->GetTracePhysVals(
1697  e, elmtToTrace[n][e], jac, auxArray1 = jacEdge);
1698 
1699  // Check the ordering of the fluxJumps and reverse
1700  // it in case of backward definition of edge 'e'
1701  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1703  {
1704  Vmath::Reverse(nEdgePts, &jacEdge[0], 1, &jacEdge[0], 1);
1705  }
1706 
1707  // Multiply the fluxJumps by the edge 'e' Jacobians
1708  // to bring the fluxJumps into the standard space
1709  for (j = 0; j < nEdgePts; j++)
1710  {
1711  fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1712  }
1713  }
1714  // Non-deformed elements
1715  else
1716  {
1717  // Multiply the fluxJumps by the edge 'e' Jacobians
1718  // to bring the fluxJumps into the standard space
1719  Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1, 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, &divCFluxE3[0], 1,
1790  &derCFlux[phys_offset], 1);
1791  }
1792  else if (direction == 1)
1793  {
1794  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE2[0], 1,
1795  &derCFlux[phys_offset], 1);
1796  }
1797  }
1798 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFR.h:64
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFR.h:65
@ eDeformed
Geometry is curved or has non-constant factors.
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1286
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:419

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().

◆ DivCFlux_2D()

void Nektar::SolverUtils::DiffusionLFR::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 
)
private

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>> &elmtToTrace =
1838  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(e, elmtToTrace[n][e],
1879  fluxX1 + phys_offset,
1880  auxArray1 = tmparrayX1);
1881 
1882  // Extract the edge values of flux-y on edge e and order
1883  // them accordingly to the order of the trace space
1884  fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1885  fluxX2 + phys_offset,
1886  auxArray1 = tmparrayX2);
1887 
1888  // Multiply the edge components of the flux by the normal
1889  for (i = 0; i < nEdgePts; ++i)
1890  {
1891  fluxN[i] = tmparrayX1[i] * m_traceNormals[0][trace_offset + i] +
1892  tmparrayX2[i] * m_traceNormals[1][trace_offset + i];
1893  }
1894 
1895  // Subtract to the Riemann flux the discontinuous flux
1896  Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1, &fluxN[0], 1,
1897  &fluxJumps[0], 1);
1898 
1899  // Check the ordering of the jump vectors
1900  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1902  {
1903  Vmath::Reverse(nEdgePts, auxArray1 = fluxJumps, 1,
1904  auxArray2 = fluxJumps, 1);
1905  }
1906 
1907  for (i = 0; i < nEdgePts; ++i)
1908  {
1909  if (m_traceNormals[0][trace_offset + i] != normals[0][i] ||
1910  m_traceNormals[1][trace_offset + i] != normals[1][i])
1911  {
1912  fluxJumps[i] = -fluxJumps[i];
1913  }
1914  }
1915 
1916  // Multiply jumps by derivatives of the correction functions
1917  switch (e)
1918  {
1919  case 0:
1920  for (i = 0; i < nquad0; ++i)
1921  {
1922  // Multiply fluxJumps by Q factors
1923  fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1924 
1925  for (j = 0; j < nquad1; ++j)
1926  {
1927  cnt = i + j * nquad0;
1928  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1929  }
1930  }
1931  break;
1932  case 1:
1933  for (i = 0; i < nquad1; ++i)
1934  {
1935  // Multiply fluxJumps by Q factors
1936  fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1937 
1938  for (j = 0; j < nquad0; ++j)
1939  {
1940  cnt = (nquad0)*i + j;
1941  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1942  }
1943  }
1944  break;
1945  case 2:
1946  for (i = 0; i < nquad0; ++i)
1947  {
1948  // Multiply fluxJumps by Q factors
1949  fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1950 
1951  for (j = 0; j < nquad1; ++j)
1952  {
1953  cnt = j * nquad0 + i;
1954  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1955  }
1956  }
1957  break;
1958  case 3:
1959  for (i = 0; i < nquad1; ++i)
1960  {
1961  // Multiply fluxJumps by Q factors
1962  fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1963  for (j = 0; j < nquad0; ++j)
1964  {
1965  cnt = j + i * nquad0;
1966  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1967  }
1968  }
1969  break;
1970 
1971  default:
1972  ASSERTL0(false, "edge value (< 3) is out of range");
1973  break;
1974  }
1975  }
1976 
1977  // Sum each edge contribution
1978  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
1979  &divCFlux[phys_offset], 1);
1980 
1981  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1982  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1983 
1984  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1985  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1986  }
1987 }
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:74
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFR.h:57
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFR.h:60
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFR.h:59
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().

◆ DivCFlux_2D_Gauss()

void Nektar::SolverUtils::DiffusionLFR::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 
)
private

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 2003 of file DiffusionLFR.cpp.

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

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().

◆ NumFluxforScalar()

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

Builds the numerical flux for the 1st order derivatives.

Definition at line 1178 of file DiffusionLFR.cpp.

1182 {
1183  int i, j;
1184  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1185  int nvariables = fields.size();
1186  int nDim = fields[0]->GetCoordim(0);
1187 
1188  Array<OneD, NekDouble> Fwd(nTracePts);
1189  Array<OneD, NekDouble> Bwd(nTracePts);
1190  Array<OneD, NekDouble> Vn(nTracePts, 0.0);
1191  Array<OneD, NekDouble> fluxtemp(nTracePts, 0.0);
1192 
1193  // Get the normal velocity Vn
1194  for (i = 0; i < nDim; ++i)
1195  {
1196  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1, Vn, 1, Vn, 1);
1197  }
1198 
1199  // Get the sign of (v \cdot n), v = an arbitrary vector
1200  // Evaluate upwind flux:
1201  // uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
1202  for (j = 0; j < nDim; ++j)
1203  {
1204  for (i = 0; i < nvariables; ++i)
1205  {
1206  // Compute Fwd and Bwd value of ufield of i direction
1207  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1208 
1209  // if Vn >= 0, flux = uFwd, i.e.,
1210  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
1211  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
1212 
1213  // else if Vn < 0, flux = uBwd, i.e.,
1214  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
1215  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
1216 
1217  fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, fluxtemp);
1218 
1219  // Imposing weak boundary condition with flux
1220  // if Vn >= 0, uflux = uBwd at Neumann, i.e.,
1221  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
1222  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
1223 
1224  // if Vn >= 0, uflux = uFwd at Neumann, i.e.,
1225  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
1226  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
1227 
1228  if (fields[0]->GetBndCondExpansions().size())
1229  {
1230  WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
1231  }
1232 
1233  // if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
1234  // i.e,
1235  // edge::eForward, uFwd \‍(\tan_{\xi}^Fwd \cdot \vec{n})
1236  // edge::eBackward, uFwd \‍(\tan_{\xi}^Bwd \cdot \vec{n})
1237 
1238  // else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}),
1239  // i.e,
1240  // edge::eForward, uBwd \‍(\tan_{\xi}^Fwd \cdot \vec{n})
1241  // edge::eBackward, uBwd \‍(\tan_{\xi}^Bwd \cdot \vec{n})
1242 
1243  Vmath::Vcopy(nTracePts, fluxtemp, 1, uflux[i][j], 1);
1244  }
1245  }
1246 }
void 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.
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:622

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

Referenced by v_Diffuse().

◆ NumFluxforVector()

void Nektar::SolverUtils::DiffusionLFR::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 
)
private

Build the numerical flux for the 2nd order derivatives.

Definition at line 1315 of file DiffusionLFR.cpp.

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

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

Referenced by v_Diffuse().

◆ SetupCFunctions()

void Nektar::SolverUtils::DiffusionLFR::SetupCFunctions ( LibUtilities::SessionReaderSharedPtr  pSession,
Array< OneD, MultiRegions::ExpListSharedPtr pFields 
)
private

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 Functions #DivCFlux_1D, DivCFlux_2D, #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 329 of file DiffusionLFR.cpp.

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

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().

◆ SetupMetrics()

void Nektar::SolverUtils::DiffusionLFR::SetupMetrics ( LibUtilities::SessionReaderSharedPtr  pSession,
Array< OneD, MultiRegions::ExpListSharedPtr pFields 
)
private

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 166 of file DiffusionLFR.cpp.

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

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_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 
)
overrideprotectedvirtual

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

Reimplemented from Nektar::SolverUtils::Diffusion.

Definition at line 874 of file DiffusionLFR.cpp.

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

References ASSERTL0, DerCFlux_1D(), DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), 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, NumFluxforScalar(), NumFluxforVector(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), and Vmath::Vmul().

◆ v_InitObject()

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

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

This routine calls the functions SetupMetrics and 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  SetupMetrics(pSession, pFields);
91  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>>>(nConvectiveFields);
101  m_DU1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
102  m_DFC1 =
103  Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
104  m_BD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
105  m_D1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
106  m_DD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
107  m_IF2 = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
108  m_DFC2 =
109  Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
110  m_divFD = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
111  m_divFC = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
112 
113  m_tmp1 =
114  Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
115  m_tmp2 =
116  Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
117 
118  for (i = 0; i < nConvectiveFields; ++i)
119  {
120  m_IF1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
121  m_DU1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
122  m_DFC1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
123  m_BD1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
124  m_D1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
125  m_DD1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
126  m_IF2[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
127  m_DFC2[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
128  m_divFD[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
129  m_divFC[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
130 
131  m_tmp1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
132  m_tmp2[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
133 
134  for (j = 0; j < nDim; ++j)
135  {
136  m_IF1[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
137  m_DU1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
138  m_DFC1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
139  m_BD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
140  m_D1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
141  m_DD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
142  m_DFC2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
143 
144  m_tmp1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
145  m_tmp2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
146  }
147  }
148 }
void 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...
LibUtilities::SessionReaderSharedPtr m_session
Definition: DiffusionLFR.h:75
void 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...

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, SetupCFunctions(), and SetupMetrics().

◆ WeakPenaltyforScalar()

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

Imposes appropriate bcs for the 1st order derivatives.

Definition at line 1252 of file DiffusionLFR.cpp.

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

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

Referenced by NumFluxforScalar().

◆ WeakPenaltyforVector()

void Nektar::SolverUtils::DiffusionLFR::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 
)
private

Imposes appropriate bcs for the 2nd order derivatives.

Definition at line 1409 of file DiffusionLFR.cpp.

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

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

Referenced by 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

◆ m_dGL_xi2

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

Definition at line 64 of file DiffusionLFR.h.

Referenced by DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and 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 SetupCFunctions().

◆ m_dGR_xi1

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

◆ m_dGR_xi2

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

Definition at line 65 of file DiffusionLFR.h.

Referenced by DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and 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 SetupCFunctions().

◆ m_diffType

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

Definition at line 91 of file DiffusionLFR.h.

Referenced by 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 SetupMetrics(), and v_Diffuse().

◆ 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 SetupMetrics(), and v_Diffuse().

◆ m_Q2D_e0

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

Definition at line 57 of file DiffusionLFR.h.

Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), and 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 DivCFlux_2D(), DivCFlux_2D_Gauss(), and 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 DivCFlux_2D(), DivCFlux_2D_Gauss(), and 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 DivCFlux_2D(), DivCFlux_2D_Gauss(), and 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

◆ type

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

Definition at line 52 of file DiffusionLFR.h.