38#include <boost/algorithm/string/predicate.hpp>
67 int shapedim =
m_fields[0]->GetShapeDimension();
69 for (
int j = 0; j < shapedim; ++j)
90 "No TESTTYPE defined in session.");
91 std::string TestTypeStr =
m_session->GetSolverInfo(
"TESTTYPE");
108 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
113 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
117 m_session->LoadParameter(
"Omega", Omegams, 7.292 * 0.00001);
118 m_Omega = Omegams * SecondToDay;
129 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
134 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
138 m_session->LoadParameter(
"Omega", Omegams, 7.292 * 0.00001);
139 m_Omega = Omegams * SecondToDay;
141 m_H0 = 133681.0 / (rad_earth * gms);
142 m_k2 = 10.0 / (rad_earth * gms);
148 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
153 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
158 m_u0 =
m_u0 * SecondToDay / rad_earth;
160 m_session->LoadParameter(
"Omega", Omegams, 7.292 * 0.00001);
161 m_Omega = Omegams * SecondToDay;
163 m_H0 = 5960.0 / rad_earth;
164 m_hs0 = 2000.0 / rad_earth;
170 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
175 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
178 m_u0 =
m_u0 * SecondToDay / rad_earth;
180 m_session->LoadParameter(
"Omega", Omegams, 7.292 * 0.00001);
181 m_Omega = Omegams * SecondToDay;
190 m_hbar = 120.0 / rad_earth;
192 std::cout <<
"m_theta0 = " <<
m_theta0
194 <<
", m_hbar = " <<
m_hbar << std::endl;
200 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
205 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
207 m_session->LoadParameter(
"Omega", Omegams, 7.292 * 0.00001);
208 m_Omega = Omegams * SecondToDay;
213 m_angfreq = 7.848 * 0.000001 * SecondToDay;
214 m_K = 7.848 * 0.000001 * SecondToDay;
237 ASSERTL0(
false,
"Implicit unsteady Advection not set up.");
259 for (i = 0; i < nfields; ++i)
263 nvariables = nfields;
275 for (i = 0; i < nvariables; ++i)
288 "Only one of IO_CheckTime and IO_CheckSteps "
292 bool doCheckTime =
false;
298 NekDouble Mass = 0.0, Energy = 0.0, Enstrophy = 0.0, Vorticity = 0.0;
301 for (
int i = 0; i < nvariables; ++i)
321 std::cout <<
"Steps: " << std::setw(8) << std::left << step + 1
323 <<
"Time: " << std::setw(12) << std::left <<
m_time;
325 std::stringstream ss;
326 ss << cpuTime <<
"s";
327 std::cout <<
" CPU Time: " << std::setw(8) << std::left << ss.str()
341 Energy = (
ComputeEnergy(fieldsprimitive[0], fieldsprimitive[1],
342 fieldsprimitive[2]) -
350 fieldsprimitive[2]) -
354 std::cout <<
"dMass = " << std::setw(8) << std::left << Mass <<
" "
355 <<
", dEnergy = " << std::setw(8) << std::left << Energy
357 <<
", dEnstrophy = " << std::setw(8) << std::left
359 <<
", dVorticity = " << std::setw(8) << std::left
360 << Vorticity << std::endl
370 for (i = 0; i < nvariables; ++i)
394 if (
m_session->GetComm()->GetRank() == 0)
400 <<
"CFL time-step : " <<
m_timestep << std::endl;
403 if (
m_session->GetSolverInfo(
"Driver") !=
"SteadyState")
405 std::cout <<
"Time-integration : " << intTime <<
"s" << std::endl;
409 for (i = 0; i < nvariables; ++i)
418 for (i = 0; i < nvariables; ++i)
430 int nvariables = inarray.size();
438 for (i = 0; i < nvariables; ++i)
451 for (i = 0; i < nvariables; ++i)
472 for (i = 0; i < nvariables; ++i)
474 m_fields[i]->MultiplyByElmtInvMass(modarray[i], modarray[i]);
475 m_fields[i]->BwdTrans(modarray[i], outarray[i]);
498 for (i = 0; i < nvariables; ++i)
500 physfield[i] = InField[i];
510 for (i = 0; i < nvariables; ++i)
521 Vmath::Vadd(ncoeffs, &tmp[0], 1, &OutField[i][0], 1,
530 for (i = 0; i < nvariables; ++i)
539 for (i = 0; i < nvariables; ++i)
542 m_fields[i]->AddFwdBwdTraceIntegral(numfluxFwd[i], numfluxBwd[i],
553 int ncoeffs = outarray[0].size();
554 int nq = physarray[0].size();
573 m_fields[0]->IProductWRTBase(tmp, tmpc);
575 Vmath::Vadd(ncoeffs, outarray[j + 1], 1, tmpc, 1, outarray[j + 1], 1);
583 int nq =
m_fields[0]->GetTotPoints();
594 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
597 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
610 Vmath::Vmul(nq, tmp, 1, physfield[1], 1, flux[1], 1);
613 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
623 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
636 Vmath::Vmul(nq, tmp, 1, physfield[2], 1, flux[0], 1);
639 Vmath::Vmul(nq, flux[0], 1, physfield[2], 1, flux[1], 1);
642 Vmath::Vmul(nq, flux[0], 1, physfield[1], 1, flux[0], 1);
689 for (i = 0; i < nvariables; ++i)
696 for (i = 0; i < nvariables; ++i)
698 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
722 for (k = 0; k < nTraceNumPoints; ++k)
746 for (k = 0; k < nTraceNumPoints; ++k)
749 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
750 Bwd[2][k], numfluxF, numfluxB);
768 denomFwd = eF1n * eF2t - eF2n * eF1t;
769 denomBwd = eB1n * eB2t - eB2n * eB1t;
771 numfluxFwd[0][k] = numfluxF[0];
772 numfluxFwd[1][k] = (1.0 / denomFwd) *
773 (eF2t * numfluxF[1] - eF2n * numfluxF[2]);
776 (-1.0 * eF1t * numfluxF[1] + eF1n * numfluxF[2]);
778 numfluxBwd[0][k] = 1.0 * numfluxB[0];
779 numfluxBwd[1][k] = (1.0 / denomBwd) *
780 (eB2t * numfluxB[1] - eB2n * numfluxB[2]);
783 (-1.0 * eB1t * numfluxB[1] + eB1n * numfluxB[2]);
795 for (k = 0; k < nTraceNumPoints; ++k)
800 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
801 Bwd[2][k], numfluxF, numfluxB);
807 Fwd[2][k], Bwd[0][k] + DepthFwd[k],
808 Bwd[1][k], Bwd[2][k], numfluxF, numfluxB);
814 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
815 Bwd[2][k], numfluxF, numfluxB);
820 for (i = 0; i < nvariables; ++i)
830 numfluxFwd[i][k] = tmpF0 + tmpF1 + numfluxF[indx + 2];
831 numfluxBwd[i][k] = tmpB0 + tmpB1 + numfluxB[indx + 2];
839 ASSERTL0(
false,
"populate switch statement for upwind flux");
866 hstarF = 0.5 * (cL + cR) + 0.25 * (uL - uRF);
870 hstarB = 0.5 * (cL + cR) + 0.25 * (uLB - uR);
881 numfluxF[0] = hfluxF;
882 numfluxF[1] = hufluxF;
883 numfluxF[2] = hvfluxF;
885 numfluxB[0] = hfluxB;
886 numfluxB[1] = hufluxB;
887 numfluxB[2] = hvfluxB;
904 SL = uL - cL *
sqrt(0.5 * ((hstar * hstar + hstar * hL) / (hL * hL)));
914 SR = uR + cR *
sqrt(0.5 * ((hstar * hstar + hstar * hR) / (hR * hR)));
921 if (fabs(hR * (uR - SR) - hL * (uL - SL)) <= 1.0e-15)
927 Sstar = (SL * hR * (uR - SR) - SR * hL * (uL - SL)) /
928 (hR * (uR - SR) - hL * (uL - SL));
934 huflux = uL * uL * hL + 0.5 * g * hL * hL;
935 hvflux = hL * uL * vL;
940 huflux = uR * uR * hR + 0.5 * g * hR * hR;
941 hvflux = hR * uR * vR;
945 if ((SL < 0) && (Sstar >= 0))
947 hC = hL * ((SL - uL) / (SL - Sstar));
951 hflux = hL * uL + SL * (hC - hL);
952 huflux = (uL * uL * hL + 0.5 * g * hL * hL) + SL * (huC - hL * uL);
953 hvflux = (uL * vL * hL) + SL * (hvC - hL * vL);
957 hC = hR * ((SR - uR) / (SR - Sstar));
961 hflux = hR * uR + SR * (hC - hR);
962 huflux = (uR * uR * hR + 0.5 * g * hR * hR) + SR * (huC - hR * uR);
963 hvflux = (uR * vR * hR) + SR * (hvC - hR * vR);
973 NekDouble MageF1, MageF2, MageB1, MageB2;
981 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
985 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
986 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
988 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
989 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
993 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
994 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
997 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
999 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
1004 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1005 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1007 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1008 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1012 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1013 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1016 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1018 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1028 NekDouble MageF1, MageF2, MageB1, MageB2;
1041 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1047 EigF[0] = velL -
sqrt(g * hL);
1049 EigF[2] = velL +
sqrt(g * hL);
1051 EigB[0] = velR -
sqrt(g * hR);
1053 EigB[2] = velR +
sqrt(g * hR);
1060 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1061 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1063 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1064 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1065 numfluxF[2] = 0.5 * lambdaF * (hL - hR);
1067 numfluxF[3] = 0.5 * (hL * uL * uL * MageF1 + hR * uRF * uRF * MageB1 +
1068 0.5 * g * (hL * hL + hR * hR));
1069 numfluxF[4] = 0.5 * (hL * uL * vL * MageF1 + hR * uRF * vRF * MageB1);
1070 numfluxF[5] = 0.5 * lambdaF * (uL * hL - uRF * hR);
1072 numfluxF[6] = 0.5 * (hL * uL * vL * MageF2 + hR * uRF * vRF * MageB2);
1073 numfluxF[7] = 0.5 * (hL * vL * vL * MageF2 + hR * vRF * vRF * MageB2 +
1074 0.5 * g * (hL * hL + hR * hR));
1075 numfluxF[8] = 0.5 * lambdaF * (vL * hL - vRF * hR);
1079 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1080 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1082 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1083 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1084 numfluxB[2] = 0.5 * lambdaB * (hL - hR);
1086 numfluxB[3] = 0.5 * (hR * uR * uR * MageB1 + hR * uLB * uLB * MageF1 +
1087 0.5 * g * (hR * hR + hL * hL));
1088 numfluxB[4] = 0.5 * (hR * uR * vR * MageB1 + hR * uLB * vLB * MageF1);
1089 numfluxB[5] = 0.5 * lambdaB * (uLB * hL - uR * hR);
1091 numfluxB[6] = 0.5 * (hR * uR * vR * MageB2 + hR * uLB * vLB * MageF2);
1092 numfluxB[7] = 0.5 * (hR * vR * vR * MageB2 + hR * vLB * vLB * MageF2 +
1093 0.5 * g * (hR * hR + hL * hL));
1094 numfluxB[8] = 0.5 * lambdaB * (vLB * hL - vR * hR);
1103 NekDouble MageF1, MageF2, MageB1, MageB2;
1116 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1123 SL = fabs(velL) +
sqrt(g * hL);
1124 SR = fabs(velR) +
sqrt(g * hR);
1138 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1139 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1141 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1142 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1143 numfluxF[2] = 0.5 * S * (hL - hR);
1146 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
1147 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1148 numfluxF[5] = 0.5 * S * (uL * hL - uRF * hR);
1150 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1152 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
1153 numfluxF[8] = 0.5 * S * (vL * hL - vRF * hR);
1157 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1158 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1160 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1161 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1162 numfluxB[2] = 0.5 * S * (hL - hR);
1165 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1166 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1167 numfluxB[5] = 0.5 * S * (uLB * hL - uR * hR);
1169 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1171 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1172 numfluxB[8] = 0.5 * S * (vLB * hL - vR * hR);
1181 NekDouble MF1x, MF1y, MF1z, MF2x, MF2y, MF2z;
1182 NekDouble MB1x, MB1y, MB1z, MB2x, MB2y, MB2z;
1201 MageF1 = MF1x * MF1x + MF1y * MF1y + MF1z * MF1z;
1202 MageF2 = MF2x * MF2x + MF2y * MF2y + MF2z * MF2z;
1203 MageB1 = MB1x * MB1x + MB1y * MB1y + MB1z * MB1z;
1204 MageB2 = MB2x * MB2x + MB2y * MB2y + MB2z * MB2z;
1206 eF1_cdot_eB1 = MF1x * MB1x + MF1y * MB1y + MF1z * MB1z;
1207 eF1_cdot_eB2 = MF1x * MB2x + MF1y * MB2y + MF1z * MB2z;
1208 eF2_cdot_eB1 = MF2x * MB1x + MF2y * MB1y + MF2z * MB1z;
1209 eF2_cdot_eB2 = MF2x * MB2x + MF2y * MB2y + MF2z * MB2z;
1216 int ncoeffs = outarray[0].size();
1217 int nq = physarray[0].size();
1252 m_fields[0]->IProductWRTBase(tmp, tmpc);
1253 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1261 int ncoeffs = outarray[0].size();
1262 int nq = physarray[0].size();
1278 m_fields[0]->IProductWRTBase(tmp, tmpc);
1280 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1291 int ncoeffs = outarray[0].size();
1292 int nq = physarray[0].size();
1309 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e0, 1, Rott1, 1);
1310 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e1, 1, Rott2, 1);
1311 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e0, 1, Rott1, 1, Rott1, 1);
1312 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e1, 1, Rott2, 1, Rott2, 1);
1315 Vmath::Vmul(nq, &h[0], 1, &Rott1[0], 1, &Rott1[0], 1);
1316 Vmath::Vmul(nq, &h[0], 1, &Rott2[0], 1, &Rott2[0], 1);
1323 m_fields[0]->IProductWRTBase(Rott1, tmpc1);
1324 m_fields[0]->IProductWRTBase(Rott2, tmpc2);
1326 Vmath::Vadd(ncoeffs, tmpc1, 1, outarray[1], 1, outarray[1], 1);
1327 Vmath::Vadd(ncoeffs, tmpc2, 1, outarray[2], 1, outarray[2], 1);
1335 const int indm,
const int indk,
1340 int nq =
m_fields[0]->GetNpoints();
1354 Vmath::Vmul(nq, &physarray[j + 1][0], 1, &tmpderiv[0], 1,
1358 &outarray[0], 1, &outarray[0], 1);
1384 ASSERTL0(
false,
"Unknown projection scheme");
1397 for (
int n = 0; n <
m_fields[0]->GetBndConditions().size(); ++n)
1401 if (
m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
"eMG")
1405 ASSERTL0(
false,
"Illegal dimension");
1414 if (
m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
"eWall")
1418 ASSERTL0(
false,
"Illegal dimension");
1427 if (
m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
1430 for (
int i = 0; i < nvariables; ++i)
1432 m_fields[i]->EvaluateBoundaryConditions(time);
1435 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
1445 int nvariables = physarray.size();
1449 for (i = 0; i < nvariables; ++i)
1452 m_fields[i]->ExtractTracePhys(physarray[i], Fwd0[i]);
1456 for (i = 0; i < nvariables; ++i)
1459 m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
1464 int e, id1, id2, npts;
1470 ->GetBndCondExpansions()[bcRegion]
1473 id1 =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
1474 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
1475 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
1493 &tmp_n[0], 1, &tmp_n[0], 1);
1498 1, &tmp_t[0], 1, &tmp_t[0], 1);
1518 &tmp_u[0], 1, &tmp_u[0], 1);
1519 Vmath::Vdiv(npts, &tmp_u[0], 1, &denom[0], 1, &tmp_u[0], 1);
1526 &tmp_v[0], 1, &tmp_v[0], 1);
1527 Vmath::Vdiv(npts, &tmp_v[0], 1, &denom[0], 1, &tmp_v[0], 1);
1534 ASSERTL0(
false,
"Illegal expansion dimension");
1538 for (i = 0; i < nvariables; ++i)
1542 ->GetBndCondExpansions()[bcRegion]
1543 ->UpdatePhys())[id1],
1561 for (
int i = 0; i <
m_fields.size(); ++i)
1593 NekDouble sin_varphi, cos_varphi, sin_theta, cos_theta;
1595 for (
int j = 0; j < nq; ++j)
1602 sin_theta, cos_theta);
1621 NekDouble phi, theta, sin_varphi, cos_varphi, sin_theta, cos_theta;
1627 thetac =
m_pi / 6.0;
1629 for (
int j = 0; j < nq; ++j)
1636 sin_theta, cos_theta);
1638 if ((
std::abs(sin(phic) - sin_varphi) +
1639 std::abs(sin(thetac) - sin_theta)) < Tol)
1641 std::cout <<
"A point " << j
1642 <<
" is coincient with the singularity "
1647 phi = atan2(sin_varphi, cos_varphi);
1648 theta = atan2(sin_theta, cos_theta);
1651 dist2 = (phi - phic) * (phi - phic) +
1652 (theta - thetac) * (theta - thetac);
1654 if (dist2 > hRad * hRad)
1656 dist2 = hRad * hRad;
1664 std::cout <<
"No point is coincident with the singularity point"
1672 for (
int j = 0; j < nq; ++j)
1703 std::cout <<
"Water Depth (m_depth) was generated with mag = "
1742 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1754 for (
int j = 0; j < nq; ++j)
1765 tmp = -1.0 * cos_varphi * cos_theta * sin(
m_alpha) +
1767 outarray[j] = 2.0 *
m_Omega * tmp;
1776 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1785 for (
int j = 0; j < nq; ++j)
1794 outarray[j] = 2.0 *
m_Omega * sin_theta;
1799 bool dumpInitialConditions,
1800 [[maybe_unused]]
const int domain)
1822 for (
int i = 0; i <
m_fields.size(); ++i)
1855 for (
int i = 0; i <
m_fields.size(); ++i)
1884 for (
int i = 0; i <
m_fields.size(); ++i)
1913 for (
int i = 0; i <
m_fields.size(); ++i)
1942 for (
int i = 0; i <
m_fields.size(); ++i)
1971 for (
int i = 0; i <
m_fields.size(); ++i)
1984 if (dumpInitialConditions)
1999 int nq =
m_fields[0]->GetNpoints();
2005 m_fields[0]->GetCoords(x0, x1, x2);
2018 for (
int i = 0; i < nq; ++i)
2020 eta0[i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2021 (1.0 / cosh(0.395 * x0[i]))) *
2022 (3.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2023 exp(-0.5 * x1[i] * x1[i]);
2024 uvec[0][i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2025 (1.0 / cosh(0.395 * x0[i]))) *
2026 (-9.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2027 exp(-0.5 * x1[i] * x1[i]);
2028 uvec[1][i] = (-2.0 * 0.395 * tanh(0.395 * x0[i])) *
2029 (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2030 (1.0 / cosh(0.395 * x0[i]))) *
2031 (2.0 * x1[i]) * exp(-0.5 * x1[i] * x1[i]);
2064 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2083 for (
int j = 0; j < nq; ++j)
2094 tmp = -1.0 * cos_varphi * cos_theta * sin(
m_alpha) +
2101 sin_theta * cos_varphi * sin(
m_alpha));
2104 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2105 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2106 uvec[2][j] = vhat * cos_theta;
2152 int nq =
m_fields[0]->GetTotPoints();
2164 int nq =
m_fields[0]->GetTotPoints();
2193 int nq =
m_fields[0]->GetTotPoints();
2224 int nq =
m_fields[0]->GetTotPoints();
2238 Vmath::Vadd(nq, tmp, 1, Vorticity, 1, Vorticity, 1);
2243 int nq =
m_fields[0]->GetNpoints();
2265 1, &velcoeff[0], 1, &velcoeff[0], 1);
2269 m_fields[0]->PhysDeriv(velcoeff, Dtmp0, Dtmp1, Dtmp2);
2279 1, &vellc[0], 1, &vellc[0], 1);
2286 1, &vellc[0], 1, &vellc[0], 1);
2293 1, &vellc[0], 1, &vellc[0], 1);
2306 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2327 for (
int j = 0; j < nq; ++j)
2340 cos_varphi * cos(
m_Omega * time) - sin_varphi * sin(
m_Omega * time);
2342 sin_varphi * cos(
m_Omega * time) + cos_varphi * sin(
m_Omega * time);
2343 tmp = -1.0 * TR * sin(
m_alpha) * cos_theta + cos(
m_alpha) * sin_theta;
2345 eta[j] = -1.0 * (
m_u0 * tmp +
m_Omega * sin_theta) *
2348 eta[j] = 0.5 * eta[j] /
m_g;
2356 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2357 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2358 uvec[2][j] = vhat * cos_theta;
2409 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2428 for (
int j = 0; j < nq; ++j)
2439 sin_theta * sin_theta;
2441 uhat =
m_u0 * cos_theta;
2444 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2445 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2446 uvec[2][j] = vhat * cos_theta;
2497 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2519 for (
int j = 0; j < nq; ++j)
2528 Ttheta = atan2(sin_theta, cos_theta);
2529 Tphi = atan2(sin_varphi, cos_varphi);
2535 dth = Ttheta / Nint;
2536 eta[j] = dth * 0.5 *
2538 for (
int i = 1; i < Nint - 1; i++)
2543 eta[j] = (-1.0 /
m_g) * eta[j];
2548 eta[j] = eta[j] +
m_hbar * cos_theta * exp(-9.0 * Tphi * Tphi) *
2549 exp(-225.0 * (
m_pi / 4.0 - Ttheta) *
2550 (
m_pi / 4.0 - Ttheta));
2553 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2554 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2555 uvec[2][j] = vhat * cos_theta;
2603 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2624 NekDouble cos2theta, cosRtheta, cos6theta, cos2Rtheta, cosRm1theta;
2625 NekDouble cos2phi, cos4phi, sin4phi, cos8phi;
2631 phi0 = 40.0 *
m_pi / 180.0;
2632 theta0 = 50.0 *
m_pi / 180.0;
2634 x0d = cos(phi0) * cos(theta0);
2635 y0d = sin(phi0) * cos(theta0);
2638 for (
int j = 0; j < nq; ++j)
2651 cos2theta = cos_theta * cos_theta;
2652 cosRm1theta = cos_theta * cos2theta;
2653 cosRtheta = cos2theta * cos2theta;
2654 cos6theta = cos2theta * cosRtheta;
2655 cos2Rtheta = cosRtheta * cosRtheta;
2658 Ath = tmp * cos2theta +
2659 0.25 *
m_K *
m_K * cos6theta *
2660 ((R + 1.0) * cosRtheta + (2 * R * R - R - 2.0) * cos2theta -
2664 Bth = tmp * cosRtheta *
2665 ((R * R + 2 * R + 2) - (R + 1.0) * (R + 1.0) * cos2theta);
2668 0.25 *
m_K *
m_K * cos2Rtheta * ((R + 1.0) * cos2theta - (R + 2.0));
2671 cos2phi = 2.0 * cos_varphi * cos_varphi - 1.0;
2672 cos4phi = 2.0 * cos2phi * cos2phi - 1.0;
2673 cos8phi = 2.0 * cos4phi * cos4phi - 1.0;
2676 sin4phi = 4.0 * sin_varphi * cos_varphi * cos2phi;
2678 eta[j] =
m_H0 + (1.0 /
m_g) * (Ath + Bth * cos4phi + Cth * cos8phi);
2684 (1.0 + (1.0 / 40.0) * (x0j * x0d + x1j * y0d + x2j * z0d));
2690 m_K * cosRm1theta * (R * sin_theta * sin_theta - cos2theta) *
2692 vhat = -1.0 *
m_K * R * cosRm1theta * sin_theta * sin4phi;
2694 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2695 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2696 uvec[2][j] = vhat * cos_theta;
2749 f = 2.0 *
m_Omega * sin(theta);
2751 dh = f * uphi + tan(theta) * uphi * uphi;
2776 int nq =
m_fields[0]->GetTotPoints();
2777 int ncoeffs =
m_fields[0]->GetNcoeffs();
2779 NekDouble rad_earth = 6.37122 * 1000000;
2784 std::vector<std::string> variables(nvariables);
2786 variables[0] =
"eta";
2787 variables[1] =
"hstar";
2788 variables[2] =
"vorticity";
2789 variables[3] =
"ux";
2790 variables[4] =
"uy";
2791 variables[5] =
"uz";
2792 variables[6] =
"null";
2796 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvariables);
2797 for (
int i = 0; i < nvariables; ++i)
2804 &fieldphys[0][0], 1);
2811 Vmath::Smul(nq, rad_earth, &fieldphys[1][0], 1, &fieldphys[1][0], 1);
2825 &fieldphys[k + indx][0], 1);
2827 &fieldphys[k + indx][0], 1, &fieldphys[k + indx][0], 1);
2830 for (
int i = 0; i < nvariables; ++i)
2832 m_fields[0]->FwdTrans(fieldphys[i], fieldcoeffs[i]);
2845 if (physin[0].get() == physout[0].get())
2849 for (
int i = 0; i < 3; ++i)
2860 Vmath::Vdiv(nq, tmp[1], 1, tmp[0], 1, physout[1], 1);
2863 Vmath::Vdiv(nq, tmp[2], 1, tmp[0], 1, physout[2], 1);
2871 Vmath::Vdiv(nq, physin[1], 1, physin[0], 1, physout[1], 1);
2874 Vmath::Vdiv(nq, physin[2], 1, physin[0], 1, physout[2], 1);
2885 if (physin[0].get() == physout[0].get())
2889 for (
int i = 0; i < 3; ++i)
2900 Vmath::Vmul(nq, physout[0], 1, tmp[1], 1, physout[1], 1);
2903 Vmath::Vmul(nq, physout[0], 1, tmp[2], 1, physout[2], 1);
2911 Vmath::Vmul(nq, physout[0], 1, physin[1], 1, physout[1], 1);
2914 Vmath::Vmul(nq, physout[0], 1, physin[2], 1, physout[2], 1);
2958 int nq =
m_fields[0]->GetTotPoints();
2963 NekDouble theta, phi, sin_theta, cos_theta, sin_varphi, cos_varphi;
2988 for (k = 0; k < nq; ++k)
2994 Re =
sqrt(xp * xp + yp * yp + zp * zp);
3001 theta = atan2(sin_theta, cos_theta);
3002 phi = atan2(sin_varphi, cos_varphi);
3004 cosntheta3 = cos(n * theta) * cos(n * theta) * cos(n * theta);
3006 beta_theta = -4.0 * n * cosntheta3 * cos(m * phi) * sin(n * theta) / Re;
3007 beta_phi = -m * cosntheta3 * sin(m * phi) / Re;
3009 thetax = -1.0 * cos_varphi * sin_theta;
3010 thetay = -1.0 * sin_varphi * sin_theta;
3013 phix = -1.0 * sin_varphi;
3017 uvec[0][k] = alpha * (beta_theta * thetax + beta_phi * phix);
3018 uvec[1][k] = alpha * (beta_theta * thetay + beta_phi * phiy);
3019 uvec[2][k] = alpha * (beta_theta * thetaz + beta_phi * phiz);
3021 vorticityexact[k] = -4.0 * n / Re / Re * cos_theta * cos_theta *
3022 cos_varphi * cos(m * phi) * sin(n * theta);
3028 std::cout <<
"chi migi1" << std::endl;
3032 Vmath::Vsub(nq, vorticityexact, 1, vorticitycompt, 1, vorticitycompt, 1);
3034 std::cout <<
"Vorticity: L2 error = " <<
AvgAbsInt(vorticitycompt)
3035 <<
", Linf error = " <<
Vmath::Vamax(nq, vorticitycompt, 1)
3043 int nq =
m_fields[field]->GetNpoints();
3048 if (
m_fields[field]->GetPhysState() ==
false)
3066 &exactsolution[0], 1, &exactsolution[0], 1);
3067 Vmath::Vabs(nq, exactsolution, 1, exactsolution, 1);
3069 L2error = (
m_fields[0]->Integral(exactsolution)) / L2exact;
3086 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3089 L2exact =
m_fields[1]->Integral(tmp);
3098 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3101 L2error = (
m_fields[1]->Integral(tmp)) / L2exact;
3115 if (Normalised ==
true)
3122 L2error =
sqrt(L2error * L2error / Vol);
3141 if (
m_fields[field]->GetPhysState() ==
false)
3147 int nq =
m_fields[field]->GetNpoints();
3171 1, &exactsolution[0], 1);
3193 Vmath::Vsub(nq, &exactu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3194 Vmath::Vsub(nq, &exactv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3196 Vmath::Vmul(nq, &tmpu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3197 Vmath::Vmul(nq, &tmpv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3199 Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &Lerr[0], 1);
3203 Vmath::Vmul(nq, &exactu[0], 1, &exactu[0], 1, &tmpu[0], 1);
3204 Vmath::Vmul(nq, &exactv[0], 1, &exactv[0], 1, &tmpv[0], 1);
3205 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)
Session reader.
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.
~MMFSWE() override
Destructor.
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).
SOLVER_UTILS_EXPORT int GetTraceNpoints()
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 GetTraceTotPoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT NekDouble LinfError(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
Linf error computation.
SOLVER_UTILS_EXPORT void EvaluateExactSolution(int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
Evaluates an exact solution.
NekDouble m_fintime
Finish time of the simulation.
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 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 GetNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
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.
int m_checksteps
Number of steps between checkpoints.
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()
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)