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

#include <DiffusionLFR.h>

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

Static Public Member Functions

static DiffusionSharedPtr create (std::string diffType)
 

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...
 
void v_InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
 Initiliase DiffusionLFR objects and store them before starting the time-stepping. More...
 
void v_Diffuse (const 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 ~Diffusion ()=default
 
virtual SOLVER_UTILS_EXPORT void v_InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)=0
 
virtual SOLVER_UTILS_EXPORT 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)=0
 
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_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 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 ()
 

Protected Attributes

std::string m_diffType
 
- Protected Attributes inherited from Nektar::SolverUtils::Diffusion
Array< OneD, NekDouble > m_divVel
 Params for Ducros sensor. More...
 
Array< OneD, NekDouble > m_divVelSquare
 
Array< OneD, NekDouble > m_curlVelSquare
 
DiffusionFluxVecCB m_fluxVector
 
DiffusionFluxVecCBNS m_fluxVectorNS
 
DiffusionFluxPenaltyNS m_fluxPenaltyNS
 
DiffusionFluxCons m_FunctorDiffusionfluxCons
 
DiffusionFluxCons m_FunctorDiffusionfluxConsTrace
 
SpecialBndTreat m_SpecialBndTreat
 
DiffusionSymmFluxCons m_FunctorSymmetricfluxCons
 
Array< OneD, Array< OneD, NekDouble > > m_gridVelocityTrace
 
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...
 

Private Attributes

Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 
LibUtilities::SessionReaderSharedPtr m_session
 
Array< OneD, NekDouble > m_jac
 
Array< OneD, Array< OneD, NekDouble > > m_gmat
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
 
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
 
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
 
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
 
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
 
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
 
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
 
DNekMatSharedPtr m_Ixm
 
DNekMatSharedPtr m_Ixp
 
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
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::SolverUtils::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 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=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, 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 SetHomoDerivs (Array< OneD, Array< OneD, NekDouble > > &deriv)
 
SOLVER_UTILS_EXPORT TensorOfArray3D< NekDouble > & GetFluxTensor ()
 
SOLVER_UTILS_EXPORT const Array< OneD, const Array< OneD, NekDouble > > & GetTraceNormal ()
 Get trace normal. More...
 
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)
 
SOLVER_UTILS_EXPORT void SetGridVelocityTrace (Array< OneD, Array< OneD, NekDouble > > &gridVelocityTrace)
 

Detailed Description

Definition at line 42 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 67 of file DiffusionLFR.cpp.

67 : m_diffType(diffType)
68{
69}

Referenced by create().

Member Function Documentation

◆ create()

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

Definition at line 45 of file DiffusionLFR.h.

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

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

