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

#include <AdvectionFR.h>

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

Static Public Member Functions

static AdvectionSharedPtr create (std::string advType)
 

Public Attributes

Array< OneD, NekDoublem_jac
 
Array< OneD, Array< OneD, NekDouble > > m_gmat
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
 
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
 
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
 
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
 
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
 
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
 
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
 
DNekMatSharedPtr m_Ixm
 
DNekMatSharedPtr m_Ixp
 

Static Public Attributes

static std::string type []
 

Protected Member Functions

 AdvectionFR (std::string advType)
 AdvectionFR uses the Flux Reconstruction (FR) approach to compute the advection term. The implementation is only for segments, quadrilaterals and hexahedra at the moment. More...
 
virtual void v_InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Initiliase AdvectionFR objects and store them before starting the time-stepping. More...
 
virtual void v_Advect (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray)
 Compute the advection term at each time-step using the Flux Reconstruction approach (FR). More...
 
virtual void v_SetupMetrics (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform the fluxes at the interfaces of each element from the physical space to the standard space). More...
 
virtual void v_SetupCFunctions (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–72. More...
 
virtual void v_DivCFlux_1D (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
 Compute the divergence of the corrective flux for 1D problems. More...
 
virtual void v_DivCFlux_2D (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
 Compute the divergence of the corrective flux for 2D problems. More...
 
virtual void v_DivCFlux_2D_Gauss (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
 Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre". More...
 
virtual void v_DivCFlux_3D (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 > &fluxX3, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
 Compute the divergence of the corrective flux for 3D problems. More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::Advection
virtual SOLVER_UTILS_EXPORT void v_AdvectVolumeFlux (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &pVolumeFlux, const NekDouble &time)
 Advects Volume Flux. More...
 
virtual SOLVER_UTILS_EXPORT void v_AdvectTraceFlux (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble >> &pTraceFlux, const NekDouble &time, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray)
 Advects Trace Flux. More...
 
virtual SOLVER_UTILS_EXPORT void v_AdvectCoeffs (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray)
 
virtual SOLVER_UTILS_EXPORT void v_SetBaseFlow (const Array< OneD, Array< OneD, NekDouble > > &inarray, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
 Overrides the base flow used during linearised advection. More...
 
virtual SOLVER_UTILS_EXPORT void v_AddVolumJacToMat (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int &nConvectiveFields, const TensorOfArray5D< NekDouble > &ElmtJacArray, Array< OneD, Array< OneD, SNekBlkMatSharedPtr >> &gmtxarray)
 

Protected Attributes

Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 
std::string m_advType
 
- Protected Attributes inherited from Nektar::SolverUtils::Advection
AdvectionFluxVecCB m_fluxVector
 Callback function to the flux vector (set when advection is in conservative form). More...
 
RiemannSolverSharedPtr m_riemann
 Riemann solver for DG-type schemes. More...
 
int m_spaceDim
 Storage for space dimension. Used for homogeneous extension. More...
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::SolverUtils::Advection
virtual SOLVER_UTILS_EXPORT ~Advection ()
 
SOLVER_UTILS_EXPORT void InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Interface function to initialise the advection object. More...
 
SOLVER_UTILS_EXPORT void Advect (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray)
 Interface function to advect the vector field. More...
 
SOLVER_UTILS_EXPORT void AdvectVolumeFlux (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, Array< OneD, NekDouble >> &pAdvVel, const Array< OneD, Array< OneD, NekDouble >> &pInarray, TensorOfArray3D< NekDouble > &pVolumeFlux, const NekDouble &pTime)
 Interface function to advect the Volume field. More...
 
SOLVER_UTILS_EXPORT void AdvectTraceFlux (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, Array< OneD, NekDouble >> &pAdvVel, const Array< OneD, Array< OneD, NekDouble >> &pInarray, Array< OneD, Array< OneD, NekDouble >> &pTraceFlux, const NekDouble &pTime, const Array< OneD, Array< OneD, NekDouble >> &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble >> &pBwd=NullNekDoubleArrayOfArray)
 Interface function to advect the Trace field. More...
 
SOLVER_UTILS_EXPORT void AdvectCoeffs (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray)
 Similar with Advection::Advect(): calculate the advection flux The difference is in the outarray: it is the coefficients of basis for AdvectCoeffs() it is the physics on quadrature points for Advect() More...
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetFluxVector (FuncPointerT func, ObjectPointerT obj)
 Set the flux vector callback function. More...
 
void SetRiemannSolver (RiemannSolverSharedPtr riemann)
 Set a Riemann solver object for this advection object. More...
 
void SetFluxVector (AdvectionFluxVecCB fluxVector)
 Set the flux vector callback function. More...
 
void SetBaseFlow (const Array< OneD, Array< OneD, NekDouble > > &inarray, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
 Set the base flow used for linearised advection objects. More...
 
template<typename DataType , typename TypeNekBlkMatSharedPtr >
SOLVER_UTILS_EXPORT void AddTraceJacToMat (const int nConvectiveFields, const int nSpaceDim, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, TypeNekBlkMatSharedPtr > &TracePntJacCons, Array< OneD, Array< OneD, TypeNekBlkMatSharedPtr > > &gmtxarray, const Array< OneD, TypeNekBlkMatSharedPtr > &TracePntJacGrad, const Array< OneD, Array< OneD, DataType > > &TracePntJacGradSign)
 
template<typename DataType , typename TypeNekBlkMatSharedPtr >
void CalcJacobTraceInteg (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int m, const int n, const Array< OneD, const TypeNekBlkMatSharedPtr > &PntJac, const Array< OneD, const Array< OneD, DataType > > &PntJacSign, Array< OneD, DNekMatSharedPtr > &TraceJacFwd, Array< OneD, DNekMatSharedPtr > &TraceJacBwd)
 
SOLVER_UTILS_EXPORT void AddVolumJacToMat (const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int &nConvectiveFields, const TensorOfArray5D< NekDouble > &ElmtJacArray, Array< OneD, Array< OneD, SNekBlkMatSharedPtr >> &gmtxarray)
 
template<typename DataType , typename TypeNekBlkMatSharedPtr >
void AddTraceJacToMat (const int nConvectiveFields, const int nSpaceDim, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, TypeNekBlkMatSharedPtr > &TracePntJacCons, Array< OneD, Array< OneD, TypeNekBlkMatSharedPtr > > &gmtxarray, const Array< OneD, TypeNekBlkMatSharedPtr > &TracePntJacGrad, const Array< OneD, Array< OneD, DataType > > &TracePntJacGradSign)
 

Detailed Description

Definition at line 46 of file AdvectionFR.h.

Constructor & Destructor Documentation

◆ AdvectionFR()

Nektar::SolverUtils::AdvectionFR::AdvectionFR ( std::string  advType)
protected

AdvectionFR uses the Flux Reconstruction (FR) approach to compute the advection 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 78 of file AdvectionFR.cpp.

78  :m_advType(advType)
79  {
80  }

Referenced by create().

Member Function Documentation

◆ create()

static AdvectionSharedPtr Nektar::SolverUtils::AdvectionFR::create ( std::string  advType)
inlinestatic

Definition at line 49 of file AdvectionFR.h.

50  {
51  return AdvectionSharedPtr(new AdvectionFR(advType));
52  }
AdvectionFR(std::string advType)
AdvectionFR uses the Flux Reconstruction (FR) approach to compute the advection term....
Definition: AdvectionFR.cpp:78
std::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
Definition: Advection.h:291

References AdvectionFR().

◆ v_Advect()

void Nektar::SolverUtils::AdvectionFR::v_Advect ( const int  nConvectiveFields,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, Array< OneD, NekDouble > > &  advVel,
const Array< OneD, Array< OneD, NekDouble > > &  inarray,
Array< OneD, Array< OneD, NekDouble > > &  outarray,
const NekDouble time,
const Array< OneD, Array< OneD, NekDouble > > &  pFwd = NullNekDoubleArrayOfArray,
const Array< OneD, Array< OneD, NekDouble > > &  pBwd = NullNekDoubleArrayOfArray 
)
protectedvirtual

Compute the advection term at each time-step using the Flux Reconstruction approach (FR).

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
advVelAdvection velocities.
inarraySolution at the previous time-step.
outarrayAdvection term to be passed at the time integration class.

Implements Nektar::SolverUtils::Advection.

Definition at line 805 of file AdvectionFR.cpp.

814  {
815  boost::ignore_unused(advVel, time, pFwd, pBwd);
816 
817  int i, j, n;
818  int phys_offset;
819 
820  Array<OneD, NekDouble> auxArray1, auxArray2;
821 
822  Array<OneD, LibUtilities::BasisSharedPtr> Basis;
823  Basis = fields[0]->GetExp(0)->GetBase();
824 
825  int nElements = fields[0]->GetExpSize();
826  int nDimensions = fields[0]->GetCoordim(0);
827  int nSolutionPts = fields[0]->GetTotPoints();
828  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
829  int nCoeffs = fields[0]->GetNcoeffs();
830 
831  // Storage for flux vector.
832  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > fluxvector
833  (nConvectiveFields);
834  // Outarray for Galerkin projection in case of primitive dealising
835  Array<OneD, Array<OneD, NekDouble> > outarrayCoeff
836  (nConvectiveFields);
837  // Store forwards/backwards space along trace space
838  Array<OneD, Array<OneD, NekDouble> > Fwd (nConvectiveFields);
839  Array<OneD, Array<OneD, NekDouble> > Bwd (nConvectiveFields);
840  Array<OneD, Array<OneD, NekDouble> > numflux(nConvectiveFields);
841 
842  // Set up storage for flux vector.
843  for (i = 0; i < nConvectiveFields; ++i)
844  {
845  fluxvector[i] =
846  Array<OneD, Array<OneD, NekDouble> >(m_spaceDim);
847 
848  for (j = 0; j < m_spaceDim; ++j)
849  {
850  fluxvector[i][j] = Array<OneD, NekDouble>(nSolutionPts);
851  }
852  }
853 
854  for (i = 0; i < nConvectiveFields; ++i)
855  {
856  outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
857  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
858  Bwd[i] = Array<OneD, NekDouble>(nTracePts);
859  numflux[i] = Array<OneD, NekDouble>(nTracePts);
860  fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
861  }
862 
863  // Computing the interface flux at each trace point
864  m_riemann->Solve(m_spaceDim, Fwd, Bwd, numflux);
865 
866  // Divergence of the flux (computing the RHS)
867  switch(nDimensions)
868  {
869  // 1D problems
870  case 1:
871  {
872  Array<OneD, NekDouble> DfluxvectorX1(nSolutionPts);
873  Array<OneD, NekDouble> divFC (nSolutionPts);
874 
875  // Get the advection flux vector (solver specific)
876  m_fluxVector(inarray, fluxvector);
877 
878  // Get the discontinuous flux divergence
879  for(i = 0; i < nConvectiveFields; ++i)
880  {
881  for (n = 0; n < nElements; n++)
882  {
883  phys_offset = fields[0]->GetPhys_Offset(n);
884 
885  fields[i]->GetExp(n)->PhysDeriv(
886  0, fluxvector[i][0] + phys_offset,
887  auxArray2 = DfluxvectorX1 + phys_offset);
888  }
889 
890  // Get the correction flux divergence
891  v_DivCFlux_1D(nConvectiveFields, fields,
892  fluxvector[i][0], numflux[i], divFC);
893 
894  // Back to the physical space using global operations
895  Vmath::Vdiv(nSolutionPts, &divFC[0], 1, &m_jac[0], 1,
896  &outarray[i][0], 1);
897 
898  // Adding the total divergence to outarray (RHS)
899  Vmath::Vadd(nSolutionPts, &outarray[i][0], 1,
900  &DfluxvectorX1[0], 1, &outarray[i][0], 1);
901 
902  // Primitive Dealiasing 1D
903  if (!(Basis[0]->Collocation()))
904  {
905  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
906  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
907  }
908  }
909  break;
910  }
911  // 2D problems
912  case 2:
913  {
914  Array<OneD, NekDouble> DfluxvectorX1(nSolutionPts);
915  Array<OneD, NekDouble> DfluxvectorX2(nSolutionPts);
916  Array<OneD, NekDouble> divFD(nSolutionPts);
917  Array<OneD, NekDouble> divFC(nSolutionPts);
918 
919  // Get the advection flux vector (solver specific)
920  m_fluxVector(inarray, fluxvector);
921 
922  for (i = 0; i < nConvectiveFields; ++i)
923  {
924  // Temporary vectors
925  Array<OneD, NekDouble> f_hat(nSolutionPts);
926  Array<OneD, NekDouble> g_hat(nSolutionPts);
927 
928  Vmath::Vvtvvtp(nSolutionPts,
929  &m_gmat[0][0], 1,
930  &fluxvector[i][0][0], 1,
931  &m_gmat[2][0], 1,
932  &fluxvector[i][1][0], 1,
933  &f_hat[0], 1);
934 
935  Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &f_hat[0], 1,
936  &f_hat[0], 1);
937 
938  Vmath::Vvtvvtp(nSolutionPts,
939  &m_gmat[1][0], 1,
940  &fluxvector[i][0][0], 1,
941  &m_gmat[3][0], 1,
942  &fluxvector[i][1][0], 1,
943  &g_hat[0], 1);
944 
945  Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &g_hat[0], 1,
946  &g_hat[0], 1);
947 
948  // Get the discontinuous flux derivatives
949  for (n = 0; n < nElements; n++)
950  {
951  phys_offset = fields[0]->GetPhys_Offset(n);
952  fields[0]->GetExp(n)->StdPhysDeriv(0,
953  auxArray1 = f_hat + phys_offset,
954  auxArray2 = DfluxvectorX1 + phys_offset);
955  fields[0]->GetExp(n)->StdPhysDeriv(1,
956  auxArray1 = g_hat + phys_offset,
957  auxArray2 = DfluxvectorX2 + phys_offset);
958  }
959 
960  // Divergence of the discontinuous flux
961  Vmath::Vadd(nSolutionPts, DfluxvectorX1, 1,
962  DfluxvectorX2, 1, divFD, 1);
963 
964  // Divergence of the correction flux
965  if (Basis[0]->GetPointsType() ==
967  Basis[1]->GetPointsType() ==
969  {
971  nConvectiveFields, fields, f_hat, g_hat,
972  numflux[i], divFC);
973  }
974  else
975  {
977  nConvectiveFields, fields,
978  fluxvector[i][0], fluxvector[i][1],
979  numflux[i], divFC);
980 
981  }
982 
983  // Divergence of the final flux
984  Vmath::Vadd(nSolutionPts, divFD, 1, divFC, 1,
985  outarray[i], 1);
986 
987  // Back to the physical space using a global operation
988  Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1,
989  &m_jac[0], 1, &outarray[i][0], 1);
990 
991  // Primitive Dealiasing 2D
992  if (!(Basis[0]->Collocation()))
993  {
994  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
995  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
996  }
997  } // close nConvectiveFields loop
998  break;
999  }
1000  // 3D problems
1001  case 3:
1002  {
1003  ASSERTL0(false,"3D FRDG case not implemented yet");
1004  break;
1005  }
1006  }
1007  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
virtual void v_DivCFlux_2D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems.
virtual void v_DivCFlux_1D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 1D problems.
Array< OneD, NekDouble > m_jac
Definition: AdvectionFR.h:56
virtual void v_DivCFlux_2D_Gauss(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre".
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: AdvectionFR.h:57
int m_spaceDim
Storage for space dimension. Used for homogeneous extension.
Definition: Advection.h:225
AdvectionFluxVecCB m_fluxVector
Callback function to the flux vector (set when advection is in conservative form).
Definition: Advection.h:221
RiemannSolverSharedPtr m_riemann
Riemann solver for DG-type schemes.
Definition: Advection.h:223
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:48
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:192
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:322
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:257
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:625

References ASSERTL0, Nektar::LibUtilities::eGaussGaussLegendre, Nektar::SolverUtils::Advection::m_fluxVector, m_gmat, m_jac, Nektar::SolverUtils::Advection::m_riemann, Nektar::SolverUtils::Advection::m_spaceDim, v_DivCFlux_1D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), Vmath::Vadd(), Vmath::Vdiv(), Vmath::Vmul(), and Vmath::Vvtvvtp().

◆ v_DivCFlux_1D()

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

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

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxX1X1-volumetric flux in physical space.
numericalFluxInterface flux in physical space.
divCFluxDivergence of the corrective flux.

Definition at line 1019 of file AdvectionFR.cpp.

1025  {
1026  boost::ignore_unused(nConvectiveFields);
1027 
1028  int n;
1029  int nLocalSolutionPts, phys_offset, t_offset;
1030 
1032  Basis = fields[0]->GetExp(0)->GetBasis(0);
1033 
1034  int nElements = fields[0]->GetExpSize();
1035  int nSolutionPts = fields[0]->GetTotPoints();
1036 
1037 
1038  vector<bool> negatedFluxNormal =
1039  std::static_pointer_cast<MultiRegions::DisContField>(
1040  fields[0])->GetNegatedFluxNormal();
1041 
1042  // Arrays to store the derivatives of the correction flux
1043  Array<OneD, NekDouble> DCL(nSolutionPts/nElements, 0.0);
1044  Array<OneD, NekDouble> DCR(nSolutionPts/nElements, 0.0);
1045 
1046  // Arrays to store the intercell numerical flux jumps
1047  Array<OneD, NekDouble> JumpL(nElements);
1048  Array<OneD, NekDouble> JumpR(nElements);
1049 
1050  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1051  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1052 
1053  for (n = 0; n < nElements; ++n)
1054  {
1055  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1056  phys_offset = fields[0]->GetPhys_Offset(n);
1057 
1058  Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1059  Array<OneD, NekDouble> tmpFluxVertex(1,0.0);
1060  Vmath::Vcopy(nLocalSolutionPts,
1061  &fluxX1[phys_offset], 1,
1062  &tmparrayX1[0], 1);
1063 
1064  fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0],
1065  tmparrayX1,
1066  tmpFluxVertex);
1067 
1068  t_offset = fields[0]->GetTrace()
1069  ->GetPhys_Offset(elmtToTrace[n][0]->GetElmtId());
1070 
1071  if(negatedFluxNormal[2*n])
1072  {
1073  JumpL[n] = numericalFlux[t_offset] - tmpFluxVertex[0];
1074  }
1075  else
1076  {
1077  JumpL[n] = -numericalFlux[t_offset] - tmpFluxVertex[0];
1078  }
1079 
1080  fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1],
1081  tmparrayX1,
1082  tmpFluxVertex);
1083 
1084  t_offset = fields[0]->GetTrace()
1085  ->GetPhys_Offset(elmtToTrace[n][1]->GetElmtId());
1086 
1087  if(negatedFluxNormal[2*n+1])
1088  {
1089  JumpR[n] = -numericalFlux[t_offset] - tmpFluxVertex[0];
1090  }
1091  else
1092  {
1093  JumpR[n] = numericalFlux[t_offset] - tmpFluxVertex[0];
1094  }
1095  }
1096 
1097  for (n = 0; n < nElements; ++n)
1098  {
1099  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1100  phys_offset = fields[0]->GetPhys_Offset(n);
1101 
1102  // Left jump multiplied by left derivative of C function
1103  Vmath::Smul(nLocalSolutionPts, JumpL[n], &m_dGL_xi1[n][0], 1,
1104  &DCL[0], 1);
1105 
1106  // Right jump multiplied by right derivative of C function
1107  Vmath::Smul(nLocalSolutionPts, JumpR[n], &m_dGR_xi1[n][0], 1,
1108  &DCR[0], 1);
1109 
1110  // Assembling divergence of the correction flux
1111  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1112  &divCFlux[phys_offset], 1);
1113  }
1114  }
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: AdvectionFR.h:64
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: AdvectionFR.h:65
std::shared_ptr< Basis > BasisSharedPtr
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199

