43 "NavierStokes equations in conservative variables.");
70 std::placeholders::_1, std::placeholders::_2,
71 std::placeholders::_3, std::placeholders::_4);
75 std::placeholders::_1, std::placeholders::_2,
76 std::placeholders::_3, std::placeholders::_4);
82 std::placeholders::_1, std::placeholders::_2,
83 std::placeholders::_3, std::placeholders::_4);
87 std::placeholders::_1, std::placeholders::_2,
88 std::placeholders::_3, std::placeholders::_4);
91 std::placeholders::_1, std::placeholders::_2,
92 std::placeholders::_3, std::placeholders::_4);
108 size_t nvariables = inarray.size();
114 for (
int i = 0; i < nvariables; ++i)
146 for (
int i = 0; i < nvariables; ++i)
148 Vmath::Vadd(ncoeffs, outarrayDiff[i], 1, outarray[i], 1,
154 ASSERTL1(
false,
"LDGNS not yet validated for implicit compressible "
160 for (
int i = 0; i < nvariables - 1; ++i)
168 m_varConv->GetTemperature(inarray, inarrayDiff[nvariables - 2]);
171 m_varConv->GetVelocityVector(inarray, inarrayDiff);
182 m_varConv->GetTemperature(pFwd, inFwd[nvariables - 2]);
183 m_varConv->GetTemperature(pBwd, inBwd[nvariables - 2]);
185 m_varConv->GetVelocityVector(pFwd, inFwd);
186 m_varConv->GetVelocityVector(pBwd, inBwd);
191 outarrayDiff, inFwd, inBwd);
193 for (
int i = 0; i < nvariables; ++i)
195 Vmath::Vadd(ncoeffs, outarrayDiff[i], 1, outarray[i], 1,
239 OneThird = 1.0 / 3.0;
240 TwoThird = 2.0 * OneThird;
241 FourThird = 4.0 * OneThird;
244 tmpArray = OutputMatrix->GetPtr();
245 int nrow = OutputMatrix->GetRows();
247 tmpArray[0 + 0 * nrow] = tmp * (-FourThird * u * nx - v * ny);
248 tmpArray[0 + 1 * nrow] = tmp * (FourThird * nx);
249 tmpArray[0 + 2 * nrow] = tmp * ny;
250 tmpArray[0 + 3 * nrow] = 0.0;
251 tmpArray[1 + 0 * nrow] = tmp * (-v * nx + TwoThird * u * ny);
252 tmpArray[1 + 1 * nrow] = tmp * (-TwoThird * ny);
253 tmpArray[1 + 2 * nrow] = tmp * nx;
254 tmpArray[1 + 3 * nrow] = 0.0;
255 tmpArray[2 + 0 * nrow] =
256 (FourThird * u * u + v * v + tmp2 * (E - q2)) * nx +
257 OneThird * u * v * ny;
258 tmpArray[2 + 0 * nrow] = -tmp * (*OutputMatrix)(2, 0);
259 tmpArray[2 + 1 * nrow] = (FourThird - tmp2) * u * nx - TwoThird * v * ny;
260 tmpArray[2 + 1 * nrow] = tmp * (*OutputMatrix)(2, 1);
261 tmpArray[2 + 2 * nrow] = (1 - tmp2) * v * nx + u * ny;
262 tmpArray[2 + 2 * nrow] = tmp * (*OutputMatrix)(2, 2);
263 tmpArray[2 + 3 * nrow] = tmp * tmp2 * nx;
297 OneThird = 1.0 / 3.0;
298 TwoThird = 2.0 * OneThird;
299 FourThird = 4.0 * OneThird;
302 tmpArray = OutputMatrix->GetPtr();
303 int nrow = OutputMatrix->GetRows();
305 tmpArray[0 + 0 * nrow] = tmp * (TwoThird * v * nx - u * ny);
306 tmpArray[0 + 1 * nrow] = tmp * ny;
307 tmpArray[0 + 2 * nrow] = tmp * (-TwoThird) * nx;
308 tmpArray[0 + 3 * nrow] = 0.0;
309 tmpArray[1 + 0 * nrow] = tmp * (-u * nx - FourThird * v * ny);
310 tmpArray[1 + 1 * nrow] = tmp * nx;
311 tmpArray[1 + 2 * nrow] = tmp * (FourThird * ny);
312 tmpArray[1 + 3 * nrow] = 0.0;
313 tmpArray[2 + 0 * nrow] = OneThird * u * v * nx +
314 (FourThird * v * v + u * u + tmp2 * (E - q2)) * ny;
315 tmpArray[2 + 0 * nrow] = -tmp * (*OutputMatrix)(2, 0);
316 tmpArray[2 + 1 * nrow] = (1 - tmp2) * u * ny + v * nx;
317 tmpArray[2 + 1 * nrow] = tmp * (*OutputMatrix)(2, 1);
318 tmpArray[2 + 2 * nrow] = (FourThird - tmp2) * v * ny - TwoThird * u * nx;
319 tmpArray[2 + 2 * nrow] = tmp * (*OutputMatrix)(2, 2);
320 tmpArray[2 + 3 * nrow] = tmp * tmp2 * ny;
362 OneThird = 1.0 / 3.0;
363 TwoThird = 2.0 * OneThird;
364 FourThird = 4.0 * OneThird;
367 tmpArray = OutputMatrix->GetPtr();
368 int nrow = OutputMatrix->GetRows();
370 tmpArray[0 + 0 * nrow] =
371 tmpx * (-FourThird * u) + tmpy * (-v) + tmpz * (-
w);
372 tmpArray[0 + 1 * nrow] = tmpx * FourThird;
373 tmpArray[0 + 2 * nrow] = tmpy;
374 tmpArray[0 + 3 * nrow] = tmpz;
375 tmpArray[0 + 4 * nrow] = 0.0;
376 tmpArray[1 + 0 * nrow] = tmpx * (-v) + tmpy * (TwoThird * u);
377 tmpArray[1 + 1 * nrow] = tmpy * (-TwoThird);
378 tmpArray[1 + 2 * nrow] = tmpx;
379 tmpArray[1 + 3 * nrow] = 0.0;
380 tmpArray[1 + 4 * nrow] = 0.0;
381 tmpArray[2 + 0 * nrow] = tmpx * (-
w) + tmpz * (TwoThird * u);
382 tmpArray[2 + 1 * nrow] = tmpz * (-TwoThird);
383 tmpArray[2 + 2 * nrow] = 0.0;
384 tmpArray[2 + 3 * nrow] = tmpx;
385 tmpArray[2 + 4 * nrow] = 0.0;
386 tmpArray[3 + 0 * nrow] =
387 -tmpx * (FourThird * u * u + v * v +
w *
w + tmp2 * (E - q2)) +
388 tmpy * (-OneThird * u * v) + tmpz * (-OneThird * u *
w);
389 tmpArray[3 + 1 * nrow] = tmpx * (FourThird - tmp2) * u +
390 tmpy * (-TwoThird * v) + tmpz * (-TwoThird *
w);
391 tmpArray[3 + 2 * nrow] = tmpx * (1.0 - tmp2) * v + tmpy * u;
392 tmpArray[3 + 3 * nrow] = tmpx * (1.0 - tmp2) *
w + tmpz * u;
393 tmpArray[3 + 4 * nrow] = tmpx * tmp2;
435 OneThird = 1.0 / 3.0;
436 TwoThird = 2.0 * OneThird;
437 FourThird = 4.0 * OneThird;
440 tmpArray = OutputMatrix->GetPtr();
441 int nrow = OutputMatrix->GetRows();
443 tmpArray[0 + 0 * nrow] = tmpx * (TwoThird * v) + tmpy * (-u);
444 tmpArray[0 + 1 * nrow] = tmpy;
445 tmpArray[0 + 2 * nrow] = tmpx * (-TwoThird);
446 tmpArray[0 + 3 * nrow] = 0.0;
447 tmpArray[0 + 4 * nrow] = 0.0;
448 tmpArray[1 + 0 * nrow] =
449 tmpx * (-u) + tmpy * (-FourThird * v) + tmpz * (-
w);
450 tmpArray[1 + 1 * nrow] = tmpx;
451 tmpArray[1 + 2 * nrow] = tmpy * FourThird;
452 tmpArray[1 + 3 * nrow] = tmpz;
453 tmpArray[1 + 4 * nrow] = 0.0;
454 tmpArray[2 + 0 * nrow] = tmpy * (-
w) + tmpz * (TwoThird * v);
455 tmpArray[2 + 1 * nrow] = 0.0;
456 tmpArray[2 + 2 * nrow] = tmpz * (-TwoThird);
457 tmpArray[2 + 3 * nrow] = tmpy;
458 tmpArray[2 + 4 * nrow] = 0.0;
459 tmpArray[3 + 0 * nrow] =
460 tmpx * (-OneThird * u * v) -
461 tmpy * (u * u + FourThird * v * v +
w *
w + tmp2 * (E - q2)) +
462 tmpz * (-OneThird * v *
w);
463 tmpArray[3 + 1 * nrow] = tmpx * v + tmpy * (1 - tmp2) * u;
464 tmpArray[3 + 2 * nrow] = tmpx * (-TwoThird * u) +
465 tmpy * (FourThird - tmp2) * v +
466 tmpz * (-TwoThird *
w);
467 tmpArray[3 + 3 * nrow] = tmpy * (1 - tmp2) *
w + tmpz * v;
468 tmpArray[3 + 4 * nrow] = tmpy * tmp2;
510 OneThird = 1.0 / 3.0;
511 TwoThird = 2.0 * OneThird;
512 FourThird = 4.0 * OneThird;
515 tmpArray = OutputMatrix->GetPtr();
516 int nrow = OutputMatrix->GetRows();
518 tmpArray[0 + 0 * nrow] = tmpx * (TwoThird *
w) + tmpz * (-u);
519 tmpArray[0 + 1 * nrow] = tmpz;
520 tmpArray[0 + 2 * nrow] = 0.0;
521 tmpArray[0 + 3 * nrow] = tmpx * (-TwoThird);
522 tmpArray[0 + 4 * nrow] = 0.0;
523 tmpArray[1 + 0 * nrow] = tmpy * (TwoThird *
w) + tmpz * (-v);
524 tmpArray[1 + 1 * nrow] = 0.0;
525 tmpArray[1 + 2 * nrow] = tmpz;
526 tmpArray[1 + 3 * nrow] = tmpy * (-TwoThird);
527 tmpArray[1 + 4 * nrow] = 0.0;
528 tmpArray[2 + 0 * nrow] =
529 tmpx * (-u) + tmpy * (-v) + tmpz * (-FourThird *
w);
530 tmpArray[2 + 1 * nrow] = tmpx;
531 tmpArray[2 + 2 * nrow] = tmpy;
532 tmpArray[2 + 3 * nrow] = tmpz * FourThird;
533 tmpArray[2 + 4 * nrow] = 0.0;
534 tmpArray[3 + 0 * nrow] =
535 tmpx * (-OneThird * u *
w) + tmpy * (-OneThird * v *
w) -
536 tmpz * (u * u + v * v + FourThird *
w *
w + tmp2 * (E - q2));
537 tmpArray[3 + 1 * nrow] = tmpx *
w + tmpz * (1 - tmp2) * u;
538 tmpArray[3 + 2 * nrow] = tmpy *
w + tmpz * (1 - tmp2) * v;
539 tmpArray[3 + 3 * nrow] = tmpx * (-TwoThird * u) + tmpy * (-TwoThird * v) +
540 tmpz * (FourThird - tmp2) *
w;
541 tmpArray[3 + 4 * nrow] = tmpz * tmp2;
561 tmpArray = OutputMatrix->GetPtr();
562 int nrow = OutputMatrix->GetRows();
586 orho2 = orho1 * orho1;
587 orho3 = orho1 * orho2;
588 orho4 = orho2 * orho2;
597 NekDouble du_dx = orho1 * (dU2_dx - u * dU1_dx);
598 NekDouble dv_dx = orho1 * (dU3_dx - v * dU1_dx);
599 NekDouble du_dy = orho1 * (dU2_dy - u * dU1_dy);
600 NekDouble dv_dy = orho1 * (dU3_dy - v * dU1_dy);
601 NekDouble s12 = FourThird * du_dx - TwoThrid * dv_dy;
604 NekDouble s23 = FourThird * dv_dy - TwoThrid * du_dx;
609 (orho1 * dU4_dx - U[3] * orho2 * dU1_dx -
610 u * (orho1 * dU2_dx - U[1] * orho2 * dU1_dx) -
611 v * (orho1 * dU3_dx - U[2] * orho2 * dU1_dx));
613 (orho1 * dU4_dy - U[3] * orho2 * dU1_dy -
614 u * (orho1 * dU2_dy - U[1] * orho2 * dU1_dy) -
615 v * (orho1 * dU3_dy - U[2] * orho2 * dU1_dy));
622 tmp[2] = snv - qn / mu;
624 dT_dU[0] = oCv * (-orho2 * U4 + orho3 * U2 * U2 + orho3 * U3 * U3);
625 dT_dU[1] = -oCv * orho2 * U2;
626 dT_dU[2] = -oCv * orho2 * U3;
627 dT_dU[3] = oCv * orho1;
628 for (
int i = 0; i < 3; i++)
630 for (
int j = 0; j < 4; j++)
632 tmpArray[i + j * nrow] = dmu_dT * dT_dU[j] * tmp[i];
649 du_dx_dU1 = -orho2 * dU2_dx + 2 * orho3 * U2 * dU1_dx;
650 du_dx_dU2 = -orho2 * dU1_dx;
651 du_dy_dU1 = -orho2 * dU2_dy + 2 * orho3 * U2 * dU1_dy;
652 du_dy_dU2 = -orho2 * dU1_dy;
653 dv_dx_dU1 = -orho2 * dU3_dx + 2 * orho3 * U3 * dU1_dx;
654 dv_dx_dU3 = du_dx_dU2;
655 dv_dy_dU1 = -orho2 * dU3_dy + 2 * orho3 * U3 * dU1_dy;
656 dv_dy_dU3 = du_dy_dU2;
657 ds12_dU1 = FourThird * du_dx_dU1 - TwoThrid * dv_dy_dU1;
658 ds12_dU2 = FourThird * du_dx_dU2;
659 ds12_dU3 = -TwoThrid * dv_dy_dU3;
660 ds13_dU1 = du_dy_dU1 + dv_dx_dU1;
661 ds13_dU2 = du_dy_dU2;
662 ds13_dU3 = dv_dx_dU3;
666 ds23_dU1 = FourThird * dv_dy_dU1 - TwoThrid * du_dx_dU1;
667 ds23_dU2 = -TwoThrid * du_dx_dU2;
668 ds23_dU3 = FourThird * dv_dy_dU3;
669 dsnx_dU1 = ds12_dU1 * nx + ds22_dU1 * ny;
670 dsnx_dU2 = ds12_dU2 * nx + ds22_dU2 * ny;
671 dsnx_dU3 = ds12_dU3 * nx + ds22_dU3 * ny;
672 dsny_dU1 = ds13_dU1 * nx + ds23_dU1 * ny;
673 dsny_dU2 = ds13_dU2 * nx + ds23_dU2 * ny;
674 dsny_dU3 = ds13_dU3 * nx + ds23_dU3 * ny;
676 u * dsnx_dU1 + v * dsny_dU1 - orho2 * U2 * snx - orho2 * U3 * sny;
677 dsnv_dU2 = u * dsnx_dU2 + v * dsny_dU2 + orho1 * snx;
678 dsnv_dU3 = u * dsnx_dU3 + v * dsny_dU3 + orho1 * sny;
679 tmpArray[0 + 0 * nrow] = tmpArray[0 + 0 * nrow] + mu * dsnx_dU1;
680 tmpArray[0 + 1 * nrow] = tmpArray[0 + 1 * nrow] + mu * dsnx_dU2;
681 tmpArray[0 + 2 * nrow] = tmpArray[0 + 2 * nrow] + mu * dsnx_dU3;
682 tmpArray[1 + 0 * nrow] = tmpArray[1 + 0 * nrow] + mu * dsny_dU1;
683 tmpArray[1 + 1 * nrow] = tmpArray[1 + 1 * nrow] + mu * dsny_dU2;
684 tmpArray[1 + 2 * nrow] = tmpArray[1 + 2 * nrow] + mu * dsny_dU3;
685 tmpArray[2 + 0 * nrow] = tmpArray[2 + 0 * nrow] + mu * dsnv_dU1;
686 tmpArray[2 + 1 * nrow] = tmpArray[2 + 1 * nrow] + mu * dsnv_dU2;
687 tmpArray[2 + 2 * nrow] = tmpArray[2 + 2 * nrow] + mu * dsnv_dU3;
690 NekDouble dqx_dU1, dqx_dU2, dqx_dU3, dqx_dU4;
691 NekDouble dqy_dU1, dqy_dU2, dqy_dU3, dqy_dU4;
693 dqx_dU1 = tmpx * (-orho2 * dU4_dx + 2 * orho3 * U4 * dU1_dx +
694 2 * orho3 * U2 * dU2_dx - 3 * orho4 * U2 * U2 * dU1_dx +
695 2 * orho3 * U3 * dU3_dx - 3 * orho4 * U3 * U3 * dU1_dx);
696 dqx_dU2 = tmpx * (-orho2 * dU2_dx + 2 * orho3 * U2 * dU1_dx);
697 dqx_dU3 = tmpx * (-orho2 * dU3_dx + 2 * orho3 * U3 * dU1_dx);
698 dqx_dU4 = -tmpx * orho2 * dU1_dx;
700 dqy_dU1 = tmpy * (-orho2 * dU4_dy + 2 * orho3 * U4 * dU1_dy +
701 2 * orho3 * U2 * dU2_dy - 3 * orho4 * U2 * U2 * dU1_dy +
702 2 * orho3 * U3 * dU3_dy - 3 * orho4 * U3 * U3 * dU1_dy);
703 dqy_dU2 = tmpy * (-orho2 * dU2_dy + 2 * orho3 * U2 * dU1_dy);
704 dqy_dU3 = tmpy * (-orho2 * dU3_dy + 2 * orho3 * U3 * dU1_dy);
705 dqy_dU4 = -tmpy * orho2 * dU1_dy;
706 tmpArray[2 + 0 * nrow] = tmpArray[2 + 0 * nrow] - dqx_dU1 - dqy_dU1;
707 tmpArray[2 + 1 * nrow] = tmpArray[2 + 1 * nrow] - dqx_dU2 - dqy_dU2;
708 tmpArray[2 + 2 * nrow] = tmpArray[2 + 2 * nrow] - dqx_dU3 - dqy_dU3;
709 tmpArray[2 + 3 * nrow] = tmpArray[2 + 3 * nrow] - dqx_dU4 - dqy_dU4;
729 tmpArray = OutputMatrix->GetPtr();
730 int nrow = OutputMatrix->GetRows();
763 orho2 = orho1 * orho1;
764 orho3 = orho1 * orho2;
765 orho4 = orho2 * orho2;
775 NekDouble du_dx = orho1 * (dU2_dx - u * dU1_dx);
776 NekDouble dv_dx = orho1 * (dU3_dx - v * dU1_dx);
777 NekDouble dw_dx = orho1 * (dU4_dx -
w * dU1_dx);
778 NekDouble du_dy = orho1 * (dU2_dy - u * dU1_dy);
779 NekDouble dv_dy = orho1 * (dU3_dy - v * dU1_dy);
780 NekDouble dw_dy = orho1 * (dU4_dy -
w * dU1_dy);
781 NekDouble du_dz = orho1 * (dU2_dz - u * dU1_dz);
782 NekDouble dv_dz = orho1 * (dU3_dz - v * dU1_dz);
783 NekDouble dw_dz = orho1 * (dU4_dz -
w * dU1_dz);
784 NekDouble s12 = FourThird * du_dx - TwoThrid * dv_dy - TwoThrid * dw_dz;
788 NekDouble s23 = FourThird * dv_dy - TwoThrid * du_dx - TwoThrid * dw_dz;
792 NekDouble s34 = FourThird * dw_dz - TwoThrid * du_dx - TwoThrid * dv_dy;
793 NekDouble snx = s12 * nx + s22 * ny + s32 * nz;
794 NekDouble sny = s13 * nx + s23 * ny + s33 * nz;
795 NekDouble snz = s14 * nz + s24 * ny + s34 * nz;
797 NekDouble qx = -tmp2 * (orho1 * dU5_dx - U5 * orho2 * dU1_dx -
798 u * (orho1 * dU2_dx - U2 * orho2 * dU1_dx) -
799 v * (orho1 * dU3_dx - U3 * orho2 * dU1_dx) -
800 w * (orho1 * dU4_dx - U4 * orho2 * dU1_dx));
801 NekDouble qy = -tmp2 * (orho1 * dU5_dy - U5 * orho2 * dU1_dy -
802 u * (orho1 * dU2_dy - U2 * orho2 * dU1_dy) -
803 v * (orho1 * dU3_dy - U3 * orho2 * dU1_dy) -
804 w * (orho1 * dU4_dy - U4 * orho2 * dU1_dy));
805 NekDouble qz = -tmp2 * (orho1 * dU5_dz - U5 * orho2 * dU1_dz -
806 u * (orho1 * dU2_dz - U2 * orho2 * dU1_dz) -
807 v * (orho1 * dU3_dz - U3 * orho2 * dU1_dz) -
808 w * (orho1 * dU4_dz - U4 * orho2 * dU1_dz));
809 NekDouble qn = qx * nx + qy * ny + qz * nz;
816 tmp[3] = snv - qn / mu;
818 dT_dU[0] = oCv * (-orho2 * U5 + orho3 * U2 * U2 + orho3 * U3 * U3 +
820 dT_dU[1] = -oCv * orho2 * U2;
821 dT_dU[2] = -oCv * orho2 * U3;
822 dT_dU[3] = -oCv * orho2 * U4;
823 dT_dU[4] = oCv * orho1;
824 for (
int i = 0; i < 4; i++)
826 for (
int j = 0; j < 5; j++)
828 tmpArray[i + j * nrow] = dmu_dT * dT_dU[j] * tmp[i];
842 NekDouble ds12_dU1, ds12_dU2, ds12_dU3, ds12_dU4;
846 NekDouble ds23_dU1, ds23_dU2, ds23_dU3, ds23_dU4;
850 NekDouble ds34_dU1, ds34_dU2, ds34_dU3, ds34_dU4;
851 NekDouble dsnx_dU1, dsnx_dU2, dsnx_dU3, dsnx_dU4;
852 NekDouble dsny_dU1, dsny_dU2, dsny_dU3, dsny_dU4;
853 NekDouble dsnz_dU1, dsnz_dU2, dsnz_dU3, dsnz_dU4;
854 NekDouble dsnv_dU1, dsnv_dU2, dsnv_dU3, dsnv_dU4;
856 du_dx_dU1 = -orho2 * dU2_dx + 2 * orho3 * U2 * dU1_dx;
857 du_dx_dU2 = -orho2 * dU1_dx;
858 du_dy_dU1 = -orho2 * dU2_dy + 2 * orho3 * U2 * dU1_dy;
859 du_dy_dU2 = -orho2 * dU1_dy;
860 du_dz_dU1 = -orho2 * dU2_dz + 2 * orho3 * U2 * dU1_dz;
861 du_dz_dU2 = -orho2 * dU1_dz;
862 dv_dx_dU1 = -orho2 * dU3_dx + 2 * orho3 * U3 * dU1_dx;
863 dv_dx_dU3 = -orho2 * dU1_dx;
864 dv_dy_dU1 = -orho2 * dU3_dy + 2 * orho3 * U3 * dU1_dy;
865 dv_dy_dU3 = -orho2 * dU1_dy;
866 dv_dz_dU1 = -orho2 * dU3_dz + 2 * orho3 * U3 * dU1_dz;
867 dv_dz_dU3 = -orho2 * dU1_dz;
868 dw_dx_dU1 = -orho2 * dU4_dx + 2 * orho3 * U4 * dU1_dx;
869 dw_dx_dU4 = -orho2 * dU1_dx;
870 dw_dy_dU1 = -orho2 * dU4_dy + 2 * orho3 * U4 * dU1_dy;
871 dw_dy_dU4 = -orho2 * dU1_dy;
872 dw_dz_dU1 = -orho2 * dU4_dz + 2 * orho3 * U4 * dU1_dz;
873 dw_dz_dU4 = -orho2 * dU1_dz;
875 FourThird * du_dx_dU1 - TwoThrid * dv_dy_dU1 - TwoThrid * dw_dz_dU1;
876 ds12_dU2 = FourThird * du_dx_dU2;
877 ds12_dU3 = -TwoThrid * dv_dy_dU3;
878 ds12_dU4 = -TwoThrid * dw_dz_dU4;
879 ds13_dU1 = du_dy_dU1 + dv_dx_dU1;
880 ds13_dU2 = du_dy_dU2;
881 ds13_dU3 = dv_dx_dU3;
882 ds14_dU1 = dw_dx_dU1 + du_dz_dU1;
883 ds14_dU2 = du_dz_dU2;
884 ds14_dU4 = dw_dx_dU4;
885 ds22_dU1 = du_dy_dU1 + dv_dx_dU1;
886 ds22_dU2 = du_dy_dU2;
887 ds22_dU3 = dv_dx_dU3;
889 FourThird * dv_dy_dU1 - TwoThrid * du_dx_dU1 - TwoThrid * dw_dz_dU1;
890 ds23_dU2 = -TwoThrid * du_dx_dU2;
891 ds23_dU3 = FourThird * dv_dy_dU3;
892 ds23_dU4 = -TwoThrid * dw_dz_dU4;
893 ds24_dU1 = dv_dz_dU1 + dw_dy_dU1;
894 ds24_dU3 = dv_dz_dU3;
895 ds24_dU4 = dw_dy_dU4;
896 ds32_dU1 = dw_dx_dU1 + du_dz_dU1;
897 ds32_dU2 = du_dz_dU2;
898 ds32_dU4 = dw_dx_dU4;
899 ds33_dU1 = dv_dz_dU1 + dw_dy_dU1;
900 ds33_dU3 = dv_dz_dU3;
901 ds33_dU4 = dw_dy_dU4;
903 FourThird * dw_dz_dU1 - TwoThrid * du_dx_dU1 - TwoThrid * dv_dy_dU1;
904 ds34_dU2 = -TwoThrid * du_dx_dU2;
905 ds34_dU3 = -TwoThrid * dv_dy_dU3;
906 ds34_dU4 = FourThird * dw_dz_dU4;
907 dsnx_dU1 = ds12_dU1 * nx + ds22_dU1 * ny + ds32_dU1 * nz;
908 dsnx_dU2 = ds12_dU2 * nx + ds22_dU2 * ny + ds32_dU2 * nz;
909 dsnx_dU3 = ds12_dU3 * nx + ds22_dU3 * ny;
910 dsnx_dU4 = ds12_dU4 * nx + ds32_dU4 * nz;
911 dsny_dU1 = ds13_dU1 * nx + ds23_dU1 * ny + ds33_dU1 * nz;
912 dsny_dU2 = ds13_dU2 * nx + ds23_dU2 * ny;
913 dsny_dU3 = ds13_dU3 * nx + ds23_dU3 * ny + ds33_dU3 * nz;
914 dsny_dU4 = ds23_dU4 * ny + ds33_dU4 * nz;
915 dsnz_dU1 = ds14_dU1 * nx + ds24_dU1 * ny + ds34_dU1 * nz;
916 dsnz_dU2 = ds14_dU2 * nx + ds34_dU2 * nz;
917 dsnz_dU3 = ds24_dU3 * ny + ds34_dU3 * nz;
919 dsnz_dU4 = ds14_dU4 * nx + ds24_dU4 * ny + ds34_dU4 * nz;
920 dsnv_dU1 = u * dsnx_dU1 + v * dsny_dU1 +
w * dsnz_dU1 - orho2 * U2 * snx -
921 orho2 * U3 * sny - orho2 * U4 * snz;
922 dsnv_dU2 = u * dsnx_dU2 + v * dsny_dU2 +
w * dsnz_dU2 + orho1 * snx;
923 dsnv_dU3 = u * dsnx_dU3 + v * dsny_dU3 +
w * dsnz_dU3 + orho1 * sny;
924 dsnv_dU4 = u * dsnx_dU4 + v * dsny_dU4 +
w * dsnz_dU4 + orho1 * snz;
925 tmpArray[0 + 0 * nrow] = tmpArray[0 + 0 * nrow] + mu * dsnx_dU1;
926 tmpArray[0 + 1 * nrow] = tmpArray[0 + 1 * nrow] + mu * dsnx_dU2;
927 tmpArray[0 + 2 * nrow] = tmpArray[0 + 2 * nrow] + mu * dsnx_dU3;
928 tmpArray[0 + 3 * nrow] = tmpArray[0 + 3 * nrow] + mu * dsnx_dU4;
929 tmpArray[1 + 0 * nrow] = tmpArray[1 + 0 * nrow] + mu * dsny_dU1;
930 tmpArray[1 + 1 * nrow] = tmpArray[1 + 1 * nrow] + mu * dsny_dU2;
931 tmpArray[1 + 2 * nrow] = tmpArray[1 + 2 * nrow] + mu * dsny_dU3;
932 tmpArray[1 + 3 * nrow] = tmpArray[1 + 3 * nrow] + mu * dsny_dU4;
933 tmpArray[2 + 0 * nrow] = tmpArray[2 + 0 * nrow] + mu * dsnz_dU1;
934 tmpArray[2 + 1 * nrow] = tmpArray[2 + 1 * nrow] + mu * dsnz_dU2;
935 tmpArray[2 + 2 * nrow] = tmpArray[2 + 2 * nrow] + mu * dsnz_dU3;
936 tmpArray[2 + 3 * nrow] = tmpArray[2 + 3 * nrow] + mu * dsnz_dU4;
937 tmpArray[3 + 0 * nrow] = tmpArray[3 + 0 * nrow] + mu * dsnv_dU1;
938 tmpArray[3 + 1 * nrow] = tmpArray[3 + 1 * nrow] + mu * dsnv_dU2;
939 tmpArray[3 + 2 * nrow] = tmpArray[3 + 2 * nrow] + mu * dsnv_dU3;
940 tmpArray[3 + 3 * nrow] = tmpArray[3 + 3 * nrow] + mu * dsnv_dU4;
943 NekDouble dqx_dU1, dqx_dU2, dqx_dU3, dqx_dU4, dqx_dU5;
944 NekDouble dqy_dU1, dqy_dU2, dqy_dU3, dqy_dU4, dqy_dU5;
945 NekDouble dqz_dU1, dqz_dU2, dqz_dU3, dqz_dU4, dqz_dU5;
947 dqx_dU1 = tmpx * (-orho2 * dU5_dx + 2 * orho3 * U5 * dU1_dx +
948 2 * orho3 * U2 * dU2_dx - 3 * orho4 * U2 * U2 * dU1_dx +
949 2 * orho3 * U3 * dU3_dx - 3 * orho4 * U3 * U3 * dU1_dx +
950 2 * orho3 * U4 * dU4_dx - 3 * orho4 * U4 * U4 * dU1_dx);
951 dqx_dU2 = tmpx * (-orho2 * dU2_dx + 2 * orho3 * U2 * dU1_dx);
952 dqx_dU3 = tmpx * (-orho2 * dU3_dx + 2 * orho3 * U3 * dU1_dx);
953 dqx_dU4 = tmpx * (-orho2 * dU4_dx + 2 * orho3 * U4 * dU1_dx);
954 dqx_dU5 = -tmpx * orho2 * dU1_dx;
956 dqy_dU1 = tmpy * (-orho2 * dU5_dy + 2 * orho3 * U5 * dU1_dy +
957 2 * orho3 * U2 * dU2_dy - 3 * orho4 * U2 * U2 * dU1_dy +
958 2 * orho3 * U3 * dU3_dy - 3 * orho4 * U3 * U3 * dU1_dy +
959 2 * orho3 * U4 * dU4_dy - 3 * orho4 * U4 * U4 * dU1_dy);
960 dqy_dU2 = tmpy * (-orho2 * dU2_dy + 2 * orho3 * U2 * dU1_dy);
961 dqy_dU3 = tmpy * (-orho2 * dU3_dy + 2 * orho3 * U3 * dU1_dy);
962 dqy_dU4 = tmpy * (-orho2 * dU4_dy + 2 * orho3 * U4 * dU1_dy);
963 dqy_dU5 = -tmpy * orho2 * dU1_dy;
965 dqz_dU1 = tmpz * (-orho2 * dU5_dz + 2 * orho3 * U5 * dU1_dz +
966 2 * orho3 * U2 * dU2_dz - 3 * orho4 * U2 * U2 * dU1_dz +
967 2 * orho3 * U3 * dU3_dz - 3 * orho4 * U3 * U3 * dU1_dz +
968 2 * orho3 * U4 * dU4_dz - 3 * orho4 * U4 * U4 * dU1_dz);
969 dqz_dU2 = tmpz * (-orho2 * dU2_dz + 2 * orho3 * U2 * dU1_dz);
970 dqz_dU3 = tmpz * (-orho2 * dU3_dz + 2 * orho3 * U3 * dU1_dz);
971 dqz_dU4 = tmpz * (-orho2 * dU4_dz + 2 * orho3 * U4 * dU1_dz);
972 dqz_dU5 = -tmpz * orho2 * dU1_dz;
973 tmpArray[3 + 0 * nrow] =
974 tmpArray[3 + 0 * nrow] - dqx_dU1 - dqy_dU1 - dqz_dU1;
975 tmpArray[3 + 1 * nrow] =
976 tmpArray[3 + 1 * nrow] - dqx_dU2 - dqy_dU2 - dqz_dU2;
977 tmpArray[3 + 2 * nrow] =
978 tmpArray[3 + 2 * nrow] - dqx_dU3 - dqy_dU3 - dqz_dU3;
979 tmpArray[3 + 3 * nrow] =
980 tmpArray[3 + 3 * nrow] - dqx_dU4 - dqy_dU4 - dqz_dU4;
981 tmpArray[3 + 4 * nrow] =
982 tmpArray[3 + 4 * nrow] - dqx_dU5 - dqy_dU5 - dqz_dU5;
986 const int nConvectiveFields,
const int nElmtPnt,
993 int nSpaceDim =
m_graph->GetSpaceDimension();
999 for (
int j = 0; j < nSpaceDim; j++)
1009 for (
int npnt = 0; npnt < nElmtPnt; npnt++)
1011 for (
int j = 0; j < nConvectiveFields; j++)
1013 pointVar[j] = locVars[j][npnt];
1015 for (
int j = 0; j < nSpaceDim; j++)
1017 for (
int k = 0; k < nConvectiveFields; k++)
1019 pointDerv[j][k] = locDerv[j][k][npnt];
1023 pointmu = locmu[npnt];
1024 pointDmuDT = locDmuDT[npnt];
1028 for (
int j = 0; j < nConvectiveFields; j++)
1030 int noffset = j * nConvectiveFields;
1033 tmp1 = PntJacArray[npnt] + (noffset + 1), 1,
1034 tmp2 = wspMatData + (noffset - j), 1,
1035 tmp3 = PntJacArray[npnt] + (noffset + 1), 1);
1049 GetdFlux_dU_2D(normals, mu, DmuDT, conservVar, conseDeriv, fluxJac);
1053 GetdFlux_dU_3D(normals, mu, DmuDT, conservVar, conseDeriv, fluxJac);
1069 int nConvectiveFields = inarray.size();
1070 std::shared_ptr<LocalRegions::ExpansionVector> expvect = explist->GetExp();
1071 int nTotElmt = (*expvect).size();
1072 int nPts = explist->GetTotPoints();
1073 int nSpaceDim =
m_graph->GetSpaceDimension();
1079 m_varConv->GetTemperature(inarray, temperature);
1090 nConvectiveFields - 1, nConvectiveFields);
1093 for (
int nelmt = 0; nelmt < nTotElmt; nelmt++)
1095 int nElmtPnt = (*expvect)[nelmt]->GetTotPoints();
1096 int noffest = explist->GetPhys_Offset(nelmt);
1098 for (
int j = 0; j < nConvectiveFields; j++)
1100 locVars[j] = inarray[j] + noffest;
1103 for (
int j = 0; j < nSpaceDim; j++)
1105 locnormal[j] = normals[j] + noffest;
1108 locmu = mu + noffest;
1109 for (
int npnt = 0; npnt < nElmtPnt; npnt++)
1111 for (
int j = 0; j < nConvectiveFields; j++)
1113 pointVar[j] = locVars[j][npnt];
1115 for (
int j = 0; j < nSpaceDim; j++)
1117 pointnormals[j] = locnormal[j][npnt];
1120 pointmu = locmu[npnt];
1124 for (
int j = 0; j < nConvectiveFields; j++)
1126 ElmtJacArray[0][j][nFluxDir][nelmt][npnt] = 0.0;
1128 for (
int j = 0; j < nConvectiveFields; j++)
1130 int noffset = j * (nConvectiveFields - 1);
1131 for (
int i = 0; i < nConvectiveFields - 1; i++)
1133 ElmtJacArray[i + 1][j][nFluxDir][nelmt][npnt] =
1134 PointFJac_data[noffset + i];
1142 const int nConvectiveFields,
const int nElmtPnt,
const int nDervDir,
1148 int nSpaceDim =
m_graph->GetSpaceDimension();
1159 for (
int npnt = 0; npnt < nElmtPnt; npnt++)
1161 for (
int j = 0; j < nConvectiveFields; j++)
1163 pointVar[j] = locVars[j][npnt];
1165 for (
int j = 0; j < nSpaceDim; j++)
1167 pointnormals[j] = locnormal[j][npnt];
1170 pointmu = locmu[npnt];
1174 Vmath::Zero(nConvectiveFields, PntJacArray[npnt], nConvectiveFields);
1175 for (
int j = 0; j < nConvectiveFields; j++)
1177 int noffset = j * (nConvectiveFields - 1);
1178 Vmath::Vcopy((nConvectiveFields - 1), tmp1 = wspMatData + noffset,
1179 1, tmp2 = PntJacArray[npnt] + noffset + j + 1, 1);
1191 int nConvectiveFields = inarray.size();
1192 std::shared_ptr<LocalRegions::ExpansionVector> expvect = explist->GetExp();
1193 int nTotElmt = (*expvect).size();
1194 int nPts = explist->GetTotPoints();
1195 int nSpaceDim =
m_graph->GetSpaceDimension();
1198 if (!ElmtJac.size())
1201 for (
int nelmt = 0; nelmt < nTotElmt; nelmt++)
1203 int nElmtPnt = (*expvect)[nelmt]->GetTotPoints();
1205 for (
int npnt = 0; npnt < nElmtPnt; npnt++)
1207 ElmtJac[nelmt][npnt] =
1209 nConvectiveFields, nConvectiveFields);
1218 m_varConv->GetTemperature(inarray, temperature);
1231 nConvectiveFields - 1, nConvectiveFields);
1235 for (
int nelmt = 0; nelmt < nTotElmt; nelmt++)
1237 int nElmtPnt = (*expvect)[nelmt]->GetTotPoints();
1238 int noffest = explist->GetPhys_Offset(nelmt);
1240 for (
int j = 0; j < nConvectiveFields; j++)
1242 locVars[j] = inarray[j] + noffest;
1245 for (
int j = 0; j < nSpaceDim; j++)
1247 locnormal[j] = normals[j] + noffest;
1250 locmu = mu + noffest;
1251 for (
int npnt = 0; npnt < nElmtPnt; npnt++)
1253 for (
int j = 0; j < nConvectiveFields; j++)
1255 pointVar[j] = locVars[j][npnt];
1257 for (
int j = 0; j < nSpaceDim; j++)
1259 pointnormals[j] = locnormal[j][npnt];
1262 pointmu = locmu[npnt];
1266 tmpMatinnData = PointFJac->GetPtr();
1267 tmpMatoutData = ElmtJac[nelmt][npnt]->GetPtr();
1269 Vmath::Fill(nConvectiveFields, 0.0, tmpMatoutData,
1271 for (
int j = 0; j < nConvectiveFields; j++)
1274 nConvectiveFields - 1,
1275 tmp1 = tmpMatinnData + (j * (nConvectiveFields - 1)), 1,
1276 tmp2 = tmpMatoutData + (1 + j * nConvectiveFields), 1);
1286 int nConvectiveFields =
m_fields.size();
1294 for (
int j = 0; j < nConvectiveFields; j++)
1311 ASSERTL1(
false,
"LDGNS not yet validated for implicit compressible "
1316 for (
int i = 0; i < nConvectiveFields - 1; ++i)
1318 inarrayDiff[i] = inarray[i];
1330 int nPts = mu.size();
1334 m_varConv->GetTemperature(inarray, temperature);
1339 if (DmuDT.size() > 0)
1341 m_varConv->GetDmuDT(temperature, mu, DmuDT);
1346 if (DmuDT.size() > 0)
1354 const std::string type)
const
1356 if (type ==
"Physical" || type ==
"Off")
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
bool m_updateShockCaptPhys
void v_InitObject(bool DeclareFields=true) override
Initialization object for CFSImplicit class.
NekDouble m_bndEvaluateTime
SolverUtils::DiffusionSharedPtr m_diffusion
VariableConverterSharedPtr m_varConv
ArtificialDiffusionSharedPtr m_artificialDiffusion
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
void GetDivCurlSquared(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &cnsVar, Array< OneD, NekDouble > &div, Array< OneD, NekDouble > &curlSquare, const Array< OneD, Array< OneD, NekDouble > > &cnsVarFwd, const Array< OneD, Array< OneD, NekDouble > > &cnsVarBwd)
Get divergence and curl squared.
void InitObject_Explicit()
bool m_is_shockCaptPhys
flag for shock capturing switch on/off an enum could be added for more options
void GetViscosityAndThermalCondFromTemp(const Array< OneD, NekDouble > &temperature, Array< OneD, NekDouble > &mu, Array< OneD, NekDouble > &thermalCond)
Update viscosity todo: add artificial viscosity here.
bool m_is_diffIP
flag to switch between IP and LDG an enum could be added for more options
std::string m_ViscosityType
void GetdFlux_dU_2D(const Array< OneD, NekDouble > &normals, const NekDouble mu, const NekDouble dmu_dT, const Array< OneD, NekDouble > &U, const Array< OneD, const Array< OneD, NekDouble > > &qfield, DNekMatSharedPtr &OutputMatrix)
return part of viscous Jacobian Input: normals:Point normals mu: dynamicviscosity dmu_dT: mu's deriva...
Array< OneD, GetdFlux_dDeriv > m_GetdFlux_dDeriv_Array
static std::string className
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void GetdFlux_dQx_3D(const Array< OneD, NekDouble > &normals, const NekDouble &mu, const Array< OneD, NekDouble > &U, DNekMatSharedPtr &OutputMatrix)
return part of viscous Jacobian derived with Qx=[drho_dx,drhou_dx,drhov_dx,drhow_dx,...
bool v_SupportsShockCaptType(const std::string type) const final
void GetdFlux_dU_3D(const Array< OneD, NekDouble > &normals, const NekDouble mu, const NekDouble dmu_dT, const Array< OneD, NekDouble > &U, const Array< OneD, const Array< OneD, NekDouble > > &qfield, DNekMatSharedPtr &OutputMatrix)
return part of viscous Jacobian Input: normals:Point normals mu: dynamicviscosity dmu_dT: mu's deriva...
virtual void v_GetDiffusionFluxJacPoint(const Array< OneD, NekDouble > &conservVar, const Array< OneD, const Array< OneD, NekDouble > > &conseDeriv, const NekDouble mu, const NekDouble DmuDT, const Array< OneD, NekDouble > &normals, DNekMatSharedPtr &fluxJac)
NavierStokesImplicitCFE(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void GetdFlux_dQx_2D(const Array< OneD, NekDouble > &normals, const NekDouble &mu, const Array< OneD, NekDouble > &U, DNekMatSharedPtr &OutputMatrix)
return part of viscous Jacobian:
void v_InitObject(bool DeclareFields=true) override
Initialization object for CompressibleFlowSystem class.
void v_DoDiffusionCoeff(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, const Array< OneD, NekDouble > > &pFwd, const Array< OneD, const Array< OneD, NekDouble > > &pBwd) override
void GetdFlux_dQy_3D(const Array< OneD, NekDouble > &normals, const NekDouble &mu, const Array< OneD, NekDouble > &U, DNekMatSharedPtr &OutputMatrix)
return part of viscous Jacobian derived with Qy=[drho_dy,drhou_dy,drhov_dy,drhow_dy,...
void v_GetFluxDerivJacDirctn(const MultiRegions::ExpListSharedPtr &explist, const Array< OneD, const Array< OneD, NekDouble > > &normals, const int nDervDir, const Array< OneD, const Array< OneD, NekDouble > > &inarray, TensorOfArray5D< NekDouble > &ElmtJacArray, const int nFluxDir) override
void v_MinusDiffusionFluxJacPoint(const int nConvectiveFields, const int nElmtPnt, const Array< OneD, const Array< OneD, NekDouble > > &locVars, const TensorOfArray3D< NekDouble > &locDerv, const Array< OneD, NekDouble > &locmu, const Array< OneD, NekDouble > &locDmuDT, const Array< OneD, NekDouble > &normals, DNekMatSharedPtr &wspMat, Array< OneD, Array< OneD, NekDouble > > &PntJacArray) override
void GetdFlux_dQz_3D(const Array< OneD, NekDouble > &normals, const NekDouble &mu, const Array< OneD, NekDouble > &U, DNekMatSharedPtr &OutputMatrix)
return part of viscous Jacobian derived with Qz=[drho_dz,drhou_dz,drhov_dz,drhow_dz,...
void GetdFlux_dQy_2D(const Array< OneD, NekDouble > &normals, const NekDouble &mu, const Array< OneD, NekDouble > &U, DNekMatSharedPtr &OutputMatrix)
return part of viscous Jacobian:
void v_GetFluxDerivJacDirctnElmt(const int nConvectiveFields, const int nElmtPnt, const int nDervDir, const Array< OneD, const Array< OneD, NekDouble > > &locVars, const Array< OneD, NekDouble > &locmu, const Array< OneD, const Array< OneD, NekDouble > > &locnormal, DNekMatSharedPtr &wspMat, Array< OneD, Array< OneD, NekDouble > > &PntJacArray) override
void v_CalcPhysDeriv(const Array< OneD, const Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &qfield) override
void v_CalcMuDmuDT(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &mu, Array< OneD, NekDouble > &DmuDT) override
int m_spacedim
Spatial dimension (>= expansion dim).
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
SOLVER_UTILS_EXPORT int GetTotPoints()
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Base class for unsteady solvers.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
EquationSystemFactory & GetEquationSystemFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::vector< double > w(NPUPPER)
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
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 Zero(int n, T *x, const int incx)
Zero vector.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
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.