Nektar++
Static Public Member Functions | Public Attributes | Static Public Attributes | Protected Member Functions | Protected Attributes | Private Member Functions | 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) override
 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) override
 Compute the advection term at each time-step using the Flux Reconstruction approach (FR). 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...
 

Private Member Functions

void SetupMetrics (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform the fluxes at the interfaces of each element from the physical space to the standard space). More...
 
void SetupCFunctions (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–72. More...
 
void 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...
 
void DivCFlux_2D (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
 Compute the divergence of the corrective flux for 2D problems. More...
 
void DivCFlux_2D_Gauss (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
 Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre". More...
 
void 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...
 

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

75  : m_advType(advType)
76 {
77 }

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:75
std::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
Definition: Advection.h:278

References AdvectionFR().

◆ DivCFlux_1D()

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

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

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

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

Referenced by v_Advect().

◆ DivCFlux_2D()

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

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

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

Definition at line 1135 of file AdvectionFR.cpp.

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

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

◆ DivCFlux_2D_Gauss()

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

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

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

Definition at line 1319 of file AdvectionFR.cpp.

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

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

Referenced by v_Advect().

◆ DivCFlux_3D()

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

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

1659 {
1660  boost::ignore_unused(nConvectiveFields, fields, fluxX1, fluxX2, fluxX3,
1661  numericalFlux, divCFlux);
1662 }

◆ SetupCFunctions()

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

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

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

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.

Definition at line 279 of file AdvectionFR.cpp.

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

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

◆ SetupMetrics()

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

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

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

117 {
118  boost::ignore_unused(pSession);
119  int i, n;
120  int nquad0, nquad1;
121  int phys_offset;
122  int nLocalSolutionPts;
123  int nElements = pFields[0]->GetExpSize();
124  int nDimensions = pFields[0]->GetCoordim(0);
125  int nSolutionPts = pFields[0]->GetTotPoints();
126  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
127 
128  Array<TwoD, const NekDouble> gmat;
129  Array<OneD, const NekDouble> jac;
130 
131  m_jac = Array<OneD, NekDouble>(nSolutionPts);
132 
133  Array<OneD, NekDouble> auxArray1;
134  Array<OneD, LibUtilities::BasisSharedPtr> base;
136 
137  switch (nDimensions)
138  {
139  case 1:
140  {
141  for (n = 0; n < nElements; ++n)
142  {
143  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
144  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
145  phys_offset = pFields[0]->GetPhys_Offset(n);
146  jac = pFields[0]
147  ->GetExp(n)
148  ->as<LocalRegions::Expansion1D>()
149  ->GetGeom1D()
150  ->GetMetricInfo()
151  ->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(0, auxArray1 =
187  m_Q2D_e0[n]);
188  pFields[0]->GetExp(n)->GetTraceQFactors(1, auxArray1 =
189  m_Q2D_e1[n]);
190  pFields[0]->GetExp(n)->GetTraceQFactors(2, auxArray1 =
191  m_Q2D_e2[n]);
192  pFields[0]->GetExp(n)->GetTraceQFactors(3, auxArray1 =
193  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]
200  ->GetExp(n)
201  ->as<LocalRegions::Expansion2D>()
202  ->GetGeom2D()
203  ->GetMetricInfo()
204  ->GetJac(ptsKeys);
205  gmat = pFields[0]
206  ->GetExp(n)
207  ->as<LocalRegions::Expansion2D>()
208  ->GetGeom2D()
209  ->GetMetricInfo()
210  ->GetDerivFactors(ptsKeys);
211 
212  if (pFields[0]
213  ->GetExp(n)
214  ->as<LocalRegions::Expansion2D>()
215  ->GetGeom2D()
216  ->GetMetricInfo()
217  ->GetGtype() == SpatialDomains::eDeformed)
218  {
219  for (i = 0; i < nLocalSolutionPts; ++i)
220  {
221  m_jac[i + phys_offset] = jac[i];
222  m_gmat[0][i + phys_offset] = gmat[0][i];
223  m_gmat[1][i + phys_offset] = gmat[1][i];
224  m_gmat[2][i + phys_offset] = gmat[2][i];
225  m_gmat[3][i + phys_offset] = gmat[3][i];
226  }
227  }
228  else
229  {
230  for (i = 0; i < nLocalSolutionPts; ++i)
231  {
232  m_jac[i + phys_offset] = jac[0];
233  m_gmat[0][i + phys_offset] = gmat[0][0];
234  m_gmat[1][i + phys_offset] = gmat[1][0];
235  m_gmat[2][i + phys_offset] = gmat[2][0];
236  m_gmat[3][i + phys_offset] = gmat[3][0];
237  }
238  }
239  }
240 
241  m_traceNormals = Array<OneD, Array<OneD, NekDouble>>(nDimensions);
242  for (i = 0; i < nDimensions; ++i)
243  {
244  m_traceNormals[i] = Array<OneD, NekDouble>(nTracePts);
245  }
246  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
247 
248  break;
249  }
250  case 3:
251  {
252  ASSERTL0(false, "3DFR Metric terms not implemented yet");
253  break;
254  }
255  default:
256  {
257  ASSERTL0(false, "Expansion dimension not recognised");
258  break;
259  }
260  }
261 }
Array< OneD, NekDouble > m_jac
Definition: AdvectionFR.h:56
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: AdvectionFR.h:57
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:250
@ 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().

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

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

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

References ASSERTL0, DivCFlux_1D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), Nektar::LibUtilities::eGaussGaussLegendre, Nektar::SolverUtils::Advection::m_fluxVector, m_gmat, m_jac, Nektar::SolverUtils::Advection::m_riemann, Nektar::SolverUtils::Advection::m_spaceDim, Vmath::Vadd(), Vmath::Vdiv(), Vmath::Vmul(), and Vmath::Vvtvvtp().

◆ v_InitObject()

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

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

This routine calls the functions SetupMetrics and SetupCFunctions to initialise the objects needed by AdvectionFR.

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.

Reimplemented from Nektar::SolverUtils::Advection.

Definition at line 89 of file AdvectionFR.cpp.

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

References SetupCFunctions(), SetupMetrics(), and Nektar::SolverUtils::Advection::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 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 DivCFlux_1D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and 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 DivCFlux_2D(), DivCFlux_2D_Gauss(), and 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 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 DivCFlux_1D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and 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 DivCFlux_2D(), DivCFlux_2D_Gauss(), and 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 SetupCFunctions().

◆ m_gmat

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

Definition at line 57 of file AdvectionFR.h.

Referenced by SetupMetrics(), and v_Advect().

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

◆ m_Q2D_e0

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

Definition at line 59 of file AdvectionFR.h.

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

◆ m_traceNormals

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

Definition at line 76 of file AdvectionFR.h.

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

◆ type

std::string Nektar::SolverUtils::AdvectionFR::type
static
Initial value:
= {
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
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.