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

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

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

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

1638{
1639}

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

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

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

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

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 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: 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 Wed Jan 31 2024 09:36:05 for Nektar++ by doxygen 1.9.4