35 #include <boost/algorithm/string.hpp>
48 string CoupledLinearNS::className =
50 "CoupledLinearisedNS", CoupledLinearNS::create);
59 CoupledLinearNS::CoupledLinearNS(
72 int expdim =
m_graph->GetMeshDimension();
82 ASSERTL0(
false,
"Last field is defined as pressure but this is not "
83 "suitable for this solver, please remove this field as "
84 "it is implicitly defined");
94 "components of velocity variables");
96 m_fields[0]->GetHomogeneousBasis()->GetBasisKey();
103 "Non binary number of planes have been specified");
127 "For single mode calculation can only have nz <= 2");
128 if (
m_session->DefinesSolverInfo(
"BetaZero"))
158 for (i = 2; i < nz; ++i)
165 else if (expdim == 3)
169 ASSERTL0(
false,
"Setup mapping aray");
173 ASSERTL0(
false,
"Exp dimension not recognised");
179 std::string vExtrapolation =
"Standard";
181 if (
m_session->DefinesSolverInfo(
"Extrapolation"))
183 vExtrapolation =
m_session->GetSolverInfo(
"Extrapolation");
200 bool IsLinearNSEquation)
212 m_session->LoadParameter(
"imagShift", lambda_imag,
218 "single mode linear NS solver");
252 for (n = 1; n < nz; ++n)
402 bool IsLinearNSEquation,
const int HomogeneousMode,
412 bool AddAdvectionTerms =
453 for (n = 0; n < nel; ++n)
455 nsize_bndry[n] = nvel *
458 nsize_bndry_p1[n] = nsize_bndry[n] + nz_loc;
462 nsize_p[n] =
m_pressure->GetExp(n)->GetNcoeffs() * nz_loc;
463 nsize_p_m1[n] = nsize_p[n] - nz_loc;
469 nsize_bndry_p1, nsize_bndry_p1, blkmatStorage);
471 nsize_bndry, nsize_int, blkmatStorage);
473 nsize_bndry, nsize_int, blkmatStorage);
475 nsize_int, nsize_int, blkmatStorage);
478 nsize_p, nsize_bndry, blkmatStorage);
481 nsize_p, nsize_int, blkmatStorage);
486 nsize_bndry_p1, nsize_p_m1, blkmatStorage);
489 nsize_p_m1, nsize_bndry_p1, blkmatStorage);
497 for (n = 0; n < nel; ++n)
499 nbndry = nsize_bndry[n];
501 k = nsize_bndry_p1[n];
517 nsize_p[n], nsize_bndry[n], zero);
519 nsize_p[n], nsize_int[n], zero);
522 locExp->GetBoundaryMap(bmap);
523 locExp->GetInteriorMap(imap);
531 int nbmap = bmap.size();
532 int nimap = imap.size();
536 int psize =
m_pressure->GetExp(n)->GetNcoeffs();
537 int pqsize =
m_pressure->GetExp(n)->GetTotPoints();
542 if (AddAdvectionTerms ==
false)
548 CondMat = locExp->GetLocStaticCondMatrix(helmkey);
550 for (k = 0; k < nvel * nz_loc; ++k)
553 rows = SubBlock.GetRows();
554 cols = SubBlock.GetColumns();
555 for (i = 0; i < rows; ++i)
557 for (j = 0; j < cols; ++j)
559 (*Ah)(i + k * rows, j + k * cols) =
565 for (k = 0; k < nvel * nz_loc; ++k)
569 rows = SubBlock.GetRows();
570 cols = SubBlock.GetColumns();
571 for (i = 0; i < rows; ++i)
573 for (j = 0; j < cols; ++j)
575 (*B)(i + k * rows, j + k * cols) = SubBlock(i, j);
576 (*C)(i + k * rows, j + k * cols) =
582 for (k = 0; k < nvel * nz_loc; ++k)
586 rows = SubBlock.GetRows();
587 cols = SubBlock.GetColumns();
588 for (i = 0; i < rows; ++i)
590 for (j = 0; j < cols; ++j)
592 (*D)(i + k * rows, j + k * cols) =
593 inv_kinvis * SubBlock(i, j);
599 for (i = 0; i < bmap.size(); ++i)
603 coeffs[bmap[i]] = 1.0;
607 for (j = 0; j < nvel; ++j)
609 if ((nz_loc == 2) && (j == 2))
615 beta, phys, 1, deriv, 1);
617 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
621 ((nz_loc * j + 1) * bmap.size() + i) *
628 ((nz_loc * j) * bmap.size() + i) *
643 for (k = 0; k < nz_loc; ++k)
646 psize, &(pcoeffs)[0], 1,
648 ((nz_loc * j + k) * bmap.size() + i) *
658 for (i = 0; i < imap.size(); ++i)
662 coeffs[imap[i]] = 1.0;
666 for (j = 0; j < nvel; ++j)
668 if ((nz_loc == 2) && (j == 2))
674 beta, phys, 1, deriv, 1);
676 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
680 ((nz_loc * j + 1) * imap.size() + i) *
687 ((nz_loc * j) * imap.size() + i) *
706 for (k = 0; k < nz_loc; ++k)
709 psize, &(pcoeffs)[0], 1,
711 ((nz_loc * j + k) * imap.size() + i) *
732 HelmMat.GetOwnedMatrix()->GetPtr();
733 NekDouble HelmMatScale = HelmMat.Scale();
734 int HelmMatRows = HelmMat.GetRows();
750 int npoints = locExp->GetTotPoints();
753 if (IsLinearNSEquation)
758 for (nv = 0; nv < nvel; ++nv)
760 for (nv1 = 0; nv1 < nvel; ++nv1)
762 if (cnt < nvel * nvel - 1)
764 AdvDeriv[cnt + 1] = AdvDeriv[cnt] + npoints;
777 Advfield[nv] + phys_offset,
778 AdvDeriv[nv * nvel + nv1]);
784 for (i = 0; i < nbmap; ++i)
789 coeffs[bmap[i]] = 1.0;
790 locExp->BwdTrans(coeffs, phys);
792 for (k = 0; k < nvel * nz_loc; ++k)
794 for (j = 0; j < nbmap; ++j)
799 Ah_data[i + k * nbmap + (j + k * nbmap) * AhRows] +=
801 HelmMat_data[bmap[i] + HelmMatRows * bmap[j]];
804 for (j = 0; j < nimap; ++j)
806 B_data[i + k * nbmap + (j + k * nimap) * nbndry] +=
808 HelmMat_data[bmap[i] + HelmMatRows * imap[j]];
815 for (k = 0; k < nvel; ++k)
817 for (j = 0; j < nbmap; ++j)
819 Ah_data[i + 2 * k * nbmap +
820 (j + (2 * k + 1) * nbmap) * AhRows] -=
821 lambda_imag * (*MassMat)(bmap[i], bmap[j]);
824 for (j = 0; j < nbmap; ++j)
826 Ah_data[i + (2 * k + 1) * nbmap +
827 (j + 2 * k * nbmap) * AhRows] +=
828 lambda_imag * (*MassMat)(bmap[i], bmap[j]);
831 for (j = 0; j < nimap; ++j)
833 B_data[i + 2 * k * nbmap +
834 (j + (2 * k + 1) * nimap) * nbndry] -=
835 lambda_imag * (*MassMat)(bmap[i], imap[j]);
838 for (j = 0; j < nimap; ++j)
840 B_data[i + (2 * k + 1) * nbmap +
841 (j + 2 * k * nimap) * nbndry] +=
842 lambda_imag * (*MassMat)(bmap[i], imap[j]);
847 for (k = 0; k < nvel; ++k)
849 if ((nz_loc == 2) && (k == 2))
856 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
859 ((nz_loc * k + 1) * bmap.size() + i) *
867 ((nz_loc * k) * bmap.size() + i) *
873 Vmath::Vmul(npoints, Advtmp = Advfield[k] + phys_offset,
874 1, deriv, 1, tmpphys, 1);
876 locExp->IProductWRTBase(tmpphys, coeffs);
879 for (nv = 0; nv < nvel; ++nv)
881 for (j = 0; j < nbmap; ++j)
883 Ah_data[j + 2 * nv * nbmap +
884 (i + (2 * nv + 1) * nbmap) * AhRows] +=
888 for (j = 0; j < nimap; ++j)
890 C_data[i + (2 * nv + 1) * nbmap +
891 (j + 2 * nv * nimap) * nbndry] +=
898 for (nv = 0; nv < nvel; ++nv)
900 for (j = 0; j < nbmap; ++j)
902 Ah_data[j + (2 * nv + 1) * nbmap +
903 (i + 2 * nv * nbmap) * AhRows] +=
907 for (j = 0; j < nimap; ++j)
909 C_data[i + 2 * nv * nbmap +
910 (j + (2 * nv + 1) * nimap) * nbndry] +=
922 Advtmp = Advfield[k] + phys_offset, 1,
923 deriv, 1, tmpphys, 1);
924 locExp->IProductWRTBase(tmpphys, coeffs);
926 for (nv = 0; nv < nvel * nz_loc; ++nv)
928 for (j = 0; j < nbmap; ++j)
930 Ah_data[j + nv * nbmap +
931 (i + nv * nbmap) * AhRows] +=
935 for (j = 0; j < nimap; ++j)
937 C_data[i + nv * nbmap +
938 (j + nv * nimap) * nbndry] +=
946 for (j = 0; j < nz_loc; ++j)
949 psize, &(pcoeffs)[0], 1,
951 ((nz_loc * k + j) * bmap.size() + i) *
959 if (IsLinearNSEquation)
961 for (nv = 0; nv < nvel; ++nv)
965 AdvDeriv[k * nvel + nv], 1, tmpphys, 1);
966 locExp->IProductWRTBase(tmpphys, coeffs);
968 for (
int n1 = 0; n1 < nz_loc; ++n1)
970 for (j = 0; j < nbmap; ++j)
972 Ah_data[j + (k * nz_loc + n1) * nbmap +
973 (i + (nv * nz_loc + n1) * nbmap) *
974 AhRows] += coeffs[bmap[j]];
977 for (j = 0; j < nimap; ++j)
979 C_data[i + (nv * nz_loc + n1) * nbmap +
980 (j + (k * nz_loc + n1) * nimap) *
981 nbndry] += coeffs[imap[j]];
989 for (i = 0; i < nimap; ++i)
993 coeffs[imap[i]] = 1.0;
994 locExp->BwdTrans(coeffs, phys);
996 for (k = 0; k < nvel * nz_loc; ++k)
998 for (j = 0; j < nbmap; ++j)
1000 C_data[j + k * nbmap + (i + k * nimap) * nbndry] +=
1002 HelmMat_data[imap[i] + HelmMatRows * bmap[j]];
1005 for (j = 0; j < nimap; ++j)
1007 D_data[i + k * nimap + (j + k * nimap) * nint] +=
1009 HelmMat_data[imap[i] + HelmMatRows * imap[j]];
1016 for (k = 0; k < nvel; ++k)
1018 for (j = 0; j < nbmap; ++j)
1020 C_data[j + 2 * k * nbmap +
1021 (i + (2 * k + 1) * nimap) * nbndry] +=
1022 lambda_imag * (*MassMat)(bmap[j], imap[i]);
1025 for (j = 0; j < nbmap; ++j)
1027 C_data[j + (2 * k + 1) * nbmap +
1028 (i + 2 * k * nimap) * nbndry] -=
1029 lambda_imag * (*MassMat)(bmap[j], imap[i]);
1032 for (j = 0; j < nimap; ++j)
1034 D_data[i + 2 * k * nimap +
1035 (j + (2 * k + 1) * nimap) * nint] -=
1036 lambda_imag * (*MassMat)(imap[i], imap[j]);
1039 for (j = 0; j < nimap; ++j)
1041 D_data[i + (2 * k + 1) * nimap +
1042 (j + 2 * k * nimap) * nint] +=
1043 lambda_imag * (*MassMat)(imap[i], imap[j]);
1048 for (k = 0; k < nvel; ++k)
1050 if ((nz_loc == 2) && (k == 2))
1056 m_pressure->GetExp(n)->IProductWRTBase(deriv, pcoeffs);
1059 ((nz_loc * k + 1) * imap.size() + i) *
1066 ((nz_loc * k) * imap.size() + i) *
1073 Vmath::Vmul(npoints, Advtmp = Advfield[k] + phys_offset,
1074 1, deriv, 1, tmpphys, 1);
1075 locExp->IProductWRTBase(tmpphys, coeffs);
1078 for (nv = 0; nv < nvel; ++nv)
1080 for (j = 0; j < nbmap; ++j)
1082 B_data[j + 2 * nv * nbmap +
1083 (i + (2 * nv + 1) * nimap) * nbndry] +=
1087 for (j = 0; j < nimap; ++j)
1089 D_data[j + 2 * nv * nimap +
1090 (i + (2 * nv + 1) * nimap) * nint] +=
1096 for (nv = 0; nv < nvel; ++nv)
1098 for (j = 0; j < nbmap; ++j)
1100 B_data[j + (2 * nv + 1) * nbmap +
1101 (i + 2 * nv * nimap) * nbndry] +=
1105 for (j = 0; j < nimap; ++j)
1107 D_data[j + (2 * nv + 1) * nimap +
1108 (i + 2 * nv * nimap) * nint] +=
1122 Advtmp = Advfield[k] + phys_offset, 1,
1123 deriv, 1, tmpphys, 1);
1124 locExp->IProductWRTBase(tmpphys, coeffs);
1126 for (nv = 0; nv < nvel * nz_loc; ++nv)
1128 for (j = 0; j < nbmap; ++j)
1130 B_data[j + nv * nbmap +
1131 (i + nv * nimap) * nbndry] +=
1135 for (j = 0; j < nimap; ++j)
1137 D_data[j + nv * nimap +
1138 (i + nv * nimap) * nint] +=
1143 m_pressure->GetExp(n)->IProductWRTBase(deriv,
1145 for (j = 0; j < nz_loc; ++j)
1148 psize, &(pcoeffs)[0], 1,
1150 ((nz_loc * k + j) * imap.size() + i) *
1158 if (IsLinearNSEquation)
1161 for (nv = 0; nv < nvel; ++nv)
1165 AdvDeriv[k * nvel + nv], 1, tmpphys, 1);
1166 locExp->IProductWRTBase(tmpphys, coeffs);
1168 for (n1 = 0; n1 < nz_loc; ++n1)
1170 for (j = 0; j < nbmap; ++j)
1172 B_data[j + (k * nz_loc + n1) * nbmap +
1173 (i + (nv * nz_loc + n1) * nimap) *
1174 nbndry] += coeffs[bmap[j]];
1177 for (j = 0; j < nimap; ++j)
1179 D_data[j + (k * nz_loc + n1) * nimap +
1180 (i + (nv * nz_loc + n1) * nimap) *
1181 nint] += coeffs[imap[j]];
1195 Blas::Dgemm(
'N',
'T', B->GetRows(), C->GetRows(), B->GetColumns(),
1196 -1.0, B->GetRawPtr(), B->GetRows(), C->GetRawPtr(),
1197 C->GetRows(), 1.0, Ah->GetRawPtr(), Ah->GetRows());
1220 DNekMat DintCinvDTint, BCinvDTint_m_DTbnd, DintCinvBTtilde_m_Dbnd;
1226 DintCinvDTint = (*Dint) * (*Cinv) *
Transpose(*Dint);
1230 DintCinvBTtilde_m_Dbnd =
1231 (*Dint) * (*Cinv) *
Transpose(*Btilde) - (*Dbnd);
1235 nsize_bndry_p1[n], nsize_p_m1[n], zero);
1237 nsize_p_m1[n], nsize_bndry_p1[n], zero);
1239 nsize_p_m1[n], nsize_p_m1[n], zero);
1244 for (n1 = 0; n1 < nz_loc; ++n1)
1246 for (i = 0; i < psize - 1; ++i)
1248 for (n2 = 0; n2 < nz_loc; ++n2)
1250 for (j = 0; j < psize - 1; ++j)
1254 Dh_data[(i + n1 * (psize - 1)) +
1255 (j + n2 * (psize - 1)) * Dh->GetRows()] =
1256 -DintCinvDTint(i + 1 + n1 * psize,
1257 j + 1 + n2 * psize);
1263 for (n1 = 0; n1 < nz_loc; ++n1)
1265 for (i = 0; i < nsize_bndry_p1[n] - nz_loc; ++i)
1267 (*Ah)(i, nsize_bndry_p1[n] - nz_loc + n1) =
1268 BCinvDTint_m_DTbnd(i, n1 * psize);
1269 (*Ah)(nsize_bndry_p1[n] - nz_loc + n1, i) =
1270 DintCinvBTtilde_m_Dbnd(n1 * psize, i);
1274 for (n1 = 0; n1 < nz_loc; ++n1)
1276 for (n2 = 0; n2 < nz_loc; ++n2)
1278 (*Ah)(nsize_bndry_p1[n] - nz_loc + n1,
1279 nsize_bndry_p1[n] - nz_loc + n2) =
1280 -DintCinvDTint(n1 * psize, n2 * psize);
1284 for (n1 = 0; n1 < nz_loc; ++n1)
1286 for (j = 0; j < psize - 1; ++j)
1288 for (i = 0; i < nsize_bndry_p1[n] - nz_loc; ++i)
1290 (*Bh)(i, j + n1 * (psize - 1)) =
1291 BCinvDTint_m_DTbnd(i, j + 1 + n1 * psize);
1292 (*Ch)(j + n1 * (psize - 1), i) =
1293 DintCinvBTtilde_m_Dbnd(j + 1 + n1 * psize, i);
1298 for (n1 = 0; n1 < nz_loc; ++n1)
1300 for (n2 = 0; n2 < nz_loc; ++n2)
1302 for (j = 0; j < psize - 1; ++j)
1304 (*Bh)(nsize_bndry_p1[n] - nz_loc + n1,
1305 j + n2 * (psize - 1)) =
1306 -DintCinvDTint(n1 * psize, j + 1 + n2 * psize);
1307 (*Ch)(j + n2 * (psize - 1),
1308 nsize_bndry_p1[n] - nz_loc + n1) =
1309 -DintCinvDTint(j + 1 + n2 * psize, n1 * psize);
1316 (*Bh) = (*Bh) * (*Dh);
1318 Blas::Dgemm(
'N',
'N', Bh->GetRows(), Ch->GetColumns(), Bh->GetColumns(),
1319 -1.0, Bh->GetRawPtr(), Bh->GetRows(), Ch->GetRawPtr(),
1320 Ch->GetRows(), 1.0, Ah->GetRawPtr(), Ah->GetRows());
1338 cout <<
"Matrix Setup Costs: " << timer.
TimePerTest(1) << endl;
1411 "Advection Velocity section must be defined in "
1414 std::vector<std::string> fieldStr;
1420 GetFunction(
"AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1439 "Restart section must be defined in session file.");
1447 std::vector<std::string> fieldStr;
1453 GetFunction(
"Restart")->Evaluate(fieldStr, Restart);
1460 cout <<
"Saving the RESTART file for m_kinvis = " <<
m_kinvis
1461 <<
" (<=> Re = " << 1 /
m_kinvis <<
")" << endl;
1472 cout <<
"Saving the Stokes Flow for m_kinvis = " <<
m_kinvis
1473 <<
" (<=> Re = " << 1 /
m_kinvis <<
")" << endl;
1489 "Advection Velocity section must be defined in "
1492 std::vector<std::string> fieldStr;
1498 GetFunction(
"AdvectionVelocity")->Evaluate(fieldStr, AdvField);
1506 "Unknown or undefined equation type for CoupledLinearNS");
1519 x->Apply(
m_fields, outarray, outarray, time);
1534 if (fabs(lambda_store - lambda) > 1e-10)
1537 lambda_store = lambda;
1566 for (
int k = 0; k < nfields; ++k)
1577 for (
int k = 0; k < nfields; ++k)
1604 Generaltimer.
Start();
1612 cout <<
"We execute INITIALLY SolveSteadyNavierStokes for m_kinvis "
1621 cout <<
"We execute SolveSteadyNavierStokes for m_kinvis = "
1633 cout <<
"We execute SolveSteadyNavierStokes for m_kinvis = "
1642 Generaltimer.
Stop();
1643 cout <<
"\nThe total calculation time is : "
1644 << Generaltimer.
TimePerTest(1) / 60 <<
" minute(s). \n\n";
1651 "Unknown or undefined equation type for CoupledLinearNS");
1667 const unsigned int ncmpt =
m_velocity.size();
1671 for (
int i = 0; i < ncmpt; ++i)
1682 x->Apply(
m_fields, forcing_phys, forcing_phys, time);
1684 for (
unsigned int i = 0; i < ncmpt; ++i)
1688 m_fields[i]->IProductWRTBase(forcing_phys[i], forcing[i]);
1689 m_fields[i]->SetWaveSpace(waveSpace);
1709 if (
m_session->DefinesFunction(
"ForcingTerm"))
1711 std::vector<std::string> fieldStr;
1726 cout <<
"'ForcingTerm' section has not been defined in the input file "
1735 Newtontimer.
Start();
1754 L2Norm(delta_velocity_Phys, L2_norm);
1757 while (max(L2_norm[0], L2_norm[1]) >
m_tol)
1807 Vmath::Vadd(Velocity_Phys[i].size(), Velocity_Phys[i], 1,
1808 delta_velocity_Phys[i], 1, Velocity_Phys[i], 1);
1812 L2Norm(delta_velocity_Phys, L2_norm);
1814 if (max(Inf_norm[0], Inf_norm[1]) > 100)
1816 cout <<
"\nThe Newton method has failed at m_kinvis = " <<
m_kinvis
1817 <<
" (<=> Re = " << 1 /
m_kinvis <<
")" << endl;
1818 ASSERTL0(0,
"The Newton method has failed... \n");
1835 cout <<
"We have done " <<
m_counter - 1 <<
" iteration(s) in "
1836 << Newtontimer.
TimePerTest(1) / 60 <<
" minute(s). \n\n";
1846 cout <<
"We apply the continuation method: " << endl;
1881 Vmath::Vadd(u_star[i].size(), u_star[i], 1, u_N[i], 1, u_star[i], 1);
1896 for (
int j = 0; j < inarray[i].size(); ++j)
1898 if (inarray[i][j] > outarray[i])
1900 outarray[i] = inarray[i][j];
1903 cout <<
"InfNorm[" << i <<
"] = " << outarray[i] << endl;
1913 for (
int j = 0; j < inarray[i].size(); ++j)
1915 outarray[i] += inarray[i][j] * inarray[i][j];
1917 outarray[i] =
sqrt(outarray[i]);
1918 cout <<
"L2Norm[" << i <<
"] = " << outarray[i] << endl;
1963 i, tmp_DerVel[i], ViscTerm[i]);
1968 Vmath::Vsub(outarray[i].size(), outarray[i], 1, AdvTerm[i], 1,
1970 Vmath::Vsub(outarray[i].size(), outarray[i], 1, ViscTerm[i], 1,
1973 Vmath::Vadd(outarray[i].size(), outarray[i], 1, Forc[i], 1, outarray[i],
1988 for (
auto &expMapIter : VelExp)
1992 for (i = 0; i < expMapIter.second->m_basisKeyVector.size(); ++i)
1997 "Velocity polynomial space not sufficiently high (>= 4)");
2001 BasisVec.push_back(newB);
2007 expMapIter.second->m_geomShPtr, BasisVec);
2008 (*returnval)[expMapIter.first] = expansionElementShPtr;
2012 m_graph->SetExpansionInfo(
"p", returnval);
2105 force[i] = forcing[i] + 2 * n * ncoeffsplane;
2122 force[i] = forcing[i];
2133 int i, j, k, n, cnt, cnt1;
2134 int nbnd, nint, offset;
2136 int nel = fields[0]->GetNumElmts();
2166 for (i = 0; i < fields.size(); ++i)
2169 tmp = fields[i]->UpdateCoeffs() + nplanecoeffs,
2178 for (k = 0; k < nvel; ++k)
2181 std::dynamic_pointer_cast<MultiRegions::ContField>(fields[k]);
2184 cfield->GetLocalToGlobalMap()->GetBndCondCoeffsToLocalCoeffsSign();
2186 cfield->GetLocalToGlobalMap()->GetBndCondCoeffsToLocalCoeffsMap();
2190 fields[k]->GetBndConditions();
2196 m_fields[k]->GetPlane(2 * mode)->GetBndCondExpansions();
2200 bndCondExp =
m_fields[k]->GetBndCondExpansions();
2203 for (n = 0; n < nz_loc; ++n)
2206 for (i = 0; i < bndCondExp.size(); ++i)
2209 bndCondExp[i]->GetCoeffs();
2212 if (bndConds[i]->GetBoundaryConditionType() ==
2214 bndConds[i]->GetBoundaryConditionType() ==
2219 for (j = 0; j < (bndCondExp[i])->
GetNcoeffs(); j++)
2221 forcing[k][n * nplanecoeffs + map[bndcnt + j]] +=
2222 sign[bndcnt + j] * bndcoeffs[j];
2227 for (j = 0; j < (bndCondExp[i])->
GetNcoeffs(); j++)
2229 forcing[k][n * nplanecoeffs + map[bndcnt + j]] +=
2235 bndcnt += bndCondExp[i]->GetNcoeffs();
2246 for (i = 0; i < nel; ++i)
2248 fields[
m_velocity[0]]->GetExp(i)->GetBoundaryMap(bmap);
2249 fields[
m_velocity[0]]->GetExp(i)->GetInteriorMap(imap);
2252 offset = fields[
m_velocity[0]]->GetCoeff_Offset(i);
2254 for (j = 0; j < nvel; ++j)
2258 for (n = 0; n < nz_loc; ++n)
2260 for (k = 0; k < nbnd; ++k)
2263 forcing[j][n * nplanecoeffs + offset + bmap[k]];
2265 bnd[bndoffset + (n + j * nz_loc) * nbnd + k] =
2266 incoeffs[n * nplanecoeffs + offset + bmap[k]];
2268 for (k = 0; k < nint; ++k)
2271 forcing[j][n * nplanecoeffs + offset + imap[k]];
2279 nvel * nz_loc * nbnd + nz_loc * (
pressure->GetExp(i)->GetNcoeffs());
2287 F_bnd = F_bnd - (*
m_mat[mode].m_BCinv) * F_int;
2288 F_p_tmp = (*
m_mat[mode].m_Cinv) * F_int;
2289 F_p = (*
m_mat[mode].m_D_int) * F_p_tmp;
2296 for (i = 0; i < nel; ++i)
2298 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2300 for (j = 0; j < nvel; ++j)
2302 for (k = 0; k < nbnd; ++k)
2304 fh_bnd[offset + j * nbnd + k] = f_bnd[cnt + k];
2309 nint =
pressure->GetExp(i)->GetNcoeffs();
2310 offset += nvel * nbnd + nint * nz_loc;
2314 for (i = 0; i < nel; ++i)
2316 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2317 nint =
pressure->GetExp(i)->GetNcoeffs();
2319 for (n = 0; n < nz_loc; ++n)
2321 for (j = 0; j < nint; ++j)
2323 fh_bnd[offset + nvel * nbnd + n * nint + j] = f_p[cnt1 + j];
2327 offset += nvel * nbnd + nz_loc * nint;
2333 int totpcoeffs =
pressure->GetNcoeffs();
2335 for (i = 0; i < nel; ++i)
2337 nbnd = nz_loc * fields[0]->GetExp(i)->NumBndryCoeffs();
2338 nint =
pressure->GetExp(i)->GetNcoeffs();
2339 for (j = 0; j < nvel; ++j)
2341 for (k = 0; k < nbnd; ++k)
2343 f_bnd[cnt + k] = bnd[offset + j * nbnd + k];
2347 offset += nvel * nbnd + nint * nz_loc;
2352 offset = cnt = cnt1 = 0;
2353 for (i = 0; i < nel; ++i)
2355 nint =
pressure->GetExp(i)->GetNcoeffs();
2356 nbnd = fields[0]->GetExp(i)->NumBndryCoeffs();
2357 cnt1 =
pressure->GetCoeff_Offset(i);
2359 for (n = 0; n < nz_loc; ++n)
2361 for (j = 0; j < nint; ++j)
2363 p_coeffs[n * totpcoeffs + cnt1 + j] = f_p[cnt + j] =
2364 bnd[offset + nvel * nz_loc * nbnd + n * nint + j];
2368 offset += (nvel * nbnd + nint) * nz_loc;
2375 F_int = (*
m_mat[mode].m_Cinv) * F_int;
2379 for (i = 0; i < nel; ++i)
2381 fields[0]->GetExp(i)->GetBoundaryMap(bmap);
2382 fields[0]->GetExp(i)->GetInteriorMap(imap);
2385 offset = fields[0]->GetCoeff_Offset(i);
2387 for (j = 0; j < nvel; ++j)
2389 for (n = 0; n < nz_loc; ++n)
2391 for (k = 0; k < nbnd; ++k)
2393 fields[j]->SetCoeff(n * nplanecoeffs + offset + bmap[k],
2397 for (k = 0; k < nint; ++k)
2399 fields[j]->SetCoeff(n * nplanecoeffs + offset + imap[k],
2408 for (j = 0; j < nvel; ++j)
2410 fields[j]->SetPhysState(
false);
2416 std::vector<Array<OneD, NekDouble>> fieldcoeffs(
m_fields.size() + 1);
2417 std::vector<std::string> variables(
m_fields.size() + 1);
2420 for (i = 0; i <
m_fields.size(); ++i)
2422 fieldcoeffs[i] =
m_fields[i]->UpdateCoeffs();
2437 m_fields[0]->GetPlane(0)->FwdTransLocalElmt(
2438 m_pressure->GetPlane(0)->GetPhys(), fieldcoeffs[i]);
2439 m_fields[0]->GetPlane(1)->FwdTransLocalElmt(
2440 m_pressure->GetPlane(1)->GetPhys(), tmpfieldcoeffs);
2441 for (
int e = 0; e <
m_fields[0]->GetNcoeffs() / 2; e++)
2443 fieldcoeffs[i][e +
m_fields[0]->GetNcoeffs() / 2] =
2461 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
virtual void v_TransPhysToCoeff(void)
Virtual function for transformation to coefficient space.
Array< OneD, CoupledLocalToGlobalC0ContMapSharedPtr > m_locToGloMap
virtual void v_Output(void)
void DefineForcingTerm(void)
void EvaluateNewtonRHS(Array< OneD, Array< OneD, NekDouble >> &Velocity, Array< OneD, Array< OneD, NekDouble >> &outarray)
void SolveSteadyNavierStokes(void)
void SetUpCoupledMatrix(const NekDouble lambda=0.0, const Array< OneD, Array< OneD, NekDouble >> &Advfield=NullNekDoubleArrayOfArray, bool IsLinearNSEquation=true)
const SpatialDomains::ExpansionInfoMap & GenPressureExp(const SpatialDomains::ExpansionInfoMap &VelExp)
virtual void v_DoSolve(void)
Solves an unsteady problem.
void EvaluateAdvection(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
void InfNorm(Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, NekDouble > &outarray)
virtual bool v_NegatedOp(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)
virtual int v_GetForceDimension(void)
bool m_zeroMode
Id to identify when single mode is mean mode (i.e. beta=0);.
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm_Coeffs
void L2Norm(Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, NekDouble > &outarray)
Array< OneD, Array< OneD, NekDouble > > m_ForcingTerm
virtual void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
virtual void v_DoInitialise(void)
Sets up initial conditions.
NekDouble m_KinvisPercentage
Array< OneD, CoupledSolverMatrices > m_mat
virtual void v_TransCoeffToPhys(void)
Virtual function for transformation to physical space.
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
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.
virtual void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
ExtrapolateSharedPtr m_extrapolation
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
EquationType m_equationType
equation type;
int m_nConvectiveFields
Number of fields to be convected;.
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
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.
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.
NekDouble m_time
Current time of simulation.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
NekDouble m_lambda
Lambda constant in real system if one required.
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.
enum HomogeneousType m_HomogeneousType
SOLVER_UTILS_EXPORT int GetNpoints()
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
SOLVER_UTILS_EXPORT int GetNcoeffs()
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
SOLVER_UTILS_EXPORT int GetTotPoints()
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
virtual SOLVER_UTILS_EXPORT void v_DoSolve()
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
The above copyright notice and this permission notice shall be included.
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
ExtrapolateFactory & GetExtrapolateFactory()
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