37 #include <boost/algorithm/string.hpp>
38 #include <boost/algorithm/string/predicate.hpp>
39 #include <boost/core/ignore_unused.hpp>
71 int shapedim =
m_fields[0]->GetShapeDimension();
124 for (
int j = 0; j < shapedim; ++j)
138 if (
m_session->DefinesSolverInfo(
"TESTMAXWELLTYPE"))
140 std::string TestMaxwellTypeStr =
141 m_session->GetSolverInfo(
"TESTMAXWELLTYPE");
159 if (
m_session->DefinesSolverInfo(
"POLTYPE"))
161 std::string PolTypeStr =
m_session->GetSolverInfo(
"POLTYPE");
177 if (
m_session->DefinesSolverInfo(
"INCTYPE"))
179 std::string IncTypeStr =
m_session->GetSolverInfo(
"INCTYPE");
195 if (
m_session->DefinesSolverInfo(
"CLOAKTYPE"))
197 std::string CloakTypeStr =
m_session->GetSolverInfo(
"CLOAKTYPE");
213 if (
m_session->DefinesSolverInfo(
"SOURCETYPE"))
215 std::string SourceTypeStr =
m_session->GetSolverInfo(
"SOURCETYPE");
258 std::cout <<
"*** rad = [ " <<
Vmath::Vmax(nq, radvec, 1) <<
" , "
259 <<
Vmath::Vmin(nq, radvec, 1) <<
" ) " << std::endl;
270 std::cout <<
"*** rad = [ " <<
Vmath::Vmax(nq, radvec, 1) <<
" , "
271 <<
Vmath::Vmin(nq, radvec, 1) <<
" ) " << std::endl;
291 NekDouble eps1min, eps1max, eps2min, eps2max, eps3min, eps3max;
292 NekDouble mu1min, mu1max, mu2min, mu2max, mu3min, mu3max;
328 std::cout <<
"*** epsvec1 = [ " << eps1min <<
" , " << eps1max
329 <<
" ], epsvec2 = [ " << eps2min <<
" , " << eps2max
330 <<
" ], epsvec3 = [ " << eps3min <<
" , " << eps3max <<
" ] "
332 std::cout <<
"*** muvec1 = [ " << mu1min <<
" , " << mu1max
333 <<
" ], muvec2 = [ " << mu2min <<
" , " << mu2max
334 <<
" ], muvec3 = [ " << mu3min <<
" , " << mu3max <<
" ] "
343 dtFactor = mu1min * eps3min;
346 dtFactor = mu2min * eps3min;
353 dtFactor = eps1min * mu3min;
354 if (eps1min > eps2min)
356 dtFactor = eps2min * mu3min;
367 std::cout <<
"*** dt factor proportional to varepsilon * mu is " << dtFactor
405 std::cout <<
"*** negepsvecminus1 = [ " << eps1min <<
" , " << eps1max
406 <<
" ], negepsvecminus1 = [ " << eps2min <<
" , " << eps2max
407 <<
" ], negepsvecminus1 = [ " << eps3min <<
" , " << eps3max
408 <<
" ] " << std::endl;
409 std::cout <<
"*** negmuvecminus1 = [ " << mu1min <<
" , " << mu1max
410 <<
" ], negmuvecminus1 = [ " << mu2min <<
" , " << mu2max
411 <<
" ], negmuvecminus1 = [ " << mu3min <<
" , " << mu3max <<
" ] "
438 ASSERTL0(
false,
"Implicit unsteady Advection not set up.");
460 for (i = 0; i < nfields; ++i)
464 nvariables = nfields;
476 for (i = 0; i < nvariables; ++i)
489 "Only one of IO_CheckTime and IO_CheckSteps "
498 bool doCheckTime =
false;
522 for (i = 0; i < nq; ++i)
527 std::cout <<
"rad" <<
rad << std::endl;
549 for (i = 0; i < nq; ++i)
558 std::cout <<
"*** Area of Planar Source = "
568 int P1indx = 0, P2indx = 0, P3indx = 0;
586 for (
int i = 0; i < nq; ++i)
588 rad =
sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i]) * (y[i]));
597 for (
int i = 0; i < nq; ++i)
600 sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i] - 1.5) * (y[i] - 1.5));
608 for (
int i = 0; i < nq; ++i)
611 sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i] - 3.0) * (y[i] - 3.0));
636 std::cout <<
"Steps: " << std::setw(8) << std::left << step + 1
638 <<
"Time: " << std::setw(12) << std::left <<
m_time;
640 std::stringstream ss;
641 ss << cpuTime / 60.0 <<
" min.";
642 std::cout <<
" CPU Time: " << std::setw(8) << std::left << ss.str()
667 for (
int i = 0; i < 3; ++i)
683 for (i = 0; i < nvariables; ++i)
698 TimeSeries[indx] =
m_time;
715 <<
", Energy = " << Energy[indx] << std::endl;
725 m_fields[0]->GetBndCondExpansions()[0]->GetExpSize();
727 ->GetBndCondExpansions()[0]
741 m_fields[0]->ExtractTracePhys(fields[0], E1Fwd);
742 m_fields[0]->ExtractTracePhys(fields[1], E2Fwd);
743 m_fields[0]->ExtractTracePhys(fields[2], H3Fwd);
754 for (
int e = 0; e < totbdryexp; ++e)
756 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
757 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
768 if (E1atPEC < E1atPECloc)
770 E1atPEC = E1atPECloc;
773 if (E2atPEC < E2atPECloc)
775 E2atPEC = E2atPECloc;
778 if (H3atPEC < H3atPECloc)
780 H3atPEC = H3atPECloc;
784 std::cout <<
"At PEC, Max. E1 = " << E1atPEC
785 <<
", E2 = " << E2atPEC <<
", H3 = " << H3atPEC
791 Ezantipod[cntap++] = fields[2][indxantipod];
796 P1[cntpml] = fields[2][P1indx];
797 P2[cntpml] = fields[2][P2indx];
798 P3[cntpml] = fields[2][P3indx];
809 if (
m_session->GetComm()->GetRank() == 0)
811 std::cout <<
"Time-integration : " << intTime <<
"s" << std::endl;
813 std::cout <<
"TimeSeries = " << std::endl;
816 std::cout << TimeSeries[i] <<
", ";
818 std::cout << std::endl << std::endl;
820 std::cout <<
"Energy Density = " << std::endl;
823 std::cout << Energy[i] <<
", ";
825 std::cout << std::endl << std::endl;
834 std::cout <<
"Ez at antipod = " << std::endl;
837 std::cout << Ezantipod[i] <<
", ";
839 std::cout << std::endl << std::endl;
844 std::cout <<
"P1 = " << std::endl;
847 std::cout << P1[i] <<
", ";
849 std::cout << std::endl << std::endl;
851 std::cout <<
"P2 = " << std::endl;
854 std::cout << P2[i] <<
", ";
856 std::cout << std::endl << std::endl;
858 std::cout <<
"P3 = " << std::endl;
861 std::cout << P3[i] <<
", ";
863 std::cout << std::endl << std::endl;
867 for (i = 0; i < nvariables; ++i)
873 for (i = 0; i < nvariables; ++i)
892 int nvar = inarray.size();
898 for (i = 0; i < nvar; ++i)
903 Vmath::Vcopy(nq, &inarray[i][0], 1, &physarray[i][0], 1);
906 for (i = 0; i < nvar; i++)
923 for (
int i = 0; i < 3; ++i)
934 for (
int i = 0; i < 3; ++i)
937 Vmath::Vcopy(ncoeffs, &tmpout[i][0], 1, &modarray[i][0], 1);
951 for (i = 0; i < nvar; ++i)
953 m_fields[i]->MultiplyByElmtInvMass(modarray[i], modarray[i]);
954 m_fields[i]->BwdTrans(modarray[i], outarray[i]);
960 for (
int j = 0; j < 2; ++j)
963 Vmath::Vadd(nq, &F[0], 1, &outarray[j][0], 1, &outarray[j][0], 1);
970 AddPML(physarray, outarray);
1004 m_fields[0]->GetCoords(x0, x1, x2);
1007 NekDouble uxx, uxy, uyy, detu, uti, uri;
1027 for (
int i = 0; i < nq; i++)
1030 tmpIN[i] * (1.0 - 2 * x0[i] * x0[i]) + (1.0 - tmpIN[i]);
1033 for (
int i = 0; i < nq; ++i)
1035 theta = atan2((x1[i] + 2.0), (x0[i] + 2.0));
1040 uxx = uti * cos(theta) * cos(theta) +
1041 uri * sin(theta) * sin(theta);
1042 uyy = uti * sin(theta) * sin(theta) +
1043 uri * cos(theta) * cos(theta);
1044 uxy = (uti - uri) * cos(theta) * sin(theta);
1046 detu = uxx * uyy - uxy * uxy;
1048 tmpx = outarray[0][i] + (1.0 - uxx) * Hxdt[i] - uxy * Hydt[i];
1049 tmpy = outarray[1][i] - uxy * Hxdt[i] + (1.0 - uyy) * Hydt[i];
1051 outarray[0][i] = (1 / detu) * (uyy * tmpx - uxy * tmpy);
1052 outarray[1][i] = (1 / detu) * (-uxy * tmpx + uxx * tmpy);
1069 outarray[0], 1, outarray[0], 1);
1073 outarray[1], 1, outarray[1], 1);
1077 outarray[2], 1, outarray[2], 1);
1084 outarray[0], 1, outarray[0], 1);
1088 outarray[1], 1, outarray[1], 1);
1092 outarray[2], 1, outarray[2], 1);
1111 outarray[0], 1, outarray[0], 1);
1115 outarray[1], 1, outarray[1], 1);
1119 outarray[2], 1, outarray[2], 1);
1126 outarray[0], 1, outarray[0], 1);
1140 outarray[1], 1, outarray[1], 1);
1145 outarray[2], 1, outarray[2], 1);
1174 int ncoeffs = outarray[0].size();
1175 int nq = physarray[0].size();
1209 m_fields[var]->IProductWRTBase(tmp, tmpc);
1210 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[var], 1, outarray[var], 1);
1217 m_fields[var]->IProductWRTBase(tmp, tmpc);
1218 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[var], 1, outarray[var], 1);
1225 &tmp[0], 1, &tmp[0], 1);
1226 m_fields[var]->IProductWRTBase(tmp, tmpc);
1227 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[var], 1, outarray[var], 1);
1263 for (i = 0; i < nvar; ++i)
1265 physfield[i] = InField[i];
1269 for (i = 0; i < nvar; ++i)
1279 fluxvector[j], tmpc);
1280 Vmath::Vadd(ncoeffs, &tmpc[0], 1, &OutField[i][0], 1,
1281 &OutField[i][0], 1);
1291 for (i = 0; i < nvar; ++i)
1302 for (i = 0; i < nvar; ++i)
1305 m_fields[i]->AddFwdBwdTraceIntegral(numfluxFwd[i], numfluxBwd[i],
1322 boost::ignore_unused(time);
1324 int var = inarray.size();
1329 for (
int i = 0; i < var; ++i)
1336 bool dumpInitialConditions,
1339 boost::ignore_unused(domain);
1375 for (
int i = 0; i < nvar; i++)
1402 for (
int i = 0; i < nvar; ++i)
1409 if (dumpInitialConditions)
1415 for (
int i = 0; i < nvar; ++i)
1417 fields[i] =
m_fields[i]->GetPhys();
1428 int nq =
m_fields[0]->GetNpoints();
1469 int nq =
m_fields[0]->GetNpoints();
1475 m_fields[0]->GetCoords(x0, x1, x2);
1493 for (
int i = 0; i < 10000; ++i)
1497 (1.0 / cos(
m_n2 * omega) / cos(
m_n2 * omega) +
1498 1.0 / cos(
m_n1 * omega) / cos(
m_n1 * omega));
1500 newomega = omega - F / Fprime;
1502 if (fabs(newomega - omega) > Tol)
1515 std::complex<double> im =
sqrt(std::complex<double>(-1));
1516 std::complex<double> A1, A2, B1, B2;
1517 std::complex<double> Ak, Bk, nk;
1518 std::complex<double> Ec, Hc;
1521 A2 = exp(-1.0 * im * omega * (
m_n1 +
m_n2));
1522 B1 = A1 * exp(-2.0 * im *
m_n1 * omega);
1523 B2 = A2 * exp(2.0 * im *
m_n2 * omega);
1525 for (
int i = 0; i < nq; ++i)
1541 Ec = (Ak * exp(im * nk * omega * x0[i]) -
1542 Bk * exp(-im * nk * omega * x0[i])) *
1543 exp(im * omega * time);
1545 (Ak * exp(im * nk * omega * x0[i]) +
1546 Bk * exp(-im * nk * omega * x0[i])) *
1547 exp(im * omega * time);
1573 const NekDouble time,
unsigned int field,
1576 int nq =
m_fields[0]->GetNpoints();
1582 m_fields[0]->GetCoords(x0, x1, x2);
1597 for (
int i = 0; i < nq; ++i)
1599 switch (Polarization)
1603 Fx = -1.0 * (npi / omega) * sin(mpi * x0[i]) *
1604 cos(npi * x1[i]) * sin(omega * time);
1605 Fy = (mpi / omega) * cos(mpi * x0[i]) * sin(npi * x1[i]) *
1612 Fz[i] = sin(mpi * x0[i]) * sin(npi * x1[i]) * cos(omega * time);
1614 dFxdt = (npi)*sin(mpi * x0[i]) * cos(npi * x1[i]) *
1616 dFydt = -1.0 * (mpi)*cos(mpi * x0[i]) * sin(npi * x1[i]) *
1623 dFzdt[i] = omega * sin(mpi * x0[i]) * sin(npi * x1[i]) *
1630 Fx = -1.0 * (npi / omega) * cos(mpi * x0[i]) *
1631 sin(npi * x1[i]) * sin(omega * time);
1632 Fy = (mpi / omega) * sin(mpi * x0[i]) * cos(npi * x1[i]) *
1639 Fz[i] = cos(mpi * x0[i]) * cos(npi * x1[i]) * cos(omega * time);
1641 dFxdt = (npi)*cos(mpi * x0[i]) * sin(npi * x1[i]) *
1643 dFydt = -1.0 * (mpi)*sin(mpi * x0[i]) * cos(npi * x1[i]) *
1650 dFzdt[i] = omega * cos(mpi * x0[i]) * cos(npi * x1[i]) *
1704 const NekDouble time,
unsigned int field,
1707 int nq =
m_fields[0]->GetNpoints();
1713 m_fields[0]->GetCoords(x0, x1, x2);
1725 for (
int i = 0; i < nq; ++i)
1727 switch (Polarization)
1731 Fx = (npi / omega) * cos(mpi * x0[i]) * sin(npi * x1[i]) *
1733 Fy = -(mpi / omega) * sin(mpi * x0[i]) * cos(npi * x1[i]) *
1740 Fz[i] = cos(mpi * x0[i]) * cos(npi * x1[i]) * cos(omega * time);
1746 Fx = (npi / omega) * sin(mpi * x0[i]) * cos(npi * x1[i]) *
1748 Fy = -(mpi / omega) * cos(mpi * x0[i]) * sin(npi * x1[i]) *
1755 Fz[i] = sin(mpi * x0[i]) * sin(npi * x1[i]) * cos(omega * time);
1793 int nq =
m_fields[0]->GetTotPoints();
1818 NekDouble xj, yj, zj, sin_varphi, cos_varphi, sin_theta, cos_theta;
1820 for (
int i = 0; i < nq; i++)
1829 vth = -4.0 * sin_varphi * cos_varphi * cos_theta * cos_theta *
1830 cos_theta * sin_theta;
1832 -1.0 * sin_varphi * sin_varphi * cos_theta * cos_theta * cos_theta;
1833 velvec[0][i] = -vth * sin_theta * cos_varphi - vphi * sin_varphi;
1834 velvec[1][i] = -vth * sin_theta * sin_varphi + vphi * cos_varphi;
1835 velvec[2][i] = vth * cos_theta;
1837 E3[i] = (-4.0 * cos_theta * cos_theta * sin_theta * cos_varphi *
1839 (1.0 / omega * sin(omega * time));
1841 Fth = -omega * vth -
1842 (8.0 / omega) * cos_theta * sin_theta * cos_varphi * sin_varphi;
1843 Fphi = -omega * vphi +
1844 (4.0 / omega) * cos_varphi * cos_varphi * cos_theta *
1845 (2.0 * sin_theta * sin_theta - cos_theta * cos_theta);
1846 Fvec[0][i] = -Fth * sin_theta * cos_varphi - Fphi * sin_varphi;
1847 Fvec[1][i] = -Fth * sin_theta * sin_varphi + Fphi * cos_varphi;
1848 Fvec[2][i] = Fth * cos_theta;
1906 int nq =
m_fields[0]->GetTotPoints();
1915 int totnpts = totbdryexp * npts;
1934 m_fields[0]->ExtractTracePhys(x, xFwd);
1935 m_fields[0]->ExtractTracePhys(y, yFwd);
1937 for (
int i = 0; i < nTraceNumPoints; ++i)
1941 phiFwd[i] = atan2(yFwd[i] / radFwd[i], -1.0 * xFwd[i] / radFwd[i]);
1949 for (
int e = 0; e < totbdryexp; ++e)
1951 id2 =
m_fields[0]->GetTrace()->GetPhys_Offset(
1952 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
1954 Vmath::Vcopy(npts, &phiFwd[id2], 1, &Jphi[e * npts], 1);
1955 Vmath::Vcopy(npts, &radFwd[id2], 1, &Jrad[e * npts], 1);
1956 Vmath::Vcopy(npts, &ntimesHFwd[id2], 1, &Jcurrent[e * npts], 1);
1963 std::cout <<
"========================================================"
1966 std::cout <<
"phi = " << std::endl;
1967 for (
int i = 0; i < totnpts; ++i)
1969 std::cout << Jphi[i] <<
", ";
1971 std::cout << std::endl << std::endl;
1973 std::cout <<
"J = " << std::endl;
1974 for (
int i = 0; i < totnpts; ++i)
1976 std::cout << Jcurrent[i] <<
", ";
1978 std::cout << std::endl << std::endl;
1984 int nq =
m_fields[0]->GetTotPoints();
1997 for (
int i = 0; i < 2; i++)
2002 m_fields[0]->ExtractTracePhys(tmp, tmpFwd);
2004 &tmpFwd[0], 1, &outfield[0], 1, &outfield[0], 1);
2015 m_fields[0]->ExtractTracePhys(tmp, outfield);
2031 int nq =
m_fields[0]->GetNpoints();
2055 Vmath::Vsub(nq, PMLRegion, 1, RecPML, 1, CurvedPML, 1);
2063 NekDouble xRlayer, xLlayer, yRlayer, yLlayer;
2071 for (
int i = 0; i < nq; i++)
2074 if (x[i] >= xRlayer)
2076 xd = (x[i] - xRlayer) / PMLthickness;
2079 else if (x[i] <= xLlayer)
2081 xd = (xLlayer - x[i]) / PMLthickness;
2089 SigmaPML[0][i] = RecPML[i] * PMLmaxsigma * (xd * xd * xd);
2092 if (y[i] >= yRlayer)
2094 yd = (y[i] - yRlayer) / PMLthickness;
2097 else if (y[i] <= yLlayer)
2099 yd = (yLlayer - y[i]) / PMLthickness;
2107 SigmaPML[1][i] = PMLRegion[i] * PMLmaxsigma * (yd * yd * yd);
2117 for (
int i = 0; i < nq; i++)
2122 if (
rad >= PMLstart)
2124 relrad = (
rad - PMLstart) / PMLthickness;
2126 PMLRegion[i] * PMLmaxsigma * pow(relrad,
m_PMLorder);
2135 NekDouble relrad, radon, radtw, radth, radfo;
2136 for (
int i = 0; i < nq; i++)
2138 radon = -1.0 * x[i] + y[i] - 7;
2139 radtw = x[i] + y[i] - 7;
2140 radth = -x[i] - y[i] - 7;
2141 radfo = x[i] - y[i] - 7;
2145 relrad = radon / PMLthickness;
2147 PMLRegion[i] * PMLmaxsigma * pow(relrad,
m_PMLorder);
2152 relrad = radtw / PMLthickness;
2154 PMLRegion[i] * PMLmaxsigma * pow(relrad,
m_PMLorder);
2159 relrad = radth / PMLthickness;
2161 PMLRegion[i] * PMLmaxsigma * pow(relrad,
m_PMLorder);
2166 relrad = radfo / PMLthickness;
2168 PMLRegion[i] * PMLmaxsigma * pow(relrad,
m_PMLorder);
2175 std::cout <<
"*** sigma1 = [ " <<
Vmath::Vmin(nq, &SigmaPML[0][0], 1)
2177 <<
" ] , sigma2 = [ " <<
Vmath::Vmin(nq, &SigmaPML[1][0], 1)
2178 <<
" , " <<
Vmath::Vmax(nq, &SigmaPML[1][0], 1) <<
" ] "
2190 for (
int i = 0; i < 3; ++i)
2192 Vmath::Vvtvp(nq, &fields[i][0], 1, &fields[i][0], 1, &tmp[0], 1,
2196 energy = 0.5 * (
m_fields[0]->Integral(tmp));
2229 m_fields[0]->GenerateElementVector(
2256 m_fields[0]->GenerateElementVector(
2258 m_fields[0]->GenerateElementVector(
2300 boost::ignore_unused(muvec, Dispersion);
2314 NekDouble boveradel = m_b / (m_b - m_adel);
2322 if (fabs(la * lb - 1.0) < 0.00001)
2324 ExactCloakArea =
m_pi * (m_b * m_b - m_a * m_a);
2329 ExactCloakArea =
m_pi * (3.0 * lb * 3.0 * la - lb * la);
2334 ExactCloakArea = ExactCloakArea - (
m_fields[0]->Integral(Cloakregion));
2335 std::cout <<
"*** Error of Cloakregion area = " << ExactCloakArea
2343 for (
int i = 0; i < nq; ++i)
2345 if (Cloakregion[i] > 0)
2349 ratio = (radvec[i] - m_adel) / radvec[i];
2351 epsvec[0][i] = boveradel * boveradel;
2356 (1.0 - boveradel * boveradel * ratio * ratio) +
2362 epsvec[1][i] = boveradel * boveradel * (ratio * ratio);
2390 Vmath::Vsub(nq, Cloakregion, 1, Vacregion, 1, Cloakregion, 1);
2393 ExactCloakArea = ExactCloakArea - (
m_fields[0]->Integral(Cloakregion));
2394 std::cout <<
"*** Error of Cloakregion area = " << ExactCloakArea
2402 for (
int i = 0; i < nq; ++i)
2404 if (Cloakregion[i] > 0)
2410 epsvec[0][i] = radvec[i] / (radvec[i] - m_adel);
2411 epsvec[1][i] = (radvec[i] - m_adel) / radvec[i];
2412 muvec[2][i] = (m_b / (m_b - m_adel)) * (m_b / (m_b - m_adel)) *
2413 (radvec[i] - m_adel) / radvec[i];
2415 muvec[0][i] = epsvec[0][i];
2416 muvec[1][i] = epsvec[1][i];
2417 epsvec[2][i] = muvec[2][i];
2426 int nq =
m_fields[0]->GetTotPoints();
2434 &Sigma1minus0[0], 1);
2436 &Sigma0minus1[0], 1);
2438 &Sigma0plus1Neg[0], 1);
2453 Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxH0][0], 1,
2454 &outarray[indxH0][0], 1, &outarray[indxH0][0], 1);
2455 Vmath::Vadd(nq, &physarray[indxQ0][0], 1, &outarray[indxH0][0], 1,
2456 &outarray[indxH0][0], 1);
2459 Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxH1][0], 1,
2460 &outarray[indxH1][0], 1, &outarray[indxH1][0], 1);
2461 Vmath::Vadd(nq, &physarray[indxQ1][0], 1, &outarray[indxH1][0], 1,
2462 &outarray[indxH1][0], 1);
2465 Vmath::Vvtvp(nq, &Sigma0plus1Neg[0], 1, &physarray[indxE2][0], 1,
2466 &outarray[indxE2][0], 1, &outarray[indxE2][0], 1);
2467 Vmath::Vadd(nq, &physarray[indxP2][0], 1, &outarray[indxE2][0], 1,
2468 &outarray[indxE2][0], 1);
2471 Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxH0][0], 1,
2472 &physarray[indxQ0][0], 1, &outarray[indxQ0][0], 1);
2474 &outarray[indxQ0][0], 1);
2478 Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxH1][0], 1,
2479 &physarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2481 &outarray[indxQ1][0], 1);
2487 &outarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2492 &outarray[indxP2][0], 1);
2493 Vmath::Vmul(nq, &physarray[indxE2][0], 1, &outarray[indxP2][0], 1,
2494 &outarray[indxP2][0], 1);
2509 Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxE0][0], 1,
2510 &outarray[indxE0][0], 1, &outarray[indxE0][0], 1);
2511 Vmath::Vsub(nq, &outarray[indxE0][0], 1, &physarray[indxQ0][0], 1,
2512 &outarray[indxE0][0], 1);
2515 Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxE1][0], 1,
2516 &outarray[indxE1][0], 1, &outarray[indxE1][0], 1);
2517 Vmath::Vsub(nq, &outarray[indxE1][0], 1, &physarray[indxQ1][0], 1,
2518 &outarray[indxE1][0], 1);
2521 Vmath::Vvtvp(nq, &Sigma0plus1Neg[0], 1, &physarray[indxH2][0], 1,
2522 &outarray[indxH2][0], 1, &outarray[indxH2][0], 1);
2523 Vmath::Vsub(nq, &outarray[indxH2][0], 1, &physarray[indxP2][0], 1,
2524 &outarray[indxH2][0], 1);
2527 Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxE0][0], 1,
2528 &physarray[indxQ0][0], 1, &outarray[indxQ0][0], 1);
2530 &outarray[indxQ0][0], 1);
2534 Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxE1][0], 1,
2535 &physarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2537 &outarray[indxQ1][0], 1);
2543 &outarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2548 &outarray[indxP2][0], 1);
2549 Vmath::Vmul(nq, &physarray[indxH2][0], 1, &outarray[indxP2][0], 1,
2550 &outarray[indxP2][0], 1);
2564 int nq =
m_fields[0]->GetTotPoints();
2565 int ncoeffs =
m_fields[0]->GetNcoeffs();
2567 std::string outname =
2568 m_sessionName +
"Tot_" + boost::lexical_cast<std::string>(n) +
".chk";
2570 std::vector<std::string> variables(nvar);
2571 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2573 for (
int i = 0; i < nvar; ++i)
2579 for (
int i = 0; i < nvar; ++i)
2582 Vmath::Vadd(nq, fieldphys[i], 1, totfield, 1, totfield, 1);
2584 m_fields[i]->FwdTrans(totfield, fieldcoeffs[i]);
2596 int nq =
m_fields[0]->GetTotPoints();
2597 int ncoeffs =
m_fields[0]->GetNcoeffs();
2599 std::string outname =
2600 m_sessionName +
"Plot_" + boost::lexical_cast<std::string>(n) +
".chk";
2602 std::vector<std::string> variables(nvar);
2603 variables[0] =
"Fx";
2604 variables[1] =
"Fy";
2605 variables[2] =
"Fz";
2607 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2608 for (
int i = 0; i < nvar; ++i)
2619 &tmp[0], 1, &tmp[0], 1);
2621 m_fields[k]->FwdTrans(tmp, fieldcoeffs[k]);
2632 int nq =
m_fields[0]->GetTotPoints();
2633 int ncoeffs =
m_fields[0]->GetNcoeffs();
2636 boost::lexical_cast<std::string>(n) +
".chk";
2638 std::vector<std::string> variables(nvar);
2639 variables[0] =
"Frx";
2640 variables[1] =
"Fry";
2641 variables[2] =
"Frz";
2642 for (
int i = 3; i < nvar; ++i)
2647 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2648 for (
int i = 0; i < nvar; ++i)
2658 Vmath::Vadd(nq, fieldphys[0], 1, totfield0, 1, totfield0, 1);
2661 Vmath::Vadd(nq, fieldphys[1], 1, totfield1, 1, totfield1, 1);
2668 &tmp[0], 1, &tmp[0], 1);
2670 m_fields[k]->FwdTrans(tmp, fieldcoeffs[k]);
2673 for (
int j = 3; j < nvar; ++j)
2675 m_fields[j]->FwdTrans(fieldphys[j], fieldcoeffs[j]);
2685 boost::ignore_unused(time);
2688 int nq =
m_fields[0]->GetTotPoints();
2689 int ncoeffs =
m_fields[0]->GetNcoeffs();
2692 boost::lexical_cast<std::string>(n) +
".chk";
2694 std::vector<std::string> variables(nvar);
2695 variables[0] =
"EDFx";
2696 variables[1] =
"EDFy";
2697 variables[2] =
"EDFz";
2698 for (
int i = 3; i < nvar; ++i)
2703 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2704 for (
int i = 0; i < nvar; ++i)
2718 &tmp[0], 1, &tmp[0], 1);
2720 Vmath::Vmul(nq, &fieldphys[2][0], 1, &tmp[0], 1, &tmp[0], 1);
2727 m_fields[k]->FwdTrans(tmp, fieldcoeffs[k]);
2737 boost::ignore_unused(time);
2740 int nq =
m_fields[0]->GetTotPoints();
2741 int ncoeffs =
m_fields[0]->GetNcoeffs();
2744 boost::lexical_cast<std::string>(n) +
".chk";
2746 std::vector<std::string> variables(nvar);
2747 variables[0] =
"Energy";
2748 variables[1] =
"EnergyFlux";
2749 variables[2] =
"Zero";
2750 for (
int i = 3; i < nvar; ++i)
2755 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2756 for (
int i = 0; i < nvar; ++i)
2770 Vmath::Vvtvp(nq, &fieldphys[k][0], 1, &fieldphys[k][0], 1, &energy[0],
2774 m_fields[0]->FwdTrans(energy, fieldcoeffs[0]);
2779 for (
int k = 0; k < 2; ++k)
2785 Vmath::Vvtvp(nq, &fieldphys[k][0], 1, &fieldphys[k][0], 1,
2786 &energyflux[0], 1, &energyflux[0], 1);
2790 Vmath::Vmul(nq, &fieldphys[2][0], 1, &energyflux[0], 1, &energyflux[0], 1);
2792 m_fields[1]->FwdTrans(energyflux, fieldcoeffs[1]);
2804 int nq =
m_fields[0]->GetTotPoints();
2805 int ncoeffs =
m_fields[0]->GetNcoeffs();
2819 1.0 / (1.0 + exp(-0.5 * (time -
m_PSduration) / SFradius));
2821 for (
int j = 0; j < nq; ++j)
2823 rad =
sqrt((x[j] - Psx) * (x[j] - Psx) + (y[j] - Psy) * (y[j] - Psy) +
2824 (z[j] - Psz) * (z[j] - Psz));
2825 outarray[j] = SmoothFactor * exp(-1.0 * (
rad / Gaussianradius) *
2826 (
rad / Gaussianradius));
2829 m_fields[0]->FwdTransLocalElmt(outarray, tmpc);
2830 m_fields[0]->BwdTrans(tmpc, outarray);
2842 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2851 for (
int j = 0; j < nq; ++j)
2860 outarray[j] = 2.0 * m_Omega * sin_theta;
2869 int nq = physarray[0].size();
2911 Vmath::Vadd(nq, tmp, 1, outarray[j], 1, outarray[j], 1);
2926 if (CloakNlayer == 8)
2929 if (fabs(la * lb - 1.0) < 0.0001)
2955 if (CloakNlayer == 16)
2958 if (fabs(la * lb - 1.0) < 0.0001)
3004 m_fields[0]->GetCoords(x0, x1, x2);
3010 m_fields[0]->GenerateElementVector(Layer[0], 1.0, 0.0, Currenttmp);
3012 for (
int i = 1; i < CloakNlayer; ++i)
3014 m_fields[0]->GenerateElementVector(Layer[i], 1.0, 0.0, Currenttmp);
3015 m_fields[0]->GenerateElementVector(Layer[i - 1], 1.0, 0.0, Formertmp);
3017 Vmath::Vsub(nq, Currenttmp, 1, Formertmp, 1, Currenttmp, 1);
3019 Vmath::Svtvp(nq, 1.0 * (i + 1), Currenttmp, 1, Cloakregion, 1,
3025 for (
int i = 0; i < nq; ++i)
3031 if ((Cloakregion[i] > 0) && (CloakNlayer > 0))
3033 radvec[i] = 1.0 + (m_b - m_a) / CloakNlayer *
3034 (Cloakregion[i] - 0.5);
3040 sqrt(x0[i] * x0[i] / la / la + x1[i] * x1[i] / lb / lb);
3047 radvec[i] =
sqrt(2.0 * x0[i] * x0[i] +
3048 x1[i] * x1[i] * x1[i] * x1[i] + x1[i] * x1[i]);
3054 radvec[i] =
sqrt(3.0 * x0[i] * x0[i] +
3055 x1[i] * x1[i] * x1[i] * x1[i] - x1[i] * x1[i]);
3071 s,
"TestMaxwellType",
3137 int Ntot = inarray.size();
3140 for (
int i = 0; i < Ntot; ++i)
3142 std::cout <<
"[" << i <<
"] = " << inarray[2][i] << std::endl;
3145 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.
Array< OneD, NekDouble > m_mu
void ComputeMaterialVector(Array< OneD, Array< OneD, NekDouble >> &epsvec, Array< OneD, Array< OneD, NekDouble >> &muvec)
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 Checkpoint_TotalFieldOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
void Checkpoint_EnergyOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
void Checkpoint_TotPlotOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
void GenerateSigmaPML(const NekDouble PMLthickness, const NekDouble PMLstart, const NekDouble PMLmaxsigma, Array< OneD, Array< OneD, NekDouble >> &SigmaPML)
void Checkpoint_PlotOutput(const int n, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
Array< OneD, NekDouble > ComputeRadCloak(const int Nlayer=5)
Array< OneD, NekDouble > TestMaxwell2DPEC(const NekDouble time, unsigned int field, const SolverUtils::PolType Polarization)
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the projection.
Array< OneD, NekDouble > m_SourceVector
virtual ~MMFMaxwell()
Destructor.
void AddPML(const Array< OneD, const Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
NekDouble m_Cloakraddelta
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the RHS.
static std::string className
Name of class.
virtual void v_SetInitialConditions(const NekDouble initialtime, bool dumpInitialConditions, const int domain)
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 .
void ComputeMaterialMicroWaveCloak(const Array< OneD, const NekDouble > &radvec, Array< OneD, Array< OneD, NekDouble >> &epsvec, Array< OneD, Array< OneD, NekDouble >> &muvec)
Array< OneD, Array< OneD, NekDouble > > m_SigmaPML
Array< OneD, NekDouble > m_varepsilon
void AddGreenDerivCompensate(const Array< OneD, const Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
MMFMaxwell(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Session reader.
void Printout_SurfaceCurrent(Array< OneD, Array< OneD, NekDouble >> &fields, const int time)
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
int m_PrintoutSurfaceCurrent
Array< OneD, NekDouble > GaussianPulse(const NekDouble time, const NekDouble Psx, const NekDouble Psy, const NekDouble Psz, const NekDouble Gaussianradius)
void AddCoriolis(Array< OneD, Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
virtual void v_DoSolve()
Solves an unsteady problem.
void print_MMF(Array< OneD, Array< OneD, NekDouble >> &inarray)
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print Summary.
Array< OneD, NekDouble > m_wp2
void Checkpoint_EDFluxOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble >> &fieldphys)
virtual void v_InitObject(bool DeclareFields=true)
Initialise the object.
virtual void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
void ComputeMaterialOpticalCloak(const Array< OneD, const NekDouble > &radvec, Array< OneD, Array< OneD, NekDouble >> &epsvec, Array< OneD, Array< OneD, NekDouble >> &muvec, const bool Dispersion=false)
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)
NekDouble ComputeEnergyDensity(Array< OneD, Array< OneD, NekDouble >> &fields)
Array< OneD, NekDouble > ComputeSurfaceCurrent(const int time, const Array< OneD, const Array< OneD, NekDouble >> &fields)
int m_spacedim
Spatial dimension (>= expansion dim).
SOLVER_UTILS_EXPORT int GetTraceNpoints()
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.
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.
std::string m_sessionName
Name of the session.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
int m_steps
Number of steps to take.
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
int m_checksteps
Number of steps between checkpoints.
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)
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 ComputeZimYim(Array< OneD, Array< OneD, NekDouble >> &epsvec, Array< OneD, Array< OneD, NekDouble >> &muvec)
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 DeriveCrossProductMF(Array< OneD, Array< OneD, NekDouble >> &CrossProductMF)
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
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble >> &Anisotropy, const int TangentXelem=-1)
Array< OneD, Array< OneD, NekDouble > > m_negepsvecminus1
SOLVER_UTILS_EXPORT void AdddedtMaxwell(const Array< OneD, const Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
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)
TestMaxwellType m_TestMaxwellType
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
SOLVER_UTILS_EXPORT void Computedemdxicdote()
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.
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 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
The above copyright notice and this permission notice shall be included.
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 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.
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)