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