Nektar++
Static Public Member Functions | Static Public Attributes | Protected Member Functions | Protected Attributes | Private Member Functions | Private 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)
 

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_InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Initialises the advection object. More...
 
virtual SOLVER_UTILS_EXPORT 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)=0
 Advects a vector field. More...
 
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

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

Private Attributes

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

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...
 
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)
 
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 }
Nektar::SolverUtils::AdvectionFR::AdvectionFR
AdvectionFR(std::string advType)
AdvectionFR uses the Flux Reconstruction (FR) approach to compute the advection term....
Definition: AdvectionFR.cpp:75
Nektar::SolverUtils::AdvectionSharedPtr
std::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
Definition: Advection.h:56

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
1042 LibUtilities::BasisSharedPtr Basis;
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}
Nektar::SolverUtils::AdvectionFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: AdvectionFR.h:87
Nektar::SolverUtils::AdvectionFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: AdvectionFR.h:88
Nektar::LibUtilities::BasisSharedPtr
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:71
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:354
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:245
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1191

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) ==
1219 StdRegions::eBackwards)
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}
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
Nektar::SolverUtils::AdvectionFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: AdvectionFR.h:90
Nektar::SolverUtils::AdvectionFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: AdvectionFR.h:82
Nektar::SolverUtils::AdvectionFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: AdvectionFR.h:85
Nektar::SolverUtils::AdvectionFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: AdvectionFR.h:83
Nektar::SolverUtils::AdvectionFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: AdvectionFR.h:89
Nektar::SolverUtils::AdvectionFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: AdvectionFR.h:77
Nektar::SolverUtils::AdvectionFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: AdvectionFR.h:84
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1222
Vmath::Vvtvvtp
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:687
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:414

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) ==
1400 StdRegions::eBackwards)
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) ==
1457 StdRegions::eBackwards)
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) ==
1516 StdRegions::eBackwards)
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) ==
1576 StdRegions::eBackwards)
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}
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:513

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}
Nektar::SolverUtils::AdvectionFR::m_dGR_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition: AdvectionFR.h:92
Nektar::SolverUtils::AdvectionFR::m_dGL_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition: AdvectionFR.h:91
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Polylib::jacobd
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition: Polylib.cpp:1318

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;
135 LibUtilities::PointsKeyVector ptsKeys;
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
172 LibUtilities::PointsKeyVector ptsKeys;
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}
Nektar::SolverUtils::AdvectionFR::m_jac
Array< OneD, NekDouble > m_jac
Definition: AdvectionFR.h:79
Nektar::SolverUtils::AdvectionFR::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: AdvectionFR.h:80
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:236
Nektar::SpatialDomains::eDeformed
@ eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:56

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() ==
982 LibUtilities::eGaussGaussLegendre &&
983 Basis[1]->GetPointsType() ==
984 LibUtilities::eGaussGaussLegendre)
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}
Nektar::SolverUtils::AdvectionFR::DivCFlux_2D_Gauss
void DivCFlux_2D_Gauss(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre".
Definition: AdvectionFR.cpp:1319
Nektar::SolverUtils::AdvectionFR::DivCFlux_2D
void DivCFlux_2D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems.
Definition: AdvectionFR.cpp:1135
Nektar::SolverUtils::AdvectionFR::DivCFlux_1D
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.
Definition: AdvectionFR.cpp:1030
Nektar::SolverUtils::Advection::m_spaceDim
int m_spaceDim
Storage for space dimension. Used for homogeneous extension.
Definition: Advection.h:174
Nektar::SolverUtils::Advection::m_fluxVector
AdvectionFluxVecCB m_fluxVector
Callback function to the flux vector (set when advection is in conservative form).
Definition: Advection.h:170
Nektar::SolverUtils::Advection::m_riemann
RiemannSolverSharedPtr m_riemann
Riemann solver for DG-type schemes.
Definition: Advection.h:172
Nektar::LibUtilities::eGaussGaussLegendre
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:48
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:207
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:280

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}
Nektar::SolverUtils::AdvectionFR::SetupMetrics
void SetupMetrics(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform...
Definition: AdvectionFR.cpp:114
Nektar::SolverUtils::AdvectionFR::SetupCFunctions
void SetupCFunctions(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–72...
Definition: AdvectionFR.cpp:279
Nektar::SolverUtils::Advection::v_InitObject
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Definition: Advection.cpp:299

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 59 of file AdvectionFR.h.

Referenced by SetupCFunctions().

◆ m_dGL_xi1

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

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

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

Definition at line 91 of file AdvectionFR.h.

Referenced by SetupCFunctions().

◆ m_dGR_xi1

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

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

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

Definition at line 92 of file AdvectionFR.h.

Referenced by SetupCFunctions().

◆ m_gmat

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

Definition at line 80 of file AdvectionFR.h.

Referenced by SetupMetrics(), and v_Advect().

◆ m_Ixm

DNekMatSharedPtr Nektar::SolverUtils::AdvectionFR::m_Ixm
private

Definition at line 93 of file AdvectionFR.h.

◆ m_Ixp

DNekMatSharedPtr Nektar::SolverUtils::AdvectionFR::m_Ixp
private

Definition at line 94 of file AdvectionFR.h.

◆ m_jac

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

Definition at line 79 of file AdvectionFR.h.

Referenced by SetupMetrics(), and v_Advect().

◆ m_Q2D_e0

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

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

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

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

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

Definition at line 77 of file AdvectionFR.h.

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

◆ type

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

Definition at line 54 of file AdvectionFR.h.

Generated on Wed Jan 31 2024 08:26:58 for Nektar++ by doxygen 1.9.4