39#include <boost/algorithm/string/predicate.hpp>
84 WARNINGL0(
false,
"LFRNS is deprecated, use LDGNS instead");
101 int nConvectiveFields = pFields.size();
102 int nScalars = nConvectiveFields - 1;
103 int nDim = pFields[0]->GetCoordim(0);
104 int nSolutionPts = pFields[0]->GetTotPoints();
105 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
108 if (pSession->DefinesSolverInfo(
"HOMOGENEOUS"))
138 for (i = 0; i < nScalars; ++i)
147 for (j = 0; j < nDim; ++j)
158 for (i = 0; i < nConvectiveFields; ++i)
166 for (j = 0; j < nDim; ++j)
177 for (j = 0; j < nScalars; ++j)
187 for (j = 0; j < nScalars + 1; ++j)
217 int nLocalSolutionPts;
218 int nElements = pFields[0]->GetExpSize();
219 int nDimensions = pFields[0]->GetCoordim(0);
220 int nSolutionPts = pFields[0]->GetTotPoints();
221 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
224 for (i = 0; i < nDimensions; ++i)
243 for (n = 0; n < nElements; ++n)
245 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
246 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
247 phys_offset = pFields[0]->GetPhys_Offset(n);
254 for (i = 0; i < nLocalSolutionPts; ++i)
256 m_jac[i + phys_offset] = jac[0];
274 for (n = 0; n < nElements; ++n)
276 base = pFields[0]->GetExp(n)->GetBase();
277 nquad0 = base[0]->GetNumPoints();
278 nquad1 = base[1]->GetNumPoints();
286 pFields[0]->GetExp(n)->GetTraceQFactors(0, auxArray1 =
288 pFields[0]->GetExp(n)->GetTraceQFactors(1, auxArray1 =
290 pFields[0]->GetExp(n)->GetTraceQFactors(2, auxArray1 =
292 pFields[0]->GetExp(n)->GetTraceQFactors(3, auxArray1 =
295 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
296 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
297 phys_offset = pFields[0]->GetPhys_Offset(n);
310 ->GetDerivFactors(ptsKeys);
319 for (i = 0; i < nLocalSolutionPts; ++i)
321 m_jac[i + phys_offset] = jac[i];
322 m_gmat[0][i + phys_offset] = gmat[0][i];
323 m_gmat[1][i + phys_offset] = gmat[1][i];
324 m_gmat[2][i + phys_offset] = gmat[2][i];
325 m_gmat[3][i + phys_offset] = gmat[3][i];
330 for (i = 0; i < nLocalSolutionPts; ++i)
332 m_jac[i + phys_offset] = jac[0];
333 m_gmat[0][i + phys_offset] = gmat[0][0];
334 m_gmat[1][i + phys_offset] = gmat[1][0];
335 m_gmat[2][i + phys_offset] = gmat[2][0];
336 m_gmat[3][i + phys_offset] = gmat[3][0];
344 ASSERTL0(
false,
"3DFR Metric terms not implemented yet");
349 ASSERTL0(
false,
"Expansion dimension not recognised");
379 int nquad0, nquad1, nquad2;
380 int nmodes0, nmodes1, nmodes2;
383 int nElements = pFields[0]->GetExpSize();
384 int nDim = pFields[0]->GetCoordim(0);
393 for (n = 0; n < nElements; ++n)
395 base = pFields[0]->GetExp(n)->GetBase();
396 nquad0 = base[0]->GetNumPoints();
397 nmodes0 = base[0]->GetNumModes();
401 base[0]->GetZW(z0, w0);
413 int p0 = nmodes0 - 1;
420 std::tgamma(2 * p0 + 1) /
421 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
431 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
432 (ap0 * std::tgamma(p0 + 1)) *
433 (ap0 * std::tgamma(p0 + 1)));
437 c0 = 2.0 * (p0 + 1.0) /
438 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
439 (ap0 * std::tgamma(p0 + 1)));
444 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
445 (ap0 * std::tgamma(p0 + 1)));
449 c0 = 10000000000000000.0;
452 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
453 (ap0 * std::tgamma(p0 + 1)) *
454 (ap0 * std::tgamma(p0 + 1));
456 NekDouble overeta0 = 1.0 / (1.0 + etap0);
469 for (i = 0; i < nquad0; ++i)
479 for (i = 0; i < nquad0; ++i)
502 for (n = 0; n < nElements; ++n)
504 base = pFields[0]->GetExp(n)->GetBase();
505 nquad0 = base[0]->GetNumPoints();
506 nquad1 = base[1]->GetNumPoints();
507 nmodes0 = base[0]->GetNumModes();
508 nmodes1 = base[1]->GetNumModes();
515 base[0]->GetZW(z0, w0);
516 base[1]->GetZW(z1, w1);
533 int p0 = nmodes0 - 1;
534 int p1 = nmodes1 - 1;
542 std::tgamma(2 * p0 + 1) /
543 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
546 std::tgamma(2 * p1 + 1) /
547 (pow(2.0, p1) * std::tgamma(p1 + 1) * std::tgamma(p1 + 1));
558 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
559 (ap0 * std::tgamma(p0 + 1)) *
560 (ap0 * std::tgamma(p0 + 1)));
563 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
564 (ap1 * std::tgamma(p1 + 1)) *
565 (ap1 * std::tgamma(p1 + 1)));
569 c0 = 2.0 * (p0 + 1.0) /
570 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
571 (ap0 * std::tgamma(p0 + 1)));
573 c1 = 2.0 * (p1 + 1.0) /
574 ((2.0 * p1 + 1.0) * p1 * (ap1 * std::tgamma(p1 + 1)) *
575 (ap1 * std::tgamma(p1 + 1)));
580 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
581 (ap0 * std::tgamma(p0 + 1)));
584 -2.0 / ((2.0 * p1 + 1.0) * (ap1 * std::tgamma(p1 + 1)) *
585 (ap1 * std::tgamma(p1 + 1)));
589 c0 = 10000000000000000.0;
590 c1 = 10000000000000000.0;
593 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
594 (ap0 * std::tgamma(p0 + 1)) *
595 (ap0 * std::tgamma(p0 + 1));
597 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
598 (ap1 * std::tgamma(p1 + 1)) *
599 (ap1 * std::tgamma(p1 + 1));
601 NekDouble overeta0 = 1.0 / (1.0 + etap0);
602 NekDouble overeta1 = 1.0 / (1.0 + etap1);
621 for (i = 0; i < nquad0; ++i)
631 for (i = 0; i < nquad1; ++i)
641 for (i = 0; i < nquad0; ++i)
651 for (i = 0; i < nquad1; ++i)
671 for (n = 0; n < nElements; ++n)
673 base = pFields[0]->GetExp(n)->GetBase();
674 nquad0 = base[0]->GetNumPoints();
675 nquad1 = base[1]->GetNumPoints();
676 nquad2 = base[2]->GetNumPoints();
677 nmodes0 = base[0]->GetNumModes();
678 nmodes1 = base[1]->GetNumModes();
679 nmodes2 = base[2]->GetNumModes();
688 base[0]->GetZW(z0, w0);
689 base[1]->GetZW(z1, w1);
690 base[1]->GetZW(z2, w2);
712 int p0 = nmodes0 - 1;
713 int p1 = nmodes1 - 1;
714 int p2 = nmodes2 - 1;
723 std::tgamma(2 * p0 + 1) /
724 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
728 std::tgamma(2 * p1 + 1) /
729 (pow(2.0, p1) * std::tgamma(p1 + 1) * std::tgamma(p1 + 1));
733 std::tgamma(2 * p2 + 1) /
734 (pow(2.0, p2) * std::tgamma(p2 + 1) * std::tgamma(p2 + 1));
746 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
747 (ap0 * std::tgamma(p0 + 1)) *
748 (ap0 * std::tgamma(p0 + 1)));
751 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
752 (ap1 * std::tgamma(p1 + 1)) *
753 (ap1 * std::tgamma(p1 + 1)));
756 ((2.0 * p2 + 1.0) * (p2 + 1.0) *
757 (ap2 * std::tgamma(p2 + 1)) *
758 (ap2 * std::tgamma(p2 + 1)));
762 c0 = 2.0 * (p0 + 1.0) /
763 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
764 (ap0 * std::tgamma(p0 + 1)));
766 c1 = 2.0 * (p1 + 1.0) /
767 ((2.0 * p1 + 1.0) * p1 * (ap1 * std::tgamma(p1 + 1)) *
768 (ap1 * std::tgamma(p1 + 1)));
770 c2 = 2.0 * (p2 + 1.0) /
771 ((2.0 * p2 + 1.0) * p2 * (ap2 * std::tgamma(p2 + 1)) *
772 (ap2 * std::tgamma(p2 + 1)));
777 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
778 (ap0 * std::tgamma(p0 + 1)));
781 -2.0 / ((2.0 * p1 + 1.0) * (ap1 * std::tgamma(p1 + 1)) *
782 (ap1 * std::tgamma(p1 + 1)));
785 -2.0 / ((2.0 * p2 + 1.0) * (ap2 * std::tgamma(p2 + 1)) *
786 (ap2 * std::tgamma(p2 + 1)));
790 c0 = 10000000000000000.0;
791 c1 = 10000000000000000.0;
792 c2 = 10000000000000000.0;
795 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
796 (ap0 * std::tgamma(p0 + 1)) *
797 (ap0 * std::tgamma(p0 + 1));
799 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
800 (ap1 * std::tgamma(p1 + 1)) *
801 (ap1 * std::tgamma(p1 + 1));
803 NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0) *
804 (ap2 * std::tgamma(p2 + 1)) *
805 (ap2 * std::tgamma(p2 + 1));
807 NekDouble overeta0 = 1.0 / (1.0 + etap0);
808 NekDouble overeta1 = 1.0 / (1.0 + etap1);
809 NekDouble overeta2 = 1.0 / (1.0 + etap2);
834 for (i = 0; i < nquad0; ++i)
844 for (i = 0; i < nquad1; ++i)
854 for (i = 0; i < nquad2; ++i)
864 for (i = 0; i < nquad0; ++i)
874 for (i = 0; i < nquad1; ++i)
884 for (i = 0; i < nquad2; ++i)
897 ASSERTL0(
false,
"Expansion dimension not recognised");
912 const std::size_t nConvectiveFields,
925 Basis = fields[0]->GetExp(0)->GetBase();
927 int nElements = fields[0]->GetExpSize();
928 int nDim = fields[0]->GetCoordim(0);
929 int nScalars = inarray.size();
930 int nSolutionPts = fields[0]->GetTotPoints();
931 int nCoeffs = fields[0]->GetNcoeffs();
934 for (i = 0; i < nConvectiveFields; ++i)
947 for (i = 0; i < nScalars; ++i)
951 for (n = 0; n < nElements; n++)
953 phys_offset = fields[0]->GetPhys_Offset(n);
955 fields[i]->GetExp(n)->PhysDeriv(
956 0, auxArray1 = inarray[i] + phys_offset,
957 auxArray2 =
m_DU1[i][0] + phys_offset);
973 1, &
m_D1[i][0][0], 1);
986 for (i = 0; i < nConvectiveFields; ++i)
990 for (n = 0; n < nElements; n++)
992 phys_offset = fields[0]->GetPhys_Offset(n);
994 fields[i]->GetExp(n)->PhysDeriv(
996 auxArray2 =
m_DD1[i][0] + phys_offset);
1012 1, &outarray[i][0], 1);
1015 if (!(
Basis[0]->Collocation()))
1017 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1018 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1026 for (i = 0; i < nScalars; ++i)
1028 for (j = 0; j < nDim; ++j)
1037 &
m_gmat[0][0], 1, &u1_hat[0], 1);
1043 &
m_gmat[1][0], 1, &u2_hat[0], 1);
1051 &
m_gmat[2][0], 1, &u1_hat[0], 1);
1057 &
m_gmat[3][0], 1, &u2_hat[0], 1);
1063 for (n = 0; n < nElements; n++)
1065 phys_offset = fields[0]->GetPhys_Offset(n);
1067 fields[i]->GetExp(n)->StdPhysDeriv(
1068 0, auxArray1 = u1_hat + phys_offset,
1069 auxArray2 =
m_tmp1[i][j] + phys_offset);
1071 fields[i]->GetExp(n)->StdPhysDeriv(
1072 1, auxArray1 = u2_hat + phys_offset,
1073 auxArray2 =
m_tmp2[i][j] + phys_offset);
1081 &
m_DU1[i][j][0], 1);
1085 DerCFlux_2D(nConvectiveFields, j, fields, inarray[i],
1091 for (j = 0; j < nSolutionPts; ++j)
1112 for (j = 0; j < nSolutionPts; j++)
1124 for (j = 0; j < nDim; ++j)
1134 for (i = 0; i < nScalars; ++i)
1149 for (i = 0; i < nConvectiveFields; ++i)
1155 for (j = 0; j < nSolutionPts; j++)
1166 for (n = 0; n < nElements; n++)
1168 phys_offset = fields[0]->GetPhys_Offset(n);
1170 fields[0]->GetExp(n)->StdPhysDeriv(
1171 0, auxArray1 = f_hat + phys_offset,
1172 auxArray2 =
m_DD1[i][0] + phys_offset);
1174 fields[0]->GetExp(n)->StdPhysDeriv(
1175 1, auxArray1 = g_hat + phys_offset,
1176 auxArray2 =
m_DD1[i][1] + phys_offset);
1184 if (
Basis[0]->GetPointsType() ==
1186 Basis[1]->GetPointsType() ==
1200 &outarray[i][0], 1);
1205 &outarray[i][0], 1);
1208 if (!(
Basis[0]->Collocation()))
1210 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1211 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1219 ASSERTL0(
false,
"3D FRDG case not implemented yet");
1235 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1236 int nScalars = inarray.size();
1237 int nDim = fields[0]->GetCoordim(0);
1243 for (i = 0; i < nDim; ++i)
1245 fields[0]->ExtractTracePhys(inarray[i],
m_traceVel[i]);
1255 for (i = 0; i < nScalars; ++i)
1260 fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
1261 fields[0]->GetTrace()->Upwind(Vn, Fwd[i], Bwd[i], numflux[i]);
1265 if (fields[0]->GetBndCondExpansions().size())
1271 for (j = 0; j < nDim; ++j)
1273 for (i = 0; i < nScalars; ++i)
1275 Vmath::Vcopy(nTracePts, numflux[i], 1, numericalFluxO1[i][j], 1);
1293 int nBndEdgePts, nBndEdges, nBndRegions;
1295 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1296 int nScalars = inarray.size();
1306 for (i = 0; i < nScalars; ++i)
1311 fields[i]->ExtractTracePhys(inarray[i], uplus[i]);
1315 for (i = 0; i < nScalars - 1; ++i)
1320 nBndRegions = fields[i + 1]->GetBndCondExpansions().size();
1321 for (j = 0; j < nBndRegions; ++j)
1324 ->GetBndConditions()[j]
1330 nBndEdges = fields[i + 1]->GetBndCondExpansions()[j]->GetExpSize();
1331 for (e = 0; e < nBndEdges; ++e)
1333 nBndEdgePts = fields[i + 1]
1334 ->GetBndCondExpansions()[j]
1339 fields[i + 1]->GetBndCondExpansions()[j]->GetPhys_Offset(e);
1341 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1342 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
1347 fields[i]->GetBndConditions()[j]->GetUserDefined(),
1350 fields[i]->GetBndConditions()[j]->GetUserDefined(),
1353 Vmath::Zero(nBndEdgePts, &scalarVariables[i][id2], 1);
1358 ->GetBndConditions()[j]
1359 ->GetBoundaryConditionType() ==
1364 ->GetBndCondExpansions()[j]
1365 ->UpdatePhys())[id1],
1368 ->GetBndCondExpansions()[j]
1369 ->UpdatePhys())[id1],
1370 1, &scalarVariables[i][id2], 1);
1375 ->GetBndConditions()[j]
1376 ->GetBoundaryConditionType() ==
1380 &penaltyfluxO1[i][id2], 1);
1384 else if ((fields[i]->GetBndConditions()[j])
1385 ->GetBoundaryConditionType() ==
1389 &penaltyfluxO1[i][id2], 1);
1393 Vmath::Vmul(nBndEdgePts, &scalarVariables[i][id2], 1,
1394 &scalarVariables[i][id2], 1, &tmp1[id2], 1);
1396 Vmath::Smul(nBndEdgePts, 0.5, &tmp1[id2], 1, &tmp1[id2], 1);
1398 Vmath::Vadd(nBndEdgePts, &tmp2[id2], 1, &tmp1[id2], 1,
1406 nBndRegions = fields[nScalars]->GetBndCondExpansions().size();
1407 for (j = 0; j < nBndRegions; ++j)
1409 nBndEdges = fields[nScalars]->GetBndCondExpansions()[j]->GetExpSize();
1411 if (fields[nScalars]
1412 ->GetBndConditions()[j]
1418 for (e = 0; e < nBndEdges; ++e)
1420 nBndEdgePts = fields[nScalars]
1421 ->GetBndCondExpansions()[j]
1426 fields[nScalars]->GetBndCondExpansions()[j]->GetPhys_Offset(e);
1428 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1429 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1433 fields[i]->GetBndConditions()[j]->GetUserDefined(),
1437 &scalarVariables[nScalars - 1][id2], 1);
1442 ->GetBndConditions()[j]
1443 ->GetBoundaryConditionType() ==
1450 ->GetBndCondExpansions()[j]
1452 1, &(fields[0]->GetBndCondExpansions()[j]->GetPhys())[id1],
1453 1, &scalarVariables[nScalars - 1][id2], 1);
1456 Vmath::Vsub(nBndEdgePts, &scalarVariables[nScalars - 1][id2], 1,
1457 &tmp2[id2], 1, &scalarVariables[nScalars - 1][id2],
1462 &scalarVariables[nScalars - 1][id2], 1,
1463 &scalarVariables[nScalars - 1][id2], 1);
1467 if (fields[nScalars]
1468 ->GetBndConditions()[j]
1469 ->GetBoundaryConditionType() ==
1472 fields[nScalars]->GetBndConditions()[j]->GetUserDefined(),
1475 Vmath::Vcopy(nBndEdgePts, &scalarVariables[nScalars - 1][id2],
1476 1, &penaltyfluxO1[nScalars - 1][id2], 1);
1480 else if (((fields[nScalars]->GetBndConditions()[j])
1481 ->GetBoundaryConditionType() ==
1483 boost::iequals(fields[nScalars]
1484 ->GetBndConditions()[j]
1488 Vmath::Vcopy(nBndEdgePts, &uplus[nScalars - 1][id2], 1,
1489 &penaltyfluxO1[nScalars - 1][id2], 1);
1506 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1507 int nVariables = fields.size();
1508 int nDim = fields[0]->GetCoordim(0);
1519 for (i = 0; i < nDim; ++i)
1521 fields[0]->ExtractTracePhys(ufield[i],
m_traceVel[i]);
1529 for (i = 1; i < nVariables; ++i)
1532 for (j = 0; j < nDim; ++j)
1535 fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
1538 fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1545 if (fields[0]->GetBndCondExpansions().size())
1551 Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1, qflux[i], 1);
1566 int nBndEdges, nBndEdgePts;
1570 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1571 int nBndRegions = fields[var]->GetBndCondExpansions().size();
1577 fields[var]->ExtractTracePhys(qfield, qtemp);
1580 for (i = 0; i < nBndRegions; ++i)
1583 nBndEdges = fields[var]->GetBndCondExpansions()[i]->GetExpSize();
1585 if (fields[var]->GetBndConditions()[i]->GetBoundaryConditionType() ==
1592 for (e = 0; e < nBndEdges; ++e)
1594 nBndEdgePts = fields[var]
1595 ->GetBndCondExpansions()[i]
1599 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1600 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1605 ->GetBndConditions()[i]
1606 ->GetBoundaryConditionType() ==
1609 fields[var]->GetBndConditions()[i]->GetUserDefined(),
1613 &qtemp[id2], 1, &penaltyflux[id2], 1);
1617 else if ((fields[var]->GetBndConditions()[i])
1618 ->GetBoundaryConditionType() ==
1621 ASSERTL0(
false,
"Neumann bcs not implemented for LFRNS");
1623 else if (boost::iequals(
1624 fields[var]->GetBndConditions()[i]->GetUserDefined(),
1635 &qtemp[id2], 1, &penaltyflux[id2], 1);
1653 [[maybe_unused]]
const int nConvectiveFields,
1659 int nLocalSolutionPts, phys_offset;
1667 Basis = fields[0]->GetExp(0)->GetBasis(0);
1669 int nElements = fields[0]->GetExpSize();
1670 int nPts = fields[0]->GetTotPoints();
1681 fields[0]->GetTraceMap()->GetElmtToTrace();
1683 for (n = 0; n < nElements; ++n)
1685 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1686 phys_offset = fields[0]->GetPhys_Offset(n);
1690 Vmath::Vcopy(nLocalSolutionPts, &flux[phys_offset], 1, &tmparrayX1[0],
1693 fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0], tmparrayX1,
1695 JumpL[n] = iFlux[n] - tmpFluxVertex[0];
1697 fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1], tmparrayX1,
1699 JumpR[n] = iFlux[n + 1] - tmpFluxVertex[0];
1702 for (n = 0; n < nElements; ++n)
1704 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1705 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1706 phys_offset = fields[0]->GetPhys_Offset(n);
1714 JumpL[n] = JumpL[n] * jac[0];
1715 JumpR[n] = JumpR[n] * jac[0];
1726 Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1727 &derCFlux[phys_offset], 1);
1747 [[maybe_unused]]
const int nConvectiveFields,
const int direction,
1752 int n, e, i, j, cnt;
1756 int nElements = fields[0]->GetExpSize();
1757 int trace_offset, phys_offset;
1758 int nLocalSolutionPts;
1766 fields[0]->GetTraceMap()->GetElmtToTrace();
1769 for (n = 0; n < nElements; ++n)
1772 phys_offset = fields[0]->GetPhys_Offset(n);
1773 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1774 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1783 base = fields[0]->GetExp(n)->GetBase();
1784 nquad0 = base[0]->GetNumPoints();
1785 nquad1 = base[1]->GetNumPoints();
1793 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1796 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1802 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1803 elmtToTrace[n][e]->GetElmtId());
1808 fields[0]->GetExp(n)->GetTracePhysVals(
1809 e, elmtToTrace[n][e], flux + phys_offset, auxArray1 = tmparray);
1814 Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1, &tmparray[0], 1,
1819 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1828 ->as<LocalRegions::Expansion2D>()
1836 fields[0]->GetExp(n)->GetTracePhysVals(
1837 e, elmtToTrace[n][e], jac, auxArray1 = jacEdge);
1841 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1849 for (j = 0; j < nEdgePts; j++)
1851 fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1859 Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1, fluxJumps, 1);
1868 for (i = 0; i < nquad0; ++i)
1870 for (j = 0; j < nquad1; ++j)
1872 cnt = i + j * nquad0;
1873 divCFluxE0[cnt] = fluxJumps[i] *
m_dGL_xi2[n][j];
1878 for (i = 0; i < nquad1; ++i)
1880 for (j = 0; j < nquad0; ++j)
1882 cnt = (nquad0)*i + j;
1883 divCFluxE1[cnt] = fluxJumps[i] *
m_dGR_xi1[n][j];
1888 for (i = 0; i < nquad0; ++i)
1890 for (j = 0; j < nquad1; ++j)
1892 cnt = j * nquad0 + i;
1893 divCFluxE2[cnt] = fluxJumps[i] *
m_dGR_xi2[n][j];
1898 for (i = 0; i < nquad1; ++i)
1900 for (j = 0; j < nquad0; ++j)
1902 cnt = j + i * nquad0;
1903 divCFluxE3[cnt] = fluxJumps[i] *
m_dGL_xi1[n][j];
1909 ASSERTL0(
false,
"edge value (< 3) is out of range");
1918 Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1, &divCFluxE3[0], 1,
1919 &derCFlux[phys_offset], 1);
1921 else if (direction == 1)
1923 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE2[0], 1,
1924 &derCFlux[phys_offset], 1);
1943 [[maybe_unused]]
const int nConvectiveFields,
1950 int n, e, i, j, cnt;
1952 int nElements = fields[0]->GetExpSize();
1954 int nLocalSolutionPts;
1965 fields[0]->GetTraceMap()->GetElmtToTrace();
1968 for (n = 0; n < nElements; ++n)
1971 phys_offset = fields[0]->GetPhys_Offset(n);
1972 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1974 base = fields[0]->GetExp(n)->GetBase();
1975 nquad0 = base[0]->GetNumPoints();
1976 nquad1 = base[1]->GetNumPoints();
1984 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1987 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1996 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1997 elmtToTrace[n][e]->GetElmtId());
2001 fields[0]->GetExp(n)->GetTraceNormal(e);
2005 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2006 fluxX1 + phys_offset,
2007 auxArray1 = tmparrayX1);
2011 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2012 fluxX2 + phys_offset,
2013 auxArray1 = tmparrayX2);
2016 for (i = 0; i < nEdgePts; ++i)
2023 Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1, &fluxN[0], 1,
2027 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2031 auxArray2 = fluxJumps, 1);
2034 for (i = 0; i < nEdgePts; ++i)
2039 fluxJumps[i] = -fluxJumps[i];
2047 for (i = 0; i < nquad0; ++i)
2050 fluxJumps[i] = -(
m_Q2D_e0[n][i]) * fluxJumps[i];
2052 for (j = 0; j < nquad1; ++j)
2054 cnt = i + j * nquad0;
2055 divCFluxE0[cnt] = fluxJumps[i] *
m_dGL_xi2[n][j];
2060 for (i = 0; i < nquad1; ++i)
2063 fluxJumps[i] = (
m_Q2D_e1[n][i]) * fluxJumps[i];
2065 for (j = 0; j < nquad0; ++j)
2067 cnt = (nquad0)*i + j;
2068 divCFluxE1[cnt] = fluxJumps[i] *
m_dGR_xi1[n][j];
2073 for (i = 0; i < nquad0; ++i)
2076 fluxJumps[i] = (
m_Q2D_e2[n][i]) * fluxJumps[i];
2078 for (j = 0; j < nquad1; ++j)
2080 cnt = j * nquad0 + i;
2081 divCFluxE2[cnt] = fluxJumps[i] *
m_dGR_xi2[n][j];
2086 for (i = 0; i < nquad1; ++i)
2089 fluxJumps[i] = -(
m_Q2D_e3[n][i]) * fluxJumps[i];
2090 for (j = 0; j < nquad0; ++j)
2092 cnt = j + i * nquad0;
2093 divCFluxE3[cnt] = fluxJumps[i] *
m_dGL_xi1[n][j];
2099 ASSERTL0(
false,
"edge value (< 3) is out of range");
2105 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
2106 &divCFlux[phys_offset], 1);
2108 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2109 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2111 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2112 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2131 [[maybe_unused]]
const int nConvectiveFields,
2138 int n, e, i, j, cnt;
2140 int nElements = fields[0]->GetExpSize();
2141 int nLocalSolutionPts;
2152 fields[0]->GetTraceMap()->GetElmtToTrace();
2155 for (n = 0; n < nElements; ++n)
2158 phys_offset = fields[0]->GetPhys_Offset(n);
2159 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
2161 base = fields[0]->GetExp(n)->GetBase();
2162 nquad0 = base[0]->GetNumPoints();
2163 nquad1 = base[1]->GetNumPoints();
2171 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
2174 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
2183 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
2184 elmtToTrace[n][e]->GetElmtId());
2188 fields[0]->GetExp(n)->GetTraceNormal(e);
2197 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2198 fluxX2 + phys_offset,
2199 auxArray1 = fluxN_D);
2204 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2208 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2212 auxArray2 = fluxN, 1);
2215 auxArray2 = fluxN_D, 1);
2219 for (i = 0; i < nquad0; ++i)
2222 fluxN_R[i] = (
m_Q2D_e0[n][i]) * fluxN[i];
2225 for (i = 0; i < nEdgePts; ++i)
2232 fluxN_R[i] = -fluxN_R[i];
2238 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2243 for (i = 0; i < nquad0; ++i)
2245 for (j = 0; j < nquad1; ++j)
2247 cnt = i + j * nquad0;
2248 divCFluxE0[cnt] = -fluxJumps[i] *
m_dGL_xi2[n][j];
2256 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2257 fluxX1 + phys_offset,
2258 auxArray1 = fluxN_D);
2261 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2265 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2269 auxArray2 = fluxN, 1);
2272 auxArray2 = fluxN_D, 1);
2276 for (i = 0; i < nquad1; ++i)
2279 fluxN_R[i] = (
m_Q2D_e1[n][i]) * fluxN[i];
2282 for (i = 0; i < nEdgePts; ++i)
2289 fluxN_R[i] = -fluxN_R[i];
2295 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2300 for (i = 0; i < nquad1; ++i)
2302 for (j = 0; j < nquad0; ++j)
2304 cnt = (nquad0)*i + j;
2305 divCFluxE1[cnt] = fluxJumps[i] *
m_dGR_xi1[n][j];
2315 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2316 fluxX2 + phys_offset,
2317 auxArray1 = fluxN_D);
2320 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2324 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2328 auxArray2 = fluxN, 1);
2331 auxArray2 = fluxN_D, 1);
2335 for (i = 0; i < nquad0; ++i)
2338 fluxN_R[i] = (
m_Q2D_e2[n][i]) * fluxN[i];
2341 for (i = 0; i < nEdgePts; ++i)
2348 fluxN_R[i] = -fluxN_R[i];
2355 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2360 for (i = 0; i < nquad0; ++i)
2362 for (j = 0; j < nquad1; ++j)
2364 cnt = j * nquad0 + i;
2365 divCFluxE2[cnt] = fluxJumps[i] *
m_dGR_xi2[n][j];
2374 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2375 fluxX1 + phys_offset,
2376 auxArray1 = fluxN_D);
2380 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2384 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2388 auxArray2 = fluxN, 1);
2391 auxArray2 = fluxN_D, 1);
2395 for (i = 0; i < nquad1; ++i)
2398 fluxN_R[i] = (
m_Q2D_e3[n][i]) * fluxN[i];
2401 for (i = 0; i < nEdgePts; ++i)
2408 fluxN_R[i] = -fluxN_R[i];
2415 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2420 for (i = 0; i < nquad1; ++i)
2422 for (j = 0; j < nquad0; ++j)
2424 cnt = j + i * nquad0;
2425 divCFluxE3[cnt] = -fluxJumps[i] *
m_dGL_xi1[n][j];
2430 ASSERTL0(
false,
"edge value (< 3) is out of range");
2436 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
2437 &divCFlux[phys_offset], 1);
2439 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2440 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2442 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2443 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
#define ASSERTL0(condition, msg)
#define WARNINGL0(condition, msg)
Represents a basis of a given type.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
const SpatialDomains::GeomFactorsSharedPtr & GetMetricInfo() const
DiffusionFluxVecCBNS m_fluxVectorNS
void v_Diffuse(const std::size_t nConvective, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd) override
Calculate FR Diffusion for the Navier-Stokes (NS) equations using an LDG interface flux.
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...
NekDouble m_thermalConductivity
static std::string type[]
Array< OneD, Array< OneD, NekDouble > > m_traceVel
Array< OneD, Array< OneD, NekDouble > > m_divFD
Array< OneD, NekDouble > m_jac
void DerCFlux_1D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, const NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux for 1D problems.
void NumericalFluxO1(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &numericalFluxO1)
Builds the numerical flux for the 1st order derivatives.
Array< OneD, Array< OneD, NekDouble > > m_viscFlux
void WeakPenaltyO1(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &penaltyfluxO1)
Imposes appropriate bcs for the 1st order derivatives.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_divFC
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Array< OneD, Array< OneD, NekDouble > > m_gmat
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_viscTensor
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...
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
void DivCFlux_2D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems.
void DerCFlux_2D(const int nConvectiveFields, const int direction, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux wrt a given coordinate for 2D problems.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
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".
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
DiffusionLFRNS(std::string diffType)
DiffusionLFRNS uses the Flux Reconstruction (FR) approach to compute the diffusion term....
LibUtilities::SessionReaderSharedPtr m_session
Array< OneD, Array< OneD, NekDouble > > m_homoDerivs
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
void NumericalFluxO2(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
Build the numerical flux for the 2nd order derivatives.
static DiffusionSharedPtr create(std::string diffType)
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
void WeakPenaltyO2(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux)
Imposes appropriate bcs for the 2nd order derivatives.
void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
Initiliase DiffusionLFRNS objects and store them before starting the time-stepping.
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
std::string m_ViscosityType
std::shared_ptr< Basis > BasisSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::vector< PointsKey > PointsKeyVector
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
DiffusionFactory & GetDiffusionFactory()
@ eDeformed
Geometry is curved or has non-constant factors.
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.
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.
void Neg(int n, T *x, const int incx)
Negate x = -x.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
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.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
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.
void Zero(int n, T *x, const int incx)
Zero vector.
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.