References m_dGL_xi1, m_dGR_xi1, Vmath::Smul(), Vmath::Vadd(), and Vmath::Vcopy().

Referenced by v_Advect().

◆ v_DivCFlux_2D()

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

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

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

Definition at line 1128 of file AdvectionFR.cpp.

1135  {
1136  boost::ignore_unused(nConvectiveFields);
1137 
1138  int n, e, i, j, cnt;
1139 
1140  int nElements = fields[0]->GetExpSize();
1141 
1142  int nLocalSolutionPts, nEdgePts;
1143  int trace_offset, phys_offset;
1144  int nquad0, nquad1;
1145 
1146  Array<OneD, NekDouble> auxArray1;
1147  Array<OneD, LibUtilities::BasisSharedPtr> base;
1148 
1149  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1150  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1151 
1152  // Loop on the elements
1153  for(n = 0; n < nElements; ++n)
1154  {
1155  // Offset of the element on the global vector
1156  phys_offset = fields[0]->GetPhys_Offset(n);
1157  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1158 
1159  base = fields[0]->GetExp(n)->GetBase();
1160  nquad0 = base[0]->GetNumPoints();
1161  nquad1 = base[1]->GetNumPoints();
1162 
1163  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1164  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1165  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1166  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1167 
1168  // Loop on the edges
1169  for(e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1170  {
1171  // Number of edge points of edge e
1172  nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1173 
1174  Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1175  Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1176  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
1177  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
1178  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
1179 
1180  // Offset of the trace space correspondent to edge e
1181  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1182  elmtToTrace[n][e]->GetElmtId());
1183 
1184  // Get the normals of edge e
1185  const Array<OneD, const Array<OneD, NekDouble> > &normals =
1186  fields[0]->GetExp(n)->GetTraceNormal(e);
1187 
1188  // Extract the edge values of flux-x on edge e and order
1189  // them accordingly to the order of the trace space
1190  fields[0]->GetExp(n)->GetTracePhysVals(
1191  e, elmtToTrace[n][e],
1192  fluxX1 + phys_offset,
1193  auxArray1 = tmparrayX1);
1194 
1195  // Extract the edge values of flux-y on edge e and order
1196  // them accordingly to the order of the trace space
1197  fields[0]->GetExp(n)->GetTracePhysVals(
1198  e, elmtToTrace[n][e],
1199  fluxX2 + phys_offset,
1200  auxArray1 = tmparrayX2);
1201 
1202  // Multiply the edge components of the flux by the normal
1203  Vmath::Vvtvvtp(nEdgePts, &tmparrayX1[0], 1,
1204  &m_traceNormals[0][trace_offset], 1,
1205  &tmparrayX2[0], 1,
1206  &m_traceNormals[1][trace_offset], 1,
1207  &fluxN[0], 1);
1208 
1209  // Subtract to the Riemann flux the discontinuous flux
1210  Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1,
1211  &fluxN[0], 1, &fluxJumps[0], 1);
1212 
1213  // Check the ordering of the jump vectors
1214  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1216  {
1217  Vmath::Reverse(nEdgePts, &fluxJumps[0], 1,
1218  &fluxJumps[0], 1);
1219  }
1220 
1221  for (i = 0; i < nEdgePts; ++i)
1222  {
1223  if (m_traceNormals[0][trace_offset+i] != normals[0][i]
1224  || m_traceNormals[1][trace_offset+i] != normals[1][i])
1225  {
1226  fluxJumps[i] = -fluxJumps[i];
1227  }
1228  }
1229 
1230  // Multiply jumps by derivatives of the correction functions
1231  switch (e)
1232  {
1233  case 0:
1234  for (i = 0; i < nquad0; ++i)
1235  {
1236  // Multiply fluxJumps by Q factors
1237  fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1238 
1239  for (j = 0; j < nquad1; ++j)
1240  {
1241  cnt = i + j*nquad0;
1242  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1243  }
1244  }
1245  break;
1246  case 1:
1247  for (i = 0; i < nquad1; ++i)
1248  {
1249  // Multiply fluxJumps by Q factors
1250  fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1251 
1252  for (j = 0; j < nquad0; ++j)
1253  {
1254  cnt = (nquad0)*i + j;
1255  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1256  }
1257  }
1258  break;
1259  case 2:
1260  for (i = 0; i < nquad0; ++i)
1261  {
1262  // Multiply fluxJumps by Q factors
1263  fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1264 
1265  for (j = 0; j < nquad1; ++j)
1266  {
1267  cnt = j*nquad0 + i;
1268  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1269  }
1270  }
1271  break;
1272  case 3:
1273  for (i = 0; i < nquad1; ++i)
1274  {
1275  // Multiply fluxJumps by Q factors
1276  fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1277  for (j = 0; j < nquad0; ++j)
1278  {
1279  cnt = j + i*nquad0;
1280  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1281  }
1282  }
1283  break;
1284  default:
1285  ASSERTL0(false,"edge value (< 3) is out of range");
1286  break;
1287  }
1288  }
1289 
1290  // Sum each edge contribution
1291  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1292  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1293 
1294  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1295  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1296 
1297  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1298  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1299  }
1300  }
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: AdvectionFR.h:67
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: AdvectionFR.h:59
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: AdvectionFR.h:62
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: AdvectionFR.h:60
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: AdvectionFR.h:66
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: AdvectionFR.h:76
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: AdvectionFR.h:61
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1226
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:372