1444{
1445 int n;
1446 int nLocalSolutionPts, phys_offset, t_offset;
1447
1448 Array<OneD, NekDouble> auxArray1, auxArray2;
1449 Array<TwoD, const NekDouble> gmat;
1450 Array<OneD, const NekDouble> jac;
1451
1452 LibUtilities::PointsKeyVector ptsKeys;
1453 LibUtilities::BasisSharedPtr Basis;
1454 Basis = fields[0]->GetExp(0)->GetBasis(0);
1455
1456 int nElements = fields[0]->GetExpSize();
1457 int nSolutionPts = fields[0]->GetTotPoints();
1458
1459 vector<bool> leftAdjacentTraces =
1460 std::static_pointer_cast<MultiRegions::DisContField>(fields[0])
1461 ->GetLeftAdjacentTraces();
1462
1463 // Arrays to store the derivatives of the correction flux
1464 Array<OneD, NekDouble> DCL(nSolutionPts / nElements, 0.0);
1465 Array<OneD, NekDouble> DCR(nSolutionPts / nElements, 0.0);
1466
1467 // Arrays to store the intercell numerical flux jumps
1468 Array<OneD, NekDouble> JumpL(nElements);
1469 Array<OneD, NekDouble> JumpR(nElements);
1470
1471 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1472 fields[0]->GetTraceMap()->GetElmtToTrace();
1473
1474 for (n = 0; n < nElements; ++n)
1475 {
1476 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1477 phys_offset = fields[0]->GetPhys_Offset(n);
1478
1479 Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1480 Array<OneD, NekDouble> tmpFluxVertex(1, 0.0);
1481 Vmath::Vcopy(nLocalSolutionPts, &flux[phys_offset], 1, &tmparrayX1[0],
1482 1);
1483
1484 fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0], tmparrayX1,
1485 tmpFluxVertex);
1486
1487 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1488 elmtToTrace[n][0]->GetElmtId());
1489
1490 if (leftAdjacentTraces[2 * n])
1491 {
1492 JumpL[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1493 }
1494 else
1495 {
1496 JumpL[n] = iFlux[t_offset] - tmpFluxVertex[0];
1497 }
1498
1499 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1500 elmtToTrace[n][1]->GetElmtId());
1501
1502 fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1], tmparrayX1,
1503 tmpFluxVertex);
1504 if (leftAdjacentTraces[2 * n + 1])
1505 {
1506 JumpR[n] = iFlux[t_offset] - tmpFluxVertex[0];
1507 }
1508 else
1509 {
1510 JumpR[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1511 }
1512 }
1513
1514 for (n = 0; n < nElements; ++n)
1515 {
1516 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1517 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1518 phys_offset = fields[0]->GetPhys_Offset(n);
1519 jac = fields[0]
1520 ->GetExp(n)
1521 ->as<LocalRegions::Expansion1D>()
1522 ->GetGeom1D()
1523 ->GetMetricInfo()
1524 ->GetJac(ptsKeys);
1525
1526 JumpL[n] = JumpL[n] * jac[0];
1527 JumpR[n] = JumpR[n] * jac[0];
1528
1529 // Left jump multiplied by left derivative of C function
1530 Vmath::Smul(nLocalSolutionPts, JumpL[n], &m_dGL_xi1[n][0], 1, &DCL[0],
1531 1);
1532
1533 // Right jump multiplied by right derivative of C function
1534 Vmath::Smul(nLocalSolutionPts, JumpR[n], &m_dGR_xi1[n][0], 1, &DCR[0],
1535 1);
1536
1537 // Assembling derivative of the correction flux
1538 Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1539 &derCFlux[phys_offset], 1);
1540 }
1541}
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:81
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:80
Nektar::LibUtilities::BasisSharedPtr
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:69
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:231
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825

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

