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_SetBaseFlow (const Array< OneD, Array< OneD, NekDouble > > &inarray, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
 Overrides the base flow used during linearised advection. More...
 

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

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.

Referenced by create().

78  :m_advType(advType)
79  {
80  }

Member Function Documentation

◆ create()

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

Definition at line 49 of file AdvectionFR.h.

References AdvectionFR().

50  {
51  return AdvectionSharedPtr(new AdvectionFR(advType));
52  }
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.
Definition: AdvectionFR.cpp:78
std::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
Definition: Advection.h:170

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

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

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
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: AdvectionFR.h:57
RiemannSolverSharedPtr m_riemann
Riemann solver for DG-type schemes.
Definition: Advection.h:143
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:244
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:48
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".
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.
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:540
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.
AdvectionFluxVecCB m_fluxVector
Callback function to the flux vector (set when advection is in conservative form).
Definition: Advection.h:141
int m_spaceDim
Storage for space dimension. Used for homogeneous extension.
Definition: Advection.h:145
Array< OneD, NekDouble > m_jac
Definition: AdvectionFR.h:56
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:302
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:186

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

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

Referenced by v_Advect().

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::DisContField1D>(
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  NekDouble tmpFluxVertex = 0;
1060  Vmath::Vcopy(nLocalSolutionPts,
1061  &fluxX1[phys_offset], 1,
1062  &tmparrayX1[0], 1);
1063 
1064  fields[0]->GetExp(n)->GetVertexPhysVals(0, tmparrayX1,
1065  tmpFluxVertex);
1066 
1067  t_offset = fields[0]->GetTrace()
1068  ->GetPhys_Offset(elmtToTrace[n][0]->GetElmtId());
1069 
1070  if(negatedFluxNormal[2*n])
1071  {
1072  JumpL[n] = numericalFlux[t_offset] - tmpFluxVertex;
1073  }
1074  else
1075  {
1076  JumpL[n] = -numericalFlux[t_offset] - tmpFluxVertex;
1077  }
1078 
1079  fields[0]->GetExp(n)->GetVertexPhysVals(1, tmparrayX1,
1080  tmpFluxVertex);
1081 
1082  t_offset = fields[0]->GetTrace()
1083  ->GetPhys_Offset(elmtToTrace[n][1]->GetElmtId());
1084 
1085  if(negatedFluxNormal[2*n+1])
1086  {
1087  JumpR[n] = -numericalFlux[t_offset] - tmpFluxVertex;
1088  }
1089  else
1090  {
1091  JumpR[n] = numericalFlux[t_offset] - tmpFluxVertex;
1092  }
1093  }
1094 
1095  for (n = 0; n < nElements; ++n)
1096  {
1097  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1098  phys_offset = fields[0]->GetPhys_Offset(n);
1099 
1100  // Left jump multiplied by left derivative of C function
1101  Vmath::Smul(nLocalSolutionPts, JumpL[n], &m_dGL_xi1[n][0], 1,
1102  &DCL[0], 1);
1103 
1104  // Right jump multiplied by right derivative of C function
1105  Vmath::Smul(nLocalSolutionPts, JumpR[n], &m_dGR_xi1[n][0], 1,
1106  &DCR[0], 1);
1107 
1108  // Assembling divergence of the correction flux
1109  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1110  &divCFlux[phys_offset], 1);
1111  }
1112  }
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*y.
Definition: Vmath.cpp:216
double NekDouble
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
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:302

◆ 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 1126 of file AdvectionFR.cpp.

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

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

◆ 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 1317 of file AdvectionFR.cpp.

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

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

◆ 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 1689 of file AdvectionFR.cpp.

1697  {
1698  boost::ignore_unused(nConvectiveFields, fields, fluxX1, fluxX2,
1699  fluxX3, numericalFlux, divCFlux);
1700  }

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

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

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 SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Definition: Advection.cpp:98
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–7...

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

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

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  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: AdvectionFR.h:64
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition: AdvectionFR.h:68
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: AdvectionFR.h:65
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition: AdvectionFR.h:69
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: AdvectionFR.h:67
std::shared_ptr< Basis > BasisSharedPtr
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:1131
double NekDouble
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: AdvectionFR.h:66

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

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

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)->GetEdgeQFactors(
187  0, auxArray1 = m_Q2D_e0[n]);
188  pFields[0]->GetExp(n)->GetEdgeQFactors(
189  1, auxArray1 = m_Q2D_e1[n]);
190  pFields[0]->GetExp(n)->GetEdgeQFactors(
191  2, auxArray1 = m_Q2D_e2[n]);
192  pFields[0]->GetExp(n)->GetEdgeQFactors(
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  }
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: AdvectionFR.h:60
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: AdvectionFR.h:59
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: AdvectionFR.h:57
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: AdvectionFR.h:62
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: AdvectionFR.h:61
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: AdvectionFR.h:76
Geometry is curved or has non-constant factors.
Array< OneD, NekDouble > m_jac
Definition: AdvectionFR.h:56

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