35#include <boost/algorithm/string.hpp>
74 int expdim =
m_graph->GetMeshDimension();
84 ASSERTL0(
false,
"Last field is defined as pressure but this is not "
85 "suitable for this solver, please remove this field as "
86 "it is implicitly defined");
96 "components of velocity variables");
98 m_fields[0]->GetHomogeneousBasis()->GetBasisKey();
105 "Non binary number of planes have been specified");
129 "For single mode calculation can only have nz <= 2");
130 if (
m_session->DefinesSolverInfo(
"BetaZero"))
160 for (i = 2; i < nz; ++i)
167 else if (expdim == 3)
171 ASSERTL0(
false,
"Setup mapping aray");
175 ASSERTL0(
false,
"Exp dimension not recognised");
181 std::string vExtrapolation =
"Standard";
183 if (
m_session->DefinesSolverInfo(
"Extrapolation"))
185 vExtrapolation =
m_session->GetSolverInfo(
"Extrapolation");
202 bool IsLinearNSEquation)
214 m_session->LoadParameter(
"imagShift", lambda_imag,
220 "single mode linear NS solver");
254 for (n = 1; n < nz; ++n)
404 bool IsLinearNSEquation,
const int HomogeneousMode,
414 bool AddAdvectionTerms =
455 for (n = 0; n < nel; ++n)
457 nsize_bndry[n] = nvel *
460 nsize_bndry_p1[n] = nsize_bndry[n] + nz_loc;
464 nsize_p[n] =
m_pressure->GetExp(n)->GetNcoeffs() * nz_loc;
465 nsize_p_m1[n] = nsize_p[n] - nz_loc;
471 nsize_bndry_p1, nsize_bndry_p1, blkmatStorage);
473 nsize_bndry, nsize_int, blkmatStorage);
475 nsize_bndry, nsize_int, blkmatStorage);
477 nsize_int, nsize_int, blkmatStorage);
480 nsize_p, nsize_bndry, blkmatStorage);
483 nsize_p, nsize_int, blkmatStorage);
488 nsize_bndry_p1, nsize_p_m1, blkmatStorage);
491 nsize_p_m1, nsize_bndry_p1, blkmatStorage);
499 for (n = 0; n < nel; ++n)
501 nbndry = nsize_bndry[n];
503 k = nsize_bndry_p1[n];
519 nsize_p[n], nsize_bndry[n], zero);
521 nsize_p[n], nsize_int[n], zero);
524 locExp->GetBoundaryMap(bmap);
525 locExp->GetInteriorMap(imap);
533 size_t nbmap = bmap.size();
534 size_t nimap = imap.size();
538 size_t psize =
m_pressure->GetExp(n)->GetNcoeffs();
539 size_t pqsize =
m_pressure->GetExp(n)->GetTotPoints();
544 if (AddAdvectionTerms ==
false)
550 CondMat = locExp->GetLocStaticCondMatrix(helmkey);
552 for (k = 0; k < nvel * nz_loc; ++k)
555 rows = SubBlock.GetRows();
556 cols = SubBlock.GetColumns();
557 for (i = 0; i < rows; ++i)
559 for (j = 0; j < cols; ++j)
561 (*Ah)(i + k * rows, j + k * cols) =
567 for (k = 0; k < nvel * nz_loc; ++k)
571 rows = SubBlock.GetRows();
572 cols = SubBlock.GetColumns();
573 for (i = 0; i < rows; ++i)
575 for (j = 0; j < cols; ++j)
577 (*B)(i + k * rows, j + k * cols) = SubBlock(i, j);
578 (*C)(i + k * rows, j + k * cols) =
584 for (k = 0; k < nvel * nz_loc; ++k)
588 rows = SubBlock.GetRows();
589 cols = SubBlock.GetColumns();
590 for (i = 0; i < rows; ++i)
592 for (j = 0; j < cols; ++j)
594 (*D)(i + k * rows, j + k * cols) =
595 inv_kinvis * SubBlock(i, j);
601 for (i = 0; i < bmap.size(); ++i)
605 coeffs[bmap[i]] = 1.0;
609 for (j = 0; j < nvel; ++j)
611 if ((nz_loc == 2) && (j == 2))
617 beta, phys, 1, deriv, 1);
619 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
623 ((nz_loc * j + 1) * bmap.size() + i) *
630 ((nz_loc * j) * bmap.size() + i) *
645 for (k = 0; k < nz_loc; ++k)
648 psize, &(pcoeffs)[0], 1,
650 ((nz_loc * j + k) * bmap.size() + i) *
660 for (i = 0; i < imap.size(); ++i)
664 coeffs[imap[i]] = 1.0;
668 for (j = 0; j < nvel; ++j)
670 if ((nz_loc == 2) && (j == 2))
676 beta, phys, 1, deriv, 1);
678 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
682 ((nz_loc * j + 1) * imap.size() + i) *
689 ((nz_loc * j) * imap.size() + i) *
708 for (k = 0; k < nz_loc; ++k)
711 psize, &(pcoeffs)[0], 1,
713 ((nz_loc * j + k) * imap.size() + i) *
734 HelmMat.GetOwnedMatrix()->GetPtr();
735 NekDouble HelmMatScale = HelmMat.Scale();
736 int HelmMatRows = HelmMat.GetRows();
752 int npoints = locExp->GetTotPoints();
755 if (IsLinearNSEquation)
760 for (nv = 0; nv < nvel; ++nv)
762 for (nv1 = 0; nv1 < nvel; ++nv1)
764 if (cnt < nvel * nvel - 1)
766 AdvDeriv[cnt + 1] = AdvDeriv[cnt] + npoints;
779 Advfield[nv] + phys_offset,
780 AdvDeriv[nv * nvel + nv1]);
786 for (i = 0; i < nbmap; ++i)
791 coeffs[bmap[i]] = 1.0;
792 locExp->BwdTrans(coeffs, phys);
794 for (k = 0; k < nvel * nz_loc; ++k)
796 for (j = 0; j < nbmap; ++j)
801 Ah_data[i + k * nbmap + (j + k * nbmap) * AhRows] +=
803 HelmMat_data[bmap[i] + HelmMatRows * bmap[j]];
806 for (j = 0; j < nimap; ++j)
808 B_data[i + k * nbmap + (j + k * nimap) * nbndry] +=
810 HelmMat_data[bmap[i] + HelmMatRows * imap[j]];
817 for (k = 0; k < nvel; ++k)
819 for (j = 0; j < nbmap; ++j)
821 Ah_data[i + 2 * k * nbmap +
822 (j + (2 * k + 1) * nbmap) * AhRows] -=
823 lambda_imag * (*MassMat)(bmap[i], bmap[j]);
826 for (j = 0; j < nbmap; ++j)
828 Ah_data[i + (2 * k + 1) * nbmap +
829 (j + 2 * k * nbmap) * AhRows] +=
830 lambda_imag * (*MassMat)(bmap[i], bmap[j]);
833 for (j = 0; j < nimap; ++j)
835 B_data[i + 2 * k * nbmap +
836 (j + (2 * k + 1) * nimap) * nbndry] -=
837 lambda_imag * (*MassMat)(bmap[i], imap[j]);
840 for (j = 0; j < nimap; ++j)
842 B_data[i + (2 * k + 1) * nbmap +
843 (j + 2 * k * nimap) * nbndry] +=
844 lambda_imag * (*MassMat)(bmap[i], imap[j]);
849 for (k = 0; k < nvel; ++k)
851 if ((nz_loc == 2) && (k == 2))
858 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
861 ((nz_loc * k + 1) * bmap.size() + i) *
869 ((nz_loc * k) * bmap.size() + i) *
875 Vmath::Vmul(npoints, Advtmp = Advfield[k] + phys_offset,
876 1, deriv, 1, tmpphys, 1);
878 locExp->IProductWRTBase(tmpphys, coeffs);
881 for (nv = 0; nv < nvel; ++nv)
883 for (j = 0; j < nbmap; ++j)
885 Ah_data[j + 2 * nv * nbmap +
886 (i + (2 * nv + 1) * nbmap) * AhRows] +=
890 for (j = 0; j < nimap; ++j)
892 C_data[i + (2 * nv + 1) * nbmap +
893 (j + 2 * nv * nimap) * nbndry] +=
900 for (nv = 0; nv < nvel; ++nv)
902 for (j = 0; j < nbmap; ++j)
904 Ah_data[j + (2 * nv + 1) * nbmap +
905 (i + 2 * nv * nbmap) * AhRows] +=
909 for (j = 0; j < nimap; ++j)
911 C_data[i + 2 * nv * nbmap +
912 (j + (2 * nv + 1) * nimap) * nbndry] +=
924 Advtmp = Advfield[k] + phys_offset, 1,
925 deriv, 1, tmpphys, 1);
926 locExp->IProductWRTBase(tmpphys, coeffs);
928 for (nv = 0; nv < nvel * nz_loc; ++nv)
930 for (j = 0; j < nbmap; ++j)
932 Ah_data[j + nv * nbmap +
933 (i + nv * nbmap) * AhRows] +=
937 for (j = 0; j < nimap; ++j)
939 C_data[i + nv * nbmap +
940 (j + nv * nimap) * nbndry] +=
948 for (j = 0; j < nz_loc; ++j)
951 psize, &(pcoeffs)[0], 1,
953 ((nz_loc * k + j) * bmap.size() + i) *
961 if (IsLinearNSEquation)
963 for (nv = 0; nv < nvel; ++nv)
967 AdvDeriv[k * nvel + nv], 1, tmpphys, 1);
968 locExp->IProductWRTBase(tmpphys, coeffs);
970 for (
size_t n1 = 0; n1 < nz_loc; ++n1)
972 for (j = 0; j < nbmap; ++j)
974 Ah_data[j + (k * nz_loc + n1) * nbmap +
975 (i + (nv * nz_loc + n1) * nbmap) *
976 AhRows] += coeffs[bmap[j]];
979 for (j = 0; j < nimap; ++j)
981 C_data[i + (nv * nz_loc + n1) * nbmap +
982 (j + (k * nz_loc + n1) * nimap) *
983 nbndry] += coeffs[imap[j]];
991 for (i = 0; i < nimap; ++i)
995 coeffs[imap[i]] = 1.0;
996 locExp->BwdTrans(coeffs, phys);
998 for (k = 0; k < nvel * nz_loc; ++k)
1000 for (j = 0; j < nbmap; ++j)
1002 C_data[j + k * nbmap + (i + k * nimap) * nbndry] +=
1004 HelmMat_data[imap[i] + HelmMatRows * bmap[j]];
1007 for (j = 0; j < nimap; ++j)
1009 D_data[i + k * nimap + (j + k * nimap) * nint] +=
1011 HelmMat_data[imap[i] + HelmMatRows * imap[j]];
1018 for (k = 0; k < nvel; ++k)
1020 for (j = 0; j < nbmap; ++j)
1022 C_data[j + 2 * k * nbmap +
1023 (i + (2 * k + 1) * nimap) * nbndry] +=
1024 lambda_imag * (*MassMat)(bmap[j], imap[i]);
1027 for (j = 0; j < nbmap; ++j)
1029 C_data[j + (2 * k + 1) * nbmap +
1030 (i + 2 * k * nimap) * nbndry] -=
1031 lambda_imag * (*MassMat)(bmap[j], imap[i]);
1034 for (j = 0; j < nimap; ++j)
1036 D_data[i + 2 * k * nimap +
1037 (j + (2 * k + 1) * nimap) * nint] -=
1038 lambda_imag * (*MassMat)(imap[i], imap[j]);
1041 for (j = 0; j < nimap; ++j)
1043 D_data[i + (2 * k + 1) * nimap +
1044 (j + 2 * k * nimap) * nint] +=
1045 lambda_imag * (*MassMat)(imap[i], imap[j]);
1050 for (k = 0; k < nvel; ++k)
1052 if ((nz_loc == 2) && (k == 2))
1058 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
1061 ((nz_loc * k + 1) * imap.size() + i) *
1068 ((nz_loc * k) * imap.size() + i) *
1075 Vmath::Vmul(npoints, Advtmp = Advfield[k] + phys_offset,
1076 1, deriv, 1, tmpphys, 1);
1077 locExp->IProductWRTBase(tmpphys, coeffs);
1080 for (nv = 0; nv < nvel; ++nv)
1082 for (j = 0; j < nbmap; ++j)
1084 B_data[j + 2 * nv * nbmap +
1085 (i + (2 * nv + 1) * nimap) * nbndry] +=
1089 for (j = 0; j < nimap; ++j)
1091 D_data[j + 2 * nv * nimap +
1092 (i + (2 * nv + 1) * nimap) * nint] +=
1098 for (nv = 0; nv < nvel; ++nv)
1100 for (j = 0; j < nbmap; ++j)
1102 B_data[j + (2 * nv + 1) * nbmap +
1103 (i + 2 * nv * nimap) * nbndry] +=
1107 for (j = 0; j < nimap; ++j)
1109 D_data[j + (2 * nv + 1) * nimap +
1110 (i + 2 * nv * nimap) * nint] +=
1124 Advtmp = Advfield[k] + phys_offset, 1,
1125 deriv, 1, tmpphys, 1);
1126 locExp->IProductWRTBase(tmpphys, coeffs);
1128 for (nv = 0; nv < nvel * nz_loc; ++nv)
1130 for (j = 0; j < nbmap; ++j)
1132 B_data[j + nv * nbmap +
1133 (i + nv * nimap) * nbndry] +=
1137 for (j = 0; j < nimap; ++j)
1139 D_data[j + nv * nimap +
1140 (i + nv * nimap) * nint] +=
1145 m_pressure->GetExp(n)->IProductWRTBase(deriv,
1147 for (j = 0; j < nz_loc; ++j)
1150 psize, &(pcoeffs)[0], 1,
1152 ((nz_loc * k + j) * imap.size() + i) *
1160 if (IsLinearNSEquation)
1163 for (nv = 0; nv < nvel; ++nv)
1167 AdvDeriv[k * nvel + nv], 1, tmpphys, 1);
1168 locExp->IProductWRTBase(tmpphys, coeffs);
1170 for (n1 = 0; n1 < nz_loc; ++n1)
1172 for (j = 0; j < nbmap; ++j)
1174 B_data[j + (k * nz_loc + n1) * nbmap +
1175 (i + (nv * nz_loc + n1) * nimap) *
1176 nbndry] += coeffs[bmap[j]];
1179 for (j = 0; j < nimap; ++j)
1181 D_data[j + (k * nz_loc + n1) * nimap +
1182 (i + (nv * nz_loc + n1) * nimap) *
1183 nint] += coeffs[imap[j]];
1197 Blas::Dgemm(
'N',
'T', B->GetRows(), C->GetRows(), B->GetColumns(),
1198 -1.0, B->GetRawPtr(), B->GetRows(), C->GetRawPtr(),
1199 C->GetRows(), 1.0, Ah->GetRawPtr(), Ah->GetRows());
1222 DNekMat DintCinvDTint, BCinvDTint_m_DTbnd, DintCinvBTtilde_m_Dbnd;
1228 DintCinvDTint = (*Dint) * (*Cinv) *
Transpose(*Dint);
1232 DintCinvBTtilde_m_Dbnd =
1233 (*Dint) * (*Cinv) *
Transpose(*Btilde) - (*Dbnd);
1237 nsize_bndry_p1[n], nsize_p_m1[n], zero);
1239 nsize_p_m1[n], nsize_bndry_p1[n], zero);
1241 nsize_p_m1[n], nsize_p_m1[n], zero);
1246 for (n1 = 0; n1 < nz_loc; ++n1)
1248 for (i = 0; i < psize - 1; ++i)
1250 for (n2 = 0; n2 < nz_loc; ++n2)
1252 for (j = 0; j < psize - 1; ++j)
1256 Dh_data[(i + n1 * (psize - 1)) +
1257 (j + n2 * (psize - 1)) * Dh->GetRows()] =
1258 -DintCinvDTint(i + 1 + n1 * psize,
1259 j + 1 + n2 * psize);
1265 for (n1 = 0; n1 < nz_loc; ++n1)
1267 for (i = 0; i < nsize_bndry_p1[n] - nz_loc; ++i)
1269 (*Ah)(i, nsize_bndry_p1[n] - nz_loc + n1) =
1270 BCinvDTint_m_DTbnd(i, n1 * psize);
1271 (*Ah)(nsize_bndry_p1[n] - nz_loc + n1, i) =
1272 DintCinvBTtilde_m_Dbnd(n1 * psize, i);
1276 for (n1 = 0; n1 < nz_loc; ++n1)
1278 for (n2 = 0; n2 < nz_loc; ++n2)
1280 (*Ah)(nsize_bndry_p1[n] - nz_loc + n1,
1281 nsize_bndry_p1[n] - nz_loc + n2) =
1282 -DintCinvDTint(n1 * psize, n2 * psize);
1286 for (n1 = 0; n1 < nz_loc; ++n1)
1288 for (j = 0; j < psize - 1; ++j)
1290 for (i = 0; i < nsize_bndry_p1[n] - nz_loc; ++i)
1292 (*Bh)(i, j + n1 * (psize - 1)) =
1293 BCinvDTint_m_DTbnd(i, j + 1 + n1 * psize);
1294 (*Ch)(j + n1 * (psize - 1), i) =
1295 DintCinvBTtilde_m_Dbnd(j + 1 + n1 * psize, i);
1300 for (n1 = 0; n1 < nz_loc; ++n1)
1302 for (n2 = 0; n2 < nz_loc; ++n2)
1304 for (j = 0; j < psize - 1; ++j)
1306 (*Bh)(nsize_bndry_p1[n] - nz_loc + n1,
1307 j + n2 * (psize - 1)) =
1308 -DintCinvDTint(n1 * psize, j + 1 + n2 * psize);
1309 (*Ch)(j + n2 * (psize - 1),
1310 nsize_bndry_p1[n] - nz_loc + n1) =
1311 -DintCinvDTint(j + 1 + n2 * psize, n1 * psize);
1318 (*Bh) = (*Bh) * (*Dh);
1320 Blas::Dgemm(
'N',
'N', Bh->GetRows(), Ch->GetColumns(), Bh->GetColumns(),
1321 -1.0, Bh->GetRawPtr(), Bh->GetRows(), Ch->GetRawPtr(),
1322 Ch->GetRows(), 1.0, Ah->GetRawPtr(), Ah->GetRows());
1340 std::cout <<
"Matrix Setup Costs: " << timer.
TimePerTest(1) << std::endl;
1407 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1414 "Advection Velocity section must be defined in "
1417 std::vector<std::string> fieldStr;
1418 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1423 GetFunction(
"AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1442 "Restart section must be defined in session file.");
1445 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1450 std::vector<std::string> fieldStr;
1451 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1456 GetFunction(
"Restart")->Evaluate(fieldStr, Restart);
1458 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1463 std::cout <<
"Saving the RESTART file for m_kinvis = "
1476 std::cout <<
"Saving the Stokes Flow for m_kinvis = "
1487 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1494 "Advection Velocity section must be defined in "
1497 std::vector<std::string> fieldStr;
1498 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1503 GetFunction(
"AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1511 "Unknown or undefined equation type for CoupledLinearNS");
1524 x->Apply(
m_fields, outarray, outarray, time);
1539 if (fabs(lambda_store - lambda) > 1e-10)
1542 lambda_store = lambda;
1571 for (
size_t k = 0; k < nfields; ++k)
1582 for (
size_t k = 0; k < nfields; ++k)
1609 Generaltimer.
Start();
1618 <<
"We execute INITIALLY SolveSteadyNavierStokes for m_kinvis "
1629 <<
"We execute SolveSteadyNavierStokes for m_kinvis = "
1642 <<
"We execute SolveSteadyNavierStokes for m_kinvis = "
1651 Generaltimer.
Stop();
1652 std::cout <<
"\nThe total calculation time is : "
1653 << Generaltimer.
TimePerTest(1) / 60 <<
" minute(s). \n\n";
1660 "Unknown or undefined equation type for CoupledLinearNS");
1680 for (
size_t i = 0; i < ncmpt; ++i)
1691 x->Apply(
m_fields, forcing_phys, forcing_phys, time);
1693 for (
size_t i = 0; i < ncmpt; ++i)
1697 m_fields[i]->IProductWRTBase(forcing_phys[i], forcing[i]);
1698 m_fields[i]->SetWaveSpace(waveSpace);
1710 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1718 if (
m_session->DefinesFunction(
"ForcingTerm"))
1720 std::vector<std::string> fieldStr;
1721 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1727 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1736 <<
"'ForcingTerm' section has not been defined in the input file "
1745 Newtontimer.
Start();
1754 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1764 L2Norm(delta_velocity_Phys, L2_norm);
1767 while (std::max(L2_norm[0], L2_norm[1]) >
m_tol)
1775 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1783 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1808 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1815 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1817 Vmath::Vadd(Velocity_Phys[i].size(), Velocity_Phys[i], 1,
1818 delta_velocity_Phys[i], 1, Velocity_Phys[i], 1);
1822 L2Norm(delta_velocity_Phys, L2_norm);
1824 if (std::max(Inf_norm[0], Inf_norm[1]) > 100)
1826 std::cout <<
"\nThe Newton method has failed at m_kinvis = "
1829 ASSERTL0(0,
"The Newton method has failed... \n");
1838 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1846 std::cout <<
"We have done " <<
m_counter - 1 <<
" iteration(s) in "
1847 << Newtontimer.
TimePerTest(1) / 60 <<
" minute(s). \n\n";
1857 std::cout <<
"We apply the continuation method: " << std::endl;
1859 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1883 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1892 Vmath::Vadd(u_star[i].size(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1904 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1907 for (
size_t j = 0; j < inarray[i].size(); ++j)
1909 if (inarray[i][j] > outarray[i])
1911 outarray[i] = inarray[i][j];
1914 std::cout <<
"InfNorm[" << i <<
"] = " << outarray[i] << std::endl;
1921 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1924 for (
size_t j = 0; j < inarray[i].size(); ++j)
1926 outarray[i] += inarray[i][j] * inarray[i][j];
1929 std::cout <<
"L2Norm[" << i <<
"] = " << outarray[i] << std::endl;
1943 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1967 for (
size_t i = 0; i <
m_velocity.size(); ++i)
1974 i, tmp_DerVel[i], ViscTerm[i]);
1979 Vmath::Vsub(outarray[i].size(), outarray[i], 1, AdvTerm[i], 1,
1981 Vmath::Vsub(outarray[i].size(), outarray[i], 1, ViscTerm[i], 1,
1984 Vmath::Vadd(outarray[i].size(), outarray[i], 1, Forc[i], 1, outarray[i],
1998 for (
auto &expMapIter : VelExp)
2002 for (
size_t i = 0; i < expMapIter.second->m_basisKeyVector.size(); ++i)
2007 "Velocity polynomial space not sufficiently high (>= 4)");
2011 BasisVec.push_back(newB);
2017 expMapIter.second->m_geomShPtr, BasisVec);
2018 (*returnval)[expMapIter.first] = expansionElementShPtr;
2022 m_graph->SetExpansionInfo(
"p", returnval);
2115 force[i] = forcing[i] + 2 * n * ncoeffsplane;
2132 force[i] = forcing[i];
2143 size_t i, j, k, n, cnt, cnt1;
2144 size_t nbnd, nint, offset;
2146 size_t nel = fields[0]->GetNumElmts();
2176 for (i = 0; i < fields.size(); ++i)
2179 tmp = fields[i]->UpdateCoeffs() + nplanecoeffs,
2188 for (k = 0; k < nvel; ++k)
2191 std::dynamic_pointer_cast<MultiRegions::ContField>(fields[k]);
2194 cfield->GetLocalToGlobalMap()->GetBndCondCoeffsToLocalCoeffsSign();
2196 cfield->GetLocalToGlobalMap()->GetBndCondCoeffsToLocalCoeffsMap();
2200 fields[k]->GetBndConditions();
2206 m_fields[k]->GetPlane(2 * mode)->GetBndCondExpansions();
2210 bndCondExp =
m_fields[k]->GetBndCondExpansions();
2213 for (n = 0; n < nz_loc; ++n)
2216 for (i = 0; i < bndCondExp.size(); ++i)
2219 bndCondExp[i]->GetCoeffs();
2220 size_t nCoeffs = (bndCondExp[i])->
GetNcoeffs();
2223 if (bndConds[i]->GetBoundaryConditionType() ==
2225 bndConds[i]->GetBoundaryConditionType() ==
2230 for (j = 0; j < nCoeffs; j++)
2232 forcing[k][n * nplanecoeffs + map[bndcnt + j]] +=
2233 sign[bndcnt + j] * bndcoeffs[j];
2238 for (j = 0; j < nCoeffs; j++)
2240 forcing[k][n * nplanecoeffs + map[bndcnt + j]] +=
2246 bndcnt += bndCondExp[i]->GetNcoeffs();
2257 for (i = 0; i < nel; ++i)
2259 fields[
m_velocity[0]]->GetExp(i)->GetBoundaryMap(bmap);
2260 fields[
m_velocity[0]]->GetExp(i)->GetInteriorMap(imap);
2263 offset = fields[
m_velocity[0]]->GetCoeff_Offset(i);
2265 for (j = 0; j < nvel; ++j)
2269 for (n = 0; n < nz_loc; ++n)
2271 for (k = 0; k < nbnd; ++k)
2274 forcing[j][n * nplanecoeffs + offset + bmap[k]];
2276 bnd[bndoffset + (n + j * nz_loc) * nbnd + k] =
2277 incoeffs[n * nplanecoeffs + offset + bmap[k]];
2279 for (k = 0; k < nint; ++k)
2282 forcing[j][n * nplanecoeffs + offset + imap[k]];
2290 nvel * nz_loc * nbnd + nz_loc * (
pressure->GetExp(i)->GetNcoeffs());
2298 F_bnd = F_bnd - (*
m_mat[mode].m_BCinv) * F_int;
2299 F_p_tmp = (*
m_mat[mode].m_Cinv) * F_int;
2300 F_p = (*
m_mat[mode].m_D_int) * F_p_tmp;
2307 for (i = 0; i < nel; ++i)
2309 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2311 for (j = 0; j < nvel; ++j)
2313 for (k = 0; k < nbnd; ++k)
2315 fh_bnd[offset + j * nbnd + k] = f_bnd[cnt + k];
2320 nint =
pressure->GetExp(i)->GetNcoeffs();
2321 offset += nvel * nbnd + nint * nz_loc;
2325 for (i = 0; i < nel; ++i)
2327 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2328 nint =
pressure->GetExp(i)->GetNcoeffs();
2330 for (n = 0; n < nz_loc; ++n)
2332 for (j = 0; j < nint; ++j)
2334 fh_bnd[offset + nvel * nbnd + n * nint + j] = f_p[cnt1 + j];
2338 offset += nvel * nbnd + nz_loc * nint;
2344 int totpcoeffs =
pressure->GetNcoeffs();
2346 for (i = 0; i < nel; ++i)
2348 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2349 nint =
pressure->GetExp(i)->GetNcoeffs();
2350 for (j = 0; j < nvel; ++j)
2352 for (k = 0; k < nbnd; ++k)
2354 f_bnd[cnt + k] = bnd[offset + j * nbnd + k];
2358 offset += nvel * nbnd + nint * nz_loc;
2363 offset = cnt = cnt1 = 0;
2364 for (i = 0; i < nel; ++i)
2366 nint =
pressure->GetExp(i)->GetNcoeffs();
2367 nbnd = fields[0]->GetExp(i)->NumBndryCoeffs();
2368 cnt1 =
pressure->GetCoeff_Offset(i);
2370 for (n = 0; n < nz_loc; ++n)
2372 for (j = 0; j < nint; ++j)
2374 p_coeffs[n * totpcoeffs + cnt1 + j] = f_p[cnt + j] =
2375 bnd[offset + nvel * nz_loc * nbnd + n * nint + j];
2379 offset += (nvel * nbnd + nint) * nz_loc;
2386 F_int = (*
m_mat[mode].m_Cinv) * F_int;
2390 for (i = 0; i < nel; ++i)
2392 fields[0]->GetExp(i)->GetBoundaryMap(bmap);
2393 fields[0]->GetExp(i)->GetInteriorMap(imap);
2396 offset = fields[0]->GetCoeff_Offset(i);
2398 for (j = 0; j < nvel; ++j)
2400 for (n = 0; n < nz_loc; ++n)
2402 for (k = 0; k < nbnd; ++k)
2404 fields[j]->SetCoeff(n * nplanecoeffs + offset + bmap[k],
2408 for (k = 0; k < nint; ++k)
2410 fields[j]->SetCoeff(n * nplanecoeffs + offset + imap[k],
2419 for (j = 0; j < nvel; ++j)
2421 fields[j]->SetPhysState(
false);
2427 std::vector<Array<OneD, NekDouble>> fieldcoeffs(
m_fields.size() + 1);
2428 std::vector<std::string> variables(
m_fields.size() + 1);
2431 for (i = 0; i <
m_fields.size(); ++i)
2433 fieldcoeffs[i] =
m_fields[i]->UpdateCoeffs();
2448 m_fields[0]->GetPlane(0)->FwdTransLocalElmt(
2449 m_pressure->GetPlane(0)->GetPhys(), fieldcoeffs[i]);
2450 m_fields[0]->GetPlane(1)->FwdTransLocalElmt(
2451 m_pressure->GetPlane(1)->GetPhys(), tmpfieldcoeffs);
2452 for (
int e = 0; e <
m_fields[0]->GetNcoeffs() / 2; e++)
2454 fieldcoeffs[i][e +
m_fields[0]->GetNcoeffs() / 2] =
2465 if (!
m_comm->IsParallelInTime())
2472 if (!fs::is_directory(newdir))
2474 fs::create_directory(newdir);
2479 m_comm->GetTimeComm()->GetRank() + 1) +
2481 m_fields[0], fieldcoeffs, variables);
2487 return m_session->GetVariables().size();
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
#define sign(a, b)
return the sign(b)*a
void L2Norm(Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
CoupledLinearNS(const LibUtilities::SessionReaderSharedPtr &pSesssion, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print a summary of time stepping parameters.
void EvaluateNewtonRHS(Array< OneD, Array< OneD, NekDouble > > &Velocity, Array< OneD, Array< OneD, NekDouble > > &outarray)
void v_DoSolve(void) override
Solves an unsteady problem.
Array< OneD, CoupledLocalToGlobalC0ContMapSharedPtr > m_locToGloMap
void DefineForcingTerm(void)
void SolveUnsteadyStokesSystem(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time, const NekDouble a_iixDt)
void SolveLinearNS(const Array< OneD, Array< OneD, NekDouble > > &forcing)
void v_Output(void) override
void SolveSteadyNavierStokes(void)
void v_DoInitialise(bool dumpInitialConditions=true) override
Sets up initial conditions.
static std::string solverTypeLookupId
const SpatialDomains::ExpansionInfoMap & GenPressureExp(const SpatialDomains::ExpansionInfoMap &VelExp)
bool v_NegatedOp(void) override
void EvaluateAdvection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
bool m_zeroMode
Id to identify when single mode is mean mode (i.e. beta=0);.
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm_Coeffs
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm
void v_TransPhysToCoeff(void) override
Virtual function for transformation to coefficient space.
void SetUpCoupledMatrix(const NekDouble lambda=0.0, const Array< OneD, Array< OneD, NekDouble > > &Advfield=NullNekDoubleArrayOfArray, bool IsLinearNSEquation=true)
static std::string className
Name of class.
int v_GetForceDimension(void) override
NekDouble m_KinvisPercentage
void InfNorm(Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
void v_TransCoeffToPhys(void) override
Virtual function for transformation to physical space.
Array< OneD, CoupledSolverMatrices > m_mat
This class is the base class for Navier Stokes problems.
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
NekDouble m_kinvis
Kinematic viscosity.
void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
ExtrapolateSharedPtr m_extrapolation
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
EquationType m_equationType
equation type;
int m_nConvectiveFields
Number of fields to be convected;.
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Describes the specification for a Basis.
BasisType GetBasisType() const
Return type of expansion basis.
PointsKey GetPointsKey() const
Return distribution of points.
int GetNumModes() const
Returns the order of the basis.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Describe a linear system.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
bool m_useFFT
Flag to determine if FFT is used for homogeneous transform.
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
LibUtilities::CommSharedPtr m_comm
Communicator.
NekDouble m_time
Current time of simulation.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
int m_windowPIT
Index of windows for parallel-in-time time iteration.
SOLVER_UTILS_EXPORT int GetNpoints()
NekDouble m_lambda
Lambda constant in real system if one required.
SOLVER_UTILS_EXPORT int GetNcoeffs()
int m_npointsZ
number of points in Z direction (if homogeneous)
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
std::string m_sessionName
Name of the session.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
bool m_singleMode
Flag to determine if single homogeneous mode is used.
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
enum HomogeneousType m_HomogeneousType
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
SOLVER_UTILS_EXPORT int GetTotPoints()
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
SOLVER_UTILS_EXPORT void v_DoSolve() override
Solves an unsteady problem.
static void Dcopy(const int &n, const double *x, const int &incx, double *y, const int &incy)
BLAS level 1: Copy x to y.
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
std::vector< BasisKey > BasisKeyVector
Name for a vector of BasisKeys.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ beta
Gauss Radau pinned at x=-1,.
MultiRegions::Direction const DirCartesianMap[]
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< ContField > ContFieldSharedPtr
static const NekDouble kNekUnsetDouble
std::vector< std::pair< std::string, std::string > > SummaryList
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
std::shared_ptr< ExpansionInfoMap > ExpansionInfoMapShPtr
std::shared_ptr< ExpansionInfo > ExpansionInfoShPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::map< int, ExpansionInfoShPtr > ExpansionInfoMap
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
@ eLinearAdvectionReaction
std::map< ConstFactorType, NekDouble > ConstFactorMap
StdRegions::ConstFactorMap factors
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
ExtrapolateFactory & GetExtrapolateFactory()
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< CoupledLocalToGlobalC0ContMap > CoupledLocalToGlobalC0ContMapSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
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 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 Zero(int n, T *x, const int incx)
Zero vector.
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.
scalarT< T > sqrt(scalarT< T > in)
DNekScalBlkMatSharedPtr m_D_int
Inner product of pressure system with divergence of the interior velocity space .
DNekScalBlkMatSharedPtr m_D_bnd
Inner product of pressure system with divergence of the boundary velocity space .
DNekScalBlkMatSharedPtr m_Cinv
Interior-Interior Laplaican plus linearised convective terms inverted, i.e. the inverse of .
DNekScalBlkMatSharedPtr m_Btilde
Interior-boundary Laplacian plus linearised convective terms .
DNekScalBlkMatSharedPtr m_BCinv
Boundary-interior Laplacian plus linearised convective terms pre-multiplying Cinv: .
MultiRegions::GlobalLinSysSharedPtr m_CoupledBndSys