1563{
1564 int n, e, i, j, cnt;
1565
1566 Array<OneD, const NekDouble> jac;
1567
1568 int nElements = fields[0]->GetExpSize();
1569 int trace_offset, phys_offset;
1570 int nLocalSolutionPts;
1571 int nquad0, nquad1;
1572 int nEdgePts;
1573
1574 LibUtilities::PointsKeyVector ptsKeys;
1575 Array<OneD, NekDouble> auxArray1, auxArray2;
1576 Array<OneD, LibUtilities::BasisSharedPtr> base;
1577 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1578 fields[0]->GetTraceMap()->GetElmtToTrace();
1579
1580 // Loop on the elements
1581 for (n = 0; n < nElements; ++n)
1582 {
1583 // Offset of the element on the global vector
1584 phys_offset = fields[0]->GetPhys_Offset(n);
1585 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1586 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1587
1588 jac = fields[0]
1589 ->GetExp(n)
1590 ->as<LocalRegions::Expansion2D>()
1591 ->GetGeom2D()
1592 ->GetMetricInfo()
1593 ->GetJac(ptsKeys);
1594
1595 base = fields[0]->GetExp(n)->GetBase();
1596 nquad0 = base[0]->GetNumPoints();
1597 nquad1 = base[1]->GetNumPoints();
1598
1599 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1600 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1601 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1602 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1603
1604 // Loop on the edges
1605 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1606 {
1607 // Number of edge points of edge 'e'
1608 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1609
1610 // Array for storing volumetric fluxes on each edge
1611 Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1612
1613 // Offset of the trace space correspondent to edge 'e'
1614 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1615 elmtToTrace[n][e]->GetElmtId());
1616
1617 // Extract the edge values of the volumetric fluxes
1618 // on edge 'e' and order them accordingly to the order
1619 // of the trace space
1620 fields[0]->GetExp(n)->GetTracePhysVals(
1621 e, elmtToTrace[n][e], flux + phys_offset, auxArray1 = tmparray);
1622
1623 // Compute the fluxJumps per each edge 'e' and each
1624 // flux point
1625 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1626 Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1, &tmparray[0], 1,
1627 &fluxJumps[0], 1);
1628
1629 // Check the ordering of the fluxJumps and reverse
1630 // it in case of backward definition of edge 'e'
1631 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1632 StdRegions::eBackwards)
1633 {
1634 Vmath::Reverse(nEdgePts, &fluxJumps[0], 1, &fluxJumps[0], 1);
1635 }
1636
1637 // Deformed elements
1638 if (fields[0]
1639 ->GetExp(n)
1640 ->as<LocalRegions::Expansion2D>()
1641 ->GetGeom2D()
1642 ->GetMetricInfo()
1643 ->GetGtype() == SpatialDomains::eDeformed)
1644 {
1645 // Extract the Jacobians along edge 'e'
1646 Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1647
1648 fields[0]->GetExp(n)->GetTracePhysVals(
1649 e, elmtToTrace[n][e], jac, auxArray1 = jacEdge);
1650
1651 // Check the ordering of the fluxJumps and reverse
1652 // it in case of backward definition of edge 'e'
1653 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1654 StdRegions::eBackwards)
1655 {
1656 Vmath::Reverse(nEdgePts, &jacEdge[0], 1, &jacEdge[0], 1);
1657 }
1658
1659 // Multiply the fluxJumps by the edge 'e' Jacobians
1660 // to bring the fluxJumps into the standard space
1661 for (j = 0; j < nEdgePts; j++)
1662 {
1663 fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1664 }
1665 }
1666 // Non-deformed elements
1667 else
1668 {
1669 // Multiply the fluxJumps by the edge 'e' Jacobians
1670 // to bring the fluxJumps into the standard space
1671 Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1, fluxJumps, 1);
1672 }
1673
1674 // Multiply jumps by derivatives of the correction functions
1675 // All the quntities at this level should be defined into
1676 // the standard space
1677 switch (e)
1678 {
1679 case 0:
1680 for (i = 0; i < nquad0; ++i)
1681 {
1682 for (j = 0; j < nquad1; ++j)
1683 {
1684 cnt = i + j * nquad0;
1685 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1686 }
1687 }
1688 break;
1689 case 1:
1690 for (i = 0; i < nquad1; ++i)
1691 {
1692 for (j = 0; j < nquad0; ++j)
1693 {
1694 cnt = (nquad0)*i + j;
1695 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1696 }
1697 }
1698 break;
1699 case 2:
1700 for (i = 0; i < nquad0; ++i)
1701 {
1702 for (j = 0; j < nquad1; ++j)
1703 {
1704 cnt = j * nquad0 + i;
1705 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1706 }
1707 }
1708 break;
1709 case 3:
1710 for (i = 0; i < nquad1; ++i)
1711 {
1712 for (j = 0; j < nquad0; ++j)
1713 {
1714 cnt = j + i * nquad0;
1715 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1716 }
1717 }
1718 break;
1719
1720 default:
1721 ASSERTL0(false, "edge value (< 3) is out of range");
1722 break;
1723 }
1724 }
1725
1726 // Summing the component of the flux for calculating the
1727 // derivatives per each direction
1728 if (direction == 0)
1729 {
1730 Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1, &divCFluxE3[0], 1,
1731 &derCFlux[phys_offset], 1);
1732 }
1733 else if (direction == 1)
1734 {
1735 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE2[0], 1,
1736 &derCFlux[phys_offset], 1);
1737 }
1738 }
1739}
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFR.h:82
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFR.h:83
Nektar::SpatialDomains::eDeformed
@ eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:54
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:844
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.hpp:220

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

