38 #include <boost/algorithm/string/predicate.hpp>
39 #include <boost/core/ignore_unused.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.");
260 for (i = 0; i < nfields; ++i)
264 nvariables = nfields;
276 for (i = 0; i < nvariables; ++i)
289 "Only one of IO_CheckTime and IO_CheckSteps "
293 bool doCheckTime =
false;
299 NekDouble Mass = 0.0, Energy = 0.0, Enstrophy = 0.0, Vorticity = 0.0;
302 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)
429 boost::ignore_unused(time);
432 int nvariables = inarray.size();
440 for (i = 0; i < nvariables; ++i)
453 for (i = 0; i < nvariables; ++i)
474 for (i = 0; i < nvariables; ++i)
476 m_fields[i]->MultiplyByElmtInvMass(modarray[i], modarray[i]);
477 m_fields[i]->BwdTrans(modarray[i], outarray[i]);
500 for (i = 0; i < nvariables; ++i)
502 physfield[i] = InField[i];
512 for (i = 0; i < nvariables; ++i)
523 Vmath::Vadd(ncoeffs, &tmp[0], 1, &OutField[i][0], 1,
532 for (i = 0; i < nvariables; ++i)
541 for (i = 0; i < nvariables; ++i)
544 m_fields[i]->AddFwdBwdTraceIntegral(numfluxFwd[i], numfluxBwd[i],
555 int ncoeffs = outarray[0].size();
556 int nq = physarray[0].size();
575 m_fields[0]->IProductWRTBase(tmp, tmpc);
577 Vmath::Vadd(ncoeffs, outarray[j + 1], 1, tmpc, 1, outarray[j + 1], 1);
585 int nq =
m_fields[0]->GetTotPoints();
596 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
599 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
612 Vmath::Vmul(nq, tmp, 1, physfield[1], 1, flux[1], 1);
615 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
625 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
638 Vmath::Vmul(nq, tmp, 1, physfield[2], 1, flux[0], 1);
641 Vmath::Vmul(nq, flux[0], 1, physfield[2], 1, flux[1], 1);
644 Vmath::Vmul(nq, flux[0], 1, physfield[1], 1, flux[0], 1);
691 for (i = 0; i < nvariables; ++i)
698 for (i = 0; i < nvariables; ++i)
700 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
724 for (k = 0; k < nTraceNumPoints; ++k)
748 for (k = 0; k < nTraceNumPoints; ++k)
751 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
752 Bwd[2][k], numfluxF, numfluxB);
770 denomFwd = eF1n * eF2t - eF2n * eF1t;
771 denomBwd = eB1n * eB2t - eB2n * eB1t;
773 numfluxFwd[0][k] = numfluxF[0];
774 numfluxFwd[1][k] = (1.0 / denomFwd) *
775 (eF2t * numfluxF[1] - eF2n * numfluxF[2]);
778 (-1.0 * eF1t * numfluxF[1] + eF1n * numfluxF[2]);
780 numfluxBwd[0][k] = 1.0 * numfluxB[0];
781 numfluxBwd[1][k] = (1.0 / denomBwd) *
782 (eB2t * numfluxB[1] - eB2n * numfluxB[2]);
785 (-1.0 * eB1t * numfluxB[1] + eB1n * numfluxB[2]);
797 for (k = 0; k < nTraceNumPoints; ++k)
802 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
803 Bwd[2][k], numfluxF, numfluxB);
809 Fwd[2][k], Bwd[0][k] + DepthFwd[k],
810 Bwd[1][k], Bwd[2][k], numfluxF, numfluxB);
816 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
817 Bwd[2][k], numfluxF, numfluxB);
822 for (i = 0; i < nvariables; ++i)
832 numfluxFwd[i][k] = tmpF0 + tmpF1 + numfluxF[indx + 2];
833 numfluxBwd[i][k] = tmpB0 + tmpB1 + numfluxB[indx + 2];
841 ASSERTL0(
false,
"populate switch statement for upwind flux");
852 boost::ignore_unused(index);
869 hstarF = 0.5 * (cL + cR) + 0.25 * (uL - uRF);
873 hstarB = 0.5 * (cL + cR) + 0.25 * (uLB - uR);
884 numfluxF[0] = hfluxF;
885 numfluxF[1] = hufluxF;
886 numfluxF[2] = hvfluxF;
888 numfluxB[0] = hfluxB;
889 numfluxB[1] = hufluxB;
890 numfluxB[2] = hvfluxB;
906 SL = uL - cL *
sqrt(0.5 * ((hstar * hstar + hstar * hL) / (hL * hL)));
912 SR = uR + cR *
sqrt(0.5 * ((hstar * hstar + hstar * hR) / (hR * hR)));
916 if (fabs(hR * (uR - SR) - hL * (uL - SL)) <= 1.0e-15)
919 Sstar = (SL * hR * (uR - SR) - SR * hL * (uL - SL)) /
920 (hR * (uR - SR) - hL * (uL - SL));
925 huflux = uL * uL * hL + 0.5 * g * hL * hL;
926 hvflux = hL * uL * vL;
931 huflux = uR * uR * hR + 0.5 * g * hR * hR;
932 hvflux = hR * uR * vR;
936 if ((SL < 0) && (Sstar >= 0))
938 hC = hL * ((SL - uL) / (SL - Sstar));
942 hflux = hL * uL + SL * (hC - hL);
943 huflux = (uL * uL * hL + 0.5 * g * hL * hL) + SL * (huC - hL * uL);
944 hvflux = (uL * vL * hL) + SL * (hvC - hL * vL);
948 hC = hR * ((SR - uR) / (SR - Sstar));
952 hflux = hR * uR + SR * (hC - hR);
953 huflux = (uR * uR * hR + 0.5 * g * hR * hR) + SR * (huC - hR * uR);
954 hvflux = (uR * vR * hR) + SR * (hvC - hR * vR);
964 NekDouble MageF1, MageF2, MageB1, MageB2;
972 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
976 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
977 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
979 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
980 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
984 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
985 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
988 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
990 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
995 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
996 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
998 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
999 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1003 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1004 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1007 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1009 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1019 NekDouble MageF1, MageF2, MageB1, MageB2;
1032 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1038 EigF[0] = velL -
sqrt(g * hL);
1040 EigF[2] = velL +
sqrt(g * hL);
1042 EigB[0] = velR -
sqrt(g * hR);
1044 EigB[2] = velR +
sqrt(g * hR);
1051 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1052 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1054 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1055 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1056 numfluxF[2] = 0.5 * lambdaF * (hL - hR);
1058 numfluxF[3] = 0.5 * (hL * uL * uL * MageF1 + hR * uRF * uRF * MageB1 +
1059 0.5 * g * (hL * hL + hR * hR));
1060 numfluxF[4] = 0.5 * (hL * uL * vL * MageF1 + hR * uRF * vRF * MageB1);
1061 numfluxF[5] = 0.5 * lambdaF * (uL * hL - uRF * hR);
1063 numfluxF[6] = 0.5 * (hL * uL * vL * MageF2 + hR * uRF * vRF * MageB2);
1064 numfluxF[7] = 0.5 * (hL * vL * vL * MageF2 + hR * vRF * vRF * MageB2 +
1065 0.5 * g * (hL * hL + hR * hR));
1066 numfluxF[8] = 0.5 * lambdaF * (vL * hL - vRF * hR);
1070 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1071 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1073 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1074 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1075 numfluxB[2] = 0.5 * lambdaB * (hL - hR);
1077 numfluxB[3] = 0.5 * (hR * uR * uR * MageB1 + hR * uLB * uLB * MageF1 +
1078 0.5 * g * (hR * hR + hL * hL));
1079 numfluxB[4] = 0.5 * (hR * uR * vR * MageB1 + hR * uLB * vLB * MageF1);
1080 numfluxB[5] = 0.5 * lambdaB * (uLB * hL - uR * hR);
1082 numfluxB[6] = 0.5 * (hR * uR * vR * MageB2 + hR * uLB * vLB * MageF2);
1083 numfluxB[7] = 0.5 * (hR * vR * vR * MageB2 + hR * vLB * vLB * MageF2 +
1084 0.5 * g * (hR * hR + hL * hL));
1085 numfluxB[8] = 0.5 * lambdaB * (vLB * hL - vR * hR);
1094 NekDouble MageF1, MageF2, MageB1, MageB2;
1107 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1114 SL = fabs(velL) +
sqrt(g * hL);
1115 SR = fabs(velR) +
sqrt(g * hR);
1125 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1126 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1128 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1129 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1130 numfluxF[2] = 0.5 * S * (hL - hR);
1133 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
1134 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1135 numfluxF[5] = 0.5 * S * (uL * hL - uRF * hR);
1137 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1139 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
1140 numfluxF[8] = 0.5 * S * (vL * hL - vRF * hR);
1144 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1145 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1147 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1148 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1149 numfluxB[2] = 0.5 * S * (hL - hR);
1152 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1153 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1154 numfluxB[5] = 0.5 * S * (uLB * hL - uR * hR);
1156 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1158 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1159 numfluxB[8] = 0.5 * S * (vLB * hL - vR * hR);
1168 NekDouble MF1x, MF1y, MF1z, MF2x, MF2y, MF2z;
1169 NekDouble MB1x, MB1y, MB1z, MB2x, MB2y, MB2z;
1188 MageF1 = MF1x * MF1x + MF1y * MF1y + MF1z * MF1z;
1189 MageF2 = MF2x * MF2x + MF2y * MF2y + MF2z * MF2z;
1190 MageB1 = MB1x * MB1x + MB1y * MB1y + MB1z * MB1z;
1191 MageB2 = MB2x * MB2x + MB2y * MB2y + MB2z * MB2z;
1193 eF1_cdot_eB1 = MF1x * MB1x + MF1y * MB1y + MF1z * MB1z;
1194 eF1_cdot_eB2 = MF1x * MB2x + MF1y * MB2y + MF1z * MB2z;
1195 eF2_cdot_eB1 = MF2x * MB1x + MF2y * MB1y + MF2z * MB1z;
1196 eF2_cdot_eB2 = MF2x * MB2x + MF2y * MB2y + MF2z * MB2z;
1203 int ncoeffs = outarray[0].size();
1204 int nq = physarray[0].size();
1239 m_fields[0]->IProductWRTBase(tmp, tmpc);
1240 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1248 int ncoeffs = outarray[0].size();
1249 int nq = physarray[0].size();
1265 m_fields[0]->IProductWRTBase(tmp, tmpc);
1267 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1278 int ncoeffs = outarray[0].size();
1279 int nq = physarray[0].size();
1296 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e0, 1, Rott1, 1);
1297 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e1, 1, Rott2, 1);
1298 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e0, 1, Rott1, 1, Rott1, 1);
1299 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e1, 1, Rott2, 1, Rott2, 1);
1302 Vmath::Vmul(nq, &h[0], 1, &Rott1[0], 1, &Rott1[0], 1);
1303 Vmath::Vmul(nq, &h[0], 1, &Rott2[0], 1, &Rott2[0], 1);
1310 m_fields[0]->IProductWRTBase(Rott1, tmpc1);
1311 m_fields[0]->IProductWRTBase(Rott2, tmpc2);
1313 Vmath::Vadd(ncoeffs, tmpc1, 1, outarray[1], 1, outarray[1], 1);
1314 Vmath::Vadd(ncoeffs, tmpc2, 1, outarray[2], 1, outarray[2], 1);
1322 const int indm,
const int indk,
1327 int nq =
m_fields[0]->GetNpoints();
1341 Vmath::Vmul(nq, &physarray[j + 1][0], 1, &tmpderiv[0], 1,
1345 &outarray[0], 1, &outarray[0], 1);
1371 ASSERTL0(
false,
"Unknown projection scheme");
1384 for (
int n = 0; n <
m_fields[0]->GetBndConditions().size(); ++n)
1388 if (
m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
"eMG")
1392 ASSERTL0(
false,
"Illegal dimension");
1401 if (
m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
"eWall")
1405 ASSERTL0(
false,
"Illegal dimension");
1414 if (
m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
1417 for (
int i = 0; i < nvariables; ++i)
1419 m_fields[i]->EvaluateBoundaryConditions(time);
1422 cnt +=
m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
1432 int nvariables = physarray.size();
1436 for (i = 0; i < nvariables; ++i)
1439 m_fields[i]->ExtractTracePhys(physarray[i], Fwd0[i]);
1443 for (i = 0; i < nvariables; ++i)
1446 m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
1451 int e, id1, id2, npts;
1457 ->GetBndCondExpansions()[bcRegion]
1460 id1 =
m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
1461 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
1462 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
1480 &tmp_n[0], 1, &tmp_n[0], 1);
1485 1, &tmp_t[0], 1, &tmp_t[0], 1);
1505 &tmp_u[0], 1, &tmp_u[0], 1);
1506 Vmath::Vdiv(npts, &tmp_u[0], 1, &denom[0], 1, &tmp_u[0], 1);
1513 &tmp_v[0], 1, &tmp_v[0], 1);
1514 Vmath::Vdiv(npts, &tmp_v[0], 1, &denom[0], 1, &tmp_v[0], 1);
1521 ASSERTL0(
false,
"Illegal expansion dimension");
1525 for (i = 0; i < nvariables; ++i)
1529 ->GetBndCondExpansions()[bcRegion]
1530 ->UpdatePhys())[id1],
1548 for (
int i = 0; i <
m_fields.size(); ++i)
1580 NekDouble sin_varphi, cos_varphi, sin_theta, cos_theta;
1582 for (
int j = 0; j < nq; ++j)
1589 sin_theta, cos_theta);
1608 NekDouble phi, theta, sin_varphi, cos_varphi, sin_theta, cos_theta;
1614 thetac =
m_pi / 6.0;
1616 for (
int j = 0; j < nq; ++j)
1623 sin_theta, cos_theta);
1625 if ((
std::abs(sin(phic) - sin_varphi) +
1626 std::abs(sin(thetac) - sin_theta)) < Tol)
1628 std::cout <<
"A point " << j
1629 <<
" is coincient with the singularity "
1634 phi = atan2(sin_varphi, cos_varphi);
1635 theta = atan2(sin_theta, cos_theta);
1638 dist2 = (phi - phic) * (phi - phic) +
1639 (theta - thetac) * (theta - thetac);
1641 if (dist2 > hRad * hRad)
1643 dist2 = hRad * hRad;
1651 std::cout <<
"No point is coincident with the singularity point"
1659 for (
int j = 0; j < nq; ++j)
1690 std::cout <<
"Water Depth (m_depth) was generated with mag = "
1729 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1741 for (
int j = 0; j < nq; ++j)
1752 tmp = -1.0 * cos_varphi * cos_theta * sin(
m_alpha) +
1754 outarray[j] = 2.0 *
m_Omega * tmp;
1763 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1772 for (
int j = 0; j < nq; ++j)
1781 outarray[j] = 2.0 *
m_Omega * sin_theta;
1786 bool dumpInitialConditions,
1789 boost::ignore_unused(domain);
1811 for (
int i = 0; i <
m_fields.size(); ++i)
1844 for (
int i = 0; i <
m_fields.size(); ++i)
1873 for (
int i = 0; i <
m_fields.size(); ++i)
1902 for (
int i = 0; i <
m_fields.size(); ++i)
1931 for (
int i = 0; i <
m_fields.size(); ++i)
1960 for (
int i = 0; i <
m_fields.size(); ++i)
1973 if (dumpInitialConditions)
1987 boost::ignore_unused(time);
1989 int nq =
m_fields[0]->GetNpoints();
1995 m_fields[0]->GetCoords(x0, x1, x2);
2008 for (
int i = 0; i < nq; ++i)
2010 eta0[i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2011 (1.0 / cosh(0.395 * x0[i]))) *
2012 (3.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2013 exp(-0.5 * x1[i] * x1[i]);
2014 uvec[0][i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2015 (1.0 / cosh(0.395 * x0[i]))) *
2016 (-9.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2017 exp(-0.5 * x1[i] * x1[i]);
2018 uvec[1][i] = (-2.0 * 0.395 * tanh(0.395 * x0[i])) *
2019 (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2020 (1.0 / cosh(0.395 * x0[i]))) *
2021 (2.0 * x1[i]) * exp(-0.5 * x1[i] * x1[i]);
2054 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2073 for (
int j = 0; j < nq; ++j)
2084 tmp = -1.0 * cos_varphi * cos_theta * sin(
m_alpha) +
2091 sin_theta * cos_varphi * sin(
m_alpha));
2094 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2095 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2096 uvec[2][j] = vhat * cos_theta;
2142 int nq =
m_fields[0]->GetTotPoints();
2154 int nq =
m_fields[0]->GetTotPoints();
2183 int nq =
m_fields[0]->GetTotPoints();
2214 int nq =
m_fields[0]->GetTotPoints();
2228 Vmath::Vadd(nq, tmp, 1, Vorticity, 1, Vorticity, 1);
2233 int nq =
m_fields[0]->GetNpoints();
2255 1, &velcoeff[0], 1, &velcoeff[0], 1);
2259 m_fields[0]->PhysDeriv(velcoeff, Dtmp0, Dtmp1, Dtmp2);
2269 1, &vellc[0], 1, &vellc[0], 1);
2276 1, &vellc[0], 1, &vellc[0], 1);
2283 1, &vellc[0], 1, &vellc[0], 1);
2296 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2317 for (
int j = 0; j < nq; ++j)
2330 cos_varphi * cos(
m_Omega * time) - sin_varphi * sin(
m_Omega * time);
2332 sin_varphi * cos(
m_Omega * time) + cos_varphi * sin(
m_Omega * time);
2333 tmp = -1.0 * TR * sin(
m_alpha) * cos_theta + cos(
m_alpha) * sin_theta;
2335 eta[j] = -1.0 * (
m_u0 * tmp +
m_Omega * sin_theta) *
2338 eta[j] = 0.5 * eta[j] /
m_g;
2346 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2347 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2348 uvec[2][j] = vhat * cos_theta;
2395 boost::ignore_unused(time);
2400 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2419 for (
int j = 0; j < nq; ++j)
2430 sin_theta * sin_theta;
2432 uhat =
m_u0 * cos_theta;
2435 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2436 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2437 uvec[2][j] = vhat * cos_theta;
2484 boost::ignore_unused(time);
2489 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2511 for (
int j = 0; j < nq; ++j)
2520 Ttheta = atan2(sin_theta, cos_theta);
2521 Tphi = atan2(sin_varphi, cos_varphi);
2527 dth = Ttheta / Nint;
2528 eta[j] = dth * 0.5 *
2530 for (
int i = 1; i < Nint - 1; i++)
2535 eta[j] = (-1.0 /
m_g) * eta[j];
2540 eta[j] = eta[j] +
m_hbar * cos_theta * exp(-9.0 * Tphi * Tphi) *
2541 exp(-225.0 * (
m_pi / 4.0 - Ttheta) *
2542 (
m_pi / 4.0 - Ttheta));
2545 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2546 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2547 uvec[2][j] = vhat * cos_theta;
2595 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2616 NekDouble cos2theta, cosRtheta, cos6theta, cos2Rtheta, cosRm1theta;
2617 NekDouble cos2phi, cos4phi, sin4phi, cos8phi;
2623 phi0 = 40.0 *
m_pi / 180.0;
2624 theta0 = 50.0 *
m_pi / 180.0;
2626 x0d = cos(phi0) * cos(theta0);
2627 y0d = sin(phi0) * cos(theta0);
2630 for (
int j = 0; j < nq; ++j)
2643 cos2theta = cos_theta * cos_theta;
2644 cosRm1theta = cos_theta * cos2theta;
2645 cosRtheta = cos2theta * cos2theta;
2646 cos6theta = cos2theta * cosRtheta;
2647 cos2Rtheta = cosRtheta * cosRtheta;
2650 Ath = tmp * cos2theta +
2651 0.25 *
m_K *
m_K * cos6theta *
2652 ((R + 1.0) * cosRtheta + (2 * R * R - R - 2.0) * cos2theta -
2656 Bth = tmp * cosRtheta *
2657 ((R * R + 2 * R + 2) - (R + 1.0) * (R + 1.0) * cos2theta);
2660 0.25 *
m_K *
m_K * cos2Rtheta * ((R + 1.0) * cos2theta - (R + 2.0));
2663 cos2phi = 2.0 * cos_varphi * cos_varphi - 1.0;
2664 cos4phi = 2.0 * cos2phi * cos2phi - 1.0;
2665 cos8phi = 2.0 * cos4phi * cos4phi - 1.0;
2668 sin4phi = 4.0 * sin_varphi * cos_varphi * cos2phi;
2670 eta[j] =
m_H0 + (1.0 /
m_g) * (Ath + Bth * cos4phi + Cth * cos8phi);
2676 (1.0 + (1.0 / 40.0) * (x0j * x0d + x1j * y0d + x2j * z0d));
2682 m_K * cosRm1theta * (R * sin_theta * sin_theta - cos2theta) *
2684 vhat = -1.0 *
m_K * R * cosRm1theta * sin_theta * sin4phi;
2686 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2687 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2688 uvec[2][j] = vhat * cos_theta;
2741 f = 2.0 *
m_Omega * sin(theta);
2743 dh = f * uphi + tan(theta) * uphi * uphi;
2768 int nq =
m_fields[0]->GetTotPoints();
2769 int ncoeffs =
m_fields[0]->GetNcoeffs();
2771 NekDouble rad_earth = 6.37122 * 1000000;
2776 std::vector<std::string> variables(nvariables);
2778 variables[0] =
"eta";
2779 variables[1] =
"hstar";
2780 variables[2] =
"vorticity";
2781 variables[3] =
"ux";
2782 variables[4] =
"uy";
2783 variables[5] =
"uz";
2784 variables[6] =
"null";
2788 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvariables);
2789 for (
int i = 0; i < nvariables; ++i)
2796 &fieldphys[0][0], 1);
2803 Vmath::Smul(nq, rad_earth, &fieldphys[1][0], 1, &fieldphys[1][0], 1);
2817 &fieldphys[k + indx][0], 1);
2819 &fieldphys[k + indx][0], 1, &fieldphys[k + indx][0], 1);
2822 for (
int i = 0; i < nvariables; ++i)
2824 m_fields[0]->FwdTrans(fieldphys[i], fieldcoeffs[i]);
2837 if (physin[0].get() == physout[0].get())
2841 for (
int i = 0; i < 3; ++i)
2852 Vmath::Vdiv(nq, tmp[1], 1, tmp[0], 1, physout[1], 1);
2855 Vmath::Vdiv(nq, tmp[2], 1, tmp[0], 1, physout[2], 1);
2863 Vmath::Vdiv(nq, physin[1], 1, physin[0], 1, physout[1], 1);
2866 Vmath::Vdiv(nq, physin[2], 1, physin[0], 1, physout[2], 1);
2877 if (physin[0].get() == physout[0].get())
2881 for (
int i = 0; i < 3; ++i)
2892 Vmath::Vmul(nq, physout[0], 1, tmp[1], 1, physout[1], 1);
2895 Vmath::Vmul(nq, physout[0], 1, tmp[2], 1, physout[2], 1);
2903 Vmath::Vmul(nq, physout[0], 1, physin[1], 1, physout[1], 1);
2906 Vmath::Vmul(nq, physout[0], 1, physin[2], 1, physout[2], 1);
2950 int nq =
m_fields[0]->GetTotPoints();
2955 NekDouble theta, phi, sin_theta, cos_theta, sin_varphi, cos_varphi;
2980 for (k = 0; k < nq; ++k)
2986 Re =
sqrt(xp * xp + yp * yp + zp * zp);
2993 theta = atan2(sin_theta, cos_theta);
2994 phi = atan2(sin_varphi, cos_varphi);
2996 cosntheta3 = cos(n * theta) * cos(n * theta) * cos(n * theta);
2998 beta_theta = -4.0 * n * cosntheta3 * cos(m * phi) * sin(n * theta) / Re;
2999 beta_phi = -m * cosntheta3 * sin(m * phi) / Re;
3001 thetax = -1.0 * cos_varphi * sin_theta;
3002 thetay = -1.0 * sin_varphi * sin_theta;
3005 phix = -1.0 * sin_varphi;
3009 uvec[0][k] = alpha * (beta_theta * thetax + beta_phi * phix);
3010 uvec[1][k] = alpha * (beta_theta * thetay + beta_phi * phiy);
3011 uvec[2][k] = alpha * (beta_theta * thetaz + beta_phi * phiz);
3013 vorticityexact[k] = -4.0 * n / Re / Re * cos_theta * cos_theta *
3014 cos_varphi * cos(m * phi) * sin(n * theta);
3020 std::cout <<
"chi migi1" << std::endl;
3031 Vmath::Vsub(nq, vorticityexact, 1, vorticitycompt, 1, vorticitycompt, 1);
3033 std::cout <<
"Vorticity: L2 error = " <<
AvgAbsInt(vorticitycompt)
3034 <<
", Linf error = " <<
Vmath::Vamax(nq, vorticitycompt, 1)
3042 boost::ignore_unused(exactsoln);
3044 int nq =
m_fields[field]->GetNpoints();
3049 if (
m_fields[field]->GetPhysState() ==
false)
3067 &exactsolution[0], 1, &exactsolution[0], 1);
3068 Vmath::Vabs(nq, exactsolution, 1, exactsolution, 1);
3070 L2error = (
m_fields[0]->Integral(exactsolution)) / L2exact;
3087 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3090 L2exact =
m_fields[1]->Integral(tmp);
3099 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3102 L2error = (
m_fields[1]->Integral(tmp)) / L2exact;
3116 if (Normalised ==
true)
3123 L2error =
sqrt(L2error * L2error / Vol);
3139 boost::ignore_unused(exactsoln);
3143 if (
m_fields[field]->GetPhysState() ==
false)
3149 int nq =
m_fields[field]->GetNpoints();
3173 1, &exactsolution[0], 1);
3195 Vmath::Vsub(nq, &exactu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3196 Vmath::Vsub(nq, &exactv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3198 Vmath::Vmul(nq, &tmpu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3199 Vmath::Vmul(nq, &tmpv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3201 Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &Lerr[0], 1);
3205 Vmath::Vmul(nq, &exactu[0], 1, &exactu[0], 1, &tmpu[0], 1);
3206 Vmath::Vmul(nq, &exactv[0], 1, &exactv[0], 1, &tmpv[0], 1);
3207 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 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)
virtual void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
NekDouble ComputeUnstableJetuphi(const NekDouble theta)
MMFSWE(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Session reader.
void AddElevationEffect(Array< OneD, Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
virtual ~MMFSWE()
Destructor.
void RossbyWave(unsigned int field, Array< OneD, NekDouble > &outfield)
void ComputeNablaCdotVelocity(Array< OneD, NekDouble > &vellc)
void GetSWEFluxVector(const int i, const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &flux)
void AddDivForGradient(Array< OneD, Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
void WeakDGSWEDirDeriv(const Array< OneD, Array< OneD, NekDouble >> &InField, Array< OneD, Array< OneD, NekDouble >> &OutField)
void EvaluateStandardCoriolis(Array< OneD, NekDouble > &outarray)
static std::string className
Name of class.
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print Summary.
void PrimitiveToConservative()
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 DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the RHS.
void Compute_demdt_cdot_ek(const int indm, const int indk, const Array< OneD, const Array< OneD, NekDouble >> &physarray, 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 Checkpoint_Output_Cartesian(std::string outname)
virtual void v_InitObject(bool DeclareFields=true)
Initialise the object.
void Computehhuhvflux(NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, NekDouble hstar, NekDouble &hflux, NekDouble &huflux, NekDouble &hvflux)
void TestVorticityComputation()
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 NumericalSWEFlux(Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd)
NekDouble ComputeMass(const Array< OneD, const NekDouble > &eta)
void EvaluateWaterDepth(void)
void UnsteadyZonalFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
virtual NekDouble v_LinfError(unsigned int field, const Array< OneD, NekDouble > &exactsoln)
Virtual function for the L_inf error computation between fields and a given exact solution.
virtual void v_DoInitialise()
Sets up initial conditions.
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)
virtual void v_DoSolve()
Solves an unsteady problem.
void EvaluateCoriolis(void)
void EvaluateCoriolisForZonalFlow(Array< OneD, NekDouble > &outarray)
virtual NekDouble v_L2Error(unsigned int field, const Array< OneD, NekDouble > &exactsoln, bool Normalised)
Virtual function for the L_2 error computation between fields and a given exact solution.
Array< OneD, Array< OneD, NekDouble > > m_velocity
Advection velocity.
virtual void v_SetInitialConditions(const NekDouble initialtime, bool dumpInitialConditions, const int domain)
void WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble >> &physarray)
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble >> &inarray, NekDouble time)
void TestSWE2Dproblem(const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield)
void AddRotation(Array< OneD, Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
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.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the projection.
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.
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.
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble >> &uvec, unsigned int field)
Array< OneD, Array< OneD, NekDouble > > m_DivMF
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFFwd
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 void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble >> &Anisotropy, const int TangentXelem=-1)
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
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
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.
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
std::vector< int > m_intVariables
int m_infosteps
Number of time steps between outputting status information.
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
The above copyright notice and this permission notice shall be included.
@ 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 plus 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)