37#include <boost/algorithm/string.hpp>
38#include <boost/algorithm/string/predicate.hpp>
71 int shapedim =
m_fields[0]->GetShapeDimension();
72 ASSERTL0(shapedim <= 2,
"This solver is only for 1D lines or 2D surfaces");
125 for (
int j = 0; j < shapedim; ++j)
139 if (
m_session->DefinesSolverInfo(
"TESTMAXWELLTYPE"))
141 std::string TestMaxwellTypeStr =
142 m_session->GetSolverInfo(
"TESTMAXWELLTYPE");
160 if (
m_session->DefinesSolverInfo(
"POLTYPE"))
162 std::string PolTypeStr =
m_session->GetSolverInfo(
"POLTYPE");
178 if (
m_session->DefinesSolverInfo(
"INCTYPE"))
180 std::string IncTypeStr =
m_session->GetSolverInfo(
"INCTYPE");
196 if (
m_session->DefinesSolverInfo(
"CLOAKTYPE"))
198 std::string CloakTypeStr =
m_session->GetSolverInfo(
"CLOAKTYPE");
214 if (
m_session->DefinesSolverInfo(
"SOURCETYPE"))
216 std::string SourceTypeStr =
m_session->GetSolverInfo(
"SOURCETYPE");
259 std::cout <<
"*** rad = [ " <<
Vmath::Vmax(nq, radvec, 1) <<
" , "
260 <<
Vmath::Vmin(nq, radvec, 1) <<
" ) " << std::endl;
271 std::cout <<
"*** rad = [ " <<
Vmath::Vmax(nq, radvec, 1) <<
" , "
272 <<
Vmath::Vmin(nq, radvec, 1) <<
" ) " << std::endl;
292 NekDouble eps1min, eps1max, eps2min, eps2max, eps3min, eps3max;
293 NekDouble mu1min, mu1max, mu2min, mu2max, mu3min, mu3max;
329 std::cout <<
"*** epsvec1 = [ " << eps1min <<
" , " << eps1max
330 <<
" ], epsvec2 = [ " << eps2min <<
" , " << eps2max
331 <<
" ], epsvec3 = [ " << eps3min <<
" , " << eps3max <<
" ] "
333 std::cout <<
"*** muvec1 = [ " << mu1min <<
" , " << mu1max
334 <<
" ], muvec2 = [ " << mu2min <<
" , " << mu2max
335 <<
" ], muvec3 = [ " << mu3min <<
" , " << mu3max <<
" ] "
344 dtFactor = mu1min * eps3min;
347 dtFactor = mu2min * eps3min;
354 dtFactor = eps1min * mu3min;
355 if (eps1min > eps2min)
357 dtFactor = eps2min * mu3min;
368 std::cout <<
"*** dt factor proportional to varepsilon * mu is " << dtFactor
406 std::cout <<
"*** negepsvecminus1 = [ " << eps1min <<
" , " << eps1max
407 <<
" ], negepsvecminus1 = [ " << eps2min <<
" , " << eps2max
408 <<
" ], negepsvecminus1 = [ " << eps3min <<
" , " << eps3max
409 <<
" ] " << std::endl;
410 std::cout <<
"*** negmuvecminus1 = [ " << mu1min <<
" , " << mu1max
411 <<
" ], negmuvecminus1 = [ " << mu2min <<
" , " << mu2max
412 <<
" ], negmuvecminus1 = [ " << mu3min <<
" , " << mu3max <<
" ] "
439 ASSERTL0(
false,
"Implicit unsteady Advection not set up.");
454 for (i = 0; i < nfields; ++i)
458 nvariables = nfields;
470 for (i = 0; i < nvariables; ++i)
483 "Only one of IO_CheckTime and IO_CheckSteps "
492 bool doCheckTime =
false;
517 for (i = 0; i < nq; ++i)
522 std::cout <<
"rad" <<
rad << std::endl;
544 for (i = 0; i < nq; ++i)
553 std::cout <<
"*** Area of Planar Source = "
563 int P1indx = 0, P2indx = 0, P3indx = 0;
581 for (
int i = 0; i < nq; ++i)
583 rad =
sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i]) * (y[i]));
592 for (
int i = 0; i < nq; ++i)
595 sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i] - 1.5) * (y[i] - 1.5));
603 for (
int i = 0; i < nq; ++i)
606 sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i] - 3.0) * (y[i] - 3.0));
632 std::cout <<
"Steps: " << std::setw(8) << std::left << step + 1
634 <<
"Time: " << std::setw(12) << std::left <<
m_time;
636 std::stringstream ss;
637 ss << cpuTime / 60.0 <<
" min.";
638 std::cout <<
" CPU Time: " << std::setw(8) << std::left << ss.str()
663 for (
int i = 0; i < 3; ++i)
679 for (i = 0; i < nvariables; ++i)
694 TimeSeries[indx] =
m_time;
711 <<
", Energy = " << Energy[indx] << std::endl;
721 m_fields[0]->GetBndCondExpansions()[0]->GetExpSize();
723 ->GetBndCondExpansions()[0]
737 m_fields[0]->ExtractTracePhys(fields[0], E1Fwd);
738 m_fields[0]->ExtractTracePhys(fields[1], E2Fwd);
739 m_fields[0]->ExtractTracePhys(fields[2], H3Fwd);
750 for (
int e = 0; e < totbdryexp; ++e)
752 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
753 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
764 if (E1atPEC < E1atPECloc)
766 E1atPEC = E1atPECloc;
769 if (E2atPEC < E2atPECloc)
771 E2atPEC = E2atPECloc;
774 if (H3atPEC < H3atPECloc)
776 H3atPEC = H3atPECloc;
780 std::cout <<
"At PEC, Max. E1 = " << E1atPEC
781 <<
", E2 = " << E2atPEC <<
", H3 = " << H3atPEC
787 Ezantipod[cntap++] = fields[2][indxantipod];
792 P1[cntpml] = fields[2][P1indx];
793 P2[cntpml] = fields[2][P2indx];
794 P3[cntpml] = fields[2][P3indx];
807 std::cout <<
"Time-integration : " << intTime <<
"s" << std::endl;
809 std::cout <<
"TimeSeries = " << std::endl;
812 std::cout << TimeSeries[i] <<
", ";
814 std::cout << std::endl << std::endl;
816 std::cout <<
"Energy Density = " << std::endl;
819 std::cout << Energy[i] <<
", ";
821 std::cout << std::endl << std::endl;
830 std::cout <<
"Ez at antipod = " << std::endl;
833 std::cout << Ezantipod[i] <<
", ";
835 std::cout << std::endl << std::endl;
840 std::cout <<
"P1 = " << std::endl;
843 std::cout << P1[i] <<
", ";
845 std::cout << std::endl << std::endl;
847 std::cout <<
"P2 = " << std::endl;
850 std::cout << P2[i] <<
", ";
852 std::cout << std::endl << std::endl;
854 std::cout <<
"P3 = " << std::endl;
857 std::cout << P3[i] <<
", ";
859 std::cout << std::endl << std::endl;
863 for (i = 0; i < nvariables; ++i)
869 for (i = 0; i < nvariables; ++i)
888 int nvar = inarray.size();
894 for (i = 0; i < nvar; ++i)
899 Vmath::Vcopy(nq, &inarray[i][0], 1, &physarray[i][0], 1);
902 for (i = 0; i < nvar; i++)
919 for (
int i = 0; i < 3; ++i)
930 for (
int i = 0; i < 3; ++i)
933 Vmath::Vcopy(ncoeffs, &tmpout[i][0], 1, &modarray[i][0], 1);
947 for (i = 0; i < nvar; ++i)
949 m_fields[i]->MultiplyByElmtInvMass(modarray[i], modarray[i]);
950 m_fields[i]->BwdTrans(modarray[i], outarray[i]);
956 for (
int j = 0; j < 2; ++j)
959 Vmath::Vadd(nq, &F[0], 1, &outarray[j][0], 1, &outarray[j][0], 1);
966 AddPML(physarray, outarray);
1000 m_fields[0]->GetCoords(x0, x1, x2);
1003 NekDouble uxx, uxy, uyy, detu, uti, uri;
1023 for (
int i = 0; i < nq; i++)
1026 tmpIN[i] * (1.0 - 2 * x0[i] * x0[i]) + (1.0 - tmpIN[i]);
1029 for (
int i = 0; i < nq; ++i)
1031 theta = atan2((x1[i] + 2.0), (x0[i] + 2.0));
1036 uxx = uti * cos(theta) * cos(theta) +
1037 uri * sin(theta) * sin(theta);
1038 uyy = uti * sin(theta) * sin(theta) +
1039 uri * cos(theta) * cos(theta);
1040 uxy = (uti - uri) * cos(theta) * sin(theta);
1042 detu = uxx * uyy - uxy * uxy;
1044 tmpx = outarray[0][i] + (1.0 - uxx) * Hxdt[i] - uxy * Hydt[i];
1045 tmpy = outarray[1][i] - uxy * Hxdt[i] + (1.0 - uyy) * Hydt[i];
1047 outarray[0][i] = (1 / detu) * (uyy * tmpx - uxy * tmpy);
1048 outarray[1][i] = (1 / detu) * (-uxy * tmpx + uxx * tmpy);
1065 outarray[0], 1, outarray[0], 1);
1069 outarray[1], 1, outarray[1], 1);
1073 outarray[2], 1, outarray[2], 1);
1080 outarray[0], 1, outarray[0], 1);
1084 outarray[1], 1, outarray[1], 1);
1088 outarray[2], 1, outarray[2], 1);
1107 outarray[0], 1, outarray[0], 1);
1111 outarray[1], 1, outarray[1], 1);
1115 outarray[2], 1, outarray[2], 1);
1122 outarray[0], 1, outarray[0], 1);
1136 outarray[1], 1, outarray[1], 1);
1141 outarray[2], 1, outarray[2], 1);
1170 int ncoeffs = outarray[0].size();
1171 int nq = physarray[0].size();
1205 m_fields[var]->IProductWRTBase(tmp, tmpc);
1206 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[var], 1, outarray[var], 1);
1213 m_fields[var]->IProductWRTBase(tmp, tmpc);
1214 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[var], 1, outarray[var], 1);
1221 &tmp[0], 1, &tmp[0], 1);
1222 m_fields[var]->IProductWRTBase(tmp, tmpc);
1223 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[var], 1, outarray[var], 1);
1259 for (i = 0; i < nvar; ++i)
1261 physfield[i] = InField[i];
1265 for (i = 0; i < nvar; ++i)
1275 fluxvector[j], tmpc);
1276 Vmath::Vadd(ncoeffs, &tmpc[0], 1, &OutField[i][0], 1,
1277 &OutField[i][0], 1);
1287 for (i = 0; i < nvar; ++i)
1298 for (i = 0; i < nvar; ++i)
1301 m_fields[i]->AddFwdBwdTraceIntegral(numfluxFwd[i], numfluxBwd[i],
1319 int var = inarray.size();
1321 if (inarray != outarray)
1324 for (
int i = 0; i < var; ++i)
1332 bool dumpInitialConditions,
1333 [[maybe_unused]]
const int domain)
1369 for (
int i = 0; i < nvar; i++)
1396 for (
int i = 0; i < nvar; ++i)
1403 if (dumpInitialConditions)
1409 for (
int i = 0; i < nvar; ++i)
1411 fields[i] =
m_fields[i]->GetPhys();
1422 int nq =
m_fields[0]->GetNpoints();
1463 int nq =
m_fields[0]->GetNpoints();
1469 m_fields[0]->GetCoords(x0, x1, x2);
1487 for (
int i = 0; i < 10000; ++i)
1491 (1.0 / cos(
m_n2 * omega) / cos(
m_n2 * omega) +
1492 1.0 / cos(
m_n1 * omega) / cos(
m_n1 * omega));
1494 newomega = omega - F / Fprime;
1496 if (fabs(newomega - omega) > Tol)
1509 std::complex<double> im =
sqrt(std::complex<double>(-1));
1510 std::complex<double> A1, A2, B1, B2;
1511 std::complex<double> Ak, Bk, nk;
1512 std::complex<double> Ec, Hc;
1515 A2 = exp(-1.0 * im * omega * (
m_n1 +
m_n2));
1516 B1 = A1 * exp(-2.0 * im *
m_n1 * omega);
1517 B2 = A2 * exp(2.0 * im *
m_n2 * omega);
1519 for (
int i = 0; i < nq; ++i)
1535 Ec = (Ak * exp(im * nk * omega * x0[i]) -
1536 Bk * exp(-im * nk * omega * x0[i])) *
1537 exp(im * omega * time);
1539 (Ak * exp(im * nk * omega * x0[i]) +
1540 Bk * exp(-im * nk * omega * x0[i])) *
1541 exp(im * omega * time);
1570 int nq =
m_fields[0]->GetNpoints();
1576 m_fields[0]->GetCoords(x0, x1, x2);
1591 for (
int i = 0; i < nq; ++i)
1593 switch (Polarization)
1597 Fx = -1.0 * (npi / omega) * sin(mpi * x0[i]) *
1598 cos(npi * x1[i]) * sin(omega * time);
1599 Fy = (mpi / omega) * cos(mpi * x0[i]) * sin(npi * x1[i]) *
1606 Fz[i] = sin(mpi * x0[i]) * sin(npi * x1[i]) * cos(omega * time);
1608 dFxdt = (npi)*sin(mpi * x0[i]) * cos(npi * x1[i]) *
1610 dFydt = -1.0 * (mpi)*cos(mpi * x0[i]) * sin(npi * x1[i]) *
1617 dFzdt[i] = omega * sin(mpi * x0[i]) * sin(npi * x1[i]) *
1624 Fx = -1.0 * (npi / omega) * cos(mpi * x0[i]) *
1625 sin(npi * x1[i]) * sin(omega * time);
1626 Fy = (mpi / omega) * sin(mpi * x0[i]) * cos(npi * x1[i]) *
1633 Fz[i] = cos(mpi * x0[i]) * cos(npi * x1[i]) * cos(omega * time);
1635 dFxdt = (npi)*cos(mpi * x0[i]) * sin(npi * x1[i]) *
1637 dFydt = -1.0 * (mpi)*sin(mpi * x0[i]) * cos(npi * x1[i]) *
1644 dFzdt[i] = omega * cos(mpi * x0[i]) * cos(npi * x1[i]) *
1701 int nq =
m_fields[0]->GetNpoints();
1707 m_fields[0]->GetCoords(x0, x1, x2);
1719 for (
int i = 0; i < nq; ++i)
1721 switch (Polarization)
1725 Fx = (npi / omega) * cos(mpi * x0[i]) * sin(npi * x1[i]) *
1727 Fy = -(mpi / omega) * sin(mpi * x0[i]) * cos(npi * x1[i]) *
1734 Fz[i] = cos(mpi * x0[i]) * cos(npi * x1[i]) * cos(omega * time);
1740 Fx = (npi / omega) * sin(mpi * x0[i]) * cos(npi * x1[i]) *
1742 Fy = -(mpi / omega) * cos(mpi * x0[i]) * sin(npi * x1[i]) *
1749 Fz[i] = sin(mpi * x0[i]) * sin(npi * x1[i]) * cos(omega * time);
1787 int nq =
m_fields[0]->GetTotPoints();
1812 NekDouble xj, yj, zj, sin_varphi, cos_varphi, sin_theta, cos_theta;
1814 for (
int i = 0; i < nq; i++)
1823 vth = -4.0 * sin_varphi * cos_varphi * cos_theta * cos_theta *
1824 cos_theta * sin_theta;
1826 -1.0 * sin_varphi * sin_varphi * cos_theta * cos_theta * cos_theta;
1827 velvec[0][i] = -vth * sin_theta * cos_varphi - vphi * sin_varphi;
1828 velvec[1][i] = -vth * sin_theta * sin_varphi + vphi * cos_varphi;
1829 velvec[2][i] = vth * cos_theta;
1831 E3[i] = (-4.0 * cos_theta * cos_theta * sin_theta * cos_varphi *
1833 (1.0 / omega * sin(omega * time));
1835 Fth = -omega * vth -
1836 (8.0 / omega) * cos_theta * sin_theta * cos_varphi * sin_varphi;
1837 Fphi = -omega * vphi +
1838 (4.0 / omega) * cos_varphi * cos_varphi * cos_theta *
1839 (2.0 * sin_theta * sin_theta - cos_theta * cos_theta);
1840 Fvec[0][i] = -Fth * sin_theta * cos_varphi - Fphi * sin_varphi;
1841 Fvec[1][i] = -Fth * sin_theta * sin_varphi + Fphi * cos_varphi;
1842 Fvec[2][i] = Fth * cos_theta;
1900 int nq =
m_fields[0]->GetTotPoints();
1909 int totnpts = totbdryexp * npts;
1928 m_fields[0]->ExtractTracePhys(x, xFwd);
1929 m_fields[0]->ExtractTracePhys(y, yFwd);
1931 for (
int i = 0; i < nTraceNumPoints; ++i)
1935 phiFwd[i] = atan2(yFwd[i] / radFwd[i], -1.0 * xFwd[i] / radFwd[i]);
1943 for (
int e = 0; e < totbdryexp; ++e)
1945 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
1946 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
1948 Vmath::Vcopy(npts, &phiFwd[id2], 1, &Jphi[e * npts], 1);
1949 Vmath::Vcopy(npts, &radFwd[id2], 1, &Jrad[e * npts], 1);
1950 Vmath::Vcopy(npts, &ntimesHFwd[id2], 1, &Jcurrent[e * npts], 1);
1957 std::cout <<
"========================================================"
1960 std::cout <<
"phi = " << std::endl;
1961 for (
int i = 0; i < totnpts; ++i)
1963 std::cout << Jphi[i] <<
", ";
1965 std::cout << std::endl << std::endl;
1967 std::cout <<
"J = " << std::endl;
1968 for (
int i = 0; i < totnpts; ++i)
1970 std::cout << Jcurrent[i] <<
", ";
1972 std::cout << std::endl << std::endl;
1978 int nq =
m_fields[0]->GetTotPoints();
1991 for (
int i = 0; i < 2; i++)
1996 m_fields[0]->ExtractTracePhys(tmp, tmpFwd);
1998 &tmpFwd[0], 1, &outfield[0], 1, &outfield[0], 1);
2009 m_fields[0]->ExtractTracePhys(tmp, outfield);
2025 int nq =
m_fields[0]->GetNpoints();
2049 Vmath::Vsub(nq, PMLRegion, 1, RecPML, 1, CurvedPML, 1);
2057 NekDouble xRlayer, xLlayer, yRlayer, yLlayer;
2065 for (
int i = 0; i < nq; i++)
2068 if (x[i] >= xRlayer)
2070 xd = (x[i] - xRlayer) / PMLthickness;
2073 else if (x[i] <= xLlayer)
2075 xd = (xLlayer - x[i]) / PMLthickness;
2083 SigmaPML[0][i] = RecPML[i] * PMLmaxsigma * (xd * xd * xd);
2086 if (y[i] >= yRlayer)
2088 yd = (y[i] - yRlayer) / PMLthickness;
2091 else if (y[i] <= yLlayer)
2093 yd = (yLlayer - y[i]) / PMLthickness;
2101 SigmaPML[1][i] = PMLRegion[i] * PMLmaxsigma * (yd * yd * yd);
2111 for (
int i = 0; i < nq; i++)
2116 if (
rad >= PMLstart)
2118 relrad = (
rad - PMLstart) / PMLthickness;
2120 PMLRegion[i] * PMLmaxsigma * pow(relrad,
m_PMLorder);
2129 NekDouble relrad, radon, radtw, radth, radfo;
2130 for (
int i = 0; i < nq; i++)
2132 radon = -1.0 * x[i] + y[i] - 7;
2133 radtw = x[i] + y[i] - 7;
2134 radth = -x[i] - y[i] - 7;
2135 radfo = x[i] - y[i] - 7;
2139 relrad = radon / PMLthickness;
2141 PMLRegion[i] * PMLmaxsigma * pow(relrad,
m_PMLorder);
2146 relrad = radtw / PMLthickness;
2148 PMLRegion[i] * PMLmaxsigma * pow(relrad,
m_PMLorder);
2153 relrad = radth / PMLthickness;
2155 PMLRegion[i] * PMLmaxsigma * pow(relrad,
m_PMLorder);
2160 relrad = radfo / PMLthickness;
2162 PMLRegion[i] * PMLmaxsigma * pow(relrad,
m_PMLorder);
2169 std::cout <<
"*** sigma1 = [ " <<
Vmath::Vmin(nq, &SigmaPML[0][0], 1)
2171 <<
" ] , sigma2 = [ " <<
Vmath::Vmin(nq, &SigmaPML[1][0], 1)
2172 <<
" , " <<
Vmath::Vmax(nq, &SigmaPML[1][0], 1) <<
" ] "
2184 for (
int i = 0; i < 3; ++i)
2186 Vmath::Vvtvp(nq, &fields[i][0], 1, &fields[i][0], 1, &tmp[0], 1,
2190 energy = 0.5 * (
m_fields[0]->Integral(tmp));
2223 m_fields[0]->GenerateElementVector(
2250 m_fields[0]->GenerateElementVector(
2252 m_fields[0]->GenerateElementVector(
2293 [[maybe_unused]]
const bool Dispersion)
2307 NekDouble boveradel = m_b / (m_b - m_adel);
2315 if (fabs(la * lb - 1.0) < 0.00001)
2317 ExactCloakArea =
m_pi * (m_b * m_b - m_a * m_a);
2322 ExactCloakArea =
m_pi * (3.0 * lb * 3.0 * la - lb * la);
2327 ExactCloakArea = ExactCloakArea - (
m_fields[0]->Integral(Cloakregion));
2328 std::cout <<
"*** Error of Cloakregion area = " << ExactCloakArea
2336 for (
int i = 0; i < nq; ++i)
2338 if (Cloakregion[i] > 0)
2342 ratio = (radvec[i] - m_adel) / radvec[i];
2344 epsvec[0][i] = boveradel * boveradel;
2349 (1.0 - boveradel * boveradel * ratio * ratio) +
2355 epsvec[1][i] = boveradel * boveradel * (ratio * ratio);
2383 Vmath::Vsub(nq, Cloakregion, 1, Vacregion, 1, Cloakregion, 1);
2386 ExactCloakArea = ExactCloakArea - (
m_fields[0]->Integral(Cloakregion));
2387 std::cout <<
"*** Error of Cloakregion area = " << ExactCloakArea
2395 for (
int i = 0; i < nq; ++i)
2397 if (Cloakregion[i] > 0)
2403 epsvec[0][i] = radvec[i] / (radvec[i] - m_adel);
2404 epsvec[1][i] = (radvec[i] - m_adel) / radvec[i];
2405 muvec[2][i] = (m_b / (m_b - m_adel)) * (m_b / (m_b - m_adel)) *
2406 (radvec[i] - m_adel) / radvec[i];
2408 muvec[0][i] = epsvec[0][i];
2409 muvec[1][i] = epsvec[1][i];
2410 epsvec[2][i] = muvec[2][i];
2419 int nq =
m_fields[0]->GetTotPoints();
2427 &Sigma1minus0[0], 1);
2429 &Sigma0minus1[0], 1);
2431 &Sigma0plus1Neg[0], 1);
2446 Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxH0][0], 1,
2447 &outarray[indxH0][0], 1, &outarray[indxH0][0], 1);
2448 Vmath::Vadd(nq, &physarray[indxQ0][0], 1, &outarray[indxH0][0], 1,
2449 &outarray[indxH0][0], 1);
2452 Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxH1][0], 1,
2453 &outarray[indxH1][0], 1, &outarray[indxH1][0], 1);
2454 Vmath::Vadd(nq, &physarray[indxQ1][0], 1, &outarray[indxH1][0], 1,
2455 &outarray[indxH1][0], 1);
2458 Vmath::Vvtvp(nq, &Sigma0plus1Neg[0], 1, &physarray[indxE2][0], 1,
2459 &outarray[indxE2][0], 1, &outarray[indxE2][0], 1);
2460 Vmath::Vadd(nq, &physarray[indxP2][0], 1, &outarray[indxE2][0], 1,
2461 &outarray[indxE2][0], 1);
2464 Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxH0][0], 1,
2465 &physarray[indxQ0][0], 1, &outarray[indxQ0][0], 1);
2467 &outarray[indxQ0][0], 1);
2471 Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxH1][0], 1,
2472 &physarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2474 &outarray[indxQ1][0], 1);
2480 &outarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2485 &outarray[indxP2][0], 1);
2486 Vmath::Vmul(nq, &physarray[indxE2][0], 1, &outarray[indxP2][0], 1,
2487 &outarray[indxP2][0], 1);
2502 Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxE0][0], 1,
2503 &outarray[indxE0][0], 1, &outarray[indxE0][0], 1);
2504 Vmath::Vsub(nq, &outarray[indxE0][0], 1, &physarray[indxQ0][0], 1,
2505 &outarray[indxE0][0], 1);
2508 Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxE1][0], 1,
2509 &outarray[indxE1][0], 1, &outarray[indxE1][0], 1);
2510 Vmath::Vsub(nq, &outarray[indxE1][0], 1, &physarray[indxQ1][0], 1,
2511 &outarray[indxE1][0], 1);
2514 Vmath::Vvtvp(nq, &Sigma0plus1Neg[0], 1, &physarray[indxH2][0], 1,
2515 &outarray[indxH2][0], 1, &outarray[indxH2][0], 1);
2516 Vmath::Vsub(nq, &outarray[indxH2][0], 1, &physarray[indxP2][0], 1,
2517 &outarray[indxH2][0], 1);
2520 Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxE0][0], 1,
2521 &physarray[indxQ0][0], 1, &outarray[indxQ0][0], 1);
2523 &outarray[indxQ0][0], 1);
2527 Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxE1][0], 1,
2528 &physarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2530 &outarray[indxQ1][0], 1);
2536 &outarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2541 &outarray[indxP2][0], 1);
2542 Vmath::Vmul(nq, &physarray[indxH2][0], 1, &outarray[indxP2][0], 1,
2543 &outarray[indxP2][0], 1);
2557 int nq =
m_fields[0]->GetTotPoints();
2558 int ncoeffs =
m_fields[0]->GetNcoeffs();
2560 std::string outname =
m_sessionName +
"Tot_" + std::to_string(n) +
".chk";
2562 std::vector<std::string> variables(nvar);
2563 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2565 for (
int i = 0; i < nvar; ++i)
2571 for (
int i = 0; i < nvar; ++i)
2574 Vmath::Vadd(nq, fieldphys[i], 1, totfield, 1, totfield, 1);
2576 m_fields[i]->FwdTrans(totfield, fieldcoeffs[i]);
2588 int nq =
m_fields[0]->GetTotPoints();
2589 int ncoeffs =
m_fields[0]->GetNcoeffs();
2591 std::string outname =
m_sessionName +
"Plot_" + std::to_string(n) +
".chk";
2593 std::vector<std::string> variables(nvar);
2594 variables[0] =
"Fx";
2595 variables[1] =
"Fy";
2596 variables[2] =
"Fz";
2598 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2599 for (
int i = 0; i < nvar; ++i)
2610 &tmp[0], 1, &tmp[0], 1);
2612 m_fields[k]->FwdTrans(tmp, fieldcoeffs[k]);
2623 int nq =
m_fields[0]->GetTotPoints();
2624 int ncoeffs =
m_fields[0]->GetNcoeffs();
2626 std::string outname =
2629 std::vector<std::string> variables(nvar);
2630 variables[0] =
"Frx";
2631 variables[1] =
"Fry";
2632 variables[2] =
"Frz";
2633 for (
int i = 3; i < nvar; ++i)
2638 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2639 for (
int i = 0; i < nvar; ++i)
2649 Vmath::Vadd(nq, fieldphys[0], 1, totfield0, 1, totfield0, 1);
2652 Vmath::Vadd(nq, fieldphys[1], 1, totfield1, 1, totfield1, 1);
2659 &tmp[0], 1, &tmp[0], 1);
2661 m_fields[k]->FwdTrans(tmp, fieldcoeffs[k]);
2664 for (
int j = 3; j < nvar; ++j)
2666 m_fields[j]->FwdTrans(fieldphys[j], fieldcoeffs[j]);
2673 const int n, [[maybe_unused]]
const NekDouble time,
2677 int nq =
m_fields[0]->GetTotPoints();
2678 int ncoeffs =
m_fields[0]->GetNcoeffs();
2680 std::string outname =
2683 std::vector<std::string> variables(nvar);
2684 variables[0] =
"EDFx";
2685 variables[1] =
"EDFy";
2686 variables[2] =
"EDFz";
2687 for (
int i = 3; i < nvar; ++i)
2692 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2693 for (
int i = 0; i < nvar; ++i)
2707 &tmp[0], 1, &tmp[0], 1);
2709 Vmath::Vmul(nq, &fieldphys[2][0], 1, &tmp[0], 1, &tmp[0], 1);
2716 m_fields[k]->FwdTrans(tmp, fieldcoeffs[k]);
2723 const int n, [[maybe_unused]]
const NekDouble time,
2727 int nq =
m_fields[0]->GetTotPoints();
2728 int ncoeffs =
m_fields[0]->GetNcoeffs();
2730 std::string outname =
2733 std::vector<std::string> variables(nvar);
2734 variables[0] =
"Energy";
2735 variables[1] =
"EnergyFlux";
2736 variables[2] =
"Zero";
2737 for (
int i = 3; i < nvar; ++i)
2742 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2743 for (
int i = 0; i < nvar; ++i)
2757 Vmath::Vvtvp(nq, &fieldphys[k][0], 1, &fieldphys[k][0], 1, &energy[0],
2761 m_fields[0]->FwdTrans(energy, fieldcoeffs[0]);
2766 for (
int k = 0; k < 2; ++k)
2772 Vmath::Vvtvp(nq, &fieldphys[k][0], 1, &fieldphys[k][0], 1,
2773 &energyflux[0], 1, &energyflux[0], 1);
2777 Vmath::Vmul(nq, &fieldphys[2][0], 1, &energyflux[0], 1, &energyflux[0], 1);
2779 m_fields[1]->FwdTrans(energyflux, fieldcoeffs[1]);
2791 int nq =
m_fields[0]->GetTotPoints();
2792 int ncoeffs =
m_fields[0]->GetNcoeffs();
2806 1.0 / (1.0 + exp(-0.5 * (time -
m_PSduration) / SFradius));
2808 for (
int j = 0; j < nq; ++j)
2810 rad =
sqrt((x[j] - Psx) * (x[j] - Psx) + (y[j] - Psy) * (y[j] - Psy) +
2811 (
z[j] - Psz) * (
z[j] - Psz));
2812 outarray[j] = SmoothFactor * exp(-1.0 * (
rad / Gaussianradius) *
2813 (
rad / Gaussianradius));
2816 m_fields[0]->FwdTransLocalElmt(outarray, tmpc);
2817 m_fields[0]->BwdTrans(tmpc, outarray);
2829 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2838 for (
int j = 0; j < nq; ++j)
2847 outarray[j] = 2.0 * m_Omega * sin_theta;
2856 int nq = physarray[0].size();
2898 Vmath::Vadd(nq, tmp, 1, outarray[j], 1, outarray[j], 1);
2913 if (CloakNlayer == 8)
2916 if (fabs(la * lb - 1.0) < 0.0001)
2942 if (CloakNlayer == 16)
2945 if (fabs(la * lb - 1.0) < 0.0001)
2991 m_fields[0]->GetCoords(x0, x1, x2);
2997 m_fields[0]->GenerateElementVector(Layer[0], 1.0, 0.0, Currenttmp);
2999 for (
int i = 1; i < CloakNlayer; ++i)
3001 m_fields[0]->GenerateElementVector(Layer[i], 1.0, 0.0, Currenttmp);
3002 m_fields[0]->GenerateElementVector(Layer[i - 1], 1.0, 0.0, Formertmp);
3004 Vmath::Vsub(nq, Currenttmp, 1, Formertmp, 1, Currenttmp, 1);
3006 Vmath::Svtvp(nq, 1.0 * (i + 1), Currenttmp, 1, Cloakregion, 1,
3012 for (
int i = 0; i < nq; ++i)
3018 if ((Cloakregion[i] > 0) && (CloakNlayer > 0))
3020 radvec[i] = 1.0 + (m_b - m_a) / CloakNlayer *
3021 (Cloakregion[i] - 0.5);
3027 sqrt(x0[i] * x0[i] / la / la + x1[i] * x1[i] / lb / lb);
3034 radvec[i] =
sqrt(2.0 * x0[i] * x0[i] +
3035 x1[i] * x1[i] * x1[i] * x1[i] + x1[i] * x1[i]);
3041 radvec[i] =
sqrt(3.0 * x0[i] * x0[i] +
3042 x1[i] * x1[i] * x1[i] * x1[i] - x1[i] * x1[i]);
3058 s,
"TestMaxwellType",
3124 int Ntot = inarray.size();
3127 for (
int i = 0; i < Ntot; ++i)
3129 std::cout <<
"[" << i <<
"] = " << inarray[2][i] << std::endl;
3132 reval =
sqrt(reval / Ntot);
#define ASSERTL0(condition, msg)
const char *const CloakTypeMap[]
const char *const SourceTypeMap[]
@ eOpticalDispersiveCloak
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 Checkpoint_EDFluxOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble > > &fieldphys)
Array< OneD, NekDouble > m_mu
void AddCoriolis(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
Array< OneD, NekDouble > TestMaxwellSphere(const NekDouble time, const NekDouble omega, unsigned int field)
Array< OneD, NekDouble > m_coriolis
Array< OneD, Array< OneD, NekDouble > > m_CrossProductMF
NekDouble m_Gaussianradius
void WeakDGMaxwellDirDeriv(const Array< OneD, const Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, const NekDouble time=0.0)
Calculate weak DG advection in the form .
Array< OneD, NekDouble > ComputeRadCloak(const int Nlayer=5)
void print_MMF(Array< OneD, Array< OneD, NekDouble > > &inarray)
Array< OneD, NekDouble > TestMaxwell2DPEC(const NekDouble time, unsigned int field, const SolverUtils::PolType Polarization)
Array< OneD, NekDouble > m_SourceVector
NekDouble ComputeEnergyDensity(Array< OneD, Array< OneD, NekDouble > > &fields)
NekDouble m_Cloakraddelta
static std::string className
Name of class.
void Checkpoint_EnergyOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble > > &fieldphys)
void Checkpoint_TotalFieldOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble > > &fieldphys)
void ComputeMaterialMicroWaveCloak(const Array< OneD, const NekDouble > &radvec, Array< OneD, Array< OneD, NekDouble > > &epsvec, Array< OneD, Array< OneD, NekDouble > > &muvec)
Array< OneD, NekDouble > ComputeSurfaceCurrent(const int time, const Array< OneD, const Array< OneD, NekDouble > > &fields)
Array< OneD, Array< OneD, NekDouble > > m_SigmaPML
void Checkpoint_TotPlotOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble > > &fieldphys)
Array< OneD, NekDouble > m_varepsilon
MMFMaxwell(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection.
void v_InitObject(bool DeclareFields=true) override
Initialise the object.
void ComputeMaterialOpticalCloak(const Array< OneD, const NekDouble > &radvec, Array< OneD, Array< OneD, NekDouble > > &epsvec, Array< OneD, Array< OneD, NekDouble > > &muvec, const bool Dispersion=false)
void AddGreenDerivCompensate(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void ComputeMaterialVector(Array< OneD, Array< OneD, NekDouble > > &epsvec, Array< OneD, Array< OneD, NekDouble > > &muvec)
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
int m_PrintoutSurfaceCurrent
void GenerateSigmaPML(const NekDouble PMLthickness, const NekDouble PMLstart, const NekDouble PMLmaxsigma, Array< OneD, Array< OneD, NekDouble > > &SigmaPML)
Array< OneD, NekDouble > GaussianPulse(const NekDouble time, const NekDouble Psx, const NekDouble Psy, const NekDouble Psz, const NekDouble Gaussianradius)
void Printout_SurfaceCurrent(Array< OneD, Array< OneD, NekDouble > > &fields, const int time)
void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time) override
void v_DoSolve() override
Solves an unsteady problem.
Array< OneD, NekDouble > m_wp2
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print Summary.
void AddPML(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void v_SetInitialConditions(const NekDouble initialtime, bool dumpInitialConditions, const int domain) override
Array< OneD, NekDouble > TestMaxwell2DPMC(const NekDouble time, unsigned int field, const SolverUtils::PolType Polarization)
Array< OneD, NekDouble > EvaluateCoriolis()
Array< OneD, NekDouble > TestMaxwell1D(const NekDouble time, unsigned int field)
void Checkpoint_PlotOutput(const int n, const Array< OneD, const Array< OneD, NekDouble > > &fieldphys)
int m_spacedim
Spatial dimension (>= expansion dim).
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.
std::string m_sessionName
Name of the session.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTotPoints()
int m_steps
Number of steps to take.
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
int m_checksteps
Number of steps between checkpoints.
A base class for PDEs which include an advection component.
SOLVER_UTILS_EXPORT void GetMaxwellFluxVector(const int var, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Virtual function for generating summary information.
SOLVER_UTILS_EXPORT void DeriveCrossProductMF(Array< OneD, Array< OneD, NekDouble > > &CrossProductMF)
SOLVER_UTILS_EXPORT void NumericalMaxwellFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxFwd, Array< OneD, Array< OneD, NekDouble > > &numfluxBwd, const NekDouble time=0.0)
SOLVER_UTILS_EXPORT void AdddedtMaxwell(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
SOLVER_UTILS_EXPORT void ComputeZimYim(Array< OneD, Array< OneD, NekDouble > > &epsvec, Array< OneD, Array< OneD, NekDouble > > &muvec)
Array< OneD, Array< OneD, NekDouble > > m_muvec
SpatialDomains::GeomMMF m_MMFdir
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Array< OneD, NekDouble > m_MMFfactors
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > GetIncidentField(const int var, const NekDouble time)
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFFwd
Array< OneD, Array< OneD, NekDouble > > m_negepsvecminus1
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble > > &uvec, unsigned int field)
Array< OneD, Array< OneD, NekDouble > > m_epsvec
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)
Array< OneD, Array< OneD, NekDouble > > m_negmuvecminus1
SOLVER_UTILS_EXPORT void ComputeNtimesMF()
SOLVER_UTILS_EXPORT NekDouble RootMeanSquare(const Array< OneD, const NekDouble > &inarray)
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble > > &Anisotropy, const int TangentXelem=-1)
TestMaxwellType m_TestMaxwellType
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
SOLVER_UTILS_EXPORT void Computedemdxicdote()
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
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 NekDouble rad(NekDouble x, NekDouble y)
std::shared_ptr< SessionReader > SessionReaderSharedPtr
static const NekDouble kNekZeroTol
const char *const IncTypeMap[]
std::vector< std::pair< std::string, std::string > > SummaryList
@ SIZE_TestMaxwellType
Length of enum list.
@ eTestMaxwell2DPECAVGFLUX
EquationSystemFactory & GetEquationSystemFactory()
const char *const TestMaxwellTypeMap[]
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
const char *const PolTypeMap[]
@ eTangentIrregular
Circular around the centre of domain.
@ eTangentCircular
Circular around the centre of domain.
@ eTangentNonconvex
Circular around the centre of domain.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::vector< double > z(NPUPPER)
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 Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvp (scalar times vector plus vector): z = alpha*x + y.
void Neg(int n, T *x, const int incx)
Negate x = -x.
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
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.
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.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
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 > sqrt(scalarT< T > in)