1761{
1762 int n, e, i, j, cnt;
1763
1764 int nElements = fields[0]->GetExpSize();
1765
1766 int nLocalSolutionPts;
1767 int nEdgePts;
1768 int trace_offset;
1769 int phys_offset;
1770 int nquad0;
1771 int nquad1;
1772
1773 Array<OneD, NekDouble> auxArray1, auxArray2;
1774 Array<OneD, LibUtilities::BasisSharedPtr> base;
1775
1776 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1777 fields[0]->GetTraceMap()->GetElmtToTrace();
1778
1779 // Loop on the elements
1780 for (n = 0; n < nElements; ++n)
1781 {
1782 // Offset of the element on the global vector
1783 phys_offset = fields[0]->GetPhys_Offset(n);
1784 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1785
1786 base = fields[0]->GetExp(n)->GetBase();
1787 nquad0 = base[0]->GetNumPoints();
1788 nquad1 = base[1]->GetNumPoints();
1789
1790 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1791 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1792 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1793 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1794
1795 // Loop on the edges
1796 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1797 {
1798 // Number of edge points of edge e
1799 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1800
1801 Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1802 Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1803 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1804 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1805 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1806
1807 // Offset of the trace space correspondent to edge e
1808 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1809 elmtToTrace[n][e]->GetElmtId());
1810
1811 // Get the normals of edge e
1812 const Array<OneD, const Array<OneD, NekDouble>> &normals =
1813 fields[0]->GetExp(n)->GetTraceNormal(e);
1814
1815 // Extract the edge values of flux-x on edge e and order
1816 // them accordingly to the order of the trace space
1817 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1818 fluxX1 + phys_offset,
1819 auxArray1 = tmparrayX1);
1820
1821 // Extract the edge values of flux-y on edge e and order
1822 // them accordingly to the order of the trace space
1823 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1824 fluxX2 + phys_offset,
1825 auxArray1 = tmparrayX2);
1826
1827 // Multiply the edge components of the flux by the normal
1828 for (i = 0; i < nEdgePts; ++i)
1829 {
1830 fluxN[i] = tmparrayX1[i] * m_traceNormals[0][trace_offset + i] +
1831 tmparrayX2[i] * m_traceNormals[1][trace_offset + i];
1832 }
1833
1834 // Subtract to the Riemann flux the discontinuous flux
1835 Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1, &fluxN[0], 1,
1836 &fluxJumps[0], 1);
1837
1838 // Check the ordering of the jump vectors
1839 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1840 StdRegions::eBackwards)
1841 {
1842 Vmath::Reverse(nEdgePts, auxArray1 = fluxJumps, 1,
1843 auxArray2 = fluxJumps, 1);
1844 }
1845
1846 for (i = 0; i < nEdgePts; ++i)
1847 {
1848 if (m_traceNormals[0][trace_offset + i] != normals[0][i] ||
1849 m_traceNormals[1][trace_offset + i] != normals[1][i])
1850 {
1851 fluxJumps[i] = -fluxJumps[i];
1852 }
1853 }
1854
1855 // Multiply jumps by derivatives of the correction functions
1856 switch (e)
1857 {
1858 case 0:
1859 for (i = 0; i < nquad0; ++i)
1860 {
1861 // Multiply fluxJumps by Q factors
1862 fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1863
1864 for (j = 0; j < nquad1; ++j)
1865 {
1866 cnt = i + j * nquad0;
1867 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1868 }
1869 }
1870 break;
1871 case 1:
1872 for (i = 0; i < nquad1; ++i)
1873 {
1874 // Multiply fluxJumps by Q factors
1875 fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1876
1877 for (j = 0; j < nquad0; ++j)
1878 {
1879 cnt = (nquad0)*i + j;
1880 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1881 }
1882 }
1883 break;
1884 case 2:
1885 for (i = 0; i < nquad0; ++i)
1886 {
1887 // Multiply fluxJumps by Q factors
1888 fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1889
1890 for (j = 0; j < nquad1; ++j)
1891 {
1892 cnt = j * nquad0 + i;
1893 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1894 }
1895 }
1896 break;
1897 case 3:
1898 for (i = 0; i < nquad1; ++i)
1899 {
1900 // Multiply fluxJumps by Q factors
1901 fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1902 for (j = 0; j < nquad0; ++j)
1903 {
1904 cnt = j + i * nquad0;
1905 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1906 }
1907 }
1908 break;
1909
1910 default:
1911 ASSERTL0(false, "edge value (< 3) is out of range");
1912 break;
1913 }
1914 }
1915
1916 // Sum each edge contribution
1917 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
1918 &divCFlux[phys_offset], 1);
1919
1920 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1921 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1922
1923 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1924 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1925 }
1926}
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:69
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFR.h:75
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFR.h:78
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFR.h:77
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFR.h:76

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

