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