References ASSERTL0, Nektar::StdRegions::eBackwards, 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(), Vmath::Vsub(), and Vmath::Vvtvvtp().

Referenced by v_Advect().

◆ v_DivCFlux_2D_Gauss()

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

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

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

Definition at line 1316 of file AdvectionFR.cpp.

1323  {
1324  boost::ignore_unused(nConvectiveFields);
1325 
1326  int n, e, i, j, cnt;
1327 
1328  int nElements = fields[0]->GetExpSize();
1329  int nLocalSolutionPts;
1330  int nEdgePts;
1331  int trace_offset;
1332  int phys_offset;
1333  int nquad0;
1334  int nquad1;
1335 
1336  Array<OneD, NekDouble> auxArray1, auxArray2;
1337  Array<OneD, LibUtilities::BasisSharedPtr> base;
1338 
1339  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1340  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1341 
1342  // Loop on the elements
1343  for(n = 0; n < nElements; ++n)
1344  {
1345  // Offset of the element on the global vector
1346  phys_offset = fields[0]->GetPhys_Offset(n);
1347  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1348 
1349  base = fields[0]->GetExp(n)->GetBase();
1350  nquad0 = base[0]->GetNumPoints();
1351  nquad1 = base[1]->GetNumPoints();
1352 
1353  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1354  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1355  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1356  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1357 
1358  // Loop on the edges
1359  for(e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1360  {
1361  // Number of edge points of edge e
1362  nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1363 
1364  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
1365  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
1366  Array<OneD, NekDouble> fluxN_R (nEdgePts, 0.0);
1367  Array<OneD, NekDouble> fluxN_D (nEdgePts, 0.0);
1368  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
1369 
1370  // Offset of the trace space correspondent to edge e
1371  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1372  elmtToTrace[n][e]->GetElmtId());
1373 
1374  // Get the normals of edge e
1375  const Array<OneD, const Array<OneD, NekDouble> > &normals =
1376  fields[0]->GetExp(n)->GetTraceNormal(e);
1377 
1378  // Extract the trasformed normal flux at each edge
1379  switch (e)
1380  {
1381  case 0:
1382  // Extract the edge values of transformed flux-y on
1383  // edge e and order them accordingly to the order of
1384  // the trace space
1385  fields[0]->GetExp(n)->GetTracePhysVals(
1386  e, elmtToTrace[n][e],
1387  fluxX2 + phys_offset,
1388  auxArray1 = fluxN_D);
1389 
1390  Vmath::Neg (nEdgePts, fluxN_D, 1);
1391 
1392  // Extract the physical Rieamann flux at each edge
1393  Vmath::Vcopy(nEdgePts,
1394  &numericalFlux[trace_offset], 1,
1395  &fluxN[0], 1);
1396 
1397  // Check the ordering of vectors
1398  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1400  {
1401  Vmath::Reverse(nEdgePts,
1402  auxArray1 = fluxN, 1,
1403  auxArray2 = fluxN, 1);
1404 
1405  Vmath::Reverse(nEdgePts,
1406  auxArray1 = fluxN_D, 1,
1407  auxArray2 = fluxN_D, 1);
1408  }
1409 
1410  // Transform Riemann Fluxes in the standard element
1411  for (i = 0; i < nquad0; ++i)
1412  {
1413  // Multiply Riemann Flux by Q factors
1414  fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
1415  }
1416 
1417  for (i = 0; i < nEdgePts; ++i)
1418  {
1419  if (m_traceNormals[0][trace_offset+i]
1420  != normals[0][i] ||
1421  m_traceNormals[1][trace_offset+i]
1422  != normals[1][i])
1423  {
1424  fluxN_R[i] = -fluxN_R[i];
1425  }
1426  }
1427 
1428  // Subtract to the Riemann flux the discontinuous
1429  // flux
1430  Vmath::Vsub(nEdgePts,
1431  &fluxN_R[0], 1,
1432  &fluxN_D[0], 1, &fluxJumps[0], 1);
1433 
1434  // Multiplicate the flux jumps for the correction
1435  // function
1436  for (i = 0; i < nquad0; ++i)
1437  {
1438  for (j = 0; j < nquad1; ++j)
1439  {
1440  cnt = i + j*nquad0;
1441  divCFluxE0[cnt] =
1442  -fluxJumps[i] * m_dGL_xi2[n][j];
1443  }
1444  }
1445  break;
1446  case 1:
1447  // Extract the edge values of transformed flux-x on
1448  // edge e and order them accordingly to the order of
1449  // the trace space
1450  fields[0]->GetExp(n)->GetTracePhysVals(
1451  e, elmtToTrace[n][e],
1452  fluxX1 + phys_offset,
1453  auxArray1 = fluxN_D);
1454 
1455  // Extract the physical Rieamann flux at each edge
1456  Vmath::Vcopy(nEdgePts,
1457  &numericalFlux[trace_offset], 1,
1458  &fluxN[0], 1);
1459 
1460  // Check the ordering of vectors
1461  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1463  {
1464  Vmath::Reverse(nEdgePts,
1465  auxArray1 = fluxN, 1,
1466  auxArray2 = fluxN, 1);
1467 
1468  Vmath::Reverse(nEdgePts,
1469  auxArray1 = fluxN_D, 1,
1470  auxArray2 = fluxN_D, 1);
1471  }
1472 
1473  // Transform Riemann Fluxes in the standard element
1474  for (i = 0; i < nquad1; ++i)
1475  {
1476  // Multiply Riemann Flux by Q factors
1477  fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
1478  }
1479 
1480  for (i = 0; i < nEdgePts; ++i)
1481  {
1482  if (m_traceNormals[0][trace_offset+i]
1483  != normals[0][i] ||
1484  m_traceNormals[1][trace_offset+i]
1485  != normals[1][i])
1486  {
1487  fluxN_R[i] = -fluxN_R[i];
1488  }
1489  }
1490 
1491  // Subtract to the Riemann flux the discontinuous
1492  // flux
1493  Vmath::Vsub(nEdgePts,
1494  &fluxN_R[0], 1,
1495  &fluxN_D[0], 1, &fluxJumps[0], 1);
1496 
1497  // Multiplicate the flux jumps for the correction
1498  // function
1499  for (i = 0; i < nquad1; ++i)
1500  {
1501  for (j = 0; j < nquad0; ++j)
1502  {
1503  cnt = (nquad0)*i + j;
1504  divCFluxE1[cnt] =
1505  fluxJumps[i] * m_dGR_xi1[n][j];
1506  }
1507  }
1508  break;
1509  case 2:
1510 
1511  // Extract the edge values of transformed flux-y on
1512  // edge e and order them accordingly to the order of
1513  // the trace space
1514 
1515  fields[0]->GetExp(n)->GetTracePhysVals(
1516  e, elmtToTrace[n][e],
1517  fluxX2 + phys_offset,
1518  auxArray1 = fluxN_D);
1519 
1520  // Extract the physical Rieamann flux at each edge
1521  Vmath::Vcopy(nEdgePts,
1522  &numericalFlux[trace_offset], 1,
1523  &fluxN[0], 1);
1524 
1525  // Check the ordering of vectors
1526  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1528  {
1529  Vmath::Reverse(nEdgePts,
1530  auxArray1 = fluxN, 1,
1531  auxArray2 = fluxN, 1);
1532 
1533  Vmath::Reverse(nEdgePts,
1534  auxArray1 = fluxN_D, 1,
1535  auxArray2 = fluxN_D, 1);
1536  }
1537 
1538  // Transform Riemann Fluxes in the standard element
1539  for (i = 0; i < nquad0; ++i)
1540  {
1541  // Multiply Riemann Flux by Q factors
1542  fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
1543  }
1544 
1545  for (i = 0; i < nEdgePts; ++i)
1546  {
1547  if (m_traceNormals[0][trace_offset+i]
1548  != normals[0][i] ||
1549  m_traceNormals[1][trace_offset+i]
1550  != normals[1][i])
1551  {
1552  fluxN_R[i] = -fluxN_R[i];
1553  }
1554  }
1555 
1556  // Subtract to the Riemann flux the discontinuous
1557  // flux
1558 
1559  Vmath::Vsub(nEdgePts,
1560  &fluxN_R[0], 1,
1561  &fluxN_D[0], 1, &fluxJumps[0], 1);
1562 
1563  // Multiplicate the flux jumps for the correction
1564  // function
1565  for (i = 0; i < nquad0; ++i)
1566  {
1567  for (j = 0; j < nquad1; ++j)
1568  {
1569  cnt = j*nquad0 + i;
1570  divCFluxE2[cnt] =
1571  fluxJumps[i] * m_dGR_xi2[n][j];
1572  }
1573  }
1574  break;
1575  case 3:
1576  // Extract the edge values of transformed flux-x on
1577  // edge e and order them accordingly to the order of
1578  // the trace space
1579 
1580  fields[0]->GetExp(n)->GetTracePhysVals(
1581  e, elmtToTrace[n][e],
1582  fluxX1 + phys_offset,
1583  auxArray1 = fluxN_D);
1584  Vmath::Neg (nEdgePts, fluxN_D, 1);
1585 
1586  // Extract the physical Rieamann flux at each edge
1587  Vmath::Vcopy(nEdgePts,
1588  &numericalFlux[trace_offset], 1,
1589  &fluxN[0], 1);
1590 
1591  // Check the ordering of vectors
1592  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1594  {
1595  Vmath::Reverse(nEdgePts,
1596  auxArray1 = fluxN, 1,
1597  auxArray2 = fluxN, 1);
1598 
1599  Vmath::Reverse(nEdgePts,
1600  auxArray1 = fluxN_D, 1,
1601  auxArray2 = fluxN_D, 1);
1602  }
1603 
1604  // Transform Riemann Fluxes in the standard element
1605  for (i = 0; i < nquad1; ++i)
1606  {
1607  // Multiply Riemann Flux by Q factors
1608  fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
1609  }
1610 
1611  for (i = 0; i < nEdgePts; ++i)
1612  {
1613  if (m_traceNormals[0][trace_offset+i]
1614  != normals[0][i] ||
1615  m_traceNormals[1][trace_offset+i]
1616  != normals[1][i])
1617  {
1618  fluxN_R[i] = -fluxN_R[i];
1619  }
1620  }
1621 
1622  // Subtract to the Riemann flux the discontinuous
1623  // flux
1624 
1625  Vmath::Vsub(nEdgePts,
1626  &fluxN_R[0], 1,
1627  &fluxN_D[0], 1, &fluxJumps[0], 1);
1628 
1629  // Multiplicate the flux jumps for the correction
1630  // function
1631  for (i = 0; i < nquad1; ++i)
1632  {
1633  for (j = 0; j < nquad0; ++j)
1634  {
1635  cnt = j + i*nquad0;
1636  divCFluxE3[cnt] =
1637  -fluxJumps[i] * m_dGL_xi1[n][j];
1638  }
1639  }
1640  break;
1641  default:
1642  ASSERTL0(false,"edge value (< 3) is out of range");
1643  break;
1644  }
1645  }
1646 
1647 
1648  // Sum each edge contribution
1649  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1650  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1651 
1652  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1653  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1654 
1655  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1656  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1657  }
1658  }
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:461

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