1949{
1950 int n, e, i, j, cnt;
1951
1952 int nElements = fields[0]->GetExpSize();
1953 int nLocalSolutionPts;
1954 int nEdgePts;
1955 int trace_offset;
1956 int phys_offset;
1957 int nquad0;
1958 int nquad1;
1959
1960 Array<OneD, NekDouble> auxArray1, auxArray2;
1961 Array<OneD, LibUtilities::BasisSharedPtr> base;
1962
1963 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1964 fields[0]->GetTraceMap()->GetElmtToTrace();
1965
1966 // Loop on the elements
1967 for (n = 0; n < nElements; ++n)
1968 {
1969 // Offset of the element on the global vector
1970 phys_offset = fields[0]->GetPhys_Offset(n);
1971 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1972
1973 base = fields[0]->GetExp(n)->GetBase();
1974 nquad0 = base[0]->GetNumPoints();
1975 nquad1 = base[1]->GetNumPoints();
1976
1977 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1978 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1979 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1980 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1981
1982 // Loop on the edges
1983 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1984 {
1985 // Number of edge points of edge e
1986 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1987
1988 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1989 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1990 Array<OneD, NekDouble> fluxN_R(nEdgePts, 0.0);
1991 Array<OneD, NekDouble> fluxN_D(nEdgePts, 0.0);
1992 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1993
1994 // Offset of the trace space correspondent to edge e
1995 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1996 elmtToTrace[n][e]->GetElmtId());
1997
1998 // Get the normals of edge e
1999 const Array<OneD, const Array<OneD, NekDouble>> &normals =
2000 fields[0]->GetExp(n)->GetTraceNormal(e);
2001
2002 // Extract the trasformed normal flux at each edge
2003 switch (e)
2004 {
2005 case 0:
2006 // Extract the edge values of transformed flux-y on
2007 // edge e and order them accordingly to the order of
2008 // the trace space
2009 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2010 fluxX2 + phys_offset,
2011 auxArray1 = fluxN_D);
2012
2013 Vmath::Neg(nEdgePts, fluxN_D, 1);
2014
2015 // Extract the physical Rieamann flux at each edge
2016 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2017 &fluxN[0], 1);
2018
2019 // Check the ordering of vectors
2020 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2021 StdRegions::eBackwards)
2022 {
2023 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2024 auxArray2 = fluxN, 1);
2025
2026 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2027 auxArray2 = fluxN_D, 1);
2028 }
2029
2030 // Transform Riemann Fluxes in the standard element
2031 for (i = 0; i < nquad0; ++i)
2032 {
2033 // Multiply Riemann Flux by Q factors
2034 fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
2035 }
2036
2037 for (i = 0; i < nEdgePts; ++i)
2038 {
2039 if (m_traceNormals[0][trace_offset + i] !=
2040 normals[0][i] ||
2041 m_traceNormals[1][trace_offset + i] !=
2042 normals[1][i])
2043 {
2044 fluxN_R[i] = -fluxN_R[i];
2045 }
2046 }
2047
2048 // Subtract to the Riemann flux the discontinuous
2049 // flux
2050 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2051 &fluxJumps[0], 1);
2052
2053 // Multiplicate the flux jumps for the correction
2054 // function
2055 for (i = 0; i < nquad0; ++i)
2056 {
2057 for (j = 0; j < nquad1; ++j)
2058 {
2059 cnt = i + j * nquad0;
2060 divCFluxE0[cnt] = -fluxJumps[i] * m_dGL_xi2[n][j];
2061 }
2062 }
2063 break;
2064 case 1:
2065 // Extract the edge values of transformed flux-x on
2066 // edge e and order them accordingly to the order of
2067 // the trace space
2068 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2069 fluxX1 + phys_offset,
2070 auxArray1 = fluxN_D);
2071
2072 // Extract the physical Rieamann flux at each edge
2073 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2074 &fluxN[0], 1);
2075
2076 // Check the ordering of vectors
2077 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2078 StdRegions::eBackwards)
2079 {
2080 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2081 auxArray2 = fluxN, 1);
2082
2083 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2084 auxArray2 = fluxN_D, 1);
2085 }
2086
2087 // Transform Riemann Fluxes in the standard element
2088 for (i = 0; i < nquad1; ++i)
2089 {
2090 // Multiply Riemann Flux by Q factors
2091 fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
2092 }
2093
2094 for (i = 0; i < nEdgePts; ++i)
2095 {
2096 if (m_traceNormals[0][trace_offset + i] !=
2097 normals[0][i] ||
2098 m_traceNormals[1][trace_offset + i] !=
2099 normals[1][i])
2100 {
2101 fluxN_R[i] = -fluxN_R[i];
2102 }
2103 }
2104
2105 // Subtract to the Riemann flux the discontinuous
2106 // flux
2107 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2108 &fluxJumps[0], 1);
2109
2110 // Multiplicate the flux jumps for the correction
2111 // function
2112 for (i = 0; i < nquad1; ++i)
2113 {
2114 for (j = 0; j < nquad0; ++j)
2115 {
2116 cnt = (nquad0)*i + j;
2117 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
2118 }
2119 }
2120 break;
2121 case 2:
2122
2123 // Extract the edge values of transformed flux-y on
2124 // edge e and order them accordingly to the order of
2125 // the trace space
2126
2127 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2128 fluxX2 + phys_offset,
2129 auxArray1 = fluxN_D);
2130
2131 // Extract the physical Rieamann flux at each edge
2132 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2133 &fluxN[0], 1);
2134
2135 // Check the ordering of vectors
2136 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2137 StdRegions::eBackwards)
2138 {
2139 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2140 auxArray2 = fluxN, 1);
2141
2142 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2143 auxArray2 = fluxN_D, 1);
2144 }
2145
2146 // Transform Riemann Fluxes in the standard element
2147 for (i = 0; i < nquad0; ++i)
2148 {
2149 // Multiply Riemann Flux by Q factors
2150 fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
2151 }
2152
2153 for (i = 0; i < nEdgePts; ++i)
2154 {
2155 if (m_traceNormals[0][trace_offset + i] !=
2156 normals[0][i] ||
2157 m_traceNormals[1][trace_offset + i] !=
2158 normals[1][i])
2159 {
2160 fluxN_R[i] = -fluxN_R[i];
2161 }
2162 }
2163
2164 // Subtract to the Riemann flux the discontinuous
2165 // flux
2166
2167 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2168 &fluxJumps[0], 1);
2169
2170 // Multiplicate the flux jumps for the correction
2171 // function
2172 for (i = 0; i < nquad0; ++i)
2173 {
2174 for (j = 0; j < nquad1; ++j)
2175 {
2176 cnt = j * nquad0 + i;
2177 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
2178 }
2179 }
2180 break;
2181 case 3:
2182 // Extract the edge values of transformed flux-x on
2183 // edge e and order them accordingly to the order of
2184 // the trace space
2185
2186 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2187 fluxX1 + phys_offset,
2188 auxArray1 = fluxN_D);
2189 Vmath::Neg(nEdgePts, fluxN_D, 1);
2190
2191 // Extract the physical Rieamann flux at each edge
2192 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2193 &fluxN[0], 1);
2194
2195 // Check the ordering of vectors
2196 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2197 StdRegions::eBackwards)
2198 {
2199 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2200 auxArray2 = fluxN, 1);
2201
2202 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2203 auxArray2 = fluxN_D, 1);
2204 }
2205
2206 // Transform Riemann Fluxes in the standard element
2207 for (i = 0; i < nquad1; ++i)
2208 {
2209 // Multiply Riemann Flux by Q factors
2210 fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
2211 }
2212
2213 for (i = 0; i < nEdgePts; ++i)
2214 {
2215 if (m_traceNormals[0][trace_offset + i] !=
2216 normals[0][i] ||
2217 m_traceNormals[1][trace_offset + i] !=
2218 normals[1][i])
2219 {
2220 fluxN_R[i] = -fluxN_R[i];
2221 }
2222 }
2223
2224 // Subtract to the Riemann flux the discontinuous
2225 // flux
2226
2227 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2228 &fluxJumps[0], 1);
2229
2230 // Multiplicate the flux jumps for the correction
2231 // function
2232 for (i = 0; i < nquad1; ++i)
2233 {
2234 for (j = 0; j < nquad0; ++j)
2235 {
2236 cnt = j + i * nquad0;
2237 divCFluxE3[cnt] = -fluxJumps[i] * m_dGL_xi1[n][j];
2238 }
2239 }
2240 break;
2241 default:
2242 ASSERTL0(false, "edge value (< 3) is out of range");
2243 break;
2244 }
2245 }
2246
2247 // Sum each edge contribution
2248 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
2249 &divCFlux[phys_offset], 1);
2250
2251 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2252 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2253
2254 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2255 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2256 }
2257}
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292

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

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

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

