38#include <boost/algorithm/string/predicate.hpp>
68 int shapedim =
m_fields[0]->GetShapeDimension();
70 for (
int j = 0; j < shapedim; ++j)
91 "No TESTTYPE defined in session.");
92 std::string TestTypeStr =
m_session->GetSolverInfo(
"TESTTYPE");
109 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
114 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
118 m_session->LoadParameter(
"Omega", Omegams, 7.292 * 0.00001);
119 m_Omega = Omegams * SecondToDay;
130 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
135 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
139 m_session->LoadParameter(
"Omega", Omegams, 7.292 * 0.00001);
140 m_Omega = Omegams * SecondToDay;
142 m_H0 = 133681.0 / (rad_earth * gms);
143 m_k2 = 10.0 / (rad_earth * gms);
149 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
154 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
159 m_u0 =
m_u0 * SecondToDay / rad_earth;
161 m_session->LoadParameter(
"Omega", Omegams, 7.292 * 0.00001);
162 m_Omega = Omegams * SecondToDay;
164 m_H0 = 5960.0 / rad_earth;
165 m_hs0 = 2000.0 / rad_earth;
171 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
176 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
179 m_u0 =
m_u0 * SecondToDay / rad_earth;
181 m_session->LoadParameter(
"Omega", Omegams, 7.292 * 0.00001);
182 m_Omega = Omegams * SecondToDay;
191 m_hbar = 120.0 / rad_earth;
193 std::cout <<
"m_theta0 = " <<
m_theta0
195 <<
", m_hbar = " <<
m_hbar << std::endl;
201 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
206 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
208 m_session->LoadParameter(
"Omega", Omegams, 7.292 * 0.00001);
209 m_Omega = Omegams * SecondToDay;
214 m_angfreq = 7.848 * 0.000001 * SecondToDay;
215 m_K = 7.848 * 0.000001 * SecondToDay;
238 ASSERTL0(
false,
"Implicit unsteady Advection not set up.");
253 for (i = 0; i < nfields; ++i)
257 nvariables = nfields;
269 for (i = 0; i < nvariables; ++i)
282 "Only one of IO_CheckTime and IO_CheckSteps "
286 bool doCheckTime =
false;
292 NekDouble Mass = 0.0, Energy = 0.0, Enstrophy = 0.0, Vorticity = 0.0;
295 for (
int i = 0; i < nvariables; ++i)
315 std::cout <<
"Steps: " << std::setw(8) << std::left << step + 1
317 <<
"Time: " << std::setw(12) << std::left <<
m_time;
319 std::stringstream ss;
320 ss << cpuTime <<
"s";
321 std::cout <<
" CPU Time: " << std::setw(8) << std::left << ss.str()
335 Energy = (
ComputeEnergy(fieldsprimitive[0], fieldsprimitive[1],
336 fieldsprimitive[2]) -
344 fieldsprimitive[2]) -
348 std::cout <<
"dMass = " << std::setw(8) << std::left << Mass <<
" "
349 <<
", dEnergy = " << std::setw(8) << std::left << Energy
351 <<
", dEnstrophy = " << std::setw(8) << std::left
353 <<
", dVorticity = " << std::setw(8) << std::left
354 << Vorticity << std::endl
364 for (i = 0; i < nvariables; ++i)
388 if (
m_session->GetComm()->GetRank() == 0)
394 <<
"CFL time-step : " <<
m_timestep << std::endl;
397 if (
m_session->GetSolverInfo(
"Driver") !=
"SteadyState")
399 std::cout <<
"Time-integration : " << intTime <<
"s" << std::endl;
403 for (i = 0; i < nvariables; ++i)
412 for (i = 0; i < nvariables; ++i)
424 int nvariables = inarray.size();
432 for (i = 0; i < nvariables; ++i)
445 for (i = 0; i < nvariables; ++i)
466 for (i = 0; i < nvariables; ++i)
468 m_fields[i]->MultiplyByElmtInvMass(modarray[i], modarray[i]);
469 m_fields[i]->BwdTrans(modarray[i], outarray[i]);
492 for (i = 0; i < nvariables; ++i)
494 physfield[i] = InField[i];
504 for (i = 0; i < nvariables; ++i)
515 Vmath::Vadd(ncoeffs, &tmp[0], 1, &OutField[i][0], 1,
524 for (i = 0; i < nvariables; ++i)
533 for (i = 0; i < nvariables; ++i)
536 m_fields[i]->AddFwdBwdTraceIntegral(numfluxFwd[i], numfluxBwd[i],
547 int ncoeffs = outarray[0].size();
548 int nq = physarray[0].size();
567 m_fields[0]->IProductWRTBase(tmp, tmpc);
569 Vmath::Vadd(ncoeffs, outarray[j + 1], 1, tmpc, 1, outarray[j + 1], 1);
577 int nq =
m_fields[0]->GetTotPoints();
588 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
591 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
604 Vmath::Vmul(nq, tmp, 1, physfield[1], 1, flux[1], 1);
607 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
617 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
630 Vmath::Vmul(nq, tmp, 1, physfield[2], 1, flux[0], 1);
633 Vmath::Vmul(nq, flux[0], 1, physfield[2], 1, flux[1], 1);
636 Vmath::Vmul(nq, flux[0], 1, physfield[1], 1, flux[0], 1);
683 for (i = 0; i < nvariables; ++i)
690 for (i = 0; i < nvariables; ++i)
692 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
716 for (k = 0; k < nTraceNumPoints; ++k)
740 for (k = 0; k < nTraceNumPoints; ++k)
743 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
744 Bwd[2][k], numfluxF, numfluxB);
762 denomFwd = eF1n * eF2t - eF2n * eF1t;
763 denomBwd = eB1n * eB2t - eB2n * eB1t;
765 numfluxFwd[0][k] = numfluxF[0];
766 numfluxFwd[1][k] = (1.0 / denomFwd) *
767 (eF2t * numfluxF[1] - eF2n * numfluxF[2]);
770 (-1.0 * eF1t * numfluxF[1] + eF1n * numfluxF[2]);
772 numfluxBwd[0][k] = 1.0 * numfluxB[0];
773 numfluxBwd[1][k] = (1.0 / denomBwd) *
774 (eB2t * numfluxB[1] - eB2n * numfluxB[2]);
777 (-1.0 * eB1t * numfluxB[1] + eB1n * numfluxB[2]);
789 for (k = 0; k < nTraceNumPoints; ++k)
794 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
795 Bwd[2][k], numfluxF, numfluxB);
801 Fwd[2][k], Bwd[0][k] + DepthFwd[k],
802 Bwd[1][k], Bwd[2][k], numfluxF, numfluxB);
808 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
809 Bwd[2][k], numfluxF, numfluxB);
814 for (i = 0; i < nvariables; ++i)
824 numfluxFwd[i][k] = tmpF0 + tmpF1 + numfluxF[indx + 2];
825 numfluxBwd[i][k] = tmpB0 + tmpB1 + numfluxB[indx + 2];
833 ASSERTL0(
false,
"populate switch statement for upwind flux");
860 hstarF = 0.5 * (cL + cR) + 0.25 * (uL - uRF);
864 hstarB = 0.5 * (cL + cR) + 0.25 * (uLB - uR);
875 numfluxF[0] = hfluxF;
876 numfluxF[1] = hufluxF;
877 numfluxF[2] = hvfluxF;
879 numfluxB[0] = hfluxB;
880 numfluxB[1] = hufluxB;
881 numfluxB[2] = hvfluxB;
898 SL = uL - cL *
sqrt(0.5 * ((hstar * hstar + hstar * hL) / (hL * hL)));
908 SR = uR + cR *
sqrt(0.5 * ((hstar * hstar + hstar * hR) / (hR * hR)));
915 if (fabs(hR * (uR - SR) - hL * (uL - SL)) <= 1.0e-15)
921 Sstar = (SL * hR * (uR - SR) - SR * hL * (uL - SL)) /
922 (hR * (uR - SR) - hL * (uL - SL));
928 huflux = uL * uL * hL + 0.5 * g * hL * hL;
929 hvflux = hL * uL * vL;
934 huflux = uR * uR * hR + 0.5 * g * hR * hR;
935 hvflux = hR * uR * vR;
939 if ((SL < 0) && (Sstar >= 0))
941 hC = hL * ((SL - uL) / (SL - Sstar));
945 hflux = hL * uL + SL * (hC - hL);
946 huflux = (uL * uL * hL + 0.5 * g * hL * hL) + SL * (huC - hL * uL);
947 hvflux = (uL * vL * hL) + SL * (hvC - hL * vL);
951 hC = hR * ((SR - uR) / (SR - Sstar));
955 hflux = hR * uR + SR * (hC - hR);
956 huflux = (uR * uR * hR + 0.5 * g * hR * hR) + SR * (huC - hR * uR);
957 hvflux = (uR * vR * hR) + SR * (hvC - hR * vR);
967 NekDouble MageF1, MageF2, MageB1, MageB2;
975 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
979 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
980 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
982 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
983 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
987 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
988 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
991 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
993 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
998 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
999 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1001 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1002 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1006 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1007 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1010 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1012 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1022 NekDouble MageF1, MageF2, MageB1, MageB2;
1035 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1041 EigF[0] = velL -
sqrt(g * hL);
1043 EigF[2] = velL +
sqrt(g * hL);
1045 EigB[0] = velR -
sqrt(g * hR);
1047 EigB[2] = velR +
sqrt(g * hR);
1054 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1055 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1057 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1058 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1059 numfluxF[2] = 0.5 * lambdaF * (hL - hR);
1061 numfluxF[3] = 0.5 * (hL * uL * uL * MageF1 + hR * uRF * uRF * MageB1 +
1062 0.5 * g * (hL * hL + hR * hR));
1063 numfluxF[4] = 0.5 * (hL * uL * vL * MageF1 + hR * uRF * vRF * MageB1);
1064 numfluxF[5] = 0.5 * lambdaF * (uL * hL - uRF * hR);
1066 numfluxF[6] = 0.5 * (hL * uL * vL * MageF2 + hR * uRF * vRF * MageB2);
1067 numfluxF[7] = 0.5 * (hL * vL * vL * MageF2 + hR * vRF * vRF * MageB2 +
1068 0.5 * g * (hL * hL + hR * hR));
1069 numfluxF[8] = 0.5 * lambdaF * (vL * hL - vRF * hR);
1073 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1074 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1076 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1077 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1078 numfluxB[2] = 0.5 * lambdaB * (hL - hR);
1080 numfluxB[3] = 0.5 * (hR * uR * uR * MageB1 + hR * uLB * uLB * MageF1 +
1081 0.5 * g * (hR * hR + hL * hL));
1082 numfluxB[4] = 0.5 * (hR * uR * vR * MageB1 + hR * uLB * vLB * MageF1);
1083 numfluxB[5] = 0.5 * lambdaB * (uLB * hL - uR * hR);
1085 numfluxB[6] = 0.5 * (hR * uR * vR * MageB2 + hR * uLB * vLB * MageF2);
1086 numfluxB[7] = 0.5 * (hR * vR * vR * MageB2 + hR * vLB * vLB * MageF2 +
1087 0.5 * g * (hR * hR + hL * hL));
1088 numfluxB[8] = 0.5 * lambdaB * (vLB * hL - vR * hR);
1097 NekDouble MageF1, MageF2, MageB1, MageB2;
1110 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1117 SL = fabs(velL) +
sqrt(g * hL);
1118 SR = fabs(velR) +
sqrt(g * hR);
1132 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1133 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1135 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1136 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1137 numfluxF[2] = 0.5 * S * (hL - hR);
1140 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
1141 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1142 numfluxF[5] = 0.5 * S * (uL * hL - uRF * hR);
1144 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1146 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
1147 numfluxF[8] = 0.5 * S * (vL * hL - vRF * hR);
1151 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1152 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1154 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1155 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1156 numfluxB[2] = 0.5 * S * (hL - hR);
1159 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1160 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1161 numfluxB[5] = 0.5 * S * (uLB * hL - uR * hR);
1163 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1165 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1166 numfluxB[8] = 0.5 * S * (vLB * hL - vR * hR);
1175 NekDouble MF1x, MF1y, MF1z, MF2x, MF2y, MF2z;
1176 NekDouble MB1x, MB1y, MB1z, MB2x, MB2y, MB2z;
1195 MageF1 = MF1x * MF1x + MF1y * MF1y + MF1z * MF1z;
1196 MageF2 = MF2x * MF2x + MF2y * MF2y + MF2z * MF2z;
1197 MageB1 = MB1x * MB1x + MB1y * MB1y + MB1z * MB1z;
1198 MageB2 = MB2x * MB2x + MB2y * MB2y + MB2z * MB2z;
1200 eF1_cdot_eB1 = MF1x * MB1x + MF1y * MB1y + MF1z * MB1z;
1201 eF1_cdot_eB2 = MF1x * MB2x + MF1y * MB2y + MF1z * MB2z;
1202 eF2_cdot_eB1 = MF2x * MB1x + MF2y * MB1y + MF2z * MB1z;
1203 eF2_cdot_eB2 = MF2x * MB2x + MF2y * MB2y + MF2z * MB2z;
1210 int ncoeffs = outarray[0].size();
1211 int nq = physarray[0].size();
1246 m_fields[0]->IProductWRTBase(tmp, tmpc);
1247 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1255 int ncoeffs = outarray[0].size();
1256 int nq = physarray[0].size();
1272 m_fields[0]->IProductWRTBase(tmp, tmpc);
1274 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1285 int ncoeffs = outarray[0].size();
1286 int nq = physarray[0].size();
1303 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e0, 1, Rott1, 1);
1304 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e1, 1, Rott2, 1);
1305 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e0, 1, Rott1, 1, Rott1, 1);
1306 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e1, 1, Rott2, 1, Rott2, 1);
1309 Vmath::Vmul(nq, &h[0], 1, &Rott1[0], 1, &Rott1[0], 1);
1310 Vmath::Vmul(nq, &h[0], 1, &Rott2[0], 1, &Rott2[0], 1);
1317 m_fields[0]->IProductWRTBase(Rott1, tmpc1);
1318 m_fields[0]->IProductWRTBase(Rott2, tmpc2);
1320 Vmath::Vadd(ncoeffs, tmpc1, 1, outarray[1], 1, outarray[1], 1);
1321 Vmath::Vadd(ncoeffs, tmpc2, 1, outarray[2], 1, outarray[2], 1);
1329 const int indm,
const int indk,
1334 int nq =
m_fields[0]->GetNpoints();
1348 Vmath::Vmul(nq, &physarray[j + 1][0], 1, &tmpderiv[0], 1,
1352 &outarray[0], 1, &outarray[0], 1);
1378 ASSERTL0(
false,
"Unknown projection scheme");
1391 for (
int n = 0; n <
m_fields[0]->GetBndConditions().size(); ++n)
1395 if (
m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
"eMG")
1399 ASSERTL0(
false,
"Illegal dimension");
1408 if (
m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
"eWall")
1412 ASSERTL0(
false,
"Illegal dimension");
1421 if (
m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
1424 for (
int i = 0; i < nvariables; ++i)
1426 m_fields[i]->EvaluateBoundaryConditions(time);
1429 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
1439 int nvariables = physarray.size();
1443 for (i = 0; i < nvariables; ++i)
1446 m_fields[i]->ExtractTracePhys(physarray[i], Fwd0[i]);
1450 for (i = 0; i < nvariables; ++i)
1453 m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
1458 int e, id1, id2, npts;
1464 ->GetBndCondExpansions()[bcRegion]
1467 id1 =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
1468 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
1469 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
1487 &tmp_n[0], 1, &tmp_n[0], 1);
1492 1, &tmp_t[0], 1, &tmp_t[0], 1);
1512 &tmp_u[0], 1, &tmp_u[0], 1);
1513 Vmath::Vdiv(npts, &tmp_u[0], 1, &denom[0], 1, &tmp_u[0], 1);
1520 &tmp_v[0], 1, &tmp_v[0], 1);
1521 Vmath::Vdiv(npts, &tmp_v[0], 1, &denom[0], 1, &tmp_v[0], 1);
1528 ASSERTL0(
false,
"Illegal expansion dimension");
1532 for (i = 0; i < nvariables; ++i)
1536 ->GetBndCondExpansions()[bcRegion]
1537 ->UpdatePhys())[id1],
1555 for (
int i = 0; i <
m_fields.size(); ++i)
1587 NekDouble sin_varphi, cos_varphi, sin_theta, cos_theta;
1589 for (
int j = 0; j < nq; ++j)
1596 sin_theta, cos_theta);
1615 NekDouble phi, theta, sin_varphi, cos_varphi, sin_theta, cos_theta;
1621 thetac =
m_pi / 6.0;
1623 for (
int j = 0; j < nq; ++j)
1630 sin_theta, cos_theta);
1632 if ((
std::abs(sin(phic) - sin_varphi) +
1633 std::abs(sin(thetac) - sin_theta)) < Tol)
1635 std::cout <<
"A point " << j
1636 <<
" is coincient with the singularity "
1641 phi = atan2(sin_varphi, cos_varphi);
1642 theta = atan2(sin_theta, cos_theta);
1645 dist2 = (phi - phic) * (phi - phic) +
1646 (theta - thetac) * (theta - thetac);
1648 if (dist2 > hRad * hRad)
1650 dist2 = hRad * hRad;
1658 std::cout <<
"No point is coincident with the singularity point"
1666 for (
int j = 0; j < nq; ++j)
1697 std::cout <<
"Water Depth (m_depth) was generated with mag = "
1736 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1748 for (
int j = 0; j < nq; ++j)
1759 tmp = -1.0 * cos_varphi * cos_theta * sin(
m_alpha) +
1761 outarray[j] = 2.0 *
m_Omega * tmp;
1770 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1779 for (
int j = 0; j < nq; ++j)
1788 outarray[j] = 2.0 *
m_Omega * sin_theta;
1793 bool dumpInitialConditions,
1794 [[maybe_unused]]
const int domain)
1816 for (
int i = 0; i <
m_fields.size(); ++i)
1849 for (
int i = 0; i <
m_fields.size(); ++i)
1878 for (
int i = 0; i <
m_fields.size(); ++i)
1907 for (
int i = 0; i <
m_fields.size(); ++i)
1936 for (
int i = 0; i <
m_fields.size(); ++i)
1965 for (
int i = 0; i <
m_fields.size(); ++i)
1978 if (dumpInitialConditions)
1993 int nq =
m_fields[0]->GetNpoints();
1999 m_fields[0]->GetCoords(x0, x1, x2);
2012 for (
int i = 0; i < nq; ++i)
2014 eta0[i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2015 (1.0 / cosh(0.395 * x0[i]))) *
2016 (3.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2017 exp(-0.5 * x1[i] * x1[i]);
2018 uvec[0][i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2019 (1.0 / cosh(0.395 * x0[i]))) *
2020 (-9.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2021 exp(-0.5 * x1[i] * x1[i]);
2022 uvec[1][i] = (-2.0 * 0.395 * tanh(0.395 * x0[i])) *
2023 (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2024 (1.0 / cosh(0.395 * x0[i]))) *
2025 (2.0 * x1[i]) * exp(-0.5 * x1[i] * x1[i]);
2058 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2077 for (
int j = 0; j < nq; ++j)
2088 tmp = -1.0 * cos_varphi * cos_theta * sin(
m_alpha) +
2095 sin_theta * cos_varphi * sin(
m_alpha));
2098 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2099 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2100 uvec[2][j] = vhat * cos_theta;
2146 int nq =
m_fields[0]->GetTotPoints();
2158 int nq =
m_fields[0]->GetTotPoints();
2187 int nq =
m_fields[0]->GetTotPoints();
2218 int nq =
m_fields[0]->GetTotPoints();
2232 Vmath::Vadd(nq, tmp, 1, Vorticity, 1, Vorticity, 1);
2237 int nq =
m_fields[0]->GetNpoints();
2259 1, &velcoeff[0], 1, &velcoeff[0], 1);
2263 m_fields[0]->PhysDeriv(velcoeff, Dtmp0, Dtmp1, Dtmp2);
2273 1, &vellc[0], 1, &vellc[0], 1);
2280 1, &vellc[0], 1, &vellc[0], 1);
2287 1, &vellc[0], 1, &vellc[0], 1);
2300 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2321 for (
int j = 0; j < nq; ++j)
2334 cos_varphi * cos(
m_Omega * time) - sin_varphi * sin(
m_Omega * time);
2336 sin_varphi * cos(
m_Omega * time) + cos_varphi * sin(
m_Omega * time);
2337 tmp = -1.0 * TR * sin(
m_alpha) * cos_theta + cos(
m_alpha) * sin_theta;
2339 eta[j] = -1.0 * (
m_u0 * tmp +
m_Omega * sin_theta) *
2342 eta[j] = 0.5 * eta[j] /
m_g;
2350 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2351 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2352 uvec[2][j] = vhat * cos_theta;
2403 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2422 for (
int j = 0; j < nq; ++j)
2433 sin_theta * sin_theta;
2435 uhat =
m_u0 * cos_theta;
2438 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2439 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2440 uvec[2][j] = vhat * cos_theta;
2491 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2513 for (
int j = 0; j < nq; ++j)
2522 Ttheta = atan2(sin_theta, cos_theta);
2523 Tphi = atan2(sin_varphi, cos_varphi);
2529 dth = Ttheta / Nint;
2530 eta[j] = dth * 0.5 *
2532 for (
int i = 1; i < Nint - 1; i++)
2537 eta[j] = (-1.0 /
m_g) * eta[j];
2542 eta[j] = eta[j] +
m_hbar * cos_theta * exp(-9.0 * Tphi * Tphi) *
2543 exp(-225.0 * (
m_pi / 4.0 - Ttheta) *
2544 (
m_pi / 4.0 - Ttheta));
2547 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2548 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2549 uvec[2][j] = vhat * cos_theta;
2597 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2618 NekDouble cos2theta, cosRtheta, cos6theta, cos2Rtheta, cosRm1theta;
2619 NekDouble cos2phi, cos4phi, sin4phi, cos8phi;
2625 phi0 = 40.0 *
m_pi / 180.0;
2626 theta0 = 50.0 *
m_pi / 180.0;
2628 x0d = cos(phi0) * cos(theta0);
2629 y0d = sin(phi0) * cos(theta0);
2632 for (
int j = 0; j < nq; ++j)
2645 cos2theta = cos_theta * cos_theta;
2646 cosRm1theta = cos_theta * cos2theta;
2647 cosRtheta = cos2theta * cos2theta;
2648 cos6theta = cos2theta * cosRtheta;
2649 cos2Rtheta = cosRtheta * cosRtheta;
2652 Ath = tmp * cos2theta +
2653 0.25 *
m_K *
m_K * cos6theta *
2654 ((R + 1.0) * cosRtheta + (2 * R * R - R - 2.0) * cos2theta -
2658 Bth = tmp * cosRtheta *
2659 ((R * R + 2 * R + 2) - (R + 1.0) * (R + 1.0) * cos2theta);
2662 0.25 *
m_K *
m_K * cos2Rtheta * ((R + 1.0) * cos2theta - (R + 2.0));
2665 cos2phi = 2.0 * cos_varphi * cos_varphi - 1.0;
2666 cos4phi = 2.0 * cos2phi * cos2phi - 1.0;
2667 cos8phi = 2.0 * cos4phi * cos4phi - 1.0;
2670 sin4phi = 4.0 * sin_varphi * cos_varphi * cos2phi;
2672 eta[j] =
m_H0 + (1.0 /
m_g) * (Ath + Bth * cos4phi + Cth * cos8phi);
2678 (1.0 + (1.0 / 40.0) * (x0j * x0d + x1j * y0d + x2j * z0d));
2684 m_K * cosRm1theta * (R * sin_theta * sin_theta - cos2theta) *
2686 vhat = -1.0 *
m_K * R * cosRm1theta * sin_theta * sin4phi;
2688 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2689 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2690 uvec[2][j] = vhat * cos_theta;
2743 f = 2.0 *
m_Omega * sin(theta);
2745 dh = f * uphi + tan(theta) * uphi * uphi;
2770 int nq =
m_fields[0]->GetTotPoints();
2771 int ncoeffs =
m_fields[0]->GetNcoeffs();
2773 NekDouble rad_earth = 6.37122 * 1000000;
2778 std::vector<std::string> variables(nvariables);
2780 variables[0] =
"eta";
2781 variables[1] =
"hstar";
2782 variables[2] =
"vorticity";
2783 variables[3] =
"ux";
2784 variables[4] =
"uy";
2785 variables[5] =
"uz";
2786 variables[6] =
"null";
2790 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvariables);
2791 for (
int i = 0; i < nvariables; ++i)
2798 &fieldphys[0][0], 1);
2805 Vmath::Smul(nq, rad_earth, &fieldphys[1][0], 1, &fieldphys[1][0], 1);
2819 &fieldphys[k + indx][0], 1);
2821 &fieldphys[k + indx][0], 1, &fieldphys[k + indx][0], 1);
2824 for (
int i = 0; i < nvariables; ++i)
2826 m_fields[0]->FwdTrans(fieldphys[i], fieldcoeffs[i]);
2839 if (physin[0].data() == physout[0].data())
2843 for (
int i = 0; i < 3; ++i)
2854 Vmath::Vdiv(nq, tmp[1], 1, tmp[0], 1, physout[1], 1);
2857 Vmath::Vdiv(nq, tmp[2], 1, tmp[0], 1, physout[2], 1);
2865 Vmath::Vdiv(nq, physin[1], 1, physin[0], 1, physout[1], 1);
2868 Vmath::Vdiv(nq, physin[2], 1, physin[0], 1, physout[2], 1);
2879 if (physin[0].data() == physout[0].data())
2883 for (
int i = 0; i < 3; ++i)
2894 Vmath::Vmul(nq, physout[0], 1, tmp[1], 1, physout[1], 1);
2897 Vmath::Vmul(nq, physout[0], 1, tmp[2], 1, physout[2], 1);
2905 Vmath::Vmul(nq, physout[0], 1, physin[1], 1, physout[1], 1);
2908 Vmath::Vmul(nq, physout[0], 1, physin[2], 1, physout[2], 1);
2952 int nq =
m_fields[0]->GetTotPoints();
2957 NekDouble theta, phi, sin_theta, cos_theta, sin_varphi, cos_varphi;
2982 for (k = 0; k < nq; ++k)
2988 Re =
sqrt(xp * xp + yp * yp + zp * zp);
2995 theta = atan2(sin_theta, cos_theta);
2996 phi = atan2(sin_varphi, cos_varphi);
2998 cosntheta3 = cos(n * theta) * cos(n * theta) * cos(n * theta);
3000 beta_theta = -4.0 * n * cosntheta3 * cos(m * phi) * sin(n * theta) / Re;
3001 beta_phi = -m * cosntheta3 * sin(m * phi) / Re;
3003 thetax = -1.0 * cos_varphi * sin_theta;
3004 thetay = -1.0 * sin_varphi * sin_theta;
3007 phix = -1.0 * sin_varphi;
3011 uvec[0][k] = alpha * (beta_theta * thetax + beta_phi * phix);
3012 uvec[1][k] = alpha * (beta_theta * thetay + beta_phi * phiy);
3013 uvec[2][k] = alpha * (beta_theta * thetaz + beta_phi * phiz);
3015 vorticityexact[k] = -4.0 * n / Re / Re * cos_theta * cos_theta *
3016 cos_varphi * cos(m * phi) * sin(n * theta);
3022 std::cout <<
"chi migi1" << std::endl;
3026 Vmath::Vsub(nq, vorticityexact, 1, vorticitycompt, 1, vorticitycompt, 1);
3028 std::cout <<
"Vorticity: L2 error = " <<
AvgAbsInt(vorticitycompt)
3029 <<
", Linf error = " <<
Vmath::Vamax(nq, vorticitycompt, 1)
3060 &exactsolution[0], 1, &exactsolution[0], 1);
3061 Vmath::Vabs(nq, exactsolution, 1, exactsolution, 1);
3063 L2error = (
m_fields[0]->Integral(exactsolution)) / L2exact;
3080 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3083 L2exact =
m_fields[1]->Integral(tmp);
3092 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3095 L2error = (
m_fields[1]->Integral(tmp)) / L2exact;
3109 if (Normalised ==
true)
3116 L2error =
sqrt(L2error * L2error / Vol);
3165 1, &exactsolution[0], 1);
3187 Vmath::Vsub(nq, &exactu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3188 Vmath::Vsub(nq, &exactv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3190 Vmath::Vmul(nq, &tmpu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3191 Vmath::Vmul(nq, &tmpv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3193 Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &Lerr[0], 1);
3197 Vmath::Vmul(nq, &exactu[0], 1, &exactu[0], 1, &tmpu[0], 1);
3198 Vmath::Vmul(nq, &exactv[0], 1, &exactv[0], 1, &tmpv[0], 1);
3199 Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &uT[0], 1);
#define ASSERTL0(condition, msg)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
void ComputeVorticity(const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v, Array< OneD, NekDouble > &Vorticity)
void AddDivForGradient(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void AddCoriolis(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void UnstableJetFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
void v_InitObject(bool DeclareFields=true) override
Initialise the object.
NekDouble ComputeUnstableJetuphi(const NekDouble theta)
MMFSWE(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void v_SetInitialConditions(const NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0) override
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
void WeakDGSWEDirDeriv(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField)
void RossbyWave(unsigned int field, Array< OneD, NekDouble > &outfield)
void ComputeNablaCdotVelocity(Array< OneD, NekDouble > &vellc)
NekDouble v_L2Error(unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised) override
Virtual function for the L_2 error computation between fields and a given exact solution.
void EvaluateStandardCoriolis(Array< OneD, NekDouble > &outarray)
static std::string className
Name of class.
void GetSWEFluxVector(const int i, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
void PrimitiveToConservative()
void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time) override
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &inarray, NekDouble time)
void IsolatedMountainFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Array< OneD, NekDouble > m_coriolis
Coriolis force.
void ComputeMagAndDot(const int index, NekDouble &MageF1, NekDouble &MageF2, NekDouble &MageB1, NekDouble &MageB2, NekDouble &eF1_cdot_eB1, NekDouble &eF1_cdot_eB2, NekDouble &eF2_cdot_eB1, NekDouble &eF2_cdot_eB2)
void AddRotation(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void AverageFlux(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
void AddElevationEffect(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void Checkpoint_Output_Cartesian(std::string outname)
void Computehhuhvflux(NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, NekDouble hstar, NekDouble &hflux, NekDouble &huflux, NekDouble &hvflux)
NekDouble v_LinfError(unsigned int field, const Array< OneD, NekDouble > &exactsoln) override
Virtual function for the L_inf error computation between fields and a given exact solution.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
void TestVorticityComputation()
void NumericalSWEFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxFwd, Array< OneD, Array< OneD, NekDouble > > &numfluxBwd)
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection.
void RiemannSolverHLLC(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
NekDouble ComputeEnstrophy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print Summary.
NekDouble ComputeMass(const Array< OneD, const NekDouble > &eta)
void EvaluateWaterDepth(void)
void UnsteadyZonalFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
NekDouble ComputeEnergy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
void LaxFriedrichFlux(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
NekDouble ComputeUnstableJetEta(const NekDouble theta)
void EvaluateCoriolis(void)
void EvaluateCoriolisForZonalFlow(Array< OneD, NekDouble > &outarray)
void v_DoInitialise(bool dumpInitialConditions=false) override
Sets up initial conditions.
Array< OneD, Array< OneD, NekDouble > > m_velocity
Advection velocity.
void v_DoSolve() override
Solves an unsteady problem.
void Compute_demdt_cdot_ek(const int indm, const int indk, const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, NekDouble > &outarray)
void WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
void TestSWE2Dproblem(const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield)
Array< OneD, Array< OneD, NekDouble > > m_Derivdepth
void ConservativeToPrimitive()
void SteadyZonalFlow(unsigned int field, Array< OneD, NekDouble > &outfield)
void RusanovFlux(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
Array< OneD, NekDouble > m_depth
Still water depth.
int m_spacedim
Spatial dimension (>= expansion dim).
int m_expdim
Expansion dimension.
LibUtilities::CommSharedPtr m_comm
Communicator.
NekDouble m_timestep
Time step size.
int m_infosteps
Number of time steps between outputting status information.
NekDouble m_time
Current time of simulation.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNpoints()
NekDouble m_fintime
Finish time of the simulation.
SOLVER_UTILS_EXPORT int GetNcoeffs()
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
NekDouble m_checktime
Time between checkpoints.
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
SOLVER_UTILS_EXPORT int GetExpSize()
std::string m_sessionName
Name of the session.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
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 Array< OneD, NekDouble > ErrorExtraPoints(unsigned int field)
Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf].
SOLVER_UTILS_EXPORT int GetTotPoints()
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
int m_steps
Number of steps to take.
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
SOLVER_UTILS_EXPORT NekDouble LinfError(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
Linf error computation.
int m_checksteps
Number of steps between checkpoints.
SOLVER_UTILS_EXPORT void EvaluateExactSolution(int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
Evaluates an exact solution.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
A base class for PDEs which include an advection component.
Array< OneD, Array< OneD, NekDouble > > m_DivMF
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFFwd
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Virtual function for generating summary information.
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFBwd
SOLVER_UTILS_EXPORT NekDouble AvgAbsInt(const Array< OneD, const NekDouble > &inarray)
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceFwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceBwd
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble > > &uvec, unsigned int field)
SOLVER_UTILS_EXPORT void CopyBoundaryTrace(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const BoundaryCopyType BDCopyType, const int var=0, const std::string btype="NoUserDefined")
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
SurfaceType m_surfaceType
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble > > &Anisotropy, const int TangentXelem=-1)
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
std::vector< int > m_intVariables
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
static const NekDouble kNekZeroTol
std::vector< std::pair< std::string, std::string > > SummaryList
EquationSystemFactory & GetEquationSystemFactory()
@ eAverage
averaged (or centred) flux
@ eLaxFriedrich
Lax-Friedrich flux.
@ eHLLC
Harten-Lax-Leer Contact wave flux.
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::vector< double > z(NPUPPER)
@ SIZE_TestType
Length of enum list.
const char *const TestTypeMap[]
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
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 Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
void Neg(int n, T *x, const int incx)
Negate x = -x.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Vvtvm(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvm (vector times vector minus vector): z = w*x - y
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
void Zero(int n, T *x, const int incx)
Zero vector.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
T Vamax(int n, const T *x, const int incx)
Return the maximum absolute element in x called vamax to avoid conflict with max.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
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.
scalarT< T > abs(scalarT< T > in)
scalarT< T > sqrt(scalarT< T > in)