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...
 
void v_InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
 Initiliase AdvectionFR objects and store them before starting the time-stepping. More...
 
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 ~Advection ()=default
 
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
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 44 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 70 of file AdvectionFR.cpp.

70 : m_advType(advType)
71{
72}

Referenced by create().

Member Function Documentation

◆ create()

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

Definition at line 47 of file AdvectionFR.h.

48 {
49 return AdvectionSharedPtr(new AdvectionFR(advType));
50 }
Nektar::SolverUtils::AdvectionFR::AdvectionFR
AdvectionFR(std::string advType)
AdvectionFR uses the Flux Reconstruction (FR) approach to compute the advection term....
Definition: AdvectionFR.cpp:70
Nektar::SolverUtils::AdvectionSharedPtr
std::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
Definition: Advection.h:54

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

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

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

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

References ASSERTL0, Nektar::StdRegions::eBackwards, 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 1302 of file AdvectionFR.cpp.

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

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

Referenced by v_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 1632 of file AdvectionFR.cpp.

1640{
1641}

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

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

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

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

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

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

References ASSERTL0, 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 84 of file AdvectionFR.cpp.

87{
88 WARNINGL0(false, "FR is deprecated, use DG instead");
89
90 Advection::v_InitObject(pSession, pFields);
91 SetupMetrics(pSession, pFields);
92 SetupCFunctions(pSession, pFields);
93}
WARNINGL0
#define WARNINGL0(condition, msg)
Definition: ErrorUtil.hpp:215
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:111
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:275
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:264

References SetupCFunctions(), SetupMetrics(), Nektar::SolverUtils::Advection::v_InitObject(), and WARNINGL0.

Member Data Documentation

◆ m_advType

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

Definition at line 57 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 85 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 87 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 89 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 86 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 88 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 90 of file AdvectionFR.h.

Referenced by SetupCFunctions().

◆ m_gmat

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

Definition at line 78 of file AdvectionFR.h.

Referenced by SetupMetrics(), and v_Advect().

◆ m_Ixm

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

Definition at line 91 of file AdvectionFR.h.

◆ m_Ixp

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

Definition at line 92 of file AdvectionFR.h.

◆ m_jac

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

Definition at line 77 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 80 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 81 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 82 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 83 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 75 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: BasicUtils/NekFactory.hpp:197
Nektar::SolverUtils::AdvectionFR::create
static AdvectionSharedPtr create(std::string advType)
Definition: AdvectionFR.h:47
Nektar::SolverUtils::GetAdvectionFactory
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43

Definition at line 52 of file AdvectionFR.h.

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