1305{
1306 int i, j;
1307 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1308 int nvariables = fields.size();
1309 int nDim = fields[0]->GetCoordim(0);
1310
1311 NekDouble C11 = 0.0;
1312 Array<OneD, NekDouble> Fwd(nTracePts);
1313 Array<OneD, NekDouble> Bwd(nTracePts);
1314 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
1315
1316 Array<OneD, NekDouble> qFwd(nTracePts);
1317 Array<OneD, NekDouble> qBwd(nTracePts);
1318 Array<OneD, NekDouble> qfluxtemp(nTracePts, 0.0);
1319 Array<OneD, NekDouble> uterm(nTracePts);
1320
1321 // Get the normal velocity Vn
1322 for (i = 0; i < nDim; ++i)
1323 {
1324 Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1, Vn, 1, Vn, 1);
1325 }
1326
1327 // Evaulate upwind flux:
1328 // qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
1329 for (i = 0; i < nvariables; ++i)
1330 {
1331 qflux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1332 for (j = 0; j < nDim; ++j)
1333 {
1334 // Compute Fwd and Bwd value of ufield of jth direction
1335 fields[i]->GetFwdBwdTracePhys(qfield[i][j], qFwd, qBwd);
1336
1337 // if Vn >= 0, flux = uFwd, i.e.,
1338 // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick
1339 // qflux = qBwd = q+
1340 // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick
1341 // qflux = qBwd = q-
1342
1343 // else if Vn < 0, flux = uBwd, i.e.,
1344 // edge::eForward, if V*n<0 <=> V*n_F<0, pick
1345 // qflux = qFwd = q-
1346 // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
1347 // qflux = qFwd = q+
1348
1349 fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1350
1351 Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
1352 qfluxtemp, 1);
1353
1354 // Imposing weak boundary condition with flux
1355 if (fields[0]->GetBndCondExpansions().size())
1356 {
1357 WeakPenaltyforVector(fields, i, j, qfield[i][j], qfluxtemp,
1358 C11);
1359 }
1360
1361 // q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
1362 // n_xi = n_x * tan_xi_x + n_y * tan_xi_y + n_z * tan_xi_z
1363 // n_xi = n_x * tan_eta_x + n_y * tan_eta_y + n_z*tan_eta_z
1364 Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1, qflux[i], 1);
1365 }
1366 }
1367}
Nektar::SolverUtils::DiffusionLFR::WeakPenaltyforVector
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.
Definition: DiffusionLFR.cpp:1373
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72

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

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

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

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

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 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 
)
overrideprotected

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

