Nektar++
MMFMaxwell.cpp
Go to the documentation of this file.
1/////////////////////////////////////////////////////////////////////////////
2//
3// File: MMFMaxwell.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: MMF Maxwell solve routines
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <MMFSolver/EquationSystems/MMFMaxwell.h>
36
37#include <boost/algorithm/string.hpp>
38#include <boost/algorithm/string/predicate.hpp>
39#include <iomanip>
40#include <iostream>
41
42#include <LibUtilities/BasicUtils/Timer.h>
43#include <LibUtilities/TimeIntegration/TimeIntegrationScheme.h>
44#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
45#include <SolverUtils/MMFSystem.h>
46
47#include <typeinfo>
48
49namespace Nektar
50{
51std::string MMFMaxwell::className =
52 SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
53 "MMFMaxwell", MMFMaxwell::create, "MMFMaxwell equation.");
54
55MMFMaxwell::MMFMaxwell(const LibUtilities::SessionReaderSharedPtr &pSession,
56 const SpatialDomains::MeshGraphSharedPtr &pGraph)
57 : UnsteadySystem(pSession, pGraph), MMFSystem(pSession, pGraph)
58{
59}
60
61/**
62 * @brief Initialisation object for the unsteady linear advection equation.
63 */
64void MMFMaxwell::v_InitObject(bool DeclareFields)
65{
66 // Call to the initialisation object
67 UnsteadySystem::v_InitObject(DeclareFields);
68
69 int nq = m_fields[0]->GetNpoints();
70 int shapedim = m_fields[0]->GetShapeDimension();
71
72 m_session->LoadParameter("ElemtGroup0", m_ElemtGroup0, 0);
73 m_session->LoadParameter("ElemtGroup1", m_ElemtGroup1, 0);
74 m_session->LoadParameter("boundaryforSF", m_boundaryforSF, 0);
75 m_session->LoadParameter("PrintoutSurfaceCurrent", m_PrintoutSurfaceCurrent,
76 0);
77
78 m_session->LoadParameter("AddRotation", m_AddRotation, 0);
79
80 m_session->LoadParameter("NoInc", m_NoInc, 0);
81
82 // PML parameters
83 m_session->LoadParameter("TestPML", m_TestPML, 0);
84 m_session->LoadParameter("PMLelement", m_PMLelement, 0);
85 m_session->LoadParameter("RecPML", m_RecPML, 0);
86
87 m_session->LoadParameter("AddPML", m_AddPML, 0);
88 if (m_AddPML == 1)
89 {
90 m_RecPML = m_PMLelement;
91 }
92 m_session->LoadParameter("PMLorder", m_PMLorder, 3);
93
94 m_session->LoadParameter("PMLthickness", m_PMLthickness, 0.0);
95 m_session->LoadParameter("PMLstart", m_PMLstart, 0.0);
96 m_session->LoadParameter("PMLmaxsigma", m_PMLmaxsigma, 100.0);
97
98 // Point Source parmaters
99 m_session->LoadParameter("Psx", m_Psx, 0.0);
100 m_session->LoadParameter("Psy", m_Psy, 0.0);
101 m_session->LoadParameter("Psz", m_Psz, 0.0);
102 m_session->LoadParameter("PSduration", m_PSduration, 1.0);
103 m_session->LoadParameter("Gaussianradius", m_Gaussianradius, 1.0);
104
105 // Cloaking parameter
106 m_session->LoadParameter("CloakNlayer", m_CloakNlayer, 5);
107 m_session->LoadParameter("Cloakraddelta", m_Cloakraddelta, 0.0);
108
109 m_varepsilon = Array<OneD, NekDouble>(m_spacedim);
110 m_session->LoadParameter("varepsilon1", m_varepsilon[0], 1.0);
111 m_session->LoadParameter("varepsilon2", m_varepsilon[1], 1.0);
112 m_session->LoadParameter("varepsilon3", m_varepsilon[2], 1.0);
113 m_n1 = sqrt(m_varepsilon[0]);
114 m_n2 = sqrt(m_varepsilon[1]);
115 m_n3 = sqrt(m_varepsilon[2]);
116
117 m_mu = Array<OneD, NekDouble>(m_spacedim);
118 m_session->LoadParameter("mu1", m_mu[0], 1.0);
119 m_session->LoadParameter("mu2", m_mu[1], 1.0);
120 m_session->LoadParameter("mu3", m_mu[2], 1.0);
121
122 Array<OneD, Array<OneD, NekDouble>> Anisotropy(shapedim);
123 for (int j = 0; j < shapedim; ++j)
124 {
125 Anisotropy[j] = Array<OneD, NekDouble>(nq, 1.0);
126 }
127
128 // Add Rectangular PML
129 MMFSystem::MMFInitObject(Anisotropy, m_RecPML);
130
131 // Compute the cross producted MF
132 DeriveCrossProductMF(m_CrossProductMF);
133
134 m_session->LoadParameter("Frequency", m_freq, 1.0);
135
136 // Define TestMaxwellType
137 if (m_session->DefinesSolverInfo("TESTMAXWELLTYPE"))
138 {
139 std::string TestMaxwellTypeStr =
140 m_session->GetSolverInfo("TESTMAXWELLTYPE");
141 for (int i = 0; i < (int)SolverUtils::SIZE_TestMaxwellType; ++i)
142 {
143 if (boost::iequals(SolverUtils::TestMaxwellTypeMap[i],
144 TestMaxwellTypeStr))
145 {
146 m_TestMaxwellType = (SolverUtils::TestMaxwellType)i;
147 break;
148 }
149 }
150 }
151
152 else
153 {
154 m_TestMaxwellType = (SolverUtils::TestMaxwellType)0;
155 }
156
157 // Define Polarization
158 if (m_session->DefinesSolverInfo("POLTYPE"))
159 {
160 std::string PolTypeStr = m_session->GetSolverInfo("POLTYPE");
161 for (int i = 0; i < (int)SolverUtils::SIZE_PolType; ++i)
162 {
163 if (boost::iequals(SolverUtils::PolTypeMap[i], PolTypeStr))
164 {
165 m_PolType = (SolverUtils::PolType)i;
166 break;
167 }
168 }
169 }
170 else
171 {
172 m_PolType = (SolverUtils::PolType)0;
173 }
174
175 // Define Incident wave Type
176 if (m_session->DefinesSolverInfo("INCTYPE"))
177 {
178 std::string IncTypeStr = m_session->GetSolverInfo("INCTYPE");
179 for (int i = 0; i < (int)SolverUtils::SIZE_IncType; ++i)
180 {
181 if (boost::iequals(SolverUtils::IncTypeMap[i], IncTypeStr))
182 {
183 m_IncType = (SolverUtils::IncType)i;
184 break;
185 }
186 }
187 }
188 else
189 {
190 m_IncType = (SolverUtils::IncType)0;
191 }
192
193 // Define Cloak Type
194 if (m_session->DefinesSolverInfo("CLOAKTYPE"))
195 {
196 std::string CloakTypeStr = m_session->GetSolverInfo("CLOAKTYPE");
197 for (int i = 0; i < (int)SIZE_CloakType; ++i)
198 {
199 if (boost::iequals(CloakTypeMap[i], CloakTypeStr))
200 {
201 m_CloakType = (CloakType)i;
202 break;
203 }
204 }
205 }
206 else
207 {
208 m_CloakType = (CloakType)0;
209 }
210
211 // Define Source Type
212 if (m_session->DefinesSolverInfo("SOURCETYPE"))
213 {
214 std::string SourceTypeStr = m_session->GetSolverInfo("SOURCETYPE");
215 for (int i = 0; i < (int)SIZE_SourceType; ++i)
216 {
217 if (boost::iequals(SourceTypeMap[i], SourceTypeStr))
218 {
219 m_SourceType = (SourceType)i;
220 break;
221 }
222 }
223 }
224 else
225 {
226 m_SourceType = (SourceType)0;
227 }
228
229 // Compute n_timesMFFwd and m_times_timesMFFwd
230 ComputeNtimesMF();
231
232 // Compute vaepsilon and mu vector (m_epsveci, m_muvec0);
233 m_epsvec = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
234 m_muvec = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
235 for (int k = 0; k < m_spacedim; ++k)
236 {
237 m_epsvec[k] = Array<OneD, NekDouble>(nq, 1.0);
238 m_muvec[k] = Array<OneD, NekDouble>(nq, 1.0);
239 }
240
241 Array<OneD, NekDouble> radvec(nq);
242 m_DispersiveCloak = false;
243 switch (m_CloakType)
244 {
245 case eOpticalCloak:
246 {
247 radvec = ComputeRadCloak();
248 ComputeMaterialOpticalCloak(radvec, m_epsvec, m_muvec, false);
249 }
250 break;
251
252 case eOpticalConstCloak:
253 {
254 radvec = ComputeRadCloak(m_CloakNlayer);
255 ComputeMaterialOpticalCloak(radvec, m_epsvec, m_muvec, false);
256
257 std::cout << "*** rad = [ " << Vmath::Vmax(nq, radvec, 1) << " , "
258 << Vmath::Vmin(nq, radvec, 1) << " ) " << std::endl;
259 }
260 break;
261
262 case eOpticalDispersiveCloak:
263 {
264 m_DispersiveCloak = true;
265 m_wp2Tol = 0.01;
266 radvec = ComputeRadCloak();
267 ComputeMaterialOpticalCloak(radvec, m_epsvec, m_muvec, true);
268
269 std::cout << "*** rad = [ " << Vmath::Vmax(nq, radvec, 1) << " , "
270 << Vmath::Vmin(nq, radvec, 1) << " ) " << std::endl;
271 std::cout << "*** wp2 = [ " << Vmath::Vmax(nq, m_wp2, 1) << " , "
272 << Vmath::Vmin(nq, m_wp2, 1) << " ) " << std::endl;
273 }
274 break;
275
276 case eMicroWaveCloak:
277 {
278 radvec = ComputeRadCloak();
279 ComputeMaterialMicroWaveCloak(radvec, m_epsvec, m_muvec);
280 }
281 break;
282
283 default:
284 {
285 ComputeMaterialVector(m_epsvec, m_muvec);
286 }
287 break;
288 }
289
290 NekDouble eps1min, eps1max, eps2min, eps2max, eps3min, eps3max;
291 NekDouble mu1min, mu1max, mu2min, mu2max, mu3min, mu3max;
292
293 eps1min = Vmath::Vmin(nq, m_epsvec[0], 1);
294 eps3min = Vmath::Vmin(nq, m_epsvec[2], 1);
295 eps1max = Vmath::Vmax(nq, m_epsvec[0], 1);
296 eps3max = Vmath::Vmax(nq, m_epsvec[2], 1);
297
298 if (m_DispersiveCloak)
299 {
300 Array<OneD, NekDouble> realepsr(nq);
301 Vmath::Sadd(nq, -m_wp2Tol, m_wp2, 1, realepsr, 1);
302 Vmath::Smul(nq, 1.0 / (m_Incfreq * m_Incfreq), realepsr, 1, realepsr,
303 1);
304 Vmath::Neg(nq, realepsr, 1);
305 Vmath::Sadd(nq, 1.0, realepsr, 1, realepsr, 1);
306
307 eps2min = Vmath::Vmin(nq, realepsr, 1);
308 eps2max = Vmath::Vmax(nq, realepsr, 1);
309 }
310
311 else
312 {
313 eps2min = Vmath::Vmin(nq, m_epsvec[1], 1);
314 eps2max = Vmath::Vmax(nq, m_epsvec[1], 1);
315 }
316
317 mu1min = Vmath::Vmin(nq, m_muvec[0], 1);
318 mu2min = Vmath::Vmin(nq, m_muvec[1], 1);
319 mu3min = Vmath::Vmin(nq, m_muvec[2], 1);
320 mu1max = Vmath::Vmax(nq, m_muvec[0], 1);
321 mu2max = Vmath::Vmax(nq, m_muvec[1], 1);
322 mu3max = Vmath::Vmax(nq, m_muvec[2], 1);
323
324 std::cout << "muvec0 = " << RootMeanSquare(m_muvec[0])
325 << ", muvec1 = " << RootMeanSquare(m_muvec[1]) << std::endl;
326
327 std::cout << "*** epsvec1 = [ " << eps1min << " , " << eps1max
328 << " ], epsvec2 = [ " << eps2min << " , " << eps2max
329 << " ], epsvec3 = [ " << eps3min << " , " << eps3max << " ] "
330 << std::endl;
331 std::cout << "*** muvec1 = [ " << mu1min << " , " << mu1max
332 << " ], muvec2 = [ " << mu2min << " , " << mu2max
333 << " ], muvec3 = [ " << mu3min << " , " << mu3max << " ] "
334 << std::endl;
335
336 NekDouble dtFactor;
337 switch (m_PolType)
338 {
339 // eTransMagnetic
340 case SolverUtils::eTransMagnetic:
341 {
342 dtFactor = mu1min * eps3min;
343 if (mu1min > mu2min)
344 {
345 dtFactor = mu2min * eps3min;
346 }
347 }
348 break;
349
350 case SolverUtils::eTransElectric:
351 {
352 dtFactor = eps1min * mu3min;
353 if (eps1min > eps2min)
354 {
355 dtFactor = eps2min * mu3min;
356 }
357 }
358 break;
359
360 default:
361 {
362 dtFactor = 1.0;
363 }
364 break;
365 }
366 std::cout << "*** dt factor proportional to varepsilon * mu is " << dtFactor
367 << std::endl;
368
369 // Compute m_Zim and m_Yim
370 ComputeZimYim(m_epsvec, m_muvec);
371
372 // Compute m_epsvecminus1 and m_muminus1
373 m_negepsvecminus1 = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
374 m_negmuvecminus1 = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
375 for (int k = 0; k < m_spacedim; ++k)
376 {
377 m_negepsvecminus1[k] = Array<OneD, NekDouble>(nq, 0.0);
378 m_negmuvecminus1[k] = Array<OneD, NekDouble>(nq, 0.0);
379
380 if (!m_NoInc)
381 {
382 Vmath::Sadd(nq, -1.0, m_muvec[k], 1, m_negmuvecminus1[k], 1);
383 Vmath::Sadd(nq, -1.0, m_epsvec[k], 1, m_negepsvecminus1[k], 1);
384
385 Vmath::Neg(nq, m_negmuvecminus1[k], 1);
386 Vmath::Neg(nq, m_negepsvecminus1[k], 1);
387 }
388 }
389
390 eps1min = Vmath::Vmin(nq, m_negepsvecminus1[0], 1);
391 eps2min = Vmath::Vmin(nq, m_negepsvecminus1[1], 1);
392 eps3min = Vmath::Vmin(nq, m_negepsvecminus1[2], 1);
393 eps1max = Vmath::Vmax(nq, m_negepsvecminus1[0], 1);
394 eps2max = Vmath::Vmax(nq, m_negepsvecminus1[1], 1);
395 eps3max = Vmath::Vmax(nq, m_negepsvecminus1[2], 1);
396
397 mu1min = Vmath::Vmin(nq, m_negmuvecminus1[0], 1);
398 mu2min = Vmath::Vmin(nq, m_negmuvecminus1[1], 1);
399 mu3min = Vmath::Vmin(nq, m_negmuvecminus1[2], 1);
400 mu1max = Vmath::Vmax(nq, m_negmuvecminus1[0], 1);
401 mu2max = Vmath::Vmax(nq, m_negmuvecminus1[1], 1);
402 mu3max = Vmath::Vmax(nq, m_negmuvecminus1[2], 1);
403
404 std::cout << "*** negepsvecminus1 = [ " << eps1min << " , " << eps1max
405 << " ], negepsvecminus1 = [ " << eps2min << " , " << eps2max
406 << " ], negepsvecminus1 = [ " << eps3min << " , " << eps3max
407 << " ] " << std::endl;
408 std::cout << "*** negmuvecminus1 = [ " << mu1min << " , " << mu1max
409 << " ], negmuvecminus1 = [ " << mu2min << " , " << mu2max
410 << " ], negmuvecminus1 = [ " << mu3min << " , " << mu3max << " ] "
411 << std::endl;
412
413 // Compute de^m/dt \cdot e^k
414 if (m_AddRotation)
415 {
416 m_coriolis = Array<OneD, NekDouble>(nq);
417 m_coriolis = EvaluateCoriolis();
418
419 Computedemdxicdote();
420 }
421
422 // Generate Sigma Block with thicknes of m_PMLthickness and m_PMLmax
423 if (m_AddPML > 0)
424 {
425 GenerateSigmaPML(m_PMLthickness, m_PMLstart, m_PMLmaxsigma, m_SigmaPML);
426 }
427
428 // If explicit it computes RHS and PROJECTION for the time integration
429 if (m_explicitAdvection)
430 {
431 m_ode.DefineOdeRhs(&MMFMaxwell::DoOdeRhs, this);
432 m_ode.DefineProjection(&MMFMaxwell::DoOdeProjection, this);
433 }
434 // Otherwise it gives an error (no implicit integration)
435 else
436 {
437 ASSERTL0(false, "Implicit unsteady Advection not set up.");
438 }
439}
440
441/**
442 * @brief Unsteady linear advection equation destructor.
443 */
444MMFMaxwell::~MMFMaxwell()
445{
446}
447
448void MMFMaxwell::v_DoSolve()
449{
450 ASSERTL0(m_intScheme != nullptr, "No time integration scheme.");
451
452 int i, nchk = 1;
453 int nq = GetTotPoints();
454 int nvariables = 0;
455 int nfields = m_fields.size();
456
457 if (m_intVariables.empty())
458 {
459 for (i = 0; i < nfields; ++i)
460 {
461 m_intVariables.push_back(i);
462 }
463 nvariables = nfields;
464 }
465 else
466 {
467 nvariables = m_intVariables.size();
468 }
469
470 // Set up wrapper to fields data storage.
471 Array<OneD, Array<OneD, NekDouble>> fields(nvariables);
472 Array<OneD, Array<OneD, NekDouble>> tmp(nvariables);
473
474 // Order storage to list time-integrated fields first.
475 for (i = 0; i < nvariables; ++i)
476 {
477 fields[i] = m_fields[m_intVariables[i]]->GetPhys();
478 m_fields[m_intVariables[i]]->SetPhysState(false);
479 }
480
481 // Initialise time integration scheme
482 m_intScheme->InitializeScheme(m_timestep, fields, m_time, m_ode);
483
484 // Check uniqueness of checkpoint output
485 ASSERTL0((m_checktime == 0.0 && m_checksteps == 0) ||
486 (m_checktime > 0.0 && m_checksteps == 0) ||
487 (m_checktime == 0.0 && m_checksteps > 0),
488 "Only one of IO_CheckTime and IO_CheckSteps "
489 "should be set!");
490
491 int Ntot = m_checksteps ? m_steps / m_checksteps + 1 : 0;
492
493 Array<OneD, NekDouble> TimeSeries(Ntot ? Ntot : 1);
494 Array<OneD, NekDouble> Energy(Ntot ? Ntot : 1);
495
496 LibUtilities::Timer timer;
497 bool doCheckTime = false;
498 int step = 0;
499 NekDouble intTime = 0.0;
500 NekDouble cpuTime = 0.0;
501 NekDouble elapsed = 0.0;
502
503 int cntap = 0;
504 Array<OneD, NekDouble> Ezantipod;
505 int indxantipod = 0;
506
507 switch (m_SourceType)
508 {
509 case ePointSource:
510 {
511 Ezantipod = Array<OneD, NekDouble>(
512 m_checksteps ? m_steps / m_checksteps : 1);
513
514 Array<OneD, NekDouble> x(nq);
515 Array<OneD, NekDouble> y(nq);
516 Array<OneD, NekDouble> z(nq);
517
518 m_fields[0]->GetCoords(x, y, z);
519
520 NekDouble Tol = 0.000001;
521 NekDouble rad;
522 for (i = 0; i < nq; ++i)
523 {
524 rad = sqrt((x[i] + m_Psx) * (x[i] + m_Psx) +
525 (y[i] + m_Psy) * (y[i] + m_Psy) +
526 (z[i] + m_Psz) * (z[i] + m_Psz));
527 std::cout << "rad" << rad << std::endl;
528 if (rad < Tol)
529 {
530 indxantipod = i;
531 break;
532 }
533 }
534 }
535 break;
536
537 case ePlanarSource:
538 {
539 m_SourceVector = Array<OneD, NekDouble>(nq, 0.0);
540
541 Array<OneD, NekDouble> x(nq);
542 Array<OneD, NekDouble> y(nq);
543 Array<OneD, NekDouble> z(nq);
544
545 m_fields[0]->GetCoords(x, y, z);
546
547 NekDouble Tol = 0.000001;
548 NekDouble rad;
549 for (i = 0; i < nq; ++i)
550 {
551 rad = sqrt((x[i] - m_Psx) * (x[i] - m_Psx));
552 if (rad < Tol)
553 {
554 m_SourceVector[i] = 1.0;
555 }
556 }
557
558 std::cout << "*** Area of Planar Source = "
559 << m_fields[0]->Integral(m_SourceVector) << std::endl;
560 }
561 break;
562
563 default:
564 break;
565 }
566
567 int cntpml = 0;
568 int P1indx = 0, P2indx = 0, P3indx = 0;
569 Array<OneD, NekDouble> P1;
570 Array<OneD, NekDouble> P2;
571 Array<OneD, NekDouble> P3;
572 if (m_TestPML)
573 {
574 P1 = Array<OneD, NekDouble>(m_checksteps ? m_steps / m_checksteps : 1);
575 P2 = Array<OneD, NekDouble>(m_checksteps ? m_steps / m_checksteps : 1);
576 P3 = Array<OneD, NekDouble>(m_checksteps ? m_steps / m_checksteps : 1);
577
578 Array<OneD, NekDouble> x(nq);
579 Array<OneD, NekDouble> y(nq);
580 Array<OneD, NekDouble> z(nq);
581
582 m_fields[0]->GetCoords(x, y, z);
583
584 NekDouble Tol = 0.000001;
585 NekDouble rad;
586 for (int i = 0; i < nq; ++i)
587 {
588 rad = sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i]) * (y[i]));
589
590 if (rad < Tol)
591 {
592 P1indx = i;
593 break;
594 }
595 }
596
597 for (int i = 0; i < nq; ++i)
598 {
599 rad =
600 sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i] - 1.5) * (y[i] - 1.5));
601 if (rad < Tol)
602 {
603 P2indx = i;
604 break;
605 }
606 }
607
608 for (int i = 0; i < nq; ++i)
609 {
610 rad =
611 sqrt((x[i] + 3.0) * (x[i] + 3.0) + (y[i] - 3.0) * (y[i] - 3.0));
612 if (rad < Tol)
613 {
614 P3indx = i;
615 break;
616 }
617 }
618 }
619
620 int indx;
621 while (step < m_steps || m_time < m_fintime - NekConstants::kNekZeroTol)
622 {
623
624 timer.Start();
625 fields = m_intScheme->TimeIntegrate(step, m_timestep);
626 timer.Stop();
627
628 m_time += m_timestep;
629 elapsed = timer.TimePerTest(1);
630 intTime += elapsed;
631 cpuTime += elapsed;
632
633 // Write out status information
634 if (m_infosteps && !((step + 1) % m_infosteps) &&
635 m_session->GetComm()->GetRank() == 0)
636 {
637 std::cout << "Steps: " << std::setw(8) << std::left << step + 1
638 << " "
639 << "Time: " << std::setw(12) << std::left << m_time;
640
641 std::stringstream ss;
642 ss << cpuTime / 60.0 << " min.";
643 std::cout << " CPU Time: " << std::setw(8) << std::left << ss.str()
644 << std::endl;
645
646 cpuTime = 0.0;
647 }
648
649 switch (m_SourceType)
650 {
651 case ePointSource:
652 {
653 if (m_time <= m_PSduration)
654 {
655 Array<OneD, NekDouble> Impulse(nq);
656 Impulse = GaussianPulse(m_time, m_Psx, m_Psy, m_Psz,
657 m_Gaussianradius);
658 Vmath::Vadd(nq, &Impulse[0], 1,
659 &fields[m_intVariables[2]][0], 1,
660 &fields[m_intVariables[2]][0], 1);
661 }
662 }
663 break;
664
665 case ePlanarSource:
666 {
667 Array<OneD, NekDouble> Impulse(nq);
668 for (int i = 0; i < 3; ++i)
669 {
670 Impulse = GetIncidentField(i, m_time);
671 Vmath::Vmul(nq, m_SourceVector, 1, Impulse, 1, Impulse, 1);
672 Vmath::Vadd(nq, &Impulse[0], 1,
673 &fields[m_intVariables[i]][0], 1,
674 &fields[m_intVariables[i]][0], 1);
675 }
676 }
677 break;
678
679 default:
680 break;
681 }
682
683 // Transform data into coefficient space
684 for (i = 0; i < nvariables; ++i)
685 {
686 m_fields[m_intVariables[i]]->SetPhys(fields[i]);
687 m_fields[m_intVariables[i]]->FwdTransLocalElmt(
688 fields[i], m_fields[m_intVariables[i]]->UpdateCoeffs());
689 m_fields[m_intVariables[i]]->SetPhysState(false);
690 }
691 // for (i = 0; i < nq; ++i)
692 // std::cout << m_fields[0][0][i] <<std::endl;
693
694 // Write out checkpoint files
695 if ((m_checksteps && step && !((step + 1) % m_checksteps)) ||
696 doCheckTime)
697 {
698 indx = (step + 1) / m_checksteps;
699 TimeSeries[indx] = m_time;
700
701 if (m_TestMaxwellType == SolverUtils::eScatField2D)
702 {
703 Checkpoint_TotalFieldOutput(nchk, m_time, fields);
704 Checkpoint_TotPlotOutput(nchk, m_time, fields);
705 }
706 Checkpoint_PlotOutput(nchk, fields);
707 Checkpoint_EDFluxOutput(nchk, m_time, fields);
708 Checkpoint_EnergyOutput(nchk, m_time, fields);
709 Checkpoint_Output(nchk++);
710
711 Energy[indx] = ComputeEnergyDensity(fields);
712
713 std::cout << "|EHr|: F1 = " << RootMeanSquare(fields[0])
714 << ", F2 = " << RootMeanSquare(fields[1])
715 << ", F3 = " << RootMeanSquare(fields[2])
716 << ", Energy = " << Energy[indx] << std::endl;
717 if (nfields > 3)
718 {
719 std::cout << "|DBr|: D1 = " << RootMeanSquare(fields[3])
720 << ", D2 = " << RootMeanSquare(fields[4])
721 << ", D3 = " << RootMeanSquare(fields[5])
722 << std::endl;
723
724 int nTraceNumPoints = GetTraceNpoints();
725 int totbdryexp =
726 m_fields[0]->GetBndCondExpansions()[0]->GetExpSize();
727 int npts = m_fields[0]
728 ->GetBndCondExpansions()[0]
729 ->GetExp(0)
730 ->GetNumPoints(0);
731
732 Array<OneD, NekDouble> x0(nq);
733 Array<OneD, NekDouble> x1(nq);
734 Array<OneD, NekDouble> x2(nq);
735
736 m_fields[0]->GetCoords(x0, x1, x2);
737
738 Array<OneD, NekDouble> E1Fwd(nTraceNumPoints);
739 Array<OneD, NekDouble> E2Fwd(nTraceNumPoints);
740 Array<OneD, NekDouble> H3Fwd(nTraceNumPoints);
741
742 m_fields[0]->ExtractTracePhys(fields[0], E1Fwd);
743 m_fields[0]->ExtractTracePhys(fields[1], E2Fwd);
744 m_fields[0]->ExtractTracePhys(fields[2], H3Fwd);
745
746 int id2, cnt = 0;
747 NekDouble E1atPECloc, E1atPEC = 0.0;
748 NekDouble E2atPECloc, E2atPEC = 0.0;
749 NekDouble H3atPECloc, H3atPEC = 0.0;
750
751 Array<OneD, NekDouble> E1Fwdloc(npts);
752 Array<OneD, NekDouble> E2Fwdloc(npts);
753 Array<OneD, NekDouble> H3Fwdloc(npts);
754
755 for (int e = 0; e < totbdryexp; ++e)
756 {
757 id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
758 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
759 cnt + e));
760
761 Vmath::Vcopy(npts, &E1Fwd[id2], 1, &E1Fwdloc[0], 1);
762 Vmath::Vcopy(npts, &E2Fwd[id2], 1, &E2Fwdloc[0], 1);
763 Vmath::Vcopy(npts, &H3Fwd[id2], 1, &H3Fwdloc[0], 1);
764
765 E1atPECloc = Vmath::Vamax(npts, E1Fwdloc, 1);
766 E2atPECloc = Vmath::Vamax(npts, E2Fwdloc, 1);
767 H3atPECloc = Vmath::Vamax(npts, H3Fwdloc, 1);
768
769 if (E1atPEC < E1atPECloc)
770 {
771 E1atPEC = E1atPECloc;
772 }
773
774 if (E2atPEC < E2atPECloc)
775 {
776 E2atPEC = E2atPECloc;
777 }
778
779 if (H3atPEC < H3atPECloc)
780 {
781 H3atPEC = H3atPECloc;
782 }
783 }
784
785 std::cout << "At PEC, Max. E1 = " << E1atPEC
786 << ", E2 = " << E2atPEC << ", H3 = " << H3atPEC
787 << std::endl;
788 }
789
790 if (m_SourceType == ePointSource)
791 {
792 Ezantipod[cntap++] = fields[2][indxantipod];
793 }
794
795 if (m_TestPML)
796 {
797 P1[cntpml] = fields[2][P1indx];
798 P2[cntpml] = fields[2][P2indx];
799 P3[cntpml] = fields[2][P3indx];
800 cntpml++;
801 }
802 doCheckTime = false;
803 }
804
805 // Step advance
806 ++step;
807 }
808
809 // Print out summary statistics
810 if (m_checksteps && m_session->GetComm()->GetRank() == 0)
811 {
812 std::cout << "Time-integration : " << intTime << "s" << std::endl;
813
814 std::cout << "TimeSeries = " << std::endl;
815 for (int i = 0; i < m_steps / m_checksteps; ++i)
816 {
817 std::cout << TimeSeries[i] << ", ";
818 }
819 std::cout << std::endl << std::endl;
820
821 std::cout << "Energy Density = " << std::endl;
822 for (int i = 0; i < m_steps / m_checksteps; ++i)
823 {
824 std::cout << Energy[i] << ", ";
825 }
826 std::cout << std::endl << std::endl;
827
828 if (m_PrintoutSurfaceCurrent)
829 {
830 Printout_SurfaceCurrent(fields, m_time);
831 }
832
833 if (m_SourceType == ePointSource)
834 {
835 std::cout << "Ez at antipod = " << std::endl;
836 for (int i = 0; i < m_steps / m_checksteps; ++i)
837 {
838 std::cout << Ezantipod[i] << ", ";
839 }
840 std::cout << std::endl << std::endl;
841 }
842
843 if (m_TestPML)
844 {
845 std::cout << "P1 = " << std::endl;
846 for (int i = 0; i < m_steps / m_checksteps; ++i)
847 {
848 std::cout << P1[i] << ", ";
849 }
850 std::cout << std::endl << std::endl;
851
852 std::cout << "P2 = " << std::endl;
853 for (int i = 0; i < m_steps / m_checksteps; ++i)
854 {
855 std::cout << P2[i] << ", ";
856 }
857 std::cout << std::endl << std::endl;
858
859 std::cout << "P3 = " << std::endl;
860 for (int i = 0; i < m_steps / m_checksteps; ++i)
861 {
862 std::cout << P3[i] << ", ";
863 }
864 std::cout << std::endl << std::endl;
865 }
866 }
867
868 for (i = 0; i < nvariables; ++i)
869 {
870 m_fields[m_intVariables[i]]->SetPhys(fields[i]);
871 m_fields[m_intVariables[i]]->SetPhysState(true);
872 }
873
874 for (i = 0; i < nvariables; ++i)
875 {
876 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
877 m_fields[i]->UpdateCoeffs());
878 }
879}
880
881/**
882 * @brief Compute the right-hand side for the linear advection equation.
883 *
884 * @param inarray Given fields.
885 * @param outarray Calculated solution.
886 * @param time Time.
887 */
888void MMFMaxwell::DoOdeRhs(
889 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
890 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
891{
892 int i;
893 int nvar = inarray.size();
894 int ncoeffs = GetNcoeffs();
895 int nq = GetTotPoints();
896
897 Array<OneD, Array<OneD, NekDouble>> physarray(nvar);
898 Array<OneD, Array<OneD, NekDouble>> modarray(nvar);
899 for (i = 0; i < nvar; ++i)
900 {
901 physarray[i] = Array<OneD, NekDouble>(nq);
902 modarray[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
903
904 Vmath::Vcopy(nq, &inarray[i][0], 1, &physarray[i][0], 1);
905 }
906
907 for (i = 0; i < nvar; i++)
908 {
909 m_fields[i]->SetPhysState(true);
910 }
911
912 // Compute Curl
913 switch (m_TestMaxwellType)
914 {
915 case SolverUtils::eScatField2D:
916 case SolverUtils::eTotField2D:
917
918 {
919
920 // Imaginary part is computed the same as Real part
921 Array<OneD, Array<OneD, NekDouble>> tmpin(3);
922 Array<OneD, Array<OneD, NekDouble>> tmpout(3);
923
924 for (int i = 0; i < 3; ++i)
925 {
926 tmpin[i] = Array<OneD, NekDouble>(nq);
927 tmpout[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
928
929 Vmath::Vcopy(nq, &physarray[i][0], 1, &tmpin[i][0], 1);
930 }
931
932 WeakDGMaxwellDirDeriv(tmpin, tmpout, time);
933 AddGreenDerivCompensate(tmpin, tmpout);
934
935 for (int i = 0; i < 3; ++i)
936 {
937 // For E and H
938 Vmath::Vcopy(ncoeffs, &tmpout[i][0], 1, &modarray[i][0], 1);
939 }
940 }
941 break;
942
943 default:
944 {
945
946 WeakDGMaxwellDirDeriv(physarray, modarray, time);
947 AddGreenDerivCompensate(physarray, modarray);
948 }
949 break;
950 }
951
952 for (i = 0; i < nvar; ++i)
953 {
954 m_fields[i]->MultiplyByElmtInvMass(modarray[i], modarray[i]);
955 m_fields[i]->BwdTrans(modarray[i], outarray[i]);
956 }
957
958 if (m_TestMaxwellType == SolverUtils::eMaxwellSphere)
959 {
960 Array<OneD, NekDouble> F(nq);
961 for (int j = 0; j < 2; ++j)
962 {
963 F = TestMaxwellSphere(time, m_freq, 3 + j);
964 Vmath::Vadd(nq, &F[0], 1, &outarray[j][0], 1, &outarray[j][0], 1);
965 }
966 }
967
968 // Add Absorbing Boundary Conditions
969 if (m_AddPML > 0)
970 {
971 AddPML(physarray, outarray);
972 }
973
974 // Add dedt component
975 if (m_AddRotation)
976 {
977 AddCoriolis(physarray, outarray);
978 AdddedtMaxwell(physarray, outarray);
979 }
980
981 // Divide it by varepsilon or mu
982 Array<OneD, NekDouble> dFdt(nq, 0.0);
983 switch (m_TestMaxwellType)
984 {
985 case SolverUtils::eMaxwell1D:
986 {
987 Vmath::Vdiv(nq, outarray[0], 1, m_epsvec[0], 1, outarray[0], 1);
988 Vmath::Vdiv(nq, outarray[1], 1, m_muvec[0], 1, outarray[1], 1);
989 }
990 break;
991
992 // TO BE CHANGED
993 case SolverUtils::eTestMaxwell2DPECAVGFLUX:
994 {
995 Array<OneD, NekDouble> Hxdt(nq, 0.0);
996 Array<OneD, NekDouble> Hydt(nq, 0.0);
997
998 Hxdt = TestMaxwell2DPEC(time, 10, m_PolType);
999 Hydt = TestMaxwell2DPEC(time, 11, m_PolType);
1000
1001 Array<OneD, NekDouble> x0(nq);
1002 Array<OneD, NekDouble> x1(nq);
1003 Array<OneD, NekDouble> x2(nq);
1004
1005 m_fields[0]->GetCoords(x0, x1, x2);
1006
1007 NekDouble theta, tmpx, tmpy;
1008 NekDouble uxx, uxy, uyy, detu, uti, uri;
1009 Array<OneD, NekDouble> utvec(nq, 1.0);
1010 Array<OneD, NekDouble> urvec(nq, 1.0);
1011 Array<OneD, NekDouble> tmpIN(nq);
1012
1013 // Case I: ut = 4.0, ur = 0.5
1014 // NekDouble ut=4.0;
1015 // NekDouble ur=0.5;
1016
1017 // m_fields[0]->GenerateElementVector(m_ElemtGroup1, ut, 1.0,
1018 // utvec);
1019 // m_fields[0]->GenerateElementVector(m_ElemtGroup1, ur, 1.0,
1020 // urvec);
1021
1022 // Case II: ut = 0.5, ur = 1 - 2x^2
1023 NekDouble ut = 0.5;
1024 m_fields[0]->GenerateElementVector(m_ElemtGroup1, ut, 1.0, utvec);
1025
1026 m_fields[0]->GenerateElementVector(m_ElemtGroup1, 1.0, 0.0, tmpIN);
1027
1028 for (int i = 0; i < nq; i++)
1029 {
1030 urvec[i] =
1031 tmpIN[i] * (1.0 - 2 * x0[i] * x0[i]) + (1.0 - tmpIN[i]);
1032 }
1033
1034 for (int i = 0; i < nq; ++i)
1035 {
1036 theta = atan2((x1[i] + 2.0), (x0[i] + 2.0));
1037
1038 uti = utvec[i];
1039 uri = urvec[i];
1040
1041 uxx = uti * cos(theta) * cos(theta) +
1042 uri * sin(theta) * sin(theta);
1043 uyy = uti * sin(theta) * sin(theta) +
1044 uri * cos(theta) * cos(theta);
1045 uxy = (uti - uri) * cos(theta) * sin(theta);
1046
1047 detu = uxx * uyy - uxy * uxy;
1048
1049 tmpx = outarray[0][i] + (1.0 - uxx) * Hxdt[i] - uxy * Hydt[i];
1050 tmpy = outarray[1][i] - uxy * Hxdt[i] + (1.0 - uyy) * Hydt[i];
1051
1052 outarray[0][i] = (1 / detu) * (uyy * tmpx - uxy * tmpy);
1053 outarray[1][i] = (1 / detu) * (-uxy * tmpx + uxx * tmpy);
1054 }
1055 }
1056 break;
1057
1058 case SolverUtils::eTestMaxwell2DPEC:
1059 case SolverUtils::eScatField2D:
1060 case SolverUtils::eTotField2D:
1061 {
1062 switch (m_PolType)
1063 {
1064 case SolverUtils::eTransMagnetic:
1065 {
1066 if (m_TestMaxwellType == SolverUtils::eTestMaxwell2DPEC)
1067 {
1068 dFdt = TestMaxwell2DPEC(time, 10, m_PolType);
1069 Vmath::Vvtvp(nq, m_negmuvecminus1[0], 1, dFdt, 1,
1070 outarray[0], 1, outarray[0], 1);
1071
1072 dFdt = TestMaxwell2DPEC(time, 11, m_PolType);
1073 Vmath::Vvtvp(nq, m_negmuvecminus1[1], 1, dFdt, 1,
1074 outarray[1], 1, outarray[1], 1);
1075
1076 dFdt = TestMaxwell2DPEC(time, 12, m_PolType);
1077 Vmath::Vvtvp(nq, m_negepsvecminus1[2], 1, dFdt, 1,
1078 outarray[2], 1, outarray[2], 1);
1079 }
1080
1081 if (m_TestMaxwellType == SolverUtils::eScatField2D)
1082 {
1083 dFdt = GetIncidentField(10, time);
1084 Vmath::Vvtvp(nq, m_negmuvecminus1[0], 1, dFdt, 1,
1085 outarray[0], 1, outarray[0], 1);
1086
1087 dFdt = GetIncidentField(11, time);
1088 Vmath::Vvtvp(nq, m_negmuvecminus1[1], 1, dFdt, 1,
1089 outarray[1], 1, outarray[1], 1);
1090
1091 dFdt = GetIncidentField(12, time);
1092 Vmath::Vvtvp(nq, m_negepsvecminus1[2], 1, dFdt, 1,
1093 outarray[2], 1, outarray[2], 1);
1094 }
1095
1096 Vmath::Vdiv(nq, outarray[0], 1, m_muvec[0], 1, outarray[0],
1097 1);
1098 Vmath::Vdiv(nq, outarray[1], 1, m_muvec[1], 1, outarray[1],
1099 1);
1100 Vmath::Vdiv(nq, outarray[2], 1, m_epsvec[2], 1, outarray[2],
1101 1);
1102 }
1103 break;
1104
1105 case SolverUtils::eTransElectric:
1106 {
1107 if (m_TestMaxwellType == SolverUtils::eTestMaxwell2DPEC)
1108 {
1109 // (I - \mu^i) d F^{inc} / dt
1110 dFdt = TestMaxwell2DPEC(time, 10, m_PolType);
1111 Vmath::Vvtvp(nq, m_negepsvecminus1[0], 1, dFdt, 1,
1112 outarray[0], 1, outarray[0], 1);
1113
1114 dFdt = TestMaxwell2DPEC(time, 11, m_PolType);
1115 Vmath::Vvtvp(nq, m_negepsvecminus1[1], 1, dFdt, 1,
1116 outarray[1], 1, outarray[1], 1);
1117
1118 dFdt = TestMaxwell2DPEC(time, 12, m_PolType);
1119 Vmath::Vvtvp(nq, m_negmuvecminus1[2], 1, dFdt, 1,
1120 outarray[2], 1, outarray[2], 1);
1121 }
1122
1123 if (m_TestMaxwellType == SolverUtils::eScatField2D)
1124 {
1125 dFdt = GetIncidentField(10, time);
1126 Vmath::Vvtvp(nq, m_negepsvecminus1[0], 1, dFdt, 1,
1127 outarray[0], 1, outarray[0], 1);
1128
1129 // Add - wp^2 \int E_2^{inc}
1130 dFdt = GetIncidentField(21, time);
1131 if (m_DispersiveCloak)
1132 {
1133 Vmath::Vmul(nq, m_wp2, 1, dFdt, 1, dFdt, 1);
1134 Vmath::Vsub(nq, outarray[1], 1, dFdt, 1,
1135 outarray[1], 1);
1136 }
1137
1138 else
1139 {
1140 Vmath::Vvtvp(nq, m_negepsvecminus1[1], 1, dFdt, 1,
1141 outarray[1], 1, outarray[1], 1);
1142 }
1143
1144 dFdt = GetIncidentField(12, time);
1145 Vmath::Vvtvp(nq, m_negmuvecminus1[2], 1, dFdt, 1,
1146 outarray[2], 1, outarray[2], 1);
1147 }
1148
1149 Vmath::Vdiv(nq, outarray[0], 1, m_epsvec[0], 1, outarray[0],
1150 1);
1151 Vmath::Vdiv(nq, outarray[1], 1, m_epsvec[1], 1, outarray[1],
1152 1);
1153 Vmath::Vdiv(nq, outarray[2], 1, m_muvec[2], 1, outarray[2],
1154 1);
1155 }
1156 break;
1157
1158 default:
1159 break;
1160 }
1161 }
1162 break; // case SolverUtils::eTestMaxwell2DPEC:
1163
1164 default:
1165 break;
1166
1167 } // switch(m_TestMaxwellType)
1168}
1169
1170void MMFMaxwell::AddGreenDerivCompensate(
1171 const Array<OneD, const Array<OneD, NekDouble>> &physarray,
1172 Array<OneD, Array<OneD, NekDouble>> &outarray)
1173{
1174 // routine works for both primitive and conservative formulations
1175 int ncoeffs = outarray[0].size();
1176 int nq = physarray[0].size();
1177
1178 Array<OneD, NekDouble> tmp(nq);
1179 Array<OneD, NekDouble> tmpc(ncoeffs);
1180
1181 Array<OneD, Array<OneD, NekDouble>> fluxvector(m_shapedim);
1182 for (int j = 0; j < m_shapedim; ++j)
1183 {
1184 fluxvector[j] = Array<OneD, NekDouble>(nq);
1185 }
1186
1187 // m_CurlMF[0][0] = e^3 \cdot (\nabla \times e^1) [ NEW m_CurlMF[0][2] ]
1188 // m_CurlMF[0][1] = 0.0
1189 // m_CurlMF[1][0] = 0.0,
1190 // m_CurlMF[1][1] = e^3 \cdot (\nabla \times e^2) [ NEW m_CurlMF[1][2] ]
1191 // m_CurlMF[2][0] = e^1 \cdot (\nabla \times e^3) [ NEW m_CurlMF[2][0] ]
1192 // m_CurlMF[2][1] = e^2 \cdot (\nabla \times e^3) [ NEW m_CurlMF[2][1] ]
1193
1194 int var;
1195
1196 switch (m_TestMaxwellType)
1197 {
1198 case SolverUtils::eTestMaxwell2DPEC:
1199 case SolverUtils::eTestMaxwell2DPECAVGFLUX:
1200 case SolverUtils::eTestMaxwell2DPMC:
1201 case SolverUtils::eScatField2D:
1202 case SolverUtils::eTotField2D:
1203 case SolverUtils::eMaxwellSphere:
1204 case SolverUtils::eELF2DSurface:
1205 {
1206 var = 0;
1207 GetMaxwellFluxVector(var, physarray, fluxvector);
1208 Vmath::Vmul(nq, &fluxvector[0][0], 1, &m_CurlMF[0][2][0], 1,
1209 &tmp[0], 1);
1210 m_fields[var]->IProductWRTBase(tmp, tmpc);
1211 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[var], 1, outarray[var], 1);
1212
1213 var = 1;
1214 GetMaxwellFluxVector(var, physarray, fluxvector);
1215 Vmath::Vmul(nq, &fluxvector[1][0], 1, &m_CurlMF[1][2][0], 1,
1216 &tmp[0], 1);
1217 Vmath::Neg(nq, tmp, 1);
1218 m_fields[var]->IProductWRTBase(tmp, tmpc);
1219 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[var], 1, outarray[var], 1);
1220
1221 var = 2;
1222 GetMaxwellFluxVector(var, physarray, fluxvector);
1223 Vmath::Vmul(nq, &fluxvector[0][0], 1, &m_CurlMF[2][0][0], 1,
1224 &tmp[0], 1);
1225 Vmath::Vvtvm(nq, &fluxvector[1][0], 1, &m_CurlMF[2][1][0], 1,
1226 &tmp[0], 1, &tmp[0], 1);
1227 m_fields[var]->IProductWRTBase(tmp, tmpc);
1228 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[var], 1, outarray[var], 1);
1229 }
1230 break;
1231
1232 default:
1233 break;
1234 }
1235}
1236
1237/**
1238 * @brief Calculate weak DG advection in the form \f$ \langle\phi,
1239 * \hat{F}\cdot n\rangle - (\nabla \phi \cdot F) \f$
1240 *
1241 * @param InField Fields.
1242 * @param OutField Storage for result.
1243 * @param NumericalFluxIncludesNormal Default: true.
1244 * @param InFieldIsPhysSpace Default: false.
1245 * @param nvariables Number of fields.
1246 */
1247void MMFMaxwell::WeakDGMaxwellDirDeriv(
1248 const Array<OneD, const Array<OneD, NekDouble>> &InField,
1249 Array<OneD, Array<OneD, NekDouble>> &OutField, const NekDouble time)
1250{
1251 int i;
1252 int nq = GetNpoints();
1253 int ncoeffs = GetNcoeffs();
1254 int nTracePointsTot = GetTraceNpoints();
1255 int nvar = 3;
1256
1257 Array<OneD, Array<OneD, NekDouble>> fluxvector(m_shapedim);
1258 for (i = 0; i < m_shapedim; ++i)
1259 {
1260 fluxvector[i] = Array<OneD, NekDouble>(nq);
1261 }
1262
1263 Array<OneD, Array<OneD, NekDouble>> physfield(nvar);
1264 for (i = 0; i < nvar; ++i)
1265 {
1266 physfield[i] = InField[i];
1267 }
1268
1269 Array<OneD, NekDouble> tmpc(ncoeffs);
1270 for (i = 0; i < nvar; ++i)
1271 {
1272 GetMaxwellFluxVector(i, physfield, fluxvector);
1273
1274 OutField[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
1275 for (int j = 0; j < m_shapedim; ++j)
1276 {
1277 // Directional derivation with respect to the j'th moving frame
1278 // tmp_j = ( \nabla \phi, fluxvector[j] \mathbf{e}^j )
1279 m_fields[i]->IProductWRTDirectionalDerivBase(m_CrossProductMF[j],
1280 fluxvector[j], tmpc);
1281 Vmath::Vadd(ncoeffs, &tmpc[0], 1, &OutField[i][0], 1,
1282 &OutField[i][0], 1);
1283 }
1284 }
1285
1286 // V the numerical flux and add to the modal coeffs
1287 // if the NumericalFlux function does not include the
1288 // normal in the output
1289 Array<OneD, Array<OneD, NekDouble>> numfluxFwd(nvar);
1290 Array<OneD, Array<OneD, NekDouble>> numfluxBwd(nvar);
1291
1292 for (i = 0; i < nvar; ++i)
1293 {
1294 numfluxFwd[i] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
1295 numfluxBwd[i] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
1296 }
1297
1298 // Evaluate numerical flux in physical space which may in
1299 // general couple all component of vectors
1300 NumericalMaxwellFlux(physfield, numfluxFwd, numfluxBwd, time);
1301
1302 // Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
1303 for (i = 0; i < nvar; ++i)
1304 {
1305 Vmath::Neg(ncoeffs, OutField[i], 1);
1306 m_fields[i]->AddFwdBwdTraceIntegral(numfluxFwd[i], numfluxBwd[i],
1307 OutField[i]);
1308 m_fields[i]->SetPhysState(false);
1309 }
1310}
1311
1312/**
1313 * @brief Compute the projection for the linear advection equation.
1314 *
1315 * @param inarray Given fields.
1316 * @param outarray Calculated solution.
1317 * @param time Time.
1318 */
1319void MMFMaxwell::DoOdeProjection(
1320 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1321 Array<OneD, Array<OneD, NekDouble>> &outarray,
1322 [[maybe_unused]] const NekDouble time)
1323{
1324 int var = inarray.size();
1325
1326 if (inarray != outarray)
1327 {
1328 int nq = GetNpoints();
1329 for (int i = 0; i < var; ++i)
1330 {
1331 Vmath::Vcopy(nq, inarray[i], 1, outarray[i], 1);
1332 }
1333 }
1334}
1335
1336void MMFMaxwell::v_SetInitialConditions(const NekDouble initialtime,
1337 bool dumpInitialConditions,
1338 [[maybe_unused]] const int domain)
1339{
1340 int nq = GetTotPoints();
1341 int nvar = m_fields.size();
1342
1343 switch (m_TestMaxwellType)
1344 {
1345 case SolverUtils::eMaxwell1D:
1346 {
1347 m_fields[0]->SetPhys(TestMaxwell1D(initialtime, 0));
1348 m_fields[1]->SetPhys(TestMaxwell1D(initialtime, 1));
1349 }
1350 break;
1351
1352 case SolverUtils::eTestMaxwell2DPEC:
1353 case SolverUtils::eTestMaxwell2DPECAVGFLUX:
1354 {
1355 m_fields[0]->SetPhys(TestMaxwell2DPEC(initialtime, 0, m_PolType));
1356 m_fields[1]->SetPhys(TestMaxwell2DPEC(initialtime, 1, m_PolType));
1357 m_fields[2]->SetPhys(TestMaxwell2DPEC(initialtime, 2, m_PolType));
1358 }
1359 break;
1360
1361 case SolverUtils::eTestMaxwell2DPMC:
1362 {
1363 m_fields[0]->SetPhys(TestMaxwell2DPMC(initialtime, 0, m_PolType));
1364 m_fields[1]->SetPhys(TestMaxwell2DPMC(initialtime, 1, m_PolType));
1365 m_fields[2]->SetPhys(TestMaxwell2DPMC(initialtime, 2, m_PolType));
1366 }
1367 break;
1368
1369 case SolverUtils::eScatField2D:
1370 case SolverUtils::eTotField2D:
1371 {
1372 Array<OneD, NekDouble> Zeros(nq, 0.0);
1373
1374 for (int i = 0; i < nvar; i++)
1375 {
1376 m_fields[i]->SetPhys(Zeros);
1377 }
1378 }
1379 break;
1380
1381 case SolverUtils::eMaxwellSphere:
1382 {
1383 m_fields[0]->SetPhys(TestMaxwellSphere(initialtime, m_freq, 0));
1384 m_fields[1]->SetPhys(TestMaxwellSphere(initialtime, m_freq, 1));
1385 m_fields[2]->SetPhys(TestMaxwellSphere(initialtime, m_freq, 2));
1386 }
1387 break;
1388
1389 case SolverUtils::eELF2DSurface:
1390 {
1391 m_fields[2]->SetPhys(GaussianPulse(initialtime, m_Psx, m_Psy, m_Psz,
1392 m_Gaussianradius));
1393 }
1394 break;
1395
1396 default:
1397 break;
1398 }
1399
1400 // forward transform to fill the modal coeffs
1401 for (int i = 0; i < nvar; ++i)
1402 {
1403 m_fields[i]->SetPhysState(true);
1404 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
1405 m_fields[i]->UpdateCoeffs());
1406 }
1407
1408 if (dumpInitialConditions)
1409 {
1410 std::string outname = m_sessionName + "_initial.chk";
1411 WriteFld(outname);
1412
1413 Array<OneD, Array<OneD, NekDouble>> fields(nvar);
1414 for (int i = 0; i < nvar; ++i)
1415 {
1416 fields[i] = m_fields[i]->GetPhys();
1417 }
1418
1419 Checkpoint_PlotOutput(0, fields);
1420 }
1421}
1422
1423void MMFMaxwell::v_EvaluateExactSolution(unsigned int field,
1424 Array<OneD, NekDouble> &outfield,
1425 const NekDouble time)
1426{
1427 int nq = m_fields[0]->GetNpoints();
1428 outfield = Array<OneD, NekDouble>(nq);
1429
1430 switch (m_TestMaxwellType)
1431 {
1432 case SolverUtils::eMaxwell1D:
1433 {
1434 outfield = TestMaxwell1D(time, field);
1435 }
1436 break;
1437
1438 case SolverUtils::eTestMaxwell2DPEC:
1439 case SolverUtils::eTestMaxwell2DPECAVGFLUX:
1440 {
1441 outfield = TestMaxwell2DPEC(time, field, m_PolType);
1442 }
1443 break;
1444
1445 case SolverUtils::eTestMaxwell2DPMC:
1446 {
1447 outfield = TestMaxwell2DPMC(time, field, m_PolType);
1448 }
1449 break;
1450
1451 case SolverUtils::eMaxwellSphere:
1452 {
1453 outfield = TestMaxwellSphere(time, m_freq, field);
1454 }
1455 break;
1456
1457 default:
1458 {
1459 outfield = Array<OneD, NekDouble>(nq, 0.0);
1460 }
1461 break;
1462 }
1463}
1464
1465Array<OneD, NekDouble> MMFMaxwell::TestMaxwell1D(const NekDouble time,
1466 unsigned int field)
1467{
1468 int nq = m_fields[0]->GetNpoints();
1469
1470 Array<OneD, NekDouble> x0(nq);
1471 Array<OneD, NekDouble> x1(nq);
1472 Array<OneD, NekDouble> x2(nq);
1473
1474 m_fields[0]->GetCoords(x0, x1, x2);
1475
1476 Array<OneD, NekDouble> E(nq);
1477 Array<OneD, NekDouble> H(nq);
1478
1479 // Derive the frequency \omega
1480 NekDouble omega;
1481 NekDouble Tol = 0.000000001;
1482 if (fabs(m_n1 - m_n2) < Tol)
1483 {
1484 omega = m_pi / m_n1;
1485 }
1486
1487 else
1488 {
1489 omega = 2.0 * m_pi / m_n2;
1490
1491 NekDouble newomega, F, Fprime;
1492 for (int i = 0; i < 10000; ++i)
1493 {
1494 F = m_n1 * tan(m_n2 * omega) + m_n2 * tan(m_n1 * omega);
1495 Fprime = m_n1 * m_n2 *
1496 (1.0 / cos(m_n2 * omega) / cos(m_n2 * omega) +
1497 1.0 / cos(m_n1 * omega) / cos(m_n1 * omega));
1498
1499 newomega = omega - F / Fprime;
1500
1501 if (fabs(newomega - omega) > Tol)
1502 {
1503 omega = newomega;
1504 }
1505
1506 else
1507 {
1508 break;
1509 }
1510 }
1511 }
1512
1513 // Generate A^k and B^k
1514 std::complex<double> im = sqrt(std::complex<double>(-1));
1515 std::complex<double> A1, A2, B1, B2;
1516 std::complex<double> Ak, Bk, nk;
1517 std::complex<double> Ec, Hc;
1518
1519 A1 = m_n2 * cos(m_n2 * omega) / (m_n1 * cos(m_n1 * omega));
1520 A2 = exp(-1.0 * im * omega * (m_n1 + m_n2));
1521 B1 = A1 * exp(-2.0 * im * m_n1 * omega);
1522 B2 = A2 * exp(2.0 * im * m_n2 * omega);
1523
1524 for (int i = 0; i < nq; ++i)
1525 {
1526 if (x0[i] > 0)
1527 {
1528 Ak = A2;
1529 Bk = B2;
1530 nk = m_n2;
1531 }
1532
1533 else
1534 {
1535 Ak = A1;
1536 Bk = B1;
1537 nk = m_n1;
1538 }
1539
1540 Ec = (Ak * exp(im * nk * omega * x0[i]) -
1541 Bk * exp(-im * nk * omega * x0[i])) *
1542 exp(im * omega * time);
1543 Hc = nk *
1544 (Ak * exp(im * nk * omega * x0[i]) +
1545 Bk * exp(-im * nk * omega * x0[i])) *
1546 exp(im * omega * time);
1547
1548 E[i] = Ec.real();
1549 H[i] = Hc.real();
1550 }
1551
1552 Array<OneD, NekDouble> outfield;
1553 switch (field)
1554 {
1555 case (0):
1556 {
1557 outfield = E;
1558 }
1559 break;
1560
1561 case (1):
1562 {
1563 outfield = H;
1564 }
1565 break;
1566 }
1567
1568 return outfield;
1569}
1570
1571Array<OneD, NekDouble> MMFMaxwell::TestMaxwell2DPEC(
1572 const NekDouble time, unsigned int field,
1573 const SolverUtils::PolType Polarization)
1574{
1575 int nq = m_fields[0]->GetNpoints();
1576
1577 Array<OneD, NekDouble> x0(nq);
1578 Array<OneD, NekDouble> x1(nq);
1579 Array<OneD, NekDouble> x2(nq);
1580
1581 m_fields[0]->GetCoords(x0, x1, x2);
1582
1583 NekDouble freqm = 1.0, freqn = 1.0;
1584 NekDouble omega = m_pi * sqrt(freqm * freqm + freqn * freqn);
1585 NekDouble mpi = freqm * m_pi;
1586 NekDouble npi = freqn * m_pi;
1587
1588 Array<OneD, NekDouble> F1(nq);
1589 Array<OneD, NekDouble> F2(nq);
1590 Array<OneD, NekDouble> Fz(nq);
1591 Array<OneD, NekDouble> dF1dt(nq);
1592 Array<OneD, NekDouble> dF2dt(nq);
1593 Array<OneD, NekDouble> dFzdt(nq);
1594 NekDouble Fx, Fy, dFxdt, dFydt;
1595
1596 for (int i = 0; i < nq; ++i)
1597 {
1598 switch (Polarization)
1599 {
1600 case SolverUtils::eTransMagnetic:
1601 {
1602 Fx = -1.0 * (npi / omega) * sin(mpi * x0[i]) *
1603 cos(npi * x1[i]) * sin(omega * time);
1604 Fy = (mpi / omega) * cos(mpi * x0[i]) * sin(npi * x1[i]) *
1605 sin(omega * time);
1606
1607 F1[i] =
1608 Fx * m_movingframes[0][i] + Fy * m_movingframes[0][i + nq];
1609 F2[i] =
1610 Fx * m_movingframes[1][i] + Fy * m_movingframes[1][i + nq];
1611 Fz[i] = sin(mpi * x0[i]) * sin(npi * x1[i]) * cos(omega * time);
1612
1613 dFxdt = (npi)*sin(mpi * x0[i]) * cos(npi * x1[i]) *
1614 cos(omega * time);
1615 dFydt = -1.0 * (mpi)*cos(mpi * x0[i]) * sin(npi * x1[i]) *
1616 cos(omega * time);
1617
1618 dF1dt[i] = dFxdt * m_movingframes[0][i] +
1619 dFydt * m_movingframes[0][i + nq];
1620 dF2dt[i] = dFxdt * m_movingframes[1][i] +
1621 dFydt * m_movingframes[1][i + nq];
1622 dFzdt[i] = omega * sin(mpi * x0[i]) * sin(npi * x1[i]) *
1623 sin(omega * time);
1624 }
1625 break;
1626
1627 case SolverUtils::eTransElectric:
1628 {
1629 Fx = -1.0 * (npi / omega) * cos(mpi * x0[i]) *
1630 sin(npi * x1[i]) * sin(omega * time);
1631 Fy = (mpi / omega) * sin(mpi * x0[i]) * cos(npi * x1[i]) *
1632 sin(omega * time);
1633
1634 F1[i] =
1635 Fx * m_movingframes[0][i] + Fy * m_movingframes[0][i + nq];
1636 F2[i] =
1637 Fx * m_movingframes[1][i] + Fy * m_movingframes[1][i + nq];
1638 Fz[i] = cos(mpi * x0[i]) * cos(npi * x1[i]) * cos(omega * time);
1639
1640 dFxdt = (npi)*cos(mpi * x0[i]) * sin(npi * x1[i]) *
1641 cos(omega * time);
1642 dFydt = -1.0 * (mpi)*sin(mpi * x0[i]) * cos(npi * x1[i]) *
1643 cos(omega * time);
1644
1645 dF1dt[i] = dFxdt * m_movingframes[0][i] +
1646 dFydt * m_movingframes[0][i + nq];
1647 dF2dt[i] = dFxdt * m_movingframes[1][i] +
1648 dFydt * m_movingframes[1][i + nq];
1649 dFzdt[i] = omega * cos(mpi * x0[i]) * cos(npi * x1[i]) *
1650 sin(omega * time);
1651 }
1652 break;
1653
1654 default:
1655 break;
1656 }
1657 }
1658
1659 Array<OneD, NekDouble> outfield;
1660 switch (field)
1661 {
1662 case (0):
1663 {
1664 outfield = F1;
1665 }
1666 break;
1667
1668 case (1):
1669 {
1670 outfield = F2;
1671 }
1672 break;
1673
1674 case (2):
1675 {
1676 outfield = Fz;
1677 }
1678 break;
1679
1680 case (10):
1681 {
1682 outfield = dF1dt;
1683 }
1684 break;
1685
1686 case (11):
1687 {
1688 outfield = dF2dt;
1689 }
1690 break;
1691
1692 case (12):
1693 {
1694 outfield = dFzdt;
1695 }
1696 break;
1697 }
1698
1699 return outfield;
1700}
1701
1702Array<OneD, NekDouble> MMFMaxwell::TestMaxwell2DPMC(
1703 const NekDouble time, unsigned int field,
1704 const SolverUtils::PolType Polarization)
1705{
1706 int nq = m_fields[0]->GetNpoints();
1707
1708 Array<OneD, NekDouble> x0(nq);
1709 Array<OneD, NekDouble> x1(nq);
1710 Array<OneD, NekDouble> x2(nq);
1711
1712 m_fields[0]->GetCoords(x0, x1, x2);
1713
1714 NekDouble freqm = 1.0, freqn = 1.0;
1715 NekDouble omega = m_pi * sqrt(freqm * freqm + freqn * freqn);
1716 NekDouble mpi = freqm * m_pi;
1717 NekDouble npi = freqn * m_pi;
1718
1719 Array<OneD, NekDouble> F1(nq);
1720 Array<OneD, NekDouble> F2(nq);
1721 Array<OneD, NekDouble> Fz(nq);
1722 NekDouble Fx, Fy;
1723
1724 for (int i = 0; i < nq; ++i)
1725 {
1726 switch (Polarization)
1727 {
1728 case SolverUtils::eTransMagnetic:
1729 {
1730 Fx = (npi / omega) * cos(mpi * x0[i]) * sin(npi * x1[i]) *
1731 sin(omega * time);
1732 Fy = -(mpi / omega) * sin(mpi * x0[i]) * cos(npi * x1[i]) *
1733 sin(omega * time);
1734
1735 F1[i] =
1736 Fx * m_movingframes[0][i] + Fy * m_movingframes[0][i + nq];
1737 F2[i] =
1738 Fx * m_movingframes[1][i] + Fy * m_movingframes[1][i + nq];
1739 Fz[i] = cos(mpi * x0[i]) * cos(npi * x1[i]) * cos(omega * time);
1740 }
1741 break;
1742
1743 case SolverUtils::eTransElectric:
1744 {
1745 Fx = (npi / omega) * sin(mpi * x0[i]) * cos(npi * x1[i]) *
1746 sin(omega * time);
1747 Fy = -(mpi / omega) * cos(mpi * x0[i]) * sin(npi * x1[i]) *
1748 sin(omega * time);
1749
1750 F1[i] =
1751 Fx * m_movingframes[0][i] + Fy * m_movingframes[0][i + nq];
1752 F2[i] =
1753 Fx * m_movingframes[1][i] + Fy * m_movingframes[1][i + nq];
1754 Fz[i] = sin(mpi * x0[i]) * sin(npi * x1[i]) * cos(omega * time);
1755 }
1756 break;
1757
1758 default:
1759 break;
1760 }
1761 }
1762
1763 Array<OneD, NekDouble> outfield;
1764 switch (field)
1765 {
1766 case (0):
1767 {
1768 outfield = F1;
1769 }
1770 break;
1771
1772 case (1):
1773 {
1774 outfield = F2;
1775 }
1776 break;
1777
1778 case (2):
1779 {
1780 outfield = Fz;
1781 }
1782 break;
1783 }
1784
1785 return outfield;
1786}
1787
1788Array<OneD, NekDouble> MMFMaxwell::TestMaxwellSphere(const NekDouble time,
1789 const NekDouble omega,
1790 unsigned int field)
1791{
1792 int nq = m_fields[0]->GetTotPoints();
1793
1794 Array<OneD, NekDouble> outfield(nq);
1795
1796 Array<OneD, NekDouble> x(nq);
1797 Array<OneD, NekDouble> y(nq);
1798 Array<OneD, NekDouble> z(nq);
1799
1800 m_fields[0]->GetCoords(x, y, z);
1801
1802 Array<OneD, NekDouble> H1(nq);
1803 Array<OneD, NekDouble> H2(nq);
1804 Array<OneD, NekDouble> E3(nq);
1805 Array<OneD, NekDouble> F1(nq);
1806 Array<OneD, NekDouble> F2(nq);
1807
1808 Array<OneD, NekDouble> curlv(nq);
1809 Array<OneD, Array<OneD, NekDouble>> velvec(m_spacedim);
1810 Array<OneD, Array<OneD, NekDouble>> Fvec(m_spacedim);
1811 for (int i = 0; i < m_spacedim; ++i)
1812 {
1813 velvec[i] = Array<OneD, NekDouble>(nq, 0.0);
1814 Fvec[i] = Array<OneD, NekDouble>(nq, 0.0);
1815 }
1816
1817 NekDouble xj, yj, zj, sin_varphi, cos_varphi, sin_theta, cos_theta;
1818 NekDouble vth, vphi, Fth, Fphi;
1819 for (int i = 0; i < nq; i++)
1820 {
1821 xj = x[i];
1822 yj = y[i];
1823 zj = z[i];
1824
1825 CartesianToSpherical(xj, yj, zj, sin_varphi, cos_varphi, sin_theta,
1826 cos_theta);
1827
1828 vth = -4.0 * sin_varphi * cos_varphi * cos_theta * cos_theta *
1829 cos_theta * sin_theta;
1830 vphi =
1831 -1.0 * sin_varphi * sin_varphi * cos_theta * cos_theta * cos_theta;
1832 velvec[0][i] = -vth * sin_theta * cos_varphi - vphi * sin_varphi;
1833 velvec[1][i] = -vth * sin_theta * sin_varphi + vphi * cos_varphi;
1834 velvec[2][i] = vth * cos_theta;
1835
1836 E3[i] = (-4.0 * cos_theta * cos_theta * sin_theta * cos_varphi *
1837 cos_varphi) *
1838 (1.0 / omega * sin(omega * time));
1839
1840 Fth = -omega * vth -
1841 (8.0 / omega) * cos_theta * sin_theta * cos_varphi * sin_varphi;
1842 Fphi = -omega * vphi +
1843 (4.0 / omega) * cos_varphi * cos_varphi * cos_theta *
1844 (2.0 * sin_theta * sin_theta - cos_theta * cos_theta);
1845 Fvec[0][i] = -Fth * sin_theta * cos_varphi - Fphi * sin_varphi;
1846 Fvec[1][i] = -Fth * sin_theta * sin_varphi + Fphi * cos_varphi;
1847 Fvec[2][i] = Fth * cos_theta;
1848 }
1849
1850 H1 = CartesianToMovingframes(velvec, 0);
1851 H2 = CartesianToMovingframes(velvec, 1);
1852
1853 Vmath::Smul(nq, cos(omega * time), H1, 1, H1, 1);
1854 Vmath::Smul(nq, cos(omega * time), H2, 1, H2, 1);
1855
1856 F1 = CartesianToMovingframes(Fvec, 0);
1857 F2 = CartesianToMovingframes(Fvec, 1);
1858
1859 Vmath::Smul(nq, sin(omega * time), F1, 1, F1, 1);
1860 Vmath::Smul(nq, sin(omega * time), F2, 1, F2, 1);
1861
1862 switch (field)
1863 {
1864 // return H1
1865 case 0:
1866 {
1867 outfield = H1;
1868 }
1869 break;
1870
1871 case 1:
1872 {
1873 outfield = H2;
1874 }
1875 break;
1876
1877 case 2:
1878 {
1879 outfield = E3;
1880 }
1881 break;
1882
1883 case 3:
1884 {
1885 outfield = F1;
1886 }
1887 break;
1888
1889 case 4:
1890 {
1891 outfield = F2;
1892 }
1893 break;
1894
1895 default:
1896 break;
1897 }
1898
1899 return outfield;
1900}
1901
1902void MMFMaxwell::Printout_SurfaceCurrent(
1903 Array<OneD, Array<OneD, NekDouble>> &fields, const int time)
1904{
1905 int nq = m_fields[0]->GetTotPoints();
1906 int nTraceNumPoints = GetTraceTotPoints();
1907
1908 int totbdryexp =
1909 m_fields[0]->GetBndCondExpansions()[m_boundaryforSF]->GetExpSize();
1910 int npts = m_fields[0]
1911 ->GetBndCondExpansions()[m_boundaryforSF]
1912 ->GetExp(0)
1913 ->GetNumPoints(0);
1914 int totnpts = totbdryexp * npts;
1915
1916 Array<OneD, NekDouble> Jphi(totnpts);
1917 Array<OneD, NekDouble> Jrad(totnpts);
1918 Array<OneD, NekDouble> Jcurrent(totnpts);
1919
1920 Array<OneD, NekDouble> phiFwd(nTraceNumPoints);
1921 Array<OneD, NekDouble> radFwd(nTraceNumPoints);
1922
1923 // Compute phiFwd = acos(-x/r) along the trace
1924 Array<OneD, NekDouble> x(nq);
1925 Array<OneD, NekDouble> y(nq);
1926 Array<OneD, NekDouble> z(nq);
1927
1928 m_fields[0]->GetCoords(x, y, z);
1929
1930 Array<OneD, NekDouble> xFwd(nTraceNumPoints);
1931 Array<OneD, NekDouble> yFwd(nTraceNumPoints);
1932
1933 m_fields[0]->ExtractTracePhys(x, xFwd);
1934 m_fields[0]->ExtractTracePhys(y, yFwd);
1935
1936 for (int i = 0; i < nTraceNumPoints; ++i)
1937 {
1938 radFwd[i] = sqrt(xFwd[i] * xFwd[i] / m_MMFfactors[0] / m_MMFfactors[0] +
1939 yFwd[i] * yFwd[i] / m_MMFfactors[1] / m_MMFfactors[1]);
1940 phiFwd[i] = atan2(yFwd[i] / radFwd[i], -1.0 * xFwd[i] / radFwd[i]);
1941 }
1942
1943 Array<OneD, NekDouble> ntimesHFwd(nTraceNumPoints);
1944 ntimesHFwd = ComputeSurfaceCurrent(time, fields);
1945
1946 // The surface for current should be the first boundary
1947 int id2, cnt = 0;
1948 for (int e = 0; e < totbdryexp; ++e)
1949 {
1950 id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
1951 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
1952
1953 Vmath::Vcopy(npts, &phiFwd[id2], 1, &Jphi[e * npts], 1);
1954 Vmath::Vcopy(npts, &radFwd[id2], 1, &Jrad[e * npts], 1);
1955 Vmath::Vcopy(npts, &ntimesHFwd[id2], 1, &Jcurrent[e * npts], 1);
1956 }
1957
1958 // Vmath::Vmul(totnpts, tmpr, 1, tmpr, 1, Jcurrent, 1);
1959 // Vmath::Vvtvp(totnpts, tmpi, 1, tmpi, 1, Jcurrent, 1, Jcurrent, 1);
1960 // Vmath::Vsqrt(totnpts, Jcurrent, 1, Jcurrent, 1);
1961
1962 std::cout << "========================================================"
1963 << std::endl;
1964
1965 std::cout << "phi = " << std::endl;
1966 for (int i = 0; i < totnpts; ++i)
1967 {
1968 std::cout << Jphi[i] << ", ";
1969 }
1970 std::cout << std::endl << std::endl;
1971
1972 std::cout << "J = " << std::endl;
1973 for (int i = 0; i < totnpts; ++i)
1974 {
1975 std::cout << Jcurrent[i] << ", ";
1976 }
1977 std::cout << std::endl << std::endl;
1978}
1979
1980Array<OneD, NekDouble> MMFMaxwell::ComputeSurfaceCurrent(
1981 const int time, const Array<OneD, const Array<OneD, NekDouble>> &fields)
1982{
1983 int nq = m_fields[0]->GetTotPoints();
1984 int nTraceNumPoints = GetTraceTotPoints();
1985
1986 Array<OneD, NekDouble> outfield(nTraceNumPoints, 0.0);
1987
1988 switch (m_PolType)
1989 {
1990 // (n \times H^r)_z = H1r (n \times e^1)_z + H2r (n \times e^2)_z,
1991 case SolverUtils::eTransMagnetic:
1992 {
1993 Array<OneD, NekDouble> tmp(nq);
1994 Array<OneD, NekDouble> tmpFwd(nTraceNumPoints);
1995
1996 for (int i = 0; i < 2; i++)
1997 {
1998 tmp = GetIncidentField(i, time);
1999 Vmath::Vadd(nq, fields[i], 1, tmp, 1, tmp, 1);
2000
2001 m_fields[0]->ExtractTracePhys(tmp, tmpFwd);
2002 Vmath::Vvtvp(nTraceNumPoints, &m_ntimesMFFwd[i][2][0], 1,
2003 &tmpFwd[0], 1, &outfield[0], 1, &outfield[0], 1);
2004 }
2005 }
2006 break;
2007
2008 case SolverUtils::eTransElectric:
2009 {
2010 Array<OneD, NekDouble> tmp(nq);
2011
2012 tmp = GetIncidentField(2, time);
2013 Vmath::Vadd(nq, fields[2], 1, tmp, 1, tmp, 1);
2014 m_fields[0]->ExtractTracePhys(tmp, outfield);
2015 }
2016 break;
2017
2018 default:
2019 break;
2020 }
2021
2022 return outfield;
2023}
2024
2025void MMFMaxwell::GenerateSigmaPML(const NekDouble PMLthickness,
2026 const NekDouble PMLstart,
2027 const NekDouble PMLmaxsigma,
2028 Array<OneD, Array<OneD, NekDouble>> &SigmaPML)
2029{
2030 int nq = m_fields[0]->GetNpoints();
2031
2032 // Construct sigmaX and sigmaY for UPML
2033 Array<OneD, NekDouble> x(nq);
2034 Array<OneD, NekDouble> y(nq);
2035 Array<OneD, NekDouble> z(nq);
2036
2037 m_fields[0]->GetCoords(x, y, z);
2038
2039 // Construction of SigmaPML
2040 SigmaPML = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
2041 for (int j = 0; j < m_shapedim; j++)
2042 {
2043 SigmaPML[j] = Array<OneD, NekDouble>(nq, 0.0);
2044 }
2045
2046 // PML region indicator: [ 0 : curvedPML : RecPML]
2047 // RecPML = [0 : 0 : 1]
2048 // PMLRegion = [0 : 1 : 1], CurvedPML = PMLRegion - RecPML = [0: 1: 0]
2049 Array<OneD, NekDouble> RecPML(nq);
2050 Array<OneD, NekDouble> CurvedPML(nq);
2051 Array<OneD, NekDouble> PMLRegion(nq);
2052 m_fields[0]->GenerateElementVector(m_RecPML, 0.0, 1.0, RecPML);
2053 m_fields[0]->GenerateElementVector(m_PMLelement, 0.0, 1.0, PMLRegion);
2054 Vmath::Vsub(nq, PMLRegion, 1, RecPML, 1, CurvedPML, 1);
2055
2056 switch (m_AddPML)
2057 {
2058 // RecPML only
2059 case 1:
2060 {
2061 // Rectangular PML
2062 NekDouble xRlayer, xLlayer, yRlayer, yLlayer;
2063
2064 xRlayer = Vmath::Vmax(nq, x, 1) - PMLthickness;
2065 xLlayer = Vmath::Vmin(nq, x, 1) + PMLthickness;
2066 yRlayer = Vmath::Vmax(nq, y, 1) - PMLthickness;
2067 yLlayer = Vmath::Vmin(nq, y, 1) + PMLthickness;
2068
2069 NekDouble xd, yd;
2070 for (int i = 0; i < nq; i++)
2071 {
2072 // SimgaPML along e^1
2073 if (x[i] >= xRlayer)
2074 {
2075 xd = (x[i] - xRlayer) / PMLthickness;
2076 }
2077
2078 else if (x[i] <= xLlayer)
2079 {
2080 xd = (xLlayer - x[i]) / PMLthickness;
2081 }
2082
2083 else
2084 {
2085 xd = 0.0;
2086 }
2087
2088 SigmaPML[0][i] = RecPML[i] * PMLmaxsigma * (xd * xd * xd);
2089
2090 // SigmaPML along e^2
2091 if (y[i] >= yRlayer)
2092 {
2093 yd = (y[i] - yRlayer) / PMLthickness;
2094 }
2095
2096 else if (y[i] <= yLlayer)
2097 {
2098 yd = (yLlayer - y[i]) / PMLthickness;
2099 }
2100
2101 else
2102 {
2103 yd = 0.0;
2104 }
2105
2106 SigmaPML[1][i] = PMLRegion[i] * PMLmaxsigma * (yd * yd * yd);
2107 }
2108 }
2109 break;
2110
2111 // CurvedPML only
2112 case 2:
2113 {
2114 // Curved PML
2115 NekDouble relrad, rad;
2116 for (int i = 0; i < nq; i++)
2117 {
2118 rad = sqrt(x[i] * x[i] / m_MMFfactors[0] / m_MMFfactors[0] +
2119 y[i] * y[i] / m_MMFfactors[1] / m_MMFfactors[1]);
2120
2121 if (rad >= PMLstart)
2122 {
2123 relrad = (rad - PMLstart) / PMLthickness;
2124 SigmaPML[1][i] =
2125 PMLRegion[i] * PMLmaxsigma * pow(relrad, m_PMLorder);
2126 }
2127 }
2128 }
2129 break;
2130
2131 // Slanted PML
2132 case 3:
2133 {
2134 NekDouble relrad, radon, radtw, radth, radfo;
2135 for (int i = 0; i < nq; i++)
2136 {
2137 radon = -1.0 * x[i] + y[i] - 7;
2138 radtw = x[i] + y[i] - 7;
2139 radth = -x[i] - y[i] - 7;
2140 radfo = x[i] - y[i] - 7;
2141
2142 if (radon >= 0.0)
2143 {
2144 relrad = radon / PMLthickness;
2145 SigmaPML[1][i] =
2146 PMLRegion[i] * PMLmaxsigma * pow(relrad, m_PMLorder);
2147 }
2148
2149 if (radtw >= 0.0)
2150 {
2151 relrad = radtw / PMLthickness;
2152 SigmaPML[0][i] =
2153 PMLRegion[i] * PMLmaxsigma * pow(relrad, m_PMLorder);
2154 }
2155
2156 if (radth >= 0.0)
2157 {
2158 relrad = radth / PMLthickness;
2159 SigmaPML[0][i] =
2160 PMLRegion[i] * PMLmaxsigma * pow(relrad, m_PMLorder);
2161 }
2162
2163 if (radfo >= 0.0)
2164 {
2165 relrad = radfo / PMLthickness;
2166 SigmaPML[1][i] =
2167 PMLRegion[i] * PMLmaxsigma * pow(relrad, m_PMLorder);
2168 }
2169 }
2170 }
2171 break;
2172 }
2173
2174 std::cout << "*** sigma1 = [ " << Vmath::Vmin(nq, &SigmaPML[0][0], 1)
2175 << " , " << Vmath::Vmax(nq, &SigmaPML[0][0], 1)
2176 << " ] , sigma2 = [ " << Vmath::Vmin(nq, &SigmaPML[1][0], 1)
2177 << " , " << Vmath::Vmax(nq, &SigmaPML[1][0], 1) << " ] "
2178 << std::endl;
2179}
2180
2181NekDouble MMFMaxwell::ComputeEnergyDensity(
2182 Array<OneD, Array<OneD, NekDouble>> &fields)
2183{
2184 int nq = GetTotPoints();
2185 NekDouble energy;
2186
2187 Array<OneD, NekDouble> tmp(nq, 0.0);
2188
2189 for (int i = 0; i < 3; ++i)
2190 {
2191 Vmath::Vvtvp(nq, &fields[i][0], 1, &fields[i][0], 1, &tmp[0], 1,
2192 &tmp[0], 1);
2193 }
2194
2195 energy = 0.5 * (m_fields[0]->Integral(tmp));
2196 return energy;
2197}
2198
2199void MMFMaxwell::ComputeMaterialVector(
2200 Array<OneD, Array<OneD, NekDouble>> &epsvec,
2201 Array<OneD, Array<OneD, NekDouble>> &muvec)
2202{
2203 switch (m_TestMaxwellType)
2204 {
2205 case SolverUtils::eMaxwell1D:
2206 {
2207 m_fields[0]->GenerateElementVector(m_ElemtGroup1, m_varepsilon[0],
2208 m_varepsilon[1], epsvec[0]);
2209 m_fields[0]->GenerateElementVector(m_ElemtGroup1, 1.0, 1.0,
2210 muvec[0]);
2211 }
2212 break;
2213
2214 case SolverUtils::eTestMaxwell2DPEC:
2215 case SolverUtils::eTestMaxwell2DPECAVGFLUX:
2216 case SolverUtils::eTestMaxwell2DPMC:
2217 case SolverUtils::eScatField2D:
2218 case SolverUtils::eTotField2D:
2219 {
2220 switch (m_PolType)
2221 {
2222 case SolverUtils::eTransMagnetic:
2223 {
2224 m_fields[0]->GenerateElementVector(m_ElemtGroup1, m_mu[0],
2225 1.0, muvec[0]);
2226 m_fields[0]->GenerateElementVector(m_ElemtGroup1, m_mu[1],
2227 1.0, muvec[1]);
2228 m_fields[0]->GenerateElementVector(
2229 m_ElemtGroup1, m_varepsilon[2], 1.0, epsvec[2]);
2230
2231 // // ONLY FOR VARIABLE ANISOTROPY TEST
2232 // int nq = GetTotPoints();
2233
2234 // Array<OneD, NekDouble> tmpIN(nq);
2235
2236 // m_fields[0]->GenerateElementVector(m_ElemtGroup1, 1.0,
2237 // 0.0, tmpIN);
2238
2239 // Array<OneD, NekDouble> x0(nq);
2240 // Array<OneD, NekDouble> x1(nq);
2241 // Array<OneD, NekDouble> x2(nq);
2242
2243 // m_fields[0]->GetCoords(x0,x1,x2);
2244
2245 // for (int i=0; i<nq; i++)
2246 // {
2247 // muvec[1][i] = tmpIN[i]*(1.0 - 2*x0[i]*x0[i]) +
2248 // (1.0-tmpIN[i]);
2249 // }
2250 }
2251 break;
2252
2253 case SolverUtils::eTransElectric:
2254 {
2255 m_fields[0]->GenerateElementVector(
2256 m_ElemtGroup1, m_varepsilon[0], 1.0, epsvec[0]);
2257 m_fields[0]->GenerateElementVector(
2258 m_ElemtGroup1, m_varepsilon[1], 1.0, epsvec[1]);
2259 m_fields[0]->GenerateElementVector(m_ElemtGroup1, m_mu[2],
2260 1.0, muvec[2]);
2261
2262 // // ONLY FOR VARIABLE ANISOTROPY TEST
2263 // int nq = GetTotPoints();
2264
2265 // Array<OneD, NekDouble> tmpIN(nq);
2266
2267 // m_fields[0]->GenerateElementVector(m_ElemtGroup1, 1.0,
2268 // 0.0, tmpIN);
2269
2270 // Array<OneD, NekDouble> x0(nq);
2271 // Array<OneD, NekDouble> x1(nq);
2272 // Array<OneD, NekDouble> x2(nq);
2273
2274 // m_fields[0]->GetCoords(x0,x1,x2);
2275
2276 // for (int i=0; i<nq; i++)
2277 // {
2278 // epsvec[1][i] = tmpIN[i]*(1.0 - 2*x0[i]*x0[i]) +
2279 // (1.0-tmpIN[i]);
2280 // }
2281 }
2282 break;
2283
2284 default:
2285 break; // Pol
2286 }
2287 }
2288
2289 default:
2290 break; // TestType
2291 }
2292}
2293
2294void MMFMaxwell::ComputeMaterialOpticalCloak(
2295 const Array<OneD, const NekDouble> &radvec,
2296 Array<OneD, Array<OneD, NekDouble>> &epsvec,
2297 [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &muvec,
2298 [[maybe_unused]] const bool Dispersion)
2299{
2300 int nq = GetNpoints();
2301
2302 // Cloaking metamaterial
2303 // \varepsilon_\theta = (b/(b-a))^2
2304 // \varepsilon_r = (b/(b-a))^2 ((r-a)/r)^2
2305 // \mu_z = 1.0m_CloakingOn
2306
2307 NekDouble m_b = 1.0 / 0.314;
2308 NekDouble m_a = 1.0;
2309
2310 NekDouble m_adel = m_a - m_Cloakraddelta;
2311
2312 NekDouble boveradel = m_b / (m_b - m_adel);
2313
2314 Array<OneD, NekDouble> Cloakregion(nq, 0.0);
2315
2316 NekDouble la = m_MMFfactors[0];
2317 NekDouble lb = m_MMFfactors[1];
2318
2319 NekDouble ExactCloakArea;
2320 if (fabs(la * lb - 1.0) < 0.00001)
2321 {
2322 ExactCloakArea = m_pi * (m_b * m_b - m_a * m_a);
2323 }
2324
2325 else
2326 {
2327 ExactCloakArea = m_pi * (3.0 * lb * 3.0 * la - lb * la);
2328 }
2329
2330 m_fields[0]->GenerateElementVector(m_ElemtGroup1, 1.0, 0.0, Cloakregion);
2331
2332 ExactCloakArea = ExactCloakArea - (m_fields[0]->Integral(Cloakregion));
2333 std::cout << "*** Error of Cloakregion area = " << ExactCloakArea
2334 << std::endl;
2335
2336 NekDouble ratio;
2337
2338 epsvec[0] = Array<OneD, NekDouble>(nq, 1.0);
2339 epsvec[1] = Array<OneD, NekDouble>(nq, 1.0);
2340 m_wp2 = Array<OneD, NekDouble>(nq, 0.0);
2341 for (int i = 0; i < nq; ++i)
2342 {
2343 if (Cloakregion[i] > 0)
2344 {
2345 // relrad = m_a +
2346 // (m_b-m_adel)*(radvec[i]-m_adel)/(Cloakradmax-m_adel);
2347 ratio = (radvec[i] - m_adel) / radvec[i];
2348
2349 epsvec[0][i] = boveradel * boveradel;
2350 if (m_DispersiveCloak)
2351 {
2352 epsvec[1][i] = 1.0;
2353 m_wp2[i] = m_Incfreq * m_Incfreq *
2354 (1.0 - boveradel * boveradel * ratio * ratio) +
2355 m_wp2Tol;
2356 }
2357
2358 else
2359 {
2360 epsvec[1][i] = boveradel * boveradel * (ratio * ratio);
2361 }
2362 }
2363 }
2364}
2365
2366void MMFMaxwell::ComputeMaterialMicroWaveCloak(
2367 const Array<OneD, const NekDouble> &radvec,
2368 Array<OneD, Array<OneD, NekDouble>> &epsvec,
2369 Array<OneD, Array<OneD, NekDouble>> &muvec)
2370{
2371 int nq = GetNpoints();
2372
2373 NekDouble m_b = 2.67;
2374 NekDouble m_a = 1.33;
2375 NekDouble m_adel;
2376
2377 m_adel = m_a - m_Cloakraddelta;
2378
2379 Array<OneD, NekDouble> Cloakregion(nq, 0.0);
2380 NekDouble ExactCloakArea = m_pi * (m_b * m_b - m_a * m_a);
2381 m_fields[0]->GenerateElementVector(m_ElemtGroup1, 1.0, 0.0, Cloakregion);
2382
2383 if (m_ElemtGroup0 > 0)
2384 {
2385 Array<OneD, NekDouble> Vacregion(nq, 0.0);
2386 m_fields[0]->GenerateElementVector(m_ElemtGroup0, 1.0, 0.0, Vacregion);
2387
2388 Vmath::Vsub(nq, Cloakregion, 1, Vacregion, 1, Cloakregion, 1);
2389 }
2390
2391 ExactCloakArea = ExactCloakArea - (m_fields[0]->Integral(Cloakregion));
2392 std::cout << "*** Error of Cloakregion area = " << ExactCloakArea
2393 << std::endl;
2394
2395 epsvec[0] = Array<OneD, NekDouble>(nq, 1.0);
2396 epsvec[1] = Array<OneD, NekDouble>(nq, 1.0);
2397
2398 muvec[0] = Array<OneD, NekDouble>(nq, 1.0);
2399 muvec[1] = Array<OneD, NekDouble>(nq, 1.0);
2400 for (int i = 0; i < nq; ++i)
2401 {
2402 if (Cloakregion[i] > 0)
2403 {
2404 // relrad = m_a +
2405 // (m_b-m_a)*(radvec[i]-Cloakradmin)/(Cloakradmax-Cloakradmin);
2406 // ratio = (relrad - m_a + m_Cloakraddelta)/relrad;
2407
2408 epsvec[0][i] = radvec[i] / (radvec[i] - m_adel);
2409 epsvec[1][i] = (radvec[i] - m_adel) / radvec[i];
2410 muvec[2][i] = (m_b / (m_b - m_adel)) * (m_b / (m_b - m_adel)) *
2411 (radvec[i] - m_adel) / radvec[i];
2412
2413 muvec[0][i] = epsvec[0][i];
2414 muvec[1][i] = epsvec[1][i];
2415 epsvec[2][i] = muvec[2][i];
2416 }
2417 }
2418}
2419
2420void MMFMaxwell::AddPML(
2421 const Array<OneD, const Array<OneD, NekDouble>> &physarray,
2422 Array<OneD, Array<OneD, NekDouble>> &outarray)
2423{
2424 int nq = m_fields[0]->GetTotPoints();
2425 Array<OneD, NekDouble> tmp(nq);
2426
2427 Array<OneD, NekDouble> Sigma0plus1Neg(nq);
2428 Array<OneD, NekDouble> Sigma0minus1(nq);
2429 Array<OneD, NekDouble> Sigma1minus0(nq);
2430
2431 Vmath::Vsub(nq, &m_SigmaPML[1][0], 1, &m_SigmaPML[0][0], 1,
2432 &Sigma1minus0[0], 1);
2433 Vmath::Vsub(nq, &m_SigmaPML[0][0], 1, &m_SigmaPML[1][0], 1,
2434 &Sigma0minus1[0], 1);
2435 Vmath::Vadd(nq, &m_SigmaPML[0][0], 1, &m_SigmaPML[1][0], 1,
2436 &Sigma0plus1Neg[0], 1);
2437 Vmath::Neg(nq, Sigma0plus1Neg, 1);
2438
2439 switch (m_PolType)
2440 {
2441 case SolverUtils::eTransMagnetic:
2442 {
2443 int indxH0 = 0;
2444 int indxH1 = 1;
2445 int indxE2 = 2;
2446 int indxQ0 = 3;
2447 int indxQ1 = 4;
2448 int indxP2 = 5;
2449
2450 // dH0/dt: Add (sigma_0 - \sigma_1) H0 + Q0
2451 Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxH0][0], 1,
2452 &outarray[indxH0][0], 1, &outarray[indxH0][0], 1);
2453 Vmath::Vadd(nq, &physarray[indxQ0][0], 1, &outarray[indxH0][0], 1,
2454 &outarray[indxH0][0], 1);
2455
2456 // dH1/dt: Add (sigma_1 - \sigma_0) H1 + Q1
2457 Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxH1][0], 1,
2458 &outarray[indxH1][0], 1, &outarray[indxH1][0], 1);
2459 Vmath::Vadd(nq, &physarray[indxQ1][0], 1, &outarray[indxH1][0], 1,
2460 &outarray[indxH1][0], 1);
2461
2462 // dHz/dt: Add -(\sigma_0 + \sigma_1) Ez + Pz
2463 Vmath::Vvtvp(nq, &Sigma0plus1Neg[0], 1, &physarray[indxE2][0], 1,
2464 &outarray[indxE2][0], 1, &outarray[indxE2][0], 1);
2465 Vmath::Vadd(nq, &physarray[indxP2][0], 1, &outarray[indxE2][0], 1,
2466 &outarray[indxE2][0], 1);
2467
2468 // dQ0/dt: Assign -\sigma_0 * ( Q0 + (\sigma_0 - \sigma_1) * H0 )
2469 Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxH0][0], 1,
2470 &physarray[indxQ0][0], 1, &outarray[indxQ0][0], 1);
2471 Vmath::Vmul(nq, &m_SigmaPML[0][0], 1, &outarray[indxQ0][0], 1,
2472 &outarray[indxQ0][0], 1);
2473 Vmath::Neg(nq, &outarray[indxQ0][0], 1);
2474
2475 // dQ1/dt: Assign -\sigma_1 * ( Q1 + (\sigma_1 - \sigma_0) * H1 )
2476 Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxH1][0], 1,
2477 &physarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2478 Vmath::Vmul(nq, &m_SigmaPML[1][0], 1, &outarray[indxQ1][0], 1,
2479 &outarray[indxQ1][0], 1);
2480 Vmath::Neg(nq, &outarray[indxQ1][0], 1);
2481
2482 if (m_DispersiveCloak)
2483 {
2484 Vmath::Vvtvp(nq, &m_wp2[0], 1, &physarray[indxH1][0], 1,
2485 &outarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2486 }
2487
2488 // dP3/dt: Assign - \sigma_1 * \sigma_2 * E_z
2489 Vmath::Vmul(nq, &m_SigmaPML[0][0], 1, &m_SigmaPML[1][0], 1,
2490 &outarray[indxP2][0], 1);
2491 Vmath::Vmul(nq, &physarray[indxE2][0], 1, &outarray[indxP2][0], 1,
2492 &outarray[indxP2][0], 1);
2493 Vmath::Neg(nq, &outarray[indxP2][0], 1);
2494 }
2495 break;
2496
2497 case SolverUtils::eTransElectric:
2498 {
2499 int indxE0 = 0;
2500 int indxE1 = 1;
2501 int indxH2 = 2;
2502 int indxQ0 = 3;
2503 int indxQ1 = 4;
2504 int indxP2 = 5;
2505
2506 // dE0/dt: Add (sigma_0 - \sigma_1) E0 - Q0
2507 Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxE0][0], 1,
2508 &outarray[indxE0][0], 1, &outarray[indxE0][0], 1);
2509 Vmath::Vsub(nq, &outarray[indxE0][0], 1, &physarray[indxQ0][0], 1,
2510 &outarray[indxE0][0], 1);
2511
2512 // dE1/dt: Add (sigma_1 - \sigma_0) E1 - Q1
2513 Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxE1][0], 1,
2514 &outarray[indxE1][0], 1, &outarray[indxE1][0], 1);
2515 Vmath::Vsub(nq, &outarray[indxE1][0], 1, &physarray[indxQ1][0], 1,
2516 &outarray[indxE1][0], 1);
2517
2518 // dHz/dt: Add -(\sigma_0 + \sigma_1) Hz - Pz
2519 Vmath::Vvtvp(nq, &Sigma0plus1Neg[0], 1, &physarray[indxH2][0], 1,
2520 &outarray[indxH2][0], 1, &outarray[indxH2][0], 1);
2521 Vmath::Vsub(nq, &outarray[indxH2][0], 1, &physarray[indxP2][0], 1,
2522 &outarray[indxH2][0], 1);
2523
2524 // dQ0/dt: Assign -\sigma_0 * ( Q0 + (\sigma_1 - \sigma_0) * E0 )
2525 Vmath::Vvtvp(nq, &Sigma1minus0[0], 1, &physarray[indxE0][0], 1,
2526 &physarray[indxQ0][0], 1, &outarray[indxQ0][0], 1);
2527 Vmath::Vmul(nq, &m_SigmaPML[0][0], 1, &outarray[indxQ0][0], 1,
2528 &outarray[indxQ0][0], 1);
2529 Vmath::Neg(nq, &outarray[indxQ0][0], 1);
2530
2531 // dQ1/dt: Assign -\sigma_1 * ( Q1 + (\sigma_0 - \sigma_1) * E1 )
2532 Vmath::Vvtvp(nq, &Sigma0minus1[0], 1, &physarray[indxE1][0], 1,
2533 &physarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2534 Vmath::Vmul(nq, &m_SigmaPML[1][0], 1, &outarray[indxQ1][0], 1,
2535 &outarray[indxQ1][0], 1);
2536 Vmath::Neg(nq, &outarray[indxQ1][0], 1);
2537
2538 if (m_DispersiveCloak)
2539 {
2540 Vmath::Vvtvp(nq, &m_wp2[0], 1, &physarray[indxE1][0], 1,
2541 &outarray[indxQ1][0], 1, &outarray[indxQ1][0], 1);
2542 }
2543
2544 // dP3/dt: Assign \sigma_1 * \sigma_2 * H_z
2545 Vmath::Vmul(nq, &m_SigmaPML[0][0], 1, &m_SigmaPML[1][0], 1,
2546 &outarray[indxP2][0], 1);
2547 Vmath::Vmul(nq, &physarray[indxH2][0], 1, &outarray[indxP2][0], 1,
2548 &outarray[indxP2][0], 1);
2549 }
2550 break;
2551
2552 default:
2553 break;
2554 }
2555}
2556
2557void MMFMaxwell::Checkpoint_TotalFieldOutput(
2558 const int n, const NekDouble time,
2559 const Array<OneD, const Array<OneD, NekDouble>> &fieldphys)
2560{
2561 int nvar = m_fields.size();
2562 int nq = m_fields[0]->GetTotPoints();
2563 int ncoeffs = m_fields[0]->GetNcoeffs();
2564
2565 std::string outname = m_sessionName + "Tot_" + std::to_string(n) + ".chk";
2566
2567 std::vector<std::string> variables(nvar);
2568 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2569
2570 for (int i = 0; i < nvar; ++i)
2571 {
2572 fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2573 }
2574
2575 Array<OneD, NekDouble> totfield(nq);
2576 for (int i = 0; i < nvar; ++i)
2577 {
2578 totfield = GetIncidentField(i, time);
2579 Vmath::Vadd(nq, fieldphys[i], 1, totfield, 1, totfield, 1);
2580
2581 m_fields[i]->FwdTrans(totfield, fieldcoeffs[i]);
2582 variables[i] = m_boundaryConditions->GetVariable(i);
2583 }
2584
2585 WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2586}
2587
2588// Only for nvar = 3
2589void MMFMaxwell::Checkpoint_PlotOutput(
2590 const int n, const Array<OneD, const Array<OneD, NekDouble>> &fieldphys)
2591{
2592 int nvar = m_fields.size();
2593 int nq = m_fields[0]->GetTotPoints();
2594 int ncoeffs = m_fields[0]->GetNcoeffs();
2595
2596 std::string outname = m_sessionName + "Plot_" + std::to_string(n) + ".chk";
2597
2598 std::vector<std::string> variables(nvar);
2599 variables[0] = "Fx";
2600 variables[1] = "Fy";
2601 variables[2] = "Fz";
2602
2603 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2604 for (int i = 0; i < nvar; ++i)
2605 {
2606 fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2607 }
2608
2609 Array<OneD, NekDouble> tmp(nq);
2610 for (int k = 0; k < m_spacedim; ++k)
2611 {
2612 Vmath::Vmul(nq, &fieldphys[0][0], 1, &m_movingframes[0][k * nq], 1,
2613 &tmp[0], 1);
2614 Vmath::Vvtvp(nq, &fieldphys[1][0], 1, &m_movingframes[1][k * nq], 1,
2615 &tmp[0], 1, &tmp[0], 1);
2616
2617 m_fields[k]->FwdTrans(tmp, fieldcoeffs[k]);
2618 }
2619
2620 WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2621}
2622
2623void MMFMaxwell::Checkpoint_TotPlotOutput(
2624 const int n, const NekDouble time,
2625 const Array<OneD, const Array<OneD, NekDouble>> &fieldphys)
2626{
2627 int nvar = m_fields.size();
2628 int nq = m_fields[0]->GetTotPoints();
2629 int ncoeffs = m_fields[0]->GetNcoeffs();
2630
2631 std::string outname =
2632 m_sessionName + "TotPlot_" + std::to_string(n) + ".chk";
2633
2634 std::vector<std::string> variables(nvar);
2635 variables[0] = "Frx";
2636 variables[1] = "Fry";
2637 variables[2] = "Frz";
2638 for (int i = 3; i < nvar; ++i)
2639 {
2640 variables[i] = m_boundaryConditions->GetVariable(i);
2641 }
2642
2643 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2644 for (int i = 0; i < nvar; ++i)
2645 {
2646 fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2647 }
2648
2649 Array<OneD, NekDouble> tmp(nq);
2650 Array<OneD, NekDouble> totfield0(nq);
2651 Array<OneD, NekDouble> totfield1(nq);
2652
2653 totfield0 = GetIncidentField(0, time);
2654 Vmath::Vadd(nq, fieldphys[0], 1, totfield0, 1, totfield0, 1);
2655
2656 totfield1 = GetIncidentField(1, time);
2657 Vmath::Vadd(nq, fieldphys[1], 1, totfield1, 1, totfield1, 1);
2658
2659 for (int k = 0; k < m_spacedim; ++k)
2660 {
2661 Vmath::Vmul(nq, &totfield0[0], 1, &m_movingframes[0][k * nq], 1,
2662 &tmp[0], 1);
2663 Vmath::Vvtvp(nq, &totfield1[0], 1, &m_movingframes[1][k * nq], 1,
2664 &tmp[0], 1, &tmp[0], 1);
2665
2666 m_fields[k]->FwdTrans(tmp, fieldcoeffs[k]);
2667 }
2668
2669 for (int j = 3; j < nvar; ++j)
2670 {
2671 m_fields[j]->FwdTrans(fieldphys[j], fieldcoeffs[j]);
2672 }
2673
2674 WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2675}
2676
2677void MMFMaxwell::Checkpoint_EDFluxOutput(
2678 const int n, [[maybe_unused]] const NekDouble time,
2679 const Array<OneD, const Array<OneD, NekDouble>> &fieldphys)
2680{
2681 int nvar = m_fields.size();
2682 int nq = m_fields[0]->GetTotPoints();
2683 int ncoeffs = m_fields[0]->GetNcoeffs();
2684
2685 std::string outname =
2686 m_sessionName + "EDFlux_" + std::to_string(n) + ".chk";
2687
2688 std::vector<std::string> variables(nvar);
2689 variables[0] = "EDFx";
2690 variables[1] = "EDFy";
2691 variables[2] = "EDFz";
2692 for (int i = 3; i < nvar; ++i)
2693 {
2694 variables[i] = m_boundaryConditions->GetVariable(i);
2695 }
2696
2697 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2698 for (int i = 0; i < nvar; ++i)
2699 {
2700 fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2701 }
2702
2703 Array<OneD, NekDouble> tmp(nq);
2704
2705 // TE: H^3 (E^2 e^1 - E^1 e^2 )
2706 // TM: -E^3 (H^2 e^1 - H^1 e^2 )
2707 for (int k = 0; k < m_spacedim; ++k)
2708 {
2709 Vmath::Vmul(nq, &fieldphys[0][0], 1, &m_movingframes[1][k * nq], 1,
2710 &tmp[0], 1);
2711 Vmath::Vvtvm(nq, &fieldphys[1][0], 1, &m_movingframes[0][k * nq], 1,
2712 &tmp[0], 1, &tmp[0], 1);
2713
2714 Vmath::Vmul(nq, &fieldphys[2][0], 1, &tmp[0], 1, &tmp[0], 1);
2715
2716 if (m_PolType == SolverUtils::eTransMagnetic)
2717 {
2718 Vmath::Neg(nq, tmp, 1);
2719 }
2720
2721 m_fields[k]->FwdTrans(tmp, fieldcoeffs[k]);
2722 }
2723
2724 WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2725}
2726
2727void MMFMaxwell::Checkpoint_EnergyOutput(
2728 const int n, [[maybe_unused]] const NekDouble time,
2729 const Array<OneD, const Array<OneD, NekDouble>> &fieldphys)
2730{
2731 int nvar = m_fields.size();
2732 int nq = m_fields[0]->GetTotPoints();
2733 int ncoeffs = m_fields[0]->GetNcoeffs();
2734
2735 std::string outname =
2736 m_sessionName + "Energy_" + std::to_string(n) + ".chk";
2737
2738 std::vector<std::string> variables(nvar);
2739 variables[0] = "Energy";
2740 variables[1] = "EnergyFlux";
2741 variables[2] = "Zero";
2742 for (int i = 3; i < nvar; ++i)
2743 {
2744 variables[i] = m_boundaryConditions->GetVariable(i);
2745 }
2746
2747 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvar);
2748 for (int i = 0; i < nvar; ++i)
2749 {
2750 fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2751 }
2752
2753 // Energy = 0.5 * ( E^2 + H^2 )
2754 Array<OneD, NekDouble> energy(nq, 0.0);
2755 Array<OneD, NekDouble> totfield(nq);
2756 for (int k = 0; k < m_spacedim; ++k)
2757 {
2758 // totfield = GetIncidentField(k,time);
2759 // Vmath::Vadd(nq, &fieldphys[k][0], 1, &totfield[0], 1, &totfield[0],
2760 // 1);
2761
2762 Vmath::Vvtvp(nq, &fieldphys[k][0], 1, &fieldphys[k][0], 1, &energy[0],
2763 1, &energy[0], 1);
2764 }
2765 Vmath::Smul(nq, 0.5, energy, 1, energy, 1);
2766 m_fields[0]->FwdTrans(energy, fieldcoeffs[0]);
2767
2768 // EnergyFlux = F3* sqrt( F1^2 + F2^2 )
2769 Array<OneD, NekDouble> energyflux(nq, 0.0);
2770 Array<OneD, NekDouble> Zero(nq, 0.0);
2771 for (int k = 0; k < 2; ++k)
2772 {
2773 // totfield = GetIncidentField(k,time);
2774 // Vmath::Vadd(nq, &fieldphys[k][0], 1, &totfield[0], 1, &totfield[0],
2775 // 1);
2776
2777 Vmath::Vvtvp(nq, &fieldphys[k][0], 1, &fieldphys[k][0], 1,
2778 &energyflux[0], 1, &energyflux[0], 1);
2779 }
2780
2781 Vmath::Vsqrt(nq, energyflux, 1, energyflux, 1);
2782 Vmath::Vmul(nq, &fieldphys[2][0], 1, &energyflux[0], 1, &energyflux[0], 1);
2783
2784 m_fields[1]->FwdTrans(energyflux, fieldcoeffs[1]);
2785 m_fields[2]->FwdTrans(Zero, fieldcoeffs[2]);
2786
2787 WriteFld(outname, m_fields[0], fieldcoeffs, variables);
2788}
2789
2790Array<OneD, NekDouble> MMFMaxwell::GaussianPulse(const NekDouble time,
2791 const NekDouble Psx,
2792 const NekDouble Psy,
2793 const NekDouble Psz,
2794 const NekDouble Gaussianradius)
2795{
2796 int nq = m_fields[0]->GetTotPoints();
2797 int ncoeffs = m_fields[0]->GetNcoeffs();
2798
2799 Array<OneD, NekDouble> x(nq);
2800 Array<OneD, NekDouble> y(nq);
2801 Array<OneD, NekDouble> z(nq);
2802
2803 m_fields[0]->GetCoords(x, y, z);
2804
2805 Array<OneD, NekDouble> outarray(nq, 0.0);
2806 Array<OneD, NekDouble> tmpc(ncoeffs);
2807 NekDouble rad;
2808
2809 NekDouble SFradius = m_PSduration * 0.1;
2810 NekDouble SmoothFactor =
2811 1.0 / (1.0 + exp(-0.5 * (time - m_PSduration) / SFradius));
2812
2813 for (int j = 0; j < nq; ++j)
2814 {
2815 rad = sqrt((x[j] - Psx) * (x[j] - Psx) + (y[j] - Psy) * (y[j] - Psy) +
2816 (z[j] - Psz) * (z[j] - Psz));
2817 outarray[j] = SmoothFactor * exp(-1.0 * (rad / Gaussianradius) *
2818 (rad / Gaussianradius));
2819 }
2820
2821 m_fields[0]->FwdTransLocalElmt(outarray, tmpc);
2822 m_fields[0]->BwdTrans(tmpc, outarray);
2823
2824 return outarray;
2825}
2826
2827Array<OneD, NekDouble> MMFMaxwell::EvaluateCoriolis()
2828{
2829 int nq = GetTotPoints();
2830
2831 NekDouble m_Omega = 1.5486 * 0.000001;
2832
2833 NekDouble x0j, x1j, x2j;
2834 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2835
2836 Array<OneD, NekDouble> x(nq);
2837 Array<OneD, NekDouble> y(nq);
2838 Array<OneD, NekDouble> z(nq);
2839
2840 m_fields[0]->GetCoords(x, y, z);
2841
2842 Array<OneD, NekDouble> outarray(nq, 0.0);
2843 for (int j = 0; j < nq; ++j)
2844 {
2845 x0j = x[j];
2846 x1j = y[j];
2847 x2j = z[j];
2848
2849 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2850 cos_theta);
2851
2852 outarray[j] = 2.0 * m_Omega * sin_theta;
2853 }
2854
2855 return outarray;
2856}
2857
2858void MMFMaxwell::AddCoriolis(Array<OneD, Array<OneD, NekDouble>> &physarray,
2859 Array<OneD, Array<OneD, NekDouble>> &outarray)
2860{
2861 int nq = physarray[0].size();
2862
2863 Array<OneD, NekDouble> tmp(nq);
2864
2865 int indx;
2866 for (int j = 0; j < m_shapedim; ++j)
2867 {
2868 if (j == 0)
2869 {
2870 indx = 2;
2871 }
2872
2873 else if (j == 1)
2874 {
2875 indx = 1;
2876 }
2877
2878 Vmath::Vmul(nq, m_coriolis, 1, physarray[indx], 1, tmp, 1);
2879
2880 switch (m_PolType)
2881 {
2882 case SolverUtils::eTransMagnetic:
2883 {
2884 Vmath::Vmul(nq, m_muvec[indx], 1, tmp, 1, tmp, 1);
2885 }
2886 break;
2887
2888 case SolverUtils::eTransElectric:
2889 {
2890 Vmath::Vmul(nq, m_epsvec[indx], 1, tmp, 1, tmp, 1);
2891 }
2892 break;
2893
2894 default:
2895 break;
2896 }
2897
2898 if (j == 1)
2899 {
2900 Vmath::Neg(nq, tmp, 1);
2901 }
2902
2903 Vmath::Vadd(nq, tmp, 1, outarray[j], 1, outarray[j], 1);
2904 }
2905}
2906
2907Array<OneD, NekDouble> MMFMaxwell::ComputeRadCloak(const int CloakNlayer)
2908{
2909 int nq = GetNpoints();
2910
2911 NekDouble m_b = 1.0 / 0.314;
2912 NekDouble m_a = 1.0;
2913
2914 Array<OneD, int> Layer(CloakNlayer);
2915
2916 NekDouble la = m_MMFfactors[0];
2917 NekDouble lb = m_MMFfactors[1];
2918 if (CloakNlayer == 8)
2919 {
2920 // Circular 8 layer
2921 if (fabs(la * lb - 1.0) < 0.0001)
2922 {
2923 Layer[0] = 119;
2924 Layer[1] = 239;
2925 Layer[2] = 323;
2926 Layer[3] = 459;
2927 Layer[4] = 575;
2928 Layer[5] = 727;
2929 Layer[6] = 891;
2930 Layer[7] = 1059;
2931 }
2932
2933 // Ellipsoid 8 layer
2934 else
2935 {
2936 Layer[0] = 115;
2937 Layer[1] = 219;
2938 Layer[2] = 335;
2939 Layer[3] = 459;
2940 Layer[4] = 607;
2941 Layer[5] = 767;
2942 Layer[6] = 951;
2943 Layer[7] = 1155;
2944 }
2945 }
2946
2947 if (CloakNlayer == 16)
2948 {
2949 // Circular 16 layer
2950 if (fabs(la * lb - 1.0) < 0.0001)
2951 {
2952 Layer[0] = 43;
2953 Layer[1] = 87;
2954 Layer[2] = 135;
2955 Layer[3] = 187;
2956 Layer[4] = 239;
2957 Layer[5] = 295;
2958 Layer[6] = 355;
2959 Layer[7] = 415;
2960 Layer[8] = 479;
2961 Layer[9] = 555;
2962 Layer[10] = 639;
2963 Layer[11] = 727;
2964 Layer[12] = 823;
2965 Layer[13] = 927;
2966 Layer[14] = 1039;
2967 Layer[15] = 1159;
2968 }
2969
2970 // Ellipsoid 8 layer
2971 else
2972 {
2973 Layer[0] = 135;
2974 Layer[1] = 259;
2975 Layer[2] = 387;
2976 Layer[3] = 523;
2977 Layer[4] = 667;
2978 Layer[5] = 803;
2979 Layer[6] = 939;
2980 Layer[7] = 1083;
2981 Layer[8] = 1235;
2982 Layer[9] = 1387;
2983 Layer[10] = 1539;
2984 Layer[11] = 1699;
2985 Layer[12] = 1867;
2986 Layer[13] = 2035;
2987 Layer[14] = 2203;
2988 Layer[15] = 2379;
2989 }
2990 }
2991
2992 Array<OneD, NekDouble> x0(nq);
2993 Array<OneD, NekDouble> x1(nq);
2994 Array<OneD, NekDouble> x2(nq);
2995
2996 m_fields[0]->GetCoords(x0, x1, x2);
2997
2998 Array<OneD, NekDouble> Formertmp(nq, 0.0);
2999 Array<OneD, NekDouble> Currenttmp(nq, 0.0);
3000 Array<OneD, NekDouble> Cloakregion(nq, 0.0);
3001
3002 m_fields[0]->GenerateElementVector(Layer[0], 1.0, 0.0, Currenttmp);
3003 Vmath::Vcopy(nq, Currenttmp, 1, Cloakregion, 1);
3004 for (int i = 1; i < CloakNlayer; ++i)
3005 {
3006 m_fields[0]->GenerateElementVector(Layer[i], 1.0, 0.0, Currenttmp);
3007 m_fields[0]->GenerateElementVector(Layer[i - 1], 1.0, 0.0, Formertmp);
3008
3009 Vmath::Vsub(nq, Currenttmp, 1, Formertmp, 1, Currenttmp, 1);
3010
3011 Vmath::Svtvp(nq, 1.0 * (i + 1), Currenttmp, 1, Cloakregion, 1,
3012 Cloakregion, 1);
3013 }
3014
3015 Array<OneD, NekDouble> radvec(nq);
3016
3017 for (int i = 0; i < nq; ++i)
3018 {
3019 switch (m_MMFdir)
3020 {
3021 case SpatialDomains::eTangentCircular:
3022 {
3023 if ((Cloakregion[i] > 0) && (CloakNlayer > 0))
3024 {
3025 radvec[i] = 1.0 + (m_b - m_a) / CloakNlayer *
3026 (Cloakregion[i] - 0.5);
3027 }
3028
3029 else
3030 {
3031 radvec[i] =
3032 sqrt(x0[i] * x0[i] / la / la + x1[i] * x1[i] / lb / lb);
3033 }
3034 }
3035 break;
3036
3037 case SpatialDomains::eTangentIrregular:
3038 {
3039 radvec[i] = sqrt(2.0 * x0[i] * x0[i] +
3040 x1[i] * x1[i] * x1[i] * x1[i] + x1[i] * x1[i]);
3041 }
3042 break;
3043
3044 case SpatialDomains::eTangentNonconvex:
3045 {
3046 radvec[i] = sqrt(3.0 * x0[i] * x0[i] +
3047 x1[i] * x1[i] * x1[i] * x1[i] - x1[i] * x1[i]);
3048 }
3049 break;
3050
3051 default:
3052 break;
3053 }
3054 }
3055
3056 return radvec;
3057}
3058
3059void MMFMaxwell::v_GenerateSummary(SolverUtils::SummaryList &s)
3060{
3061 MMFSystem::v_GenerateSummary(s);
3062 SolverUtils::AddSummaryItem(
3063 s, "TestMaxwellType",
3064 SolverUtils::TestMaxwellTypeMap[m_TestMaxwellType]);
3065 SolverUtils::AddSummaryItem(s, "PolType",
3066 SolverUtils::PolTypeMap[m_PolType]);
3067 SolverUtils::AddSummaryItem(s, "IncType",
3068 SolverUtils::IncTypeMap[m_IncType]);
3069
3070 if (m_varepsilon[0] * m_varepsilon[1] * m_varepsilon[2] > 1.0)
3071 {
3072 SolverUtils::AddSummaryItem(s, "varepsilon1", m_varepsilon[0]);
3073 SolverUtils::AddSummaryItem(s, "varepsilon2", m_varepsilon[1]);
3074 SolverUtils::AddSummaryItem(s, "varepsilon3", m_varepsilon[2]);
3075 }
3076
3077 if (m_mu[0] * m_mu[1] * m_mu[2] > 1.0)
3078 {
3079 SolverUtils::AddSummaryItem(s, "mu1", m_mu[0]);
3080 SolverUtils::AddSummaryItem(s, "mu2", m_mu[1]);
3081 SolverUtils::AddSummaryItem(s, "mu3", m_mu[2]);
3082 }
3083
3084 if (m_boundaryforSF > 0)
3085 {
3086 SolverUtils::AddSummaryItem(s, "boundarySF", m_boundaryforSF);
3087 }
3088
3089 if (m_ElemtGroup1 > 0)
3090 {
3091 SolverUtils::AddSummaryItem(s, "CloakNlayer", m_CloakNlayer);
3092 SolverUtils::AddSummaryItem(s, "ElemtGroup1", m_ElemtGroup1);
3093 }
3094
3095 SolverUtils::AddSummaryItem(s, "AddRotation", m_AddRotation);
3096
3097 if (m_AddPML > 0)
3098 {
3099 SolverUtils::AddSummaryItem(s, "AddPML", m_AddPML);
3100 SolverUtils::AddSummaryItem(s, "PMLelement", m_PMLelement);
3101 SolverUtils::AddSummaryItem(s, "RecPML", m_RecPML);
3102 SolverUtils::AddSummaryItem(s, "PMLorder", m_PMLorder);
3103 SolverUtils::AddSummaryItem(s, "PMLthickness", m_PMLthickness);
3104 SolverUtils::AddSummaryItem(s, "PMLstart", m_PMLstart);
3105 SolverUtils::AddSummaryItem(s, "PMLmaxsigma", m_PMLmaxsigma);
3106 }
3107
3108 if (m_SourceType)
3109 {
3110 SolverUtils::AddSummaryItem(s, "SourceType",
3111 SourceTypeMap[m_SourceType]);
3112 SolverUtils::AddSummaryItem(s, "Psx", m_Psx);
3113 SolverUtils::AddSummaryItem(s, "Psy", m_Psy);
3114 SolverUtils::AddSummaryItem(s, "Psz", m_Psz);
3115 SolverUtils::AddSummaryItem(s, "PSduration", m_PSduration);
3116 SolverUtils::AddSummaryItem(s, "Gaussianradius", m_Gaussianradius);
3117 }
3118
3119 if (m_CloakType)
3120 {
3121 SolverUtils::AddSummaryItem(s, "CloakType", CloakTypeMap[m_CloakType]);
3122 SolverUtils::AddSummaryItem(s, "DispersiveCloak", m_DispersiveCloak);
3123 SolverUtils::AddSummaryItem(s, "CloakNlayer", m_CloakNlayer);
3124 SolverUtils::AddSummaryItem(s, "Cloakraddelta", m_Cloakraddelta);
3125 }
3126}
3127void MMFMaxwell::print_MMF(Array<OneD, Array<OneD, NekDouble>> &inarray)
3128{
3129 int Ntot = inarray.size();
3130
3131 NekDouble reval = 0.0;
3132 for (int i = 0; i < Ntot; ++i)
3133 {
3134 std::cout << "[" << i << "] = " << inarray[2][i] << std::endl;
3135 // reval = reval + inarray[i]*inarray[i];
3136 }
3137 reval = sqrt(reval / Ntot);
3138}
3139} // namespace Nektar
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
CloakTypeMap
const char *const CloakTypeMap[]
Definition: MMFMaxwell.h:51
SourceType
SourceType
Definition: MMFMaxwell.h:58
ePointSource
@ ePointSource
Definition: MMFMaxwell.h:60
ePlanarSource
@ ePlanarSource
Definition: MMFMaxwell.h:61
SIZE_SourceType
@ SIZE_SourceType
Definition: MMFMaxwell.h:62
SourceTypeMap
const char *const SourceTypeMap[]
Definition: MMFMaxwell.h:65
CloakType
CloakType
Definition: MMFMaxwell.h:42
eOpticalCloak
@ eOpticalCloak
Definition: MMFMaxwell.h:44
eMicroWaveCloak
@ eMicroWaveCloak
Definition: MMFMaxwell.h:47
eOpticalDispersiveCloak
@ eOpticalDispersiveCloak
Definition: MMFMaxwell.h:46
eOpticalConstCloak
@ eOpticalConstCloak
Definition: MMFMaxwell.h:45
SIZE_CloakType
@ SIZE_CloakType
Definition: MMFMaxwell.h:48
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: BasicUtils/NekFactory.hpp:197
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineProjection
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationSchemeOperators.h:92
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineOdeRhs
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationSchemeOperators.h:84
Nektar::LibUtilities::Timer::TimePerTest
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition: Timer.cpp:65
Nektar::MMFMaxwell::Checkpoint_EDFluxOutput
void Checkpoint_EDFluxOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble > > &fieldphys)
Definition: MMFMaxwell.cpp:2677
Nektar::MMFMaxwell::m_CloakNlayer
NekDouble m_CloakNlayer
Definition: MMFMaxwell.h:115
Nektar::MMFMaxwell::m_PSduration
NekDouble m_PSduration
Definition: MMFMaxwell.h:123
Nektar::MMFMaxwell::m_CloakType
CloakType m_CloakType
Definition: MMFMaxwell.h:78
Nektar::MMFMaxwell::m_mu
Array< OneD, NekDouble > m_mu
Definition: MMFMaxwell.h:135
Nektar::MMFMaxwell::m_PMLmaxsigma
NekDouble m_PMLmaxsigma
Definition: MMFMaxwell.h:139
Nektar::MMFMaxwell::AddCoriolis
void AddCoriolis(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFMaxwell.cpp:2858
Nektar::MMFMaxwell::DoOdeRhs
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
Definition: MMFMaxwell.cpp:888
Nektar::MMFMaxwell::TestMaxwellSphere
Array< OneD, NekDouble > TestMaxwellSphere(const NekDouble time, const NekDouble omega, unsigned int field)
Definition: MMFMaxwell.cpp:1788
Nektar::MMFMaxwell::m_coriolis
Array< OneD, NekDouble > m_coriolis
Definition: MMFMaxwell.h:144
Nektar::MMFMaxwell::m_CrossProductMF
Array< OneD, Array< OneD, NekDouble > > m_CrossProductMF
Definition: MMFMaxwell.h:125
Nektar::MMFMaxwell::m_Gaussianradius
NekDouble m_Gaussianradius
Definition: MMFMaxwell.h:123
Nektar::MMFMaxwell::WeakDGMaxwellDirDeriv
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 .
Definition: MMFMaxwell.cpp:1247
Nektar::MMFMaxwell::m_PMLstart
NekDouble m_PMLstart
Definition: MMFMaxwell.h:139
Nektar::MMFMaxwell::ComputeRadCloak
Array< OneD, NekDouble > ComputeRadCloak(const int Nlayer=5)
Definition: MMFMaxwell.cpp:2907
Nektar::MMFMaxwell::print_MMF
void print_MMF(Array< OneD, Array< OneD, NekDouble > > &inarray)
Definition: MMFMaxwell.cpp:3127
Nektar::MMFMaxwell::TestMaxwell2DPEC
Array< OneD, NekDouble > TestMaxwell2DPEC(const NekDouble time, unsigned int field, const SolverUtils::PolType Polarization)
Definition: MMFMaxwell.cpp:1571
Nektar::MMFMaxwell::m_SourceVector
Array< OneD, NekDouble > m_SourceVector
Definition: MMFMaxwell.h:121
Nektar::MMFMaxwell::ComputeEnergyDensity
NekDouble ComputeEnergyDensity(Array< OneD, Array< OneD, NekDouble > > &fields)
Definition: MMFMaxwell.cpp:2181
Nektar::MMFMaxwell::m_Cloakraddelta
NekDouble m_Cloakraddelta
Definition: MMFMaxwell.h:116
Nektar::MMFMaxwell::className
static std::string className
Name of class.
Definition: MMFMaxwell.h:93
Nektar::MMFMaxwell::Checkpoint_EnergyOutput
void Checkpoint_EnergyOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble > > &fieldphys)
Definition: MMFMaxwell.cpp:2727
Nektar::MMFMaxwell::~MMFMaxwell
~MMFMaxwell() override
Destructor.
Definition: MMFMaxwell.cpp:444
Nektar::MMFMaxwell::Checkpoint_TotalFieldOutput
void Checkpoint_TotalFieldOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble > > &fieldphys)
Definition: MMFMaxwell.cpp:2557
Nektar::MMFMaxwell::ComputeMaterialMicroWaveCloak
void ComputeMaterialMicroWaveCloak(const Array< OneD, const NekDouble > &radvec, Array< OneD, Array< OneD, NekDouble > > &epsvec, Array< OneD, Array< OneD, NekDouble > > &muvec)
Definition: MMFMaxwell.cpp:2366
Nektar::MMFMaxwell::ComputeSurfaceCurrent
Array< OneD, NekDouble > ComputeSurfaceCurrent(const int time, const Array< OneD, const Array< OneD, NekDouble > > &fields)
Definition: MMFMaxwell.cpp:1980
Nektar::MMFMaxwell::m_SigmaPML
Array< OneD, Array< OneD, NekDouble > > m_SigmaPML
Definition: MMFMaxwell.h:140
Nektar::MMFMaxwell::Checkpoint_TotPlotOutput
void Checkpoint_TotPlotOutput(const int n, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble > > &fieldphys)
Definition: MMFMaxwell.cpp:2623
Nektar::MMFMaxwell::m_varepsilon
Array< OneD, NekDouble > m_varepsilon
Definition: MMFMaxwell.h:134
Nektar::MMFMaxwell::MMFMaxwell
MMFMaxwell(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Session reader.
Definition: MMFMaxwell.cpp:55
Nektar::MMFMaxwell::DoOdeProjection
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection.
Definition: MMFMaxwell.cpp:1319
Nektar::MMFMaxwell::v_InitObject
void v_InitObject(bool DeclareFields=true) override
Initialise the object.
Definition: MMFMaxwell.cpp:64
Nektar::MMFMaxwell::ComputeMaterialOpticalCloak
void ComputeMaterialOpticalCloak(const Array< OneD, const NekDouble > &radvec, Array< OneD, Array< OneD, NekDouble > > &epsvec, Array< OneD, Array< OneD, NekDouble > > &muvec, const bool Dispersion=false)
Definition: MMFMaxwell.cpp:2294
Nektar::MMFMaxwell::AddGreenDerivCompensate
void AddGreenDerivCompensate(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFMaxwell.cpp:1170
Nektar::MMFMaxwell::m_SourceType
SourceType m_SourceType
Definition: MMFMaxwell.h:79
Nektar::MMFMaxwell::ComputeMaterialVector
void ComputeMaterialVector(Array< OneD, Array< OneD, NekDouble > > &epsvec, Array< OneD, Array< OneD, NekDouble > > &muvec)
Definition: MMFMaxwell.cpp:2199
Nektar::MMFMaxwell::create
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: MMFMaxwell.h:83
Nektar::MMFMaxwell::m_freq
NekDouble m_freq
Definition: MMFMaxwell.h:131
Nektar::MMFMaxwell::m_PrintoutSurfaceCurrent
int m_PrintoutSurfaceCurrent
Definition: MMFMaxwell.h:107
Nektar::MMFMaxwell::GenerateSigmaPML
void GenerateSigmaPML(const NekDouble PMLthickness, const NekDouble PMLstart, const NekDouble PMLmaxsigma, Array< OneD, Array< OneD, NekDouble > > &SigmaPML)
Definition: MMFMaxwell.cpp:2025
Nektar::MMFMaxwell::m_PMLthickness
NekDouble m_PMLthickness
Definition: MMFMaxwell.h:139
Nektar::MMFMaxwell::GaussianPulse
Array< OneD, NekDouble > GaussianPulse(const NekDouble time, const NekDouble Psx, const NekDouble Psy, const NekDouble Psz, const NekDouble Gaussianradius)
Definition: MMFMaxwell.cpp:2790
Nektar::MMFMaxwell::Printout_SurfaceCurrent
void Printout_SurfaceCurrent(Array< OneD, Array< OneD, NekDouble > > &fields, const int time)
Definition: MMFMaxwell.cpp:1902
Nektar::MMFMaxwell::v_EvaluateExactSolution
void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time) override
Definition: MMFMaxwell.cpp:1423
Nektar::MMFMaxwell::m_wp2Tol
NekDouble m_wp2Tol
Definition: MMFMaxwell.h:118
Nektar::MMFMaxwell::v_DoSolve
void v_DoSolve() override
Solves an unsteady problem.
Definition: MMFMaxwell.cpp:448
Nektar::MMFMaxwell::m_wp2
Array< OneD, NekDouble > m_wp2
Definition: MMFMaxwell.h:119
Nektar::MMFMaxwell::v_GenerateSummary
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print Summary.
Definition: MMFMaxwell.cpp:3059
Nektar::MMFMaxwell::AddPML
void AddPML(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFMaxwell.cpp:2420
Nektar::MMFMaxwell::v_SetInitialConditions
void v_SetInitialConditions(const NekDouble initialtime, bool dumpInitialConditions, const int domain) override
Definition: MMFMaxwell.cpp:1336
Nektar::MMFMaxwell::TestMaxwell2DPMC
Array< OneD, NekDouble > TestMaxwell2DPMC(const NekDouble time, unsigned int field, const SolverUtils::PolType Polarization)
Definition: MMFMaxwell.cpp:1702
Nektar::MMFMaxwell::EvaluateCoriolis
Array< OneD, NekDouble > EvaluateCoriolis()
Definition: MMFMaxwell.cpp:2827
Nektar::MMFMaxwell::TestMaxwell1D
Array< OneD, NekDouble > TestMaxwell1D(const NekDouble time, unsigned int field)
Definition: MMFMaxwell.cpp:1465
Nektar::MMFMaxwell::Checkpoint_PlotOutput
void Checkpoint_PlotOutput(const int n, const Array< OneD, const Array< OneD, NekDouble > > &fieldphys)
Definition: MMFMaxwell.cpp:2589
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:396
Nektar::SolverUtils::EquationSystem::GetTraceNpoints
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Definition: EquationSystem.h:744
Nektar::SolverUtils::EquationSystem::m_timestep
NekDouble m_timestep
Time step size.
Definition: EquationSystem.h:373
Nektar::SolverUtils::EquationSystem::m_infosteps
int m_infosteps
Number of time steps between outputting status information.
Definition: EquationSystem.h:390
Nektar::SolverUtils::EquationSystem::m_time
NekDouble m_time
Current time of simulation.
Definition: EquationSystem.h:367
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:739
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:359
Nektar::SolverUtils::EquationSystem::m_fintime
NekDouble m_fintime
Finish time of the simulation.
Definition: EquationSystem.h:371
Nektar::SolverUtils::EquationSystem::Checkpoint_Output
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
Definition: EquationSystem.cpp:1176
Nektar::SolverUtils::EquationSystem::m_checktime
NekDouble m_checktime
Time between checkpoints.
Definition: EquationSystem.h:378
Nektar::SolverUtils::EquationSystem::WriteFld
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition: EquationSystem.cpp:1238
Nektar::SolverUtils::EquationSystem::m_sessionName
std::string m_sessionName
Name of the session.
Definition: EquationSystem.h:365
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:352
Nektar::SolverUtils::EquationSystem::GetNpoints
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:774
Nektar::SolverUtils::EquationSystem::GetNcoeffs
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:706
Nektar::SolverUtils::EquationSystem::m_steps
int m_steps
Number of steps to take.
Definition: EquationSystem.h:386
Nektar::SolverUtils::EquationSystem::m_boundaryConditions
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
Pointer to boundary conditions object.
Definition: EquationSystem.h:361
Nektar::SolverUtils::EquationSystem::m_checksteps
int m_checksteps
Number of steps between checkpoints.
Definition: EquationSystem.h:388
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:764
Nektar::SolverUtils::MMFSystem
A base class for PDEs which include an advection component.
Definition: MMFSystem.h:144
Nektar::SolverUtils::MMFSystem::GetMaxwellFluxVector
SOLVER_UTILS_EXPORT void GetMaxwellFluxVector(const int var, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
Definition: MMFSystem.cpp:1603
Nektar::SolverUtils::MMFSystem::v_GenerateSummary
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Virtual function for generating summary information.
Definition: MMFSystem.cpp:2463
Nektar::SolverUtils::MMFSystem::DeriveCrossProductMF
SOLVER_UTILS_EXPORT void DeriveCrossProductMF(Array< OneD, Array< OneD, NekDouble > > &CrossProductMF)
Definition: MMFSystem.cpp:568
Nektar::SolverUtils::MMFSystem::NumericalMaxwellFlux
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)
Definition: MMFSystem.cpp:1724
Nektar::SolverUtils::MMFSystem::AdddedtMaxwell
SOLVER_UTILS_EXPORT void AdddedtMaxwell(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSystem.cpp:1364
Nektar::SolverUtils::MMFSystem::ComputeZimYim
SOLVER_UTILS_EXPORT void ComputeZimYim(Array< OneD, Array< OneD, NekDouble > > &epsvec, Array< OneD, Array< OneD, NekDouble > > &muvec)
Definition: MMFSystem.cpp:1102
Nektar::SolverUtils::MMFSystem::m_muvec
Array< OneD, Array< OneD, NekDouble > > m_muvec
Definition: MMFSystem.h:210
Nektar::SolverUtils::MMFSystem::m_MMFdir
SpatialDomains::GeomMMF m_MMFdir
Definition: MMFSystem.h:219
Nektar::SolverUtils::MMFSystem::m_movingframes
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:183
Nektar::SolverUtils::MMFSystem::m_MMFfactors
Array< OneD, NekDouble > m_MMFfactors
Definition: MMFSystem.h:157
Nektar::SolverUtils::MMFSystem::GetIncidentField
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > GetIncidentField(const int var, const NekDouble time)
Definition: MMFSystem.cpp:2005
Nektar::SolverUtils::MMFSystem::m_ntimesMFFwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFFwd
Definition: MMFSystem.h:199
Nektar::SolverUtils::MMFSystem::m_negepsvecminus1
Array< OneD, Array< OneD, NekDouble > > m_negepsvecminus1
Definition: MMFSystem.h:212
Nektar::SolverUtils::MMFSystem::CartesianToMovingframes
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble > > &uvec, unsigned int field)
Definition: MMFSystem.cpp:771
Nektar::SolverUtils::MMFSystem::m_epsvec
Array< OneD, Array< OneD, NekDouble > > m_epsvec
Definition: MMFSystem.h:209
Nektar::SolverUtils::MMFSystem::CartesianToSpherical
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)
Definition: MMFSystem.cpp:792
Nektar::SolverUtils::MMFSystem::m_negmuvecminus1
Array< OneD, Array< OneD, NekDouble > > m_negmuvecminus1
Definition: MMFSystem.h:213
Nektar::SolverUtils::MMFSystem::ComputeNtimesMF
SOLVER_UTILS_EXPORT void ComputeNtimesMF()
Definition: MMFSystem.cpp:615
Nektar::SolverUtils::MMFSystem::RootMeanSquare
SOLVER_UTILS_EXPORT NekDouble RootMeanSquare(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2338
Nektar::SolverUtils::MMFSystem::MMFInitObject
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble > > &Anisotropy, const int TangentXelem=-1)
Definition: MMFSystem.cpp:51
Nektar::SolverUtils::MMFSystem::m_TestMaxwellType
TestMaxwellType m_TestMaxwellType
Definition: MMFSystem.h:153
Nektar::SolverUtils::MMFSystem::m_CurlMF
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
Definition: MMFSystem.h:193
Nektar::SolverUtils::MMFSystem::Computedemdxicdote
SOLVER_UTILS_EXPORT void Computedemdxicdote()
Definition: MMFSystem.cpp:1274
Nektar::SolverUtils::UnsteadySystem
Base class for unsteady solvers.
Definition: UnsteadySystem.h:47
Nektar::SolverUtils::UnsteadySystem::m_ode
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:89
Nektar::SolverUtils::UnsteadySystem::m_intScheme
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
Definition: UnsteadySystem.h:86
Nektar::SolverUtils::UnsteadySystem::m_explicitAdvection
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
Definition: UnsteadySystem.h:109
Nektar::SolverUtils::UnsteadySystem::v_InitObject
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
Definition: UnsteadySystem.cpp:85
FilterPython_Function.field
field
Definition: FilterPython_Function.py:4
Nektar::LibUtilities::rad
static NekDouble rad(NekDouble x, NekDouble y)
Definition: Interpreter/Interpreter.cpp:76
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:113
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:49
Nektar::SolverUtils::IncTypeMap
const char *const IncTypeMap[]
Definition: MMFSystem.h:136
Nektar::SolverUtils::SummaryList
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
Nektar::SolverUtils::SIZE_TestMaxwellType
@ SIZE_TestMaxwellType
Length of enum list.
Definition: MMFSystem.h:105
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:99
Nektar::SolverUtils::TestMaxwellTypeMap
const char *const TestMaxwellTypeMap[]
Definition: MMFSystem.h:108
Nektar::SolverUtils::AddSummaryItem
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:47
Nektar::SolverUtils::PolTypeMap
const char *const PolTypeMap[]
Definition: MMFSystem.h:123
Nektar::SpatialDomains::eTangentIrregular
@ eTangentIrregular
Circular around the centre of domain.
Definition: SpatialDomains.hpp:78
Nektar::SpatialDomains::eTangentCircular
@ eTangentCircular
Circular around the centre of domain.
Definition: SpatialDomains.hpp:77
Nektar::SpatialDomains::eTangentNonconvex
@ eTangentNonconvex
Circular around the centre of domain.
Definition: SpatialDomains.hpp:79
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
Nektar::UnitTests::z
std::vector< double > z(NPUPPER)
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
Vmath::Vmul
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.
Definition: Vmath.hpp:72
Vmath::Svtvp
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.
Definition: Vmath.hpp:396
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
Vmath::Vmin
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.hpp:725
Vmath::Vvtvp
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
Definition: Vmath.hpp:366
Vmath::Vadd
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.
Definition: Vmath.hpp:180
Vmath::Vvtvm
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
Definition: Vmath.hpp:381
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
Vmath::Vdiv
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.
Definition: Vmath.hpp:126
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
Vmath::Vamax
T Vamax(int n, const T *x, const int incx)
Return the maximum absolute element in x called vamax to avoid conflict with max.
Definition: Vmath.hpp:685
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.hpp:194
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.hpp:644
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
Vmath::Vsub
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.
Definition: Vmath.hpp:220
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294
Generated on Wed May 15 2024 07:37:00 for Nektar++ by doxygen 1.9.4