Referenced by v_Advect().

◆ v_DivCFlux_3D()

void Nektar::SolverUtils::AdvectionFR::v_DivCFlux_3D ( 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 > &  fluxX3,
const Array< OneD, const NekDouble > &  numericalFlux,
Array< OneD, NekDouble > &  divCFlux 
)
protectedvirtual

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

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

Definition at line 1674 of file AdvectionFR.cpp.

1682  {
1683  boost::ignore_unused(nConvectiveFields, fields, fluxX1, fluxX2,
1684  fluxX3, numericalFlux, divCFlux);
1685  }

◆ v_InitObject()

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

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

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

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.

Reimplemented from Nektar::SolverUtils::Advection.

Definition at line 92 of file AdvectionFR.cpp.

95  {
96  Advection::v_InitObject(pSession, pFields);
97  v_SetupMetrics (pSession, pFields);
98  v_SetupCFunctions (pSession, pFields);
99  }
virtual void v_SetupMetrics(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform...
virtual void v_SetupCFunctions(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–72...
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Definition: Advection.cpp:366

References Nektar::SolverUtils::Advection::v_InitObject(), v_SetupCFunctions(), and v_SetupMetrics().

◆ v_SetupCFunctions()

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

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

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

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.

Definition at line 272 of file AdvectionFR.cpp.

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

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

Referenced by v_InitObject().

◆ v_SetupMetrics()

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

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

This routine calls the function #GetEdgeQFactors 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 117 of file AdvectionFR.cpp.

120  {
121  boost::ignore_unused(pSession);
122  int i, n;
123  int nquad0, nquad1;
124  int phys_offset;
125  int nLocalSolutionPts;
126  int nElements = pFields[0]->GetExpSize();
127  int nDimensions = pFields[0]->GetCoordim(0);
128  int nSolutionPts = pFields[0]->GetTotPoints();
129  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
130 
131  Array<TwoD, const NekDouble> gmat;
132  Array<OneD, const NekDouble> jac;
133 
134  m_jac = Array<OneD, NekDouble>(nSolutionPts);
135 
136  Array<OneD, NekDouble> auxArray1;
137  Array<OneD, LibUtilities::BasisSharedPtr> base;
139 
140  switch (nDimensions)
141  {
142  case 1:
143  {
144  for (n = 0; n < nElements; ++n)
145  {
146  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
147  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
148  phys_offset = pFields[0]->GetPhys_Offset(n);
149  jac = pFields[0]->GetExp(n)
150  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
151  ->GetMetricInfo()->GetJac(ptsKeys);
152  for (i = 0; i < nLocalSolutionPts; ++i)
153  {
154  m_jac[i+phys_offset] = jac[0];
155  }
156  }
157  break;
158  }
159  case 2:
160  {
161  m_gmat = Array<OneD, Array<OneD, NekDouble> >(4);
162  m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
163  m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
164  m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
165  m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
166 
167  m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble> >(nElements);
168  m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
169  m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
170  m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
171 
173 
174  for (n = 0; n < nElements; ++n)
175  {
176  base = pFields[0]->GetExp(n)->GetBase();
177  nquad0 = base[0]->GetNumPoints();
178  nquad1 = base[1]->GetNumPoints();
179 
180  m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
181  m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
182  m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
183  m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
184 
185  // Extract the Q factors at each edge point
186  pFields[0]->GetExp(n)->GetTraceQFactors(
187  0, auxArray1 = m_Q2D_e0[n]);
188  pFields[0]->GetExp(n)->GetTraceQFactors(
189  1, auxArray1 = m_Q2D_e1[n]);
190  pFields[0]->GetExp(n)->GetTraceQFactors(
191  2, auxArray1 = m_Q2D_e2[n]);
192  pFields[0]->GetExp(n)->GetTraceQFactors(
193  3, auxArray1 = m_Q2D_e3[n]);
194 
195  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
196  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
197  phys_offset = pFields[0]->GetPhys_Offset(n);
198 
199  jac = pFields[0]->GetExp(n)
200  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
201  ->GetMetricInfo()->GetJac(ptsKeys);
202  gmat = pFields[0]->GetExp(n)
203  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
204  ->GetMetricInfo()->GetDerivFactors(ptsKeys);
205 
206  if (pFields[0]->GetExp(n)
207  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
208  ->GetMetricInfo()->GetGtype()
210  {
211  for (i = 0; i < nLocalSolutionPts; ++i)
212  {
213  m_jac[i+phys_offset] = jac[i];
214  m_gmat[0][i+phys_offset] = gmat[0][i];
215  m_gmat[1][i+phys_offset] = gmat[1][i];
216  m_gmat[2][i+phys_offset] = gmat[2][i];
217  m_gmat[3][i+phys_offset] = gmat[3][i];
218  }
219  }
220  else
221  {
222  for (i = 0; i < nLocalSolutionPts; ++i)
223  {
224  m_jac[i+phys_offset] = jac[0];
225  m_gmat[0][i+phys_offset] = gmat[0][0];
226  m_gmat[1][i+phys_offset] = gmat[1][0];
227  m_gmat[2][i+phys_offset] = gmat[2][0];
228  m_gmat[3][i+phys_offset] = gmat[3][0];
229  }
230  }
231  }
232 
233  m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(
234  nDimensions);
235  for(i = 0; i < nDimensions; ++i)
236  {
237  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
238  }
239  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
240 
241  break;
242  }
243  case 3:
244  {
245  ASSERTL0(false,"3DFR Metric terms not implemented yet");
246  break;
247  }
248  default:
249  {
250  ASSERTL0(false, "Expansion dimension not recognised");
251  break;
252  }
253  }
254  }
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
@ eDeformed
Geometry is curved or has non-constant factors.

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

Member Data Documentation

◆ m_advType

std::string Nektar::SolverUtils::AdvectionFR::m_advType
protected

Definition at line 78 of file AdvectionFR.h.

Referenced by v_SetupCFunctions().

◆ m_dGL_xi1

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

Definition at line 64 of file AdvectionFR.h.

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

◆ m_dGL_xi2

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

Definition at line 66 of file AdvectionFR.h.

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

◆ m_dGL_xi3

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

Definition at line 68 of file AdvectionFR.h.

Referenced by v_SetupCFunctions().

◆ m_dGR_xi1

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

Definition at line 65 of file AdvectionFR.h.

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

◆ m_dGR_xi2

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

Definition at line 67 of file AdvectionFR.h.

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

◆ m_dGR_xi3

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

Definition at line 69 of file AdvectionFR.h.

Referenced by v_SetupCFunctions().

◆ m_gmat

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

Definition at line 57 of file AdvectionFR.h.

Referenced by v_Advect(), and v_SetupMetrics().

◆ m_Ixm

DNekMatSharedPtr Nektar::SolverUtils::AdvectionFR::m_Ixm

Definition at line 70 of file AdvectionFR.h.

◆ m_Ixp

DNekMatSharedPtr Nektar::SolverUtils::AdvectionFR::m_Ixp

Definition at line 71 of file AdvectionFR.h.

◆ m_jac

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

Definition at line 56 of file AdvectionFR.h.

Referenced by v_Advect(), and v_SetupMetrics().

◆ m_Q2D_e0

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

Definition at line 59 of file AdvectionFR.h.

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

◆ m_Q2D_e1

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

Definition at line 60 of file AdvectionFR.h.

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

◆ m_Q2D_e2

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

Definition at line 61 of file AdvectionFR.h.

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

◆ m_Q2D_e3

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

Definition at line 62 of file AdvectionFR.h.

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

◆ m_traceNormals

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

Definition at line 76 of file AdvectionFR.h.

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

◆ type

std::string Nektar::SolverUtils::AdvectionFR::type
static
Initial value:
= {
"FRcmin", AdvectionFR::create),
"FRcinf", AdvectionFR::create)}
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:200
static AdvectionSharedPtr create(std::string advType)
Definition: AdvectionFR.h:49
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:47

Definition at line 54 of file AdvectionFR.h.