Definition at line 863 of file DiffusionLFR.cpp.

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

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.

Implements Nektar::SolverUtils::Diffusion.

Definition at line 81 of file DiffusionLFR.cpp.

84{
85 WARNINGL0(false, "LFR is deprecated, use LDG instead");
86
87 m_session = pSession;
88 SetupMetrics(pSession, pFields);
89 SetupCFunctions(pSession, pFields);
90
91 // Initialising arrays
92 int i, j;
93 int nConvectiveFields = pFields.size();
94 int nDim = pFields[0]->GetCoordim(0);
95 int nSolutionPts = pFields[0]->GetTotPoints();
96 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
97
98 m_IF1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
99 m_DU1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
100 m_DFC1 =
101 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
102 m_BD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
103 m_D1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
104 m_DD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
105 m_IF2 = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
106 m_DFC2 =
107 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
108 m_divFD = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
109 m_divFC = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
110
111 m_tmp1 =
112 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
113 m_tmp2 =
114 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
115
116 for (i = 0; i < nConvectiveFields; ++i)
117 {
118 m_IF1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
119 m_DU1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
120 m_DFC1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
121 m_BD1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
122 m_D1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
123 m_DD1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
124 m_IF2[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
125 m_DFC2[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
126 m_divFD[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
127 m_divFC[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
128
129 m_tmp1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
130 m_tmp2[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
131
132 for (j = 0; j < nDim; ++j)
133 {
134 m_IF1[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
135 m_DU1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
136 m_DFC1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
137 m_BD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
138 m_D1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
139 m_DD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
140 m_DFC2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
141
142 m_tmp1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
143 m_tmp2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
144 }
145 }
146}
WARNINGL0
#define WARNINGL0(condition, msg)
Definition: ErrorUtil.hpp:215
Nektar::SolverUtils::DiffusionLFR::SetupMetrics
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...
Definition: DiffusionLFR.cpp:164
Nektar::SolverUtils::DiffusionLFR::m_session
LibUtilities::SessionReaderSharedPtr m_session
Definition: DiffusionLFR.h:70
Nektar::SolverUtils::DiffusionLFR::SetupCFunctions
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...
Definition: DiffusionLFR.cpp:326

References m_BD1, m_D1, m_DD1, m_DFC1, m_DFC2, m_divFC, m_divFD, m_DU1, m_IF1, m_IF2, m_session, m_tmp1, m_tmp2, SetupCFunctions(), SetupMetrics(), and WARNINGL0.

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

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

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

1377{
1378 int i, e, id1, id2;
1379 int nBndEdges, nBndEdgePts;
1380 int nBndRegions = fields[var]->GetBndCondExpansions().size();
1381 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1382
1383 Array<OneD, NekDouble> uterm(nTracePts);
1384 Array<OneD, NekDouble> qtemp(nTracePts);
1385 int cnt = 0;
1386
1387 fields[var]->ExtractTracePhys(qfield, qtemp);
1388
1389 for (i = 0; i < nBndRegions; ++i)
1390 {
1391 nBndEdges = fields[var]->GetBndCondExpansions()[i]->GetExpSize();
1392
1393 // Weakly impose boundary conditions by modifying flux values
1394 for (e = 0; e < nBndEdges; ++e)
1395 {
1396 nBndEdgePts = fields[var]
1397 ->GetBndCondExpansions()[i]
1398 ->GetExp(e)
1399 ->GetTotPoints();
1400
1401 id1 = fields[var]->GetBndCondExpansions()[i]->GetPhys_Offset(e);
1402
1403 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1404 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1405
1406 // For Dirichlet boundary condition:
1407 // qflux = q+ - C_11 (u+ - g_D) (nx, ny)
1408 if (fields[var]
1409 ->GetBndConditions()[i]
1410 ->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1411 {
1412 Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1413 &qtemp[id2], 1, &penaltyflux[id2], 1);
1414 }
1415 // For Neumann boundary condition: qflux = g_N
1416 else if ((fields[var]->GetBndConditions()[i])
1417 ->GetBoundaryConditionType() ==
1418 SpatialDomains::eNeumann)
1419 {
1420 Vmath::Vmul(
1421 nBndEdgePts, &m_traceNormals[dir][id2], 1,
1422 &(fields[var]->GetBndCondExpansions()[i]->GetPhys())[id1],
1423 1, &penaltyflux[id2], 1);
1424 }
1425 }
1426 }
1427}

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
private

Definition at line 92 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
private

Definition at line 93 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
private

Definition at line 94 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
private

Definition at line 91 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
private

Definition at line 96 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
private

Definition at line 80 of file DiffusionLFR.h.

Referenced by DerCFlux_1D(), DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupCFunctions().

◆ m_dGL_xi2

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

Definition at line 82 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
private

Definition at line 84 of file DiffusionLFR.h.

Referenced by SetupCFunctions().

◆ m_dGR_xi1

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

Definition at line 81 of file DiffusionLFR.h.

Referenced by DerCFlux_1D(), DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupCFunctions().

◆ m_dGR_xi2

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

Definition at line 83 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
private

Definition at line 85 of file DiffusionLFR.h.

Referenced by SetupCFunctions().

◆ m_diffType

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

Definition at line 55 of file DiffusionLFR.h.

Referenced by SetupCFunctions().

◆ m_divFC

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

Definition at line 98 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_divFD

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

Definition at line 97 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
private

Definition at line 90 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_gmat

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

Definition at line 73 of file DiffusionLFR.h.

Referenced by SetupMetrics(), and v_Diffuse().

◆ m_IF1

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

Definition at line 89 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_IF2

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

Definition at line 95 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_Ixm

DNekMatSharedPtr Nektar::SolverUtils::DiffusionLFR::m_Ixm
private

Definition at line 86 of file DiffusionLFR.h.

◆ m_Ixp

DNekMatSharedPtr Nektar::SolverUtils::DiffusionLFR::m_Ixp
private

Definition at line 87 of file DiffusionLFR.h.

◆ m_jac

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

Definition at line 72 of file DiffusionLFR.h.

Referenced by SetupMetrics(), and v_Diffuse().

◆ m_Q2D_e0

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

Definition at line 75 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
private

Definition at line 76 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
private

Definition at line 77 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
private

Definition at line 78 of file DiffusionLFR.h.

Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupMetrics().

◆ m_session

LibUtilities::SessionReaderSharedPtr Nektar::SolverUtils::DiffusionLFR::m_session
private

Definition at line 70 of file DiffusionLFR.h.

Referenced by v_InitObject().

◆ m_tmp1

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

Definition at line 100 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
private

Definition at line 101 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_traceNormals

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

Definition at line 69 of file DiffusionLFR.h.

Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), NumFluxforScalar(), NumFluxforVector(), SetupMetrics(), and WeakPenaltyforVector().

◆ type

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

Definition at line 50 of file DiffusionLFR.h.

Generated on Sun Jun 15 2025 07:52:42 for Nektar++ by doxygen 1.9.4