Nektar++
MMFSystem.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: MMFSystem.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: Base class for MMF systems.
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
36#include <SolverUtils/MMFSystem.h>
37
38namespace Nektar::SolverUtils
39{
40
41MMFSystem::MMFSystem(const LibUtilities::SessionReaderSharedPtr &pSession,
42 const SpatialDomains::MeshGraphSharedPtr &pGraph)
43 : UnsteadySystem(pSession, pGraph)
44{
45}
46
47MMFSystem::~MMFSystem()
48{
49}
50
51void MMFSystem::MMFInitObject(
52 const Array<OneD, const Array<OneD, NekDouble>> &Anisotropy,
53 const int TangentXelem)
54{
55 m_pi = 3.14159265358979323846;
56 m_shapedim = m_fields[0]->GetShapeDimension();
57
58 ASSERTL0(m_spacedim == 3, "Space Dimension should be 3.");
59
60 // Factor for Numerical Flux
61 m_session->LoadParameter("alpha", m_alpha, 1.0);
62
63 // Factor for Numerical Flux
64 m_session->LoadParameter("Incfreq", m_Incfreq, 1.0);
65
66 // SmoothFactor
67 m_session->LoadParameter("SmoothFactor", m_SmoothFactor, 1);
68 m_session->LoadParameter("SFinit", m_SFinit, 0.0);
69
70 // Define SurfaceType
71 if (m_session->DefinesSolverInfo("SURFACETYPE"))
72 {
73 std::string SurfaceTypeStr = m_session->GetSolverInfo("SURFACETYPE");
74 int i;
75 for (i = 0; i < (int)SIZE_SurfaceType; ++i)
76 {
77 if (boost::iequals(SurfaceTypeMap[i], SurfaceTypeStr))
78 {
79 m_surfaceType = (SurfaceType)i;
80 break;
81 }
82 }
83 }
84 else
85 {
86 m_surfaceType = (SurfaceType)0;
87 }
88
89 // if discontinuous Galerkin determine numerical flux to use
90 for (int i = 0; i < (int)SIZE_UpwindType; ++i)
91 {
92 bool match;
93 m_session->MatchSolverInfo("UPWINDTYPE", UpwindTypeMap[i], match,
94 false);
95 if (match)
96 {
97 m_upwindType = (UpwindType)i;
98 break;
99 }
100 }
101
102 // SetUpMovingFrames: To generate m_movingframes
103 SetUpMovingFrames(Anisotropy, TangentXelem);
104
105 // Derive movingfrmaes at interfaces
106 // Get: m_MFtraceFwd and m_MFtraceBwd
107 ComputeMFtrace();
108
109 // Check Movingframes and surfraceNormal
110 // Get: m_ncdotMFFwd,m_ncdotMFBwd,m_nperpcdotMFFwd,m_nperpcdotMFBwd
111 ComputencdotMF();
112
113 // Compute nabla cdot m_movingframes
114 // Get: m_DivMF, m_CurlMF
115 ComputeDivCurlMF();
116}
117
118void MMFSystem::SetUpMovingFrames(
119 const Array<OneD, const Array<OneD, NekDouble>> &Anisotropy,
120 const int TangentXelem)
121{
122
123 int nq = m_fields[0]->GetNpoints();
124 int MFdim = 3;
125
126 // Construct The Moving Frames
127 m_movingframes = Array<OneD, Array<OneD, NekDouble>>(MFdim);
128
129 for (int j = 0; j < MFdim; ++j)
130 {
131 m_movingframes[j] = Array<OneD, NekDouble>(m_spacedim * nq, 0.0);
132 }
133
134 // Read MMF Geom Info
135 std::string conn = "TangentX";
136 m_MMFfactors = Array<OneD, NekDouble>(4);
137
138 m_session->LoadSolverInfo("MMFDir", conn, "LOCAL");
139
140 // (x-x_0)^2/a^2 + (y-y_0)^2/b^2 = 1
141 // factors[0] = a
142 // factors[1] = b
143 // factors[2] = x_0
144 // factors[3] = y_0
145 m_session->LoadParameter("MMFCircAxisX", m_MMFfactors[0], 1.0);
146 m_session->LoadParameter("MMFCircAxisY", m_MMFfactors[1], 1.0);
147 m_session->LoadParameter("MMFCircCentreX", m_MMFfactors[2], 0.0);
148 m_session->LoadParameter("MMFCircCentreY", m_MMFfactors[3], 0.0);
149
150 if (conn == "TangentX")
151 {
152 m_MMFdir = SpatialDomains::eTangentX;
153 }
154 if (conn == "TangentY")
155 {
156 m_MMFdir = SpatialDomains::eTangentY;
157 }
158 if (conn == "TangentXY")
159 {
160 m_MMFdir = SpatialDomains::eTangentXY;
161 }
162 if (conn == "TangentZ")
163 {
164 m_MMFdir = SpatialDomains::eTangentZ;
165 }
166 if (conn == "TangentCircular")
167 {
168 m_MMFdir = SpatialDomains::eTangentCircular;
169 }
170 if (conn == "TangentIrregular")
171 {
172 m_MMFdir = SpatialDomains::eTangentIrregular;
173 }
174 if (conn == "TangentNonconvex")
175 {
176 m_MMFdir = SpatialDomains::eTangentNonconvex;
177 }
178 if (conn == "LOCAL")
179 {
180 m_MMFdir = SpatialDomains::eLOCAL;
181 }
182
183 // Get Tangetn vectors from GeomFactors2D, Orthonormalized = true
184 m_fields[0]->GetMovingFrames(m_MMFdir, m_MMFfactors, m_movingframes);
185
186 // Align the tangentX direction after TangentXelem
187 if (TangentXelem > 0)
188 {
189 Array<OneD, NekDouble> tmp(nq);
190 m_fields[0]->GenerateElementVector(TangentXelem, 1.0, 0.0, tmp);
191 for (int j = 0; j < m_shapedim; j++)
192 {
193 for (int k = 0; k < m_spacedim; k++)
194 {
195 Vmath::Vmul(nq, &tmp[0], 1, &m_movingframes[j][k * nq], 1,
196 &m_movingframes[j][k * nq], 1);
197 }
198 }
199
200 m_fields[0]->GenerateElementVector(TangentXelem, 0.0, 1.0, tmp);
201 Vmath::Vadd(nq, &tmp[0], 1, &m_movingframes[0][0 * nq], 1,
202 &m_movingframes[0][0 * nq], 1);
203 Vmath::Vadd(nq, &tmp[0], 1, &m_movingframes[1][1 * nq], 1,
204 &m_movingframes[1][1 * nq], 1);
205
206 int indxtmp = Vmath::Imax(nq, tmp, 1);
207 std::cout << "*** MF in PML Region is aligned as MF1 = ( "
208 << m_movingframes[0][indxtmp] << " , "
209 << m_movingframes[0][nq + indxtmp] << " , "
210 << m_movingframes[0][2 * nq + indxtmp] << " ) "
211 << ", MF2 = ( " << m_movingframes[1][indxtmp] << " , "
212 << m_movingframes[1][nq + indxtmp] << " , "
213 << m_movingframes[1][2 * nq + indxtmp] << " ) " << std::endl;
214 }
215
216 // Multiply Anisotropy to movingframes
217 for (int j = 0; j < m_shapedim; ++j)
218 {
219 for (int k = 0; k < m_spacedim; ++k)
220 {
221 Vmath::Vmul(nq, &Anisotropy[j][0], 1, &m_movingframes[j][k * nq], 1,
222 &m_movingframes[j][k * nq], 1);
223 }
224 }
225
226 // Test the moving frames
227 CheckMovingFrames(m_movingframes);
228}
229
230void MMFSystem::CheckMovingFrames(
231 const Array<OneD, const Array<OneD, NekDouble>> &movingframes)
232{
233 NekDouble t1x, t1y, t1z, t2x, t2y, t2z, t3x, t3y, t3z;
234 NekDouble dot12 = 0.0, dot23 = 0.0, dot31 = 0.0;
235 NekDouble Tol = 0.0001;
236
237 int nq = m_fields[0]->GetNpoints();
238
239 for (int i = 0; i < nq; ++i)
240 {
241 t1x = movingframes[0][i];
242 t1y = movingframes[0][i + nq];
243 t1z = movingframes[0][i + 2 * nq];
244
245 t2x = movingframes[1][i];
246 t2y = movingframes[1][i + nq];
247 t2z = movingframes[1][i + 2 * nq];
248
249 t3x = movingframes[2][i];
250 t3y = movingframes[2][i + nq];
251 t3z = movingframes[2][i + 2 * nq];
252
253 dot12 = t1x * t2x + t1y * t2y + t1z * t2z;
254 dot23 = t2x * t3x + t2y * t3y + t2z * t3z;
255 dot31 = t3x * t1x + t3y * t1y + t3z * t1z;
256 }
257
258 std::cout << "======================================================"
259 << std::endl;
260 std::cout << "======================================================"
261 << std::endl;
262 std::cout << "*** The first moving frame is alinged along"
263 << SpatialDomains::GeomMMFMap[m_MMFdir] << std::endl;
264
265 Array<OneD, NekDouble> tmpx(nq), tmpy(nq), tmpz(nq);
266
267 Vmath::Vcopy(nq, &movingframes[0][0], 1, &tmpx[0], 1);
268 Vmath::Vcopy(nq, &movingframes[0][nq], 1, &tmpy[0], 1);
269 Vmath::Vcopy(nq, &movingframes[0][2 * nq], 1, &tmpz[0], 1);
270 std::cout << nq << " , "
271 << "*** Avg MF1 = ( " << AvgAbsInt(tmpx) << " , "
272 << AvgAbsInt(tmpy) << " , " << AvgAbsInt(tmpz) << " ) "
273 << std::endl;
274
275 Vmath::Vcopy(nq, &movingframes[1][0], 1, &tmpx[0], 1);
276 Vmath::Vcopy(nq, &movingframes[1][nq], 1, &tmpy[0], 1);
277 Vmath::Vcopy(nq, &movingframes[1][2 * nq], 1, &tmpz[0], 1);
278 std::cout << "*** Avg MF2 = ( " << AvgAbsInt(tmpx) << " , "
279 << AvgAbsInt(tmpy) << " , " << AvgAbsInt(tmpz) << " ) "
280 << std::endl;
281
282 if (m_shapedim == 3)
283 {
284 Vmath::Vcopy(nq, &movingframes[2][0], 1, &tmpx[0], 1);
285 Vmath::Vcopy(nq, &movingframes[2][nq], 1, &tmpy[0], 1);
286 Vmath::Vcopy(nq, &movingframes[2][2 * nq], 1, &tmpz[0], 1);
287 std::cout << "*** Avg MF3 = ( " << AvgAbsInt(tmpx) << " , "
288 << AvgAbsInt(tmpy) << " , " << AvgAbsInt(tmpz) << " ) "
289 << std::endl;
290 }
291
292 if ((fabs(dot12) + fabs(dot23) + fabs(dot31)) < Tol)
293 {
294 std::cout << "*** Moving frames are Orthogonal" << std::endl;
295 }
296
297 else
298 {
299 std::cout << "*** Moving frames are NOT Orthogonal" << std::endl;
300 }
301
302 Array<OneD, NekDouble> tmp;
303 for (int j = 0; j < m_shapedim; ++j)
304 {
305 tmp = Array<OneD, NekDouble>(nq, 0.0);
306 for (int k = 0; k < m_spacedim; ++k)
307 {
308 Vmath::Vvtvp(nq, &movingframes[j][k * nq], 1,
309 &movingframes[j][k * nq], 1, &tmp[0], 1, &tmp[0], 1);
310 }
311 Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
312 std::cout << "*** Avg. Magnitude of MF" << j << " = " << AvgAbsInt(tmp)
313 << std::endl;
314 }
315}
316
317void MMFSystem::ComputencdotMF()
318{
319 int nq = m_fields[0]->GetNpoints();
320 int nTracePointsTot = GetTraceNpoints();
321
322 // Compute MFjFwd and MFjBwd
323 Array<OneD, NekDouble> tmp(nq);
324
325 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> MFtraceFwd(m_shapedim);
326 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> MFtraceBwd(m_shapedim);
327 Array<OneD, Array<OneD, NekDouble>> SurfaceNormalFwd;
328 Array<OneD, Array<OneD, NekDouble>> SurfaceNormalBwd;
329
330 for (int j = 0; j < m_shapedim; ++j)
331 {
332 MFtraceFwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
333 MFtraceBwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
334
335 SurfaceNormalFwd = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
336 SurfaceNormalBwd = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
337
338 for (int k = 0; k < m_spacedim; ++k)
339 {
340 MFtraceFwd[j][k] = Array<OneD, NekDouble>(nTracePointsTot);
341 MFtraceBwd[j][k] = Array<OneD, NekDouble>(nTracePointsTot);
342
343 SurfaceNormalFwd[k] = Array<OneD, NekDouble>(nTracePointsTot);
344 SurfaceNormalBwd[k] = Array<OneD, NekDouble>(nTracePointsTot);
345
346 Vmath::Vcopy(nq, &m_movingframes[j][k * nq], 1, &tmp[0], 1);
347
348 m_fields[0]->GetFwdBwdTracePhys(tmp, MFtraceFwd[j][k],
349 MFtraceBwd[j][k]);
350
351 CopyBoundaryTrace(MFtraceFwd[j][k], MFtraceBwd[j][k], eFwdEQBwd);
352 }
353 }
354
355 VectorCrossProd(MFtraceFwd[0], MFtraceFwd[1], SurfaceNormalFwd);
356 VectorCrossProd(MFtraceBwd[0], MFtraceBwd[1], SurfaceNormalBwd);
357
358 // Compute n \times e^i
359 m_ncdotMFFwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
360 m_ncdotMFBwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
361 for (int j = 0; j < m_shapedim; ++j)
362 {
363 m_ncdotMFFwd[j] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
364 m_ncdotMFBwd[j] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
365
366 VectorDotProd(m_traceNormals, MFtraceFwd[j], m_ncdotMFFwd[j]);
367 VectorDotProd(m_traceNormals, MFtraceBwd[j], m_ncdotMFBwd[j]);
368 }
369
370 // Compute n^{\perp} \times e^i
371 Array<OneD, Array<OneD, NekDouble>> SurfaceNormal(m_spacedim);
372 Array<OneD, Array<OneD, NekDouble>> Tracevector(m_spacedim);
373 for (int k = 0; k < m_spacedim; k++)
374 {
375 SurfaceNormal[k] = Array<OneD, NekDouble>(nTracePointsTot);
376 Tracevector[k] = Array<OneD, NekDouble>(nTracePointsTot);
377
378 Vmath::Vcopy(nTracePointsTot, &m_MFtraceFwd[2][k][0], 1,
379 &SurfaceNormal[k][0], 1);
380 }
381
382 VectorCrossProd(m_traceNormals, SurfaceNormal, Tracevector);
383
384 m_nperpcdotMFFwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
385 m_nperpcdotMFBwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
386 for (int j = 0; j < m_shapedim; ++j)
387 {
388 m_nperpcdotMFFwd[j] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
389 m_nperpcdotMFBwd[j] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
390
391 VectorDotProd(Tracevector, MFtraceFwd[j], m_nperpcdotMFFwd[j]);
392 VectorDotProd(Tracevector, MFtraceBwd[j], m_nperpcdotMFBwd[j]);
393 }
394
395 if (m_shapedim == 2)
396 {
397 std::cout << "*** m_ncdotMFFwd = ( " << RootMeanSquare(m_ncdotMFFwd[0])
398 << " , " << RootMeanSquare(m_ncdotMFFwd[1]) << " ) "
399 << std::endl;
400 std::cout << "*** m_ncdotMFBwd = ( " << RootMeanSquare(m_ncdotMFBwd[0])
401 << " , " << RootMeanSquare(m_ncdotMFBwd[1]) << " ) "
402 << std::endl;
403
404 std::cout << "*** m_nperpcdotMFFwd = ( "
405 << RootMeanSquare(m_nperpcdotMFFwd[0]) << " , "
406 << RootMeanSquare(m_nperpcdotMFFwd[1]) << " ) " << std::endl;
407 std::cout << "*** m_nperpcdotMFBwd = ( "
408 << RootMeanSquare(m_nperpcdotMFBwd[0]) << " , "
409 << RootMeanSquare(m_nperpcdotMFBwd[1]) << " ) " << std::endl;
410 }
411
412 else if (m_shapedim == 3)
413 {
414 std::cout << "*** m_ncdotMFFwd = ( "
415 << Vmath::Vsum(nTracePointsTot, m_ncdotMFFwd[0], 1) << " , "
416 << Vmath::Vsum(nTracePointsTot, m_ncdotMFFwd[1], 1) << " , "
417 << Vmath::Vsum(nTracePointsTot, m_ncdotMFFwd[2], 1) << " ) "
418 << std::endl;
419 std::cout << "*** m_ncdotMFBwd = ( "
420 << Vmath::Vsum(nTracePointsTot, m_ncdotMFBwd[0], 1) << " , "
421 << Vmath::Vsum(nTracePointsTot, m_ncdotMFBwd[1], 1) << " , "
422 << Vmath::Vsum(nTracePointsTot, m_ncdotMFBwd[2], 1) << " ) "
423 << std::endl;
424
425 std::cout << "*** m_nperpcdotMFFwd = ( "
426 << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFFwd[0], 1)
427 << " , "
428 << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFFwd[1], 1)
429 << " , "
430 << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFFwd[2], 1)
431 << " ) " << std::endl;
432 std::cout << "*** m_nperpcdotMFBwd = ( "
433 << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFBwd[0], 1)
434 << " , "
435 << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFBwd[1], 1)
436 << " , "
437 << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFBwd[2], 1)
438 << " ) " << std::endl;
439 }
440}
441
442void MMFSystem::ComputeDivCurlMF()
443{
444 int nq = m_fields[0]->GetNpoints();
445 int MMFdim = 3;
446
447 Array<OneD, NekDouble> tmp(nq);
448 Array<OneD, NekDouble> Dtmp(nq);
449
450 m_DivMF = Array<OneD, Array<OneD, NekDouble>>(MMFdim);
451 for (int j = 0; j < MMFdim; ++j)
452 {
453 m_DivMF[j] = Array<OneD, NekDouble>(nq, 0.0);
454 for (int k = 0; k < m_spacedim; ++k)
455 {
456 Vmath::Vcopy(nq, &m_movingframes[j][k * nq], 1, &tmp[0], 1);
457
458 m_fields[0]->PhysDeriv(k, tmp, Dtmp);
459 Vmath::Vadd(nq, &Dtmp[0], 1, &m_DivMF[j][0], 1, &m_DivMF[j][0], 1);
460 }
461 }
462
463 std::cout << "*** Divergence of MF1 = " << AvgInt(m_DivMF[0])
464 << ", MF2 = " << AvgInt(m_DivMF[1])
465 << ", MF3 = " << AvgInt(m_DivMF[2]) << std::endl;
466
467 // Compute Curl of MF: CurlMF[i][j] = (\nabla \times e^i) cdot e^j
468 m_CurlMF = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MMFdim);
469 for (int i = 0; i < MMFdim; ++i)
470 {
471 m_CurlMF[i] = Array<OneD, Array<OneD, NekDouble>>(MMFdim);
472
473 for (int j = 0; j < MMFdim; ++j)
474 {
475 m_CurlMF[i][j] = Array<OneD, NekDouble>(nq, 0.0);
476 }
477 }
478
479 Array<OneD, Array<OneD, NekDouble>> MFtmp(m_spacedim);
480 Array<OneD, Array<OneD, NekDouble>> CurlMFtmp(m_spacedim);
481
482 for (int i = 0; i < m_spacedim; ++i)
483 {
484 MFtmp[i] = Array<OneD, NekDouble>(nq);
485 CurlMFtmp[i] = Array<OneD, NekDouble>(nq);
486 }
487
488 for (int dir = 0; dir < MMFdim; dir++)
489 {
490 for (int i = 0; i < m_spacedim; ++i)
491 {
492 Vmath::Vcopy(nq, &m_movingframes[dir][i * nq], 1, &MFtmp[i][0], 1);
493 }
494
495 ComputeCurl(MFtmp, CurlMFtmp);
496
497 for (int j = 0; j < MMFdim; ++j)
498 {
499 for (int i = 0; i < m_spacedim; ++i)
500 {
501 Vmath::Vvtvp(nq, &m_movingframes[j][i * nq], 1,
502 &CurlMFtmp[i][0], 1, &m_CurlMF[dir][j][0], 1,
503 &m_CurlMF[dir][j][0], 1);
504 }
505 }
506 }
507
508 std::cout << "*** Curl of MF1 = ( " << AvgInt(m_CurlMF[0][0]) << " , "
509 << AvgInt(m_CurlMF[0][1]) << " , " << AvgInt(m_CurlMF[0][2])
510 << " ) " << std::endl;
511 std::cout << "*** Curl of MF2 = ( " << AvgInt(m_CurlMF[1][0]) << " , "
512 << AvgInt(m_CurlMF[1][1]) << " , " << AvgInt(m_CurlMF[1][2])
513 << " ) " << std::endl;
514 std::cout << "*** Curl of MF3 = ( " << AvgInt(m_CurlMF[2][0]) << " , "
515 << AvgInt(m_CurlMF[2][1]) << " , " << AvgInt(m_CurlMF[2][2])
516 << " ) " << std::endl;
517}
518
519void MMFSystem::ComputeMFtrace()
520{
521 int MFdim = 3;
522
523 int nq = m_fields[0]->GetNpoints();
524 int nTraceNumPoints = GetTraceTotPoints();
525
526 Array<OneD, NekDouble> tmp(nq);
527 Array<OneD, NekDouble> Fwdtmp(nq);
528 Array<OneD, NekDouble> Bwdtmp(nq);
529
530 m_MFtraceFwd = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
531 m_MFtraceBwd = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
532
533 for (int j = 0; j < MFdim; ++j)
534 {
535 m_MFtraceFwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
536 m_MFtraceBwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
537 }
538
539 // m_MFtraceFwd[0] = e^1_{Fwd}, m_MFtraceFwd[1] = e^2_{Fwd}
540 for (int j = 0; j < MFdim; ++j)
541 {
542 for (int i = 0; i < m_spacedim; ++i)
543 {
544 m_MFtraceFwd[j][i] = Array<OneD, NekDouble>(nTraceNumPoints);
545 m_MFtraceBwd[j][i] = Array<OneD, NekDouble>(nTraceNumPoints);
546
547 Vmath::Vcopy(nq, &m_movingframes[j][i * nq], 1, &tmp[0], 1);
548
549 m_fields[0]->GetFwdBwdTracePhys(tmp, Fwdtmp, Bwdtmp);
550
551 CopyBoundaryTrace(Fwdtmp, Bwdtmp, eFwdEQBwd);
552
553 Vmath::Vcopy(nTraceNumPoints, &Fwdtmp[0], 1, &m_MFtraceFwd[j][i][0],
554 1);
555 Vmath::Vcopy(nTraceNumPoints, &Bwdtmp[0], 1, &m_MFtraceBwd[j][i][0],
556 1);
557 }
558 }
559
560 std::cout << "*** MFtraceFwd = ( " << VectorAvgMagnitude(m_MFtraceFwd[0])
561 << " , " << VectorAvgMagnitude(m_MFtraceFwd[1]) << " , "
562 << VectorAvgMagnitude(m_MFtraceFwd[2]) << " ) " << std::endl;
563 std::cout << "*** MFtraceBwd = ( " << VectorAvgMagnitude(m_MFtraceBwd[0])
564 << " , " << VectorAvgMagnitude(m_MFtraceBwd[1]) << " , "
565 << VectorAvgMagnitude(m_MFtraceBwd[2]) << " ) " << std::endl;
566}
567
568void MMFSystem::DeriveCrossProductMF(
569 Array<OneD, Array<OneD, NekDouble>> &CrossProductMF)
570{
571 int MFdim = 3;
572 int nq = GetTotPoints();
573
574 CrossProductMF = Array<OneD, Array<OneD, NekDouble>>(MFdim);
575
576 Array<OneD, Array<OneD, NekDouble>> MF1tmp(m_spacedim);
577 Array<OneD, Array<OneD, NekDouble>> MF2tmp(m_spacedim);
578 Array<OneD, Array<OneD, NekDouble>> MF3tmp(m_spacedim);
579 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> MFtmpCurl(MFdim);
580 for (int j = 0; j < MFdim; ++j)
581 {
582 CrossProductMF[j] = Array<OneD, NekDouble>(nq * m_spacedim);
583 MFtmpCurl[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
584 for (int k = 0; k < m_spacedim; ++k)
585 {
586 MFtmpCurl[j][k] = Array<OneD, NekDouble>(nq);
587 }
588 }
589
590 for (int k = 0; k < m_spacedim; ++k)
591 {
592 MF1tmp[k] = Array<OneD, NekDouble>(nq);
593 MF2tmp[k] = Array<OneD, NekDouble>(nq);
594 MF3tmp[k] = Array<OneD, NekDouble>(nq);
595
596 Vmath::Vcopy(nq, &m_movingframes[0][k * nq], 1, &MF1tmp[k][0], 1);
597 Vmath::Vcopy(nq, &m_movingframes[1][k * nq], 1, &MF2tmp[k][0], 1);
598 Vmath::Vcopy(nq, &m_movingframes[2][k * nq], 1, &MF3tmp[k][0], 1);
599 }
600
601 VectorCrossProd(MF3tmp, MF1tmp, MFtmpCurl[0]);
602 VectorCrossProd(MF2tmp, MF3tmp, MFtmpCurl[1]);
603 VectorCrossProd(MF1tmp, MF2tmp, MFtmpCurl[2]);
604
605 for (int j = 0; j < MFdim; ++j)
606 {
607 for (int k = 0; k < m_spacedim; ++k)
608 {
609 Vmath::Vcopy(nq, &MFtmpCurl[j][k][0], 1, &CrossProductMF[j][k * nq],
610 1);
611 }
612 }
613}
614
615void MMFSystem::ComputeNtimesMF()
616{
617 int MFdim = 3;
618 int nTracePointsTot = GetTraceNpoints();
619
620 m_ntimesMFFwd = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
621 m_ntimesMFBwd = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
622 m_ntimes_ntimesMFFwd =
623 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
624 m_ntimes_ntimesMFBwd =
625 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
626 for (int j = 0; j < MFdim; ++j)
627 {
628 m_ntimesMFFwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
629 m_ntimesMFBwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
630 m_ntimes_ntimesMFFwd[j] =
631 Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
632 m_ntimes_ntimesMFBwd[j] =
633 Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
634
635 for (int k = 0; k < m_spacedim; ++k)
636 {
637 m_ntimesMFFwd[j][k] = Array<OneD, NekDouble>(nTracePointsTot);
638 m_ntimesMFBwd[j][k] = Array<OneD, NekDouble>(nTracePointsTot);
639 m_ntimes_ntimesMFFwd[j][k] =
640 Array<OneD, NekDouble>(nTracePointsTot);
641 m_ntimes_ntimesMFBwd[j][k] =
642 Array<OneD, NekDouble>(nTracePointsTot);
643 }
644
645 VectorCrossProd(m_traceNormals, m_MFtraceFwd[j], m_ntimesMFFwd[j]);
646 VectorCrossProd(m_traceNormals, m_MFtraceBwd[j], m_ntimesMFBwd[j]);
647 VectorCrossProd(m_traceNormals, m_ntimesMFFwd[j],
648 m_ntimes_ntimesMFFwd[j]);
649 VectorCrossProd(m_traceNormals, m_ntimesMFBwd[j],
650 m_ntimes_ntimesMFBwd[j]);
651 }
652
653 std::cout << "*** m_ntimesMFFwd = ( "
654 << VectorAvgMagnitude(m_ntimesMFFwd[0]) << " , "
655 << VectorAvgMagnitude(m_ntimesMFFwd[1]) << " , "
656 << VectorAvgMagnitude(m_ntimesMFFwd[2]) << " ) " << std::endl;
657 std::cout << "*** m_ntimesMFBwd = ( "
658 << VectorAvgMagnitude(m_ntimesMFBwd[0]) << " , "
659 << VectorAvgMagnitude(m_ntimesMFBwd[1]) << " , "
660 << VectorAvgMagnitude(m_ntimesMFBwd[2]) << " ) " << std::endl;
661 std::cout << "*** m_ntimes_ntimesMFFwd = ( "
662 << VectorAvgMagnitude(m_ntimes_ntimesMFFwd[0]) << " , "
663 << VectorAvgMagnitude(m_ntimes_ntimesMFFwd[1]) << " , "
664 << VectorAvgMagnitude(m_ntimes_ntimesMFFwd[2]) << " ) "
665 << std::endl;
666 std::cout << "*** m_ntimes_ntimesMFBwd = ( "
667 << VectorAvgMagnitude(m_ntimes_ntimesMFBwd[0]) << " , "
668 << VectorAvgMagnitude(m_ntimes_ntimesMFBwd[1]) << " , "
669 << VectorAvgMagnitude(m_ntimes_ntimesMFBwd[2]) << " ) "
670 << std::endl;
671}
672
673void MMFSystem::VectorDotProd(
674 const Array<OneD, const Array<OneD, NekDouble>> &v1,
675 const Array<OneD, const Array<OneD, NekDouble>> &v2,
676 Array<OneD, NekDouble> &v3)
677{
678 int coordim = v1.size();
679 int nq = v1[0].size();
680
681 v3 = Array<OneD, NekDouble>(nq, 0.0);
682 for (int i = 0; i < coordim; ++i)
683 {
684 Vmath::Vvtvp(nq, &v1[i][0], 1, &v2[i][0], 1, &v3[0], 1, &v3[0], 1);
685 }
686}
687
688/**
689 * Computes the vector cross-product in 3D of \a v1 and \a v2, storing
690 * the result in \a v3.
691 * @param v1 First input vector.
692 * @param v2 Second input vector.
693 * @param v3 Output vector computed to be orthogonal to
694 * both \a v1 and \a v2.
695 */
696void MMFSystem::VectorCrossProd(
697 const Array<OneD, const Array<OneD, NekDouble>> &v1,
698 const Array<OneD, const Array<OneD, NekDouble>> &v2,
699 Array<OneD, Array<OneD, NekDouble>> &v3)
700{
701 ASSERTL0(v1.size() == 3, "Input 1 has dimension not equal to 3.");
702 ASSERTL0(v2.size() == 3, "Input 2 has dimension not equal to 3.");
703 ASSERTL0(v3.size() == 3, "Output vector has dimension not equal to 3.");
704
705 int nq = v1[0].size();
706 Array<OneD, NekDouble> temp(nq);
707
708 Vmath::Vmul(nq, v1[2], 1, v2[1], 1, temp, 1);
709 Vmath::Vvtvm(nq, v1[1], 1, v2[2], 1, temp, 1, v3[0], 1);
710
711 Vmath::Vmul(nq, v1[0], 1, v2[2], 1, temp, 1);
712 Vmath::Vvtvm(nq, v1[2], 1, v2[0], 1, temp, 1, v3[1], 1);
713
714 Vmath::Vmul(nq, v1[1], 1, v2[0], 1, temp, 1);
715 Vmath::Vvtvm(nq, v1[0], 1, v2[1], 1, temp, 1, v3[2], 1);
716}
717
718void MMFSystem::VectorCrossProd(const Array<OneD, NekDouble> &v1,
719 const Array<OneD, NekDouble> &v2,
720 Array<OneD, NekDouble> &v3)
721{
722 ASSERTL0(v1.size() == 3, "Input 1 has dimension not equal to 3.");
723 ASSERTL0(v2.size() == 3, "Input 2 has dimension not equal to 3.");
724 ASSERTL0(v3.size() == 3, "Output vector has dimension not equal to 3.");
725
726 v3[0] = v1[1] * v2[2] - v1[2] * v2[1];
727 v3[1] = v1[2] * v2[0] - v1[0] * v2[2];
728 v3[2] = v1[0] * v2[1] - v1[1] * v2[0];
729}
730
731void MMFSystem::ComputeCurl(
732 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
733 Array<OneD, Array<OneD, NekDouble>> &outarray)
734
735{
736 int nq = inarray[0].size();
737
738 Array<OneD, NekDouble> tmpx, tmpy, tmpz;
739 Array<OneD, NekDouble> Dtmpzdx, Dtmpydx, Dtmpxdy, Dtmpzdy, Dtmpxdz, Dtmpydz;
740
741 tmpx = Array<OneD, NekDouble>(nq);
742 tmpy = Array<OneD, NekDouble>(nq);
743 tmpz = Array<OneD, NekDouble>(nq);
744
745 Dtmpzdx = Array<OneD, NekDouble>(nq);
746 Dtmpydx = Array<OneD, NekDouble>(nq);
747 Dtmpxdy = Array<OneD, NekDouble>(nq);
748 Dtmpzdy = Array<OneD, NekDouble>(nq);
749 Dtmpxdz = Array<OneD, NekDouble>(nq);
750 Dtmpydz = Array<OneD, NekDouble>(nq);
751
752 for (int k = 0; k < m_spacedim; ++k)
753 {
754 Vmath::Vcopy(nq, &inarray[0][0], 1, &tmpx[0], 1);
755 Vmath::Vcopy(nq, &inarray[1][0], 1, &tmpy[0], 1);
756 Vmath::Vcopy(nq, &inarray[2][0], 1, &tmpz[0], 1);
757
758 m_fields[0]->PhysDeriv(0, tmpz, Dtmpzdx);
759 m_fields[0]->PhysDeriv(0, tmpy, Dtmpydx);
760 m_fields[0]->PhysDeriv(1, tmpx, Dtmpxdy);
761 m_fields[0]->PhysDeriv(1, tmpz, Dtmpzdy);
762 m_fields[0]->PhysDeriv(2, tmpx, Dtmpxdz);
763 m_fields[0]->PhysDeriv(2, tmpy, Dtmpydz);
764
765 Vmath::Vsub(nq, &Dtmpzdy[0], 1, &Dtmpydz[0], 1, &outarray[0][0], 1);
766 Vmath::Vsub(nq, &Dtmpxdz[0], 1, &Dtmpzdx[0], 1, &outarray[1][0], 1);
767 Vmath::Vsub(nq, &Dtmpydx[0], 1, &Dtmpxdy[0], 1, &outarray[2][0], 1);
768 }
769}
770
771Array<OneD, NekDouble> MMFSystem::CartesianToMovingframes(
772 const Array<OneD, const Array<OneD, NekDouble>> &uvec, unsigned int field)
773{
774 int nq = m_fields[0]->GetNpoints();
775
776 Array<OneD, NekDouble> outarray(nq, 0.0);
777
778 // new u0 = ( [u v] \cdot e^1 )/ |e^1|^2
779 Vmath::Vmul(nq, &m_movingframes[field][0], 1, &uvec[0][0], 1, &outarray[0],
780 1);
781 Vmath::Vvtvp(nq, &m_movingframes[field][nq], 1, &uvec[1][0], 1,
782 &outarray[0], 1, &outarray[0], 1);
783 Vmath::Vvtvp(nq, &m_movingframes[field][2 * nq], 1, &uvec[2][0], 1,
784 &outarray[0], 1, &outarray[0], 1);
785
786 return outarray;
787}
788
789// x = r \cos \theta \cos \varphi
790// y = r \cos \theta \sin \varphi
791// z = r \sin \theta
792void MMFSystem::CartesianToSpherical(const NekDouble x0j, const NekDouble x1j,
793 const NekDouble x2j, NekDouble &sin_varphi,
794 NekDouble &cos_varphi,
795 NekDouble &sin_theta, NekDouble &cos_theta)
796{
797 NekDouble radius;
798 NekDouble radxy;
799 NekDouble Tol = 0.0000001;
800
801 NekDouble m_Xscale = 1.0;
802 NekDouble m_Yscale = 1.0;
803 NekDouble m_Zscale = 1.0;
804
805 radius = sqrt(x0j * x0j / (m_Xscale * m_Xscale) +
806 x1j * x1j / (m_Yscale * m_Yscale) +
807 x2j * x2j / (m_Zscale * m_Zscale));
808 radxy = sqrt(x0j * x0j / (m_Xscale * m_Xscale) +
809 x1j * x1j / (m_Yscale * m_Yscale));
810
811 if (radxy > Tol)
812 {
813 sin_varphi = x1j / (radxy * m_Yscale);
814 cos_varphi = x0j / (radxy * m_Xscale);
815 }
816
817 else
818 {
819 sin_varphi = 0.0;
820 if (x2j > 0)
821 {
822 cos_varphi = 1.0;
823 }
824
825 else
826 {
827 cos_varphi = -1.0;
828 }
829 }
830
831 sin_theta = x2j / (radius * m_Zscale);
832 cos_theta = radxy / radius;
833}
834
835void MMFSystem::CopyBoundaryTrace(const Array<OneD, const NekDouble> &Fwd,
836 Array<OneD, NekDouble> &Bwd,
837 const BoundaryCopyType BDCopyType,
838 const int var, const std::string BDtype)
839{
840 int id1, id2, npts, nptselem, cnt = 0;
841 Array<OneD, NekDouble> Dirichlet, x0, x1, x2;
842
843 // loop over Boundary Regions
844 for (int n = 0; n < m_fields[var]->GetBndConditions().size(); ++n)
845 {
846 nptselem = m_fields[var]->GetBndCondExpansions()[n]->GetNpoints();
847
848 Dirichlet = Array<OneD, NekDouble>(nptselem);
849 x0 = Array<OneD, NekDouble>(nptselem);
850 x1 = Array<OneD, NekDouble>(nptselem);
851 x2 = Array<OneD, NekDouble>(nptselem);
852
853 if (BDCopyType == eDirichlet)
854 {
855 m_fields[var]->GetBndCondExpansions()[n]->GetCoords(x0, x1, x2);
856 LibUtilities::EquationSharedPtr ifunc =
857 m_session->GetFunction("BoundaryConditions", 0);
858 ifunc->Evaluate(x0, x1, x2, 0.0, Dirichlet);
859 }
860
861 for (int e = 0;
862 e < m_fields[var]->GetBndCondExpansions()[n]->GetExpSize(); ++e)
863 {
864 npts = m_fields[var]
865 ->GetBndCondExpansions()[n]
866 ->GetExp(e)
867 ->GetNumPoints(0);
868 id1 = m_fields[var]->GetBndCondExpansions()[n]->GetPhys_Offset(e);
869 id2 = m_fields[var]->GetTrace()->GetPhys_Offset(
870 m_fields[var]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt +
871 e));
872
873 if (m_fields[var]->GetBndConditions()[n]->GetUserDefined() ==
874 BDtype ||
875 BDtype == "NoUserDefined")
876 {
877 switch (BDCopyType)
878 {
879 case eDirichlet:
880 {
881 Vmath::Vcopy(npts, &Dirichlet[id1], 1, &Bwd[id2], 1);
882 }
883 break;
884
885 case eFwdEQBwd:
886 {
887 Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
888 }
889 break;
890
891 case eFwdEQNegBwd:
892 {
893 Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
894 Vmath::Neg(npts, &Bwd[id2], 1);
895 }
896 break;
897
898 default:
899 break;
900 }
901 }
902 }
903
904 cnt += m_fields[var]->GetBndCondExpansions()[n]->GetExpSize();
905 }
906}
907
908// Compute (e^[dir] \cdot ( n \times e^3)) * (imFwd * Fwd^3 + imBwd * Bwd^3 )
909void MMFSystem::ComputeNtimesFz(const int dir,
910 const Array<OneD, Array<OneD, NekDouble>> &Fwd,
911 const Array<OneD, Array<OneD, NekDouble>> &Bwd,
912 const Array<OneD, const NekDouble> &imFwd,
913 const Array<OneD, const NekDouble> &imBwd,
914 Array<OneD, NekDouble> &outarrayFwd,
915 Array<OneD, NekDouble> &outarrayBwd)
916
917{
918 int nTraceNumPoints = GetTraceTotPoints();
919
920 NekDouble tmpFwd, tmpBwd;
921 NekDouble Aver, ntimesz;
922
923 for (int i = 0; i < nTraceNumPoints; ++i)
924 {
925 Aver = 0.5 * (imFwd[i] + imBwd[i]);
926
927 tmpFwd = 0.0;
928 tmpBwd = 0.0;
929 for (int k = 0; k < m_spacedim; ++k)
930 {
931 ntimesz = 0.5 * (imFwd[i] * Fwd[2][i] + imBwd[i] * Bwd[2][i]);
932
933 tmpFwd +=
934 m_MFtraceFwd[dir][k][i] * m_ntimesMFFwd[2][k][i] * ntimesz;
935 tmpBwd +=
936 m_MFtraceBwd[dir][k][i] * m_ntimesMFBwd[2][k][i] * ntimesz;
937 }
938
939 outarrayFwd[i] = tmpFwd / Aver;
940 outarrayBwd[i] = tmpBwd / Aver;
941 }
942}
943
944// Compute e^3 \cdot ( n1e1 \times ( imFwd EFwd + imBwd EBwd ) ) / 2{{Y_i}}
945void MMFSystem::ComputeNtimesF12(const Array<OneD, Array<OneD, NekDouble>> &Fwd,
946 const Array<OneD, Array<OneD, NekDouble>> &Bwd,
947 const Array<OneD, const NekDouble> &im1Fwd,
948 const Array<OneD, const NekDouble> &im1Bwd,
949 const Array<OneD, const NekDouble> &im2Fwd,
950 const Array<OneD, const NekDouble> &im2Bwd,
951 Array<OneD, NekDouble> &outarrayFwd,
952 Array<OneD, NekDouble> &outarrayBwd)
953{
954 int nTraceNumPoints = GetTraceTotPoints();
955
956 NekDouble tmpFwd, tmpBwd, Aver1, Aver2, HFwdk, HBwdk;
957
958 Array<OneD, NekDouble> z1HAver(m_spacedim);
959 Array<OneD, NekDouble> z2HAver(m_spacedim);
960 Array<OneD, NekDouble> n1e1(m_spacedim);
961 Array<OneD, NekDouble> n2e2(m_spacedim);
962
963 Array<OneD, NekDouble> n1e1_times_z1HAver(m_spacedim);
964 Array<OneD, NekDouble> n2e2_times_z2HAver(m_spacedim);
965
966 for (int i = 0; i < nTraceNumPoints; ++i)
967 {
968 Aver1 = 0.5 * (im1Fwd[i] + im1Bwd[i]);
969 Aver2 = 0.5 * (im2Fwd[i] + im2Bwd[i]);
970
971 for (int k = 0; k < m_spacedim; k++)
972 {
973 // Compute \vec{HFwd} and \vec{HBwd}
974 HFwdk = Fwd[0][i] * m_MFtraceFwd[0][k][i] +
975 Fwd[1][i] * m_MFtraceFwd[1][k][i];
976 HBwdk = Bwd[0][i] * m_MFtraceBwd[0][k][i] +
977 Bwd[1][i] * m_MFtraceBwd[1][k][i];
978
979 // Compute z_i {{ \vec{H} }}
980 z1HAver[k] = 0.5 * (im1Fwd[i] * HFwdk + im1Bwd[i] * HBwdk);
981 z2HAver[k] = 0.5 * (im2Fwd[i] * HFwdk + im2Bwd[i] * HBwdk);
982
983 // Choose e^i for the one in anisotropy region
984 n1e1[k] = m_ncdotMFFwd[0][i] * m_MFtraceFwd[0][k][i];
985 n2e2[k] = m_ncdotMFFwd[1][i] * m_MFtraceFwd[1][k][i];
986 }
987
988 // Compute n1e1 \times z1HAver and n2e2 \times z2HAver
989 VectorCrossProd(n1e1, z1HAver, n1e1_times_z1HAver);
990 VectorCrossProd(n2e2, z2HAver, n2e2_times_z2HAver);
991
992 // e^3 \cdot ( n1e1 \times z1HAver + n2e2 \times z2HAver)
993 tmpFwd = 0.0;
994 tmpBwd = 0.0;
995 for (int k = 0; k < m_spacedim; k++)
996 {
997 tmpFwd += m_MFtraceFwd[2][k][i] * (n1e1_times_z1HAver[k] / Aver1 +
998 n2e2_times_z2HAver[k] / Aver2);
999 tmpBwd += m_MFtraceBwd[2][k][i] * (n1e1_times_z1HAver[k] / Aver1 +
1000 n2e2_times_z2HAver[k] / Aver2);
1001 }
1002
1003 outarrayFwd[i] = tmpFwd;
1004 outarrayBwd[i] = tmpBwd;
1005 }
1006}
1007
1008void MMFSystem::ComputeNtimestimesdFz(
1009 const int dir, const Array<OneD, Array<OneD, NekDouble>> &Fwd,
1010 const Array<OneD, Array<OneD, NekDouble>> &Bwd,
1011 const Array<OneD, const NekDouble> &imFwd,
1012 const Array<OneD, const NekDouble> &imBwd,
1013 Array<OneD, NekDouble> &outarrayFwd, Array<OneD, NekDouble> &outarrayBwd)
1014{
1015 int nTraceNumPoints = GetTraceTotPoints();
1016
1017 Array<OneD, NekDouble> dH(m_spacedim);
1018 Array<OneD, NekDouble> nFwd(m_spacedim);
1019
1020 Array<OneD, NekDouble> eiFwd(m_spacedim);
1021 Array<OneD, NekDouble> eiBwd(m_spacedim);
1022
1023 Array<OneD, NekDouble> eitimesdHFwd(m_spacedim);
1024 Array<OneD, NekDouble> eitimesdHBwd(m_spacedim);
1025
1026 Array<OneD, NekDouble> ntimeseitimesdHFwd(m_spacedim);
1027 Array<OneD, NekDouble> ntimeseitimesdHBwd(m_spacedim);
1028
1029 NekDouble Aver, HFwdk, HBwdk, tmpFwd, tmpBwd;
1030 for (int i = 0; i < nTraceNumPoints; ++i)
1031 {
1032 Aver = 0.5 * (imFwd[i] + imBwd[i]);
1033
1034 // Get [H]
1035 for (int k = 0; k < m_spacedim; k++)
1036 {
1037 HFwdk = Fwd[0][i] * m_MFtraceFwd[0][k][i] +
1038 Fwd[1][i] * m_MFtraceFwd[1][k][i];
1039 HBwdk = Bwd[0][i] * m_MFtraceBwd[0][k][i] +
1040 Bwd[1][i] * m_MFtraceBwd[1][k][i];
1041 dH[k] = HFwdk - HBwdk;
1042
1043 eiFwd[k] = m_MFtraceFwd[dir][k][i];
1044 eiBwd[k] = m_MFtraceBwd[dir][k][i];
1045
1046 nFwd[k] = m_traceNormals[k][i];
1047 }
1048
1049 // MFtraceFwd (MFtraceBwd) \times [H]
1050 // VectorCrossProd(eiFwd, dH, eitimesdHFwd);
1051 // VectorCrossProd(eiBwd, dH, eitimesdHBwd);
1052 VectorCrossProd(nFwd, dH, eitimesdHFwd);
1053 VectorCrossProd(nFwd, dH, eitimesdHBwd);
1054
1055 // n times eitimesdH
1056 VectorCrossProd(nFwd, eitimesdHFwd, ntimeseitimesdHFwd);
1057 VectorCrossProd(nFwd, eitimesdHBwd, ntimeseitimesdHBwd);
1058
1059 // MFtraceFwd \cdot ntimeseitimesdH
1060 tmpFwd = 0.0;
1061 tmpBwd = 0.0;
1062 for (int k = 0; k < m_spacedim; k++)
1063 {
1064 tmpFwd += eiFwd[k] * ntimeseitimesdHFwd[k];
1065 tmpBwd += eiBwd[k] * ntimeseitimesdHBwd[k];
1066 }
1067
1068 outarrayFwd[i] = 0.5 * m_alpha * tmpFwd / Aver;
1069 outarrayBwd[i] = 0.5 * m_alpha * tmpBwd / Aver;
1070 }
1071}
1072
1073// Compute - \alpha [E3] / ( {{im1}} + {{im2}} )
1074void MMFSystem::ComputeNtimestimesdF12(
1075 const Array<OneD, Array<OneD, NekDouble>> &Fwd,
1076 const Array<OneD, Array<OneD, NekDouble>> &Bwd,
1077 const Array<OneD, const NekDouble> &im1Fwd,
1078 const Array<OneD, const NekDouble> &im1Bwd,
1079 const Array<OneD, const NekDouble> &im2Fwd,
1080 const Array<OneD, const NekDouble> &im2Bwd,
1081 Array<OneD, NekDouble> &outarrayFwd, Array<OneD, NekDouble> &outarrayBwd)
1082{
1083 int nTraceNumPoints = GetTraceTotPoints();
1084
1085 Array<OneD, NekDouble> directFwd(nTraceNumPoints);
1086 Array<OneD, NekDouble> directBwd(nTraceNumPoints);
1087
1088 NekDouble Aver1, Aver2;
1089 for (int i = 0; i < nTraceNumPoints; ++i)
1090 {
1091 Aver1 = im1Fwd[i] + im1Bwd[i];
1092 Aver2 = im2Fwd[i] + im2Bwd[i];
1093
1094 outarrayFwd[i] =
1095 -m_alpha * (Fwd[2][i] - Bwd[2][i]) * (1.0 / Aver1 + 1.0 / Aver2);
1096 outarrayBwd[i] =
1097 -m_alpha * (Fwd[2][i] - Bwd[2][i]) * (1.0 / Aver1 + 1.0 / Aver2);
1098 }
1099}
1100
1101// dim = dimension, pol = polarization, 0 = TM, 1 = TE.
1102void MMFSystem::ComputeZimYim(Array<OneD, Array<OneD, NekDouble>> &epsvec,
1103 Array<OneD, Array<OneD, NekDouble>> &muvec)
1104
1105{
1106 int nTraceNumPoints = GetTraceNpoints();
1107
1108 switch (m_TestMaxwellType)
1109 {
1110 case eMaxwell1D:
1111 case eScatField1D:
1112 {
1113 Array<OneD, NekDouble> Fwdeps(nTraceNumPoints, 1.0);
1114 Array<OneD, NekDouble> Bwdeps(nTraceNumPoints, 1.0);
1115 Array<OneD, NekDouble> Fwdmu(nTraceNumPoints, 1.0);
1116 Array<OneD, NekDouble> Bwdmu(nTraceNumPoints, 1.0);
1117 m_fields[0]->GetFwdBwdTracePhys(epsvec[0], Fwdeps, Bwdeps);
1118 m_fields[0]->GetFwdBwdTracePhys(muvec[0], Fwdeps, Bwdeps);
1119
1120 CopyBoundaryTrace(Fwdeps, Bwdeps, eFwdEQBwd, 0);
1121 CopyBoundaryTrace(Fwdmu, Bwdmu, eFwdEQBwd, 1);
1122
1123 m_ZimFwd = Array<OneD, Array<OneD, NekDouble>>(1);
1124 m_ZimBwd = Array<OneD, Array<OneD, NekDouble>>(1);
1125 m_YimFwd = Array<OneD, Array<OneD, NekDouble>>(1);
1126 m_YimBwd = Array<OneD, Array<OneD, NekDouble>>(1);
1127
1128 m_ZimFwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1129 m_ZimBwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1130 m_YimFwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1131 m_YimBwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1132
1133 // ZimFwd = sqrt( muFwd / epsFwd), ZimBwd = sqrt( muBwd / epsBwd)
1134 for (int i = 0; i < nTraceNumPoints; ++i)
1135 {
1136 m_ZimFwd[0][i] = sqrt(Fwdmu[i] / Fwdeps[i]);
1137 m_ZimBwd[0][i] = sqrt(Bwdmu[i] / Bwdeps[i]);
1138
1139 m_YimFwd[0][i] = 1.0 / m_ZimFwd[0][i];
1140 m_YimBwd[0][i] = 1.0 / m_ZimBwd[0][i];
1141 }
1142
1143 std::cout << "*** ZimFwd = " << RootMeanSquare(m_ZimFwd[0])
1144 << ", ZimBwd = " << RootMeanSquare(m_ZimBwd[0])
1145 << ", YimFwd = " << RootMeanSquare(m_YimFwd[0])
1146 << ", YimBwd = " << RootMeanSquare(m_YimBwd[0])
1147 << std::endl;
1148 }
1149 break;
1150
1151 case eTestMaxwell2DPEC:
1152 case eTestMaxwell2DPECAVGFLUX:
1153 case eTestMaxwell2DPMC:
1154 case eScatField2D:
1155 case eTotField2D:
1156 case eMaxwellSphere:
1157 case eELF2DSurface:
1158 {
1159 m_ZimFwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1160 m_ZimBwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1161 m_YimFwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1162 m_YimBwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1163
1164 for (int j = 0; j < m_shapedim; ++j)
1165 {
1166 m_ZimFwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1167 m_ZimBwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1168 m_YimFwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1169 m_YimBwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1170 }
1171
1172 switch (m_PolType)
1173 {
1174 case eTransMagnetic:
1175 {
1176 Array<OneD, NekDouble> Fwdmu1(nTraceNumPoints);
1177 Array<OneD, NekDouble> Bwdmu1(nTraceNumPoints);
1178 Array<OneD, NekDouble> Fwdmu2(nTraceNumPoints);
1179 Array<OneD, NekDouble> Bwdmu2(nTraceNumPoints);
1180 Array<OneD, NekDouble> Fwdeps3(nTraceNumPoints);
1181 Array<OneD, NekDouble> Bwdeps3(nTraceNumPoints);
1182
1183 m_fields[0]->GetFwdBwdTracePhys(muvec[0], Fwdmu1, Bwdmu1);
1184 m_fields[0]->GetFwdBwdTracePhys(muvec[1], Fwdmu2, Bwdmu2);
1185 m_fields[0]->GetFwdBwdTracePhys(epsvec[2], Fwdeps3,
1186 Bwdeps3);
1187
1188 CopyBoundaryTrace(Fwdmu1, Bwdmu1, eFwdEQBwd, 0);
1189 CopyBoundaryTrace(Fwdmu2, Bwdmu2, eFwdEQBwd, 1);
1190 CopyBoundaryTrace(Fwdeps3, Bwdeps3, eFwdEQBwd, 2);
1191
1192 // ZimFwd = sqrt( muFwd / epsFwd), ZimBwd = sqrt( muBwd /
1193 // epsBwd)
1194 for (int i = 0; i < nTraceNumPoints; ++i)
1195 {
1196 m_ZimFwd[0][i] = sqrt(Fwdmu2[i] / Fwdeps3[i]);
1197 m_ZimBwd[0][i] = sqrt(Bwdmu2[i] / Bwdeps3[i]);
1198
1199 m_YimFwd[0][i] = 1.0 / m_ZimFwd[0][i];
1200 m_YimBwd[0][i] = 1.0 / m_ZimBwd[0][i];
1201
1202 m_ZimFwd[1][i] = sqrt(Fwdmu1[i] / Fwdeps3[i]);
1203 m_ZimBwd[1][i] = sqrt(Bwdmu1[i] / Bwdeps3[i]);
1204
1205 m_YimFwd[1][i] = 1.0 / m_ZimFwd[1][i];
1206 m_YimBwd[1][i] = 1.0 / m_ZimBwd[1][i];
1207 }
1208 }
1209 break; // eTransMagnetic
1210
1211 case eTransElectric:
1212 {
1213 Array<OneD, NekDouble> Fwdeps1(nTraceNumPoints);
1214 Array<OneD, NekDouble> Bwdeps1(nTraceNumPoints);
1215 Array<OneD, NekDouble> Fwdeps2(nTraceNumPoints);
1216 Array<OneD, NekDouble> Bwdeps2(nTraceNumPoints);
1217 Array<OneD, NekDouble> Fwdmu3(nTraceNumPoints);
1218 Array<OneD, NekDouble> Bwdmu3(nTraceNumPoints);
1219
1220 m_fields[0]->GetFwdBwdTracePhys(epsvec[0], Fwdeps1,
1221 Bwdeps1);
1222 m_fields[0]->GetFwdBwdTracePhys(epsvec[1], Fwdeps2,
1223 Bwdeps2);
1224 m_fields[0]->GetFwdBwdTracePhys(muvec[2], Fwdmu3, Bwdmu3);
1225
1226 CopyBoundaryTrace(Fwdeps1, Bwdeps1, eFwdEQBwd, 0);
1227 CopyBoundaryTrace(Fwdeps2, Bwdeps2, eFwdEQBwd, 1);
1228 CopyBoundaryTrace(Fwdmu3, Bwdmu3, eFwdEQBwd, 2);
1229
1230 for (int i = 0; i < nTraceNumPoints; ++i)
1231 {
1232 m_ZimFwd[0][i] = sqrt(Fwdmu3[i] / Fwdeps2[i]);
1233 m_ZimBwd[0][i] = sqrt(Bwdmu3[i] / Bwdeps2[i]);
1234
1235 m_YimFwd[0][i] = 1.0 / m_ZimFwd[0][i];
1236 m_YimBwd[0][i] = 1.0 / m_ZimBwd[0][i];
1237
1238 m_ZimFwd[1][i] = sqrt(Fwdmu3[i] / Fwdeps1[i]);
1239 m_ZimBwd[1][i] = sqrt(Bwdmu3[i] / Bwdeps1[i]);
1240
1241 m_YimFwd[1][i] = 1.0 / m_ZimFwd[1][i];
1242 m_YimBwd[1][i] = 1.0 / m_ZimBwd[1][i];
1243 }
1244 }
1245 break; // eTransELectric
1246
1247 default:
1248 break;
1249 } // PolType
1250
1251 std::cout << "*** ZimFwd0 = [ "
1252 << Vmath::Vmin(nTraceNumPoints, m_ZimFwd[0], 1) << " , "
1253 << Vmath::Vmax(nTraceNumPoints, m_ZimFwd[0], 1)
1254 << " ], ZimBwd0 = [ "
1255 << Vmath::Vmin(nTraceNumPoints, m_ZimBwd[0], 1) << " , "
1256 << Vmath::Vmax(nTraceNumPoints, m_ZimBwd[0], 1) << " ] "
1257 << std::endl;
1258 std::cout << "*** ZimFwd1 = [ "
1259 << Vmath::Vmin(nTraceNumPoints, m_ZimFwd[1], 1) << " , "
1260 << Vmath::Vmax(nTraceNumPoints, m_ZimFwd[1], 1)
1261 << " ], ZimBwd1 = [ "
1262 << Vmath::Vmin(nTraceNumPoints, m_ZimBwd[1], 1) << " , "
1263 << Vmath::Vmax(nTraceNumPoints, m_ZimBwd[1], 1) << " ] "
1264 << std::endl;
1265 }
1266 break; // eMaxwell2D
1267
1268 default:
1269 break;
1270 } // TestMaxwellType
1271}
1272
1273// m_dedxi_cdot_e[m][j][k][] = de^m / d \xi^j \cdot e^n
1274void MMFSystem::Computedemdxicdote()
1275{
1276 int MFdim = 3;
1277 int nq = GetTotPoints();
1278
1279 // Initialization
1280 m_dedxi_cdot_e =
1281 Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble>>>>(MFdim);
1282 for (int indm = 0; indm < MFdim; ++indm)
1283 {
1284 m_dedxi_cdot_e[indm] =
1285 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
1286 for (int indj = 0; indj < MFdim; ++indj)
1287 {
1288 m_dedxi_cdot_e[indm][indj] =
1289 Array<OneD, Array<OneD, NekDouble>>(MFdim);
1290 for (int indn = 0; indn < MFdim; ++indn)
1291 {
1292 m_dedxi_cdot_e[indm][indj][indn] =
1293 Array<OneD, NekDouble>(nq, 0.0);
1294 }
1295 }
1296 }
1297
1298 Array<OneD, NekDouble> tmp(nq);
1299 Array<OneD, NekDouble> tmpderiv(nq);
1300 Array<OneD, NekDouble> dedt;
1301 for (int indm = 0; indm < MFdim; ++indm)
1302 {
1303 for (int indj = 0; indj < MFdim; ++indj)
1304 {
1305 for (int indn = 0; indn < MFdim; ++indn)
1306 {
1307 dedt = Array<OneD, NekDouble>(nq, 0.0);
1308 for (int k = 0; k < m_spacedim; ++k)
1309 {
1310 // Compute d e^m / d \xi_j cdot e^n
1311 Vmath::Vcopy(nq, &m_movingframes[indm][k * nq], 1, &tmp[0],
1312 1);
1313 m_fields[0]->PhysDirectionalDeriv(m_movingframes[indj], tmp,
1314 tmpderiv);
1315
1316 Vmath::Vvtvp(nq, &tmpderiv[0], 1,
1317 &m_movingframes[indn][k * nq], 1, &dedt[0], 1,
1318 &dedt[0], 1);
1319 }
1320
1321 Vmath::Vcopy(nq, &dedt[0], 1,
1322 &m_dedxi_cdot_e[indm][indj][indn][0], 1);
1323 }
1324 }
1325 }
1326
1327 int indx = 0;
1328 std::cout << "*** m_dedxi_cdot_e[0]/dxi1 = ( "
1329 << RootMeanSquare(m_dedxi_cdot_e[indx][0][0]) << " , "
1330 << RootMeanSquare(m_dedxi_cdot_e[indx][0][1]) << " , "
1331 << RootMeanSquare(m_dedxi_cdot_e[indx][0][2]) << " )_1, "
1332 << std::endl;
1333 std::cout << "*** m_dedxi_cdot_e[0]/dxi2 = ( "
1334 << RootMeanSquare(m_dedxi_cdot_e[indx][1][0]) << " , "
1335 << RootMeanSquare(m_dedxi_cdot_e[indx][1][1]) << " , "
1336 << RootMeanSquare(m_dedxi_cdot_e[indx][1][2]) << " )_2 "
1337 << std::endl;
1338
1339 indx = 1;
1340 std::cout << "*** m_dedxi_cdot_e[1]/dxi1 = ( "
1341 << RootMeanSquare(m_dedxi_cdot_e[indx][0][0]) << " , "
1342 << RootMeanSquare(m_dedxi_cdot_e[indx][0][1]) << " , "
1343 << RootMeanSquare(m_dedxi_cdot_e[indx][0][2]) << " )_1, "
1344 << std::endl;
1345 std::cout << "*** m_dedxi_cdot_e[1]/dxi2 = ( "
1346 << RootMeanSquare(m_dedxi_cdot_e[indx][1][0]) << " , "
1347 << RootMeanSquare(m_dedxi_cdot_e[indx][1][1]) << " , "
1348 << RootMeanSquare(m_dedxi_cdot_e[indx][1][2]) << " )_2 "
1349 << std::endl;
1350
1351 indx = 2;
1352 std::cout << "*** m_dedxi_cdot_e[2]/dxi1 = ( "
1353 << RootMeanSquare(m_dedxi_cdot_e[indx][0][0]) << " , "
1354 << RootMeanSquare(m_dedxi_cdot_e[indx][0][1]) << " , "
1355 << RootMeanSquare(m_dedxi_cdot_e[indx][0][2]) << " )_1, "
1356 << std::endl;
1357 std::cout << "*** m_dedxi_cdot_e[2]/dxi2 = ( "
1358 << RootMeanSquare(m_dedxi_cdot_e[indx][1][0]) << " , "
1359 << RootMeanSquare(m_dedxi_cdot_e[indx][1][1]) << " , "
1360 << RootMeanSquare(m_dedxi_cdot_e[indx][1][2]) << " )_2 "
1361 << std::endl;
1362}
1363
1364void MMFSystem::AdddedtMaxwell(
1365 const Array<OneD, const Array<OneD, NekDouble>> &physarray,
1366 Array<OneD, Array<OneD, NekDouble>> &outarray)
1367{
1368 int nq = GetTotPoints();
1369
1370 // m_dedxi_cdot_e[m][j][n][] = de^m / d \xi^j \cdot e^n
1371 Array<OneD, NekDouble> dedtej(nq);
1372 Array<OneD, NekDouble> de1dtcdotej(nq);
1373 Array<OneD, NekDouble> de2dtcdotej(nq);
1374
1375 Array<OneD, NekDouble> normH(nq);
1376 Array<OneD, NekDouble> NormedH1(nq, 0.0);
1377 Array<OneD, NekDouble> NormednegH2(nq, 0.0);
1378
1379 Vmath::Vmul(nq, &physarray[0][0], 1, &physarray[0][0], 1, &normH[0], 1);
1380 Vmath::Vvtvp(nq, &physarray[1][0], 1, &physarray[1][0], 1, &normH[0], 1,
1381 &normH[0], 1);
1382 Vmath::Vsqrt(nq, normH, 1, normH, 1);
1383
1384 NekDouble Tol = 0.001;
1385 for (int i = 0; i < nq; ++i)
1386 {
1387 if (normH[i] > Tol)
1388 {
1389 NormedH1[i] = physarray[0][i] / normH[i];
1390 NormednegH2[i] = -1.0 * physarray[1][i] / normH[i];
1391 }
1392 }
1393
1394 for (int j = 0; j < m_shapedim; ++j)
1395 {
1396 // Compute de1 / dt \cdot ej = (-H2 de^1/d\xi1 \cdot e^j + H1 de^1/d\xi2
1397 // \cdot e^j) / sqrt{ H1^2 + H2^2 }
1398 Vmath::Vmul(nq, &NormednegH2[0], 1, &m_dedxi_cdot_e[0][0][j][0], 1,
1399 &de1dtcdotej[0], 1);
1400 Vmath::Vvtvp(nq, &NormedH1[0], 1, &m_dedxi_cdot_e[0][1][j][0], 1,
1401 &de1dtcdotej[0], 1, &de1dtcdotej[0], 1);
1402
1403 // Compute de2 / dt \cdot ej = (-H2 de2/d\xi1 \cdot e^j + H1 de2/d\xi2
1404 // \cdot e^j) / sqrt{ H1^2 + H2^2 }
1405 Vmath::Vmul(nq, &NormednegH2[0], 1, &m_dedxi_cdot_e[1][0][j][0], 1,
1406 &de2dtcdotej[0], 1);
1407 Vmath::Vvtvp(nq, &NormedH1[0], 1, &m_dedxi_cdot_e[1][1][j][0], 1,
1408 &de2dtcdotej[0], 1, &de2dtcdotej[0], 1);
1409
1410 // Add dedt component: (H1 (de1/dt) + H2 (de2/dt) ) \cdot ej
1411 Vmath::Vmul(nq, &physarray[0][0], 1, &de1dtcdotej[0], 1, &dedtej[0], 1);
1412 Vmath::Vvtvp(nq, &physarray[1][0], 1, &de2dtcdotej[0], 1, &dedtej[0], 1,
1413 &dedtej[0], 1);
1414
1415 Vmath::Neg(nq, dedtej, 1);
1416
1417 switch (m_PolType)
1418 {
1419 case SolverUtils::eTransMagnetic:
1420 {
1421 if (j == 0)
1422 {
1423 Vmath::Vmul(nq, m_muvec[0], 1, dedtej, 1, dedtej, 1);
1424 }
1425
1426 else if (j == 1)
1427 {
1428 Vmath::Vmul(nq, m_muvec[1], 1, dedtej, 1, dedtej, 1);
1429 }
1430 }
1431 break;
1432
1433 case SolverUtils::eTransElectric:
1434 {
1435 if (j == 0)
1436 {
1437 Vmath::Vmul(nq, m_epsvec[0], 1, dedtej, 1, dedtej, 1);
1438 }
1439
1440 else if (j == 1)
1441 {
1442 Vmath::Vmul(nq, m_epsvec[1], 1, dedtej, 1, dedtej, 1);
1443 }
1444 }
1445 break;
1446
1447 default:
1448 break;
1449 }
1450
1451 Vmath::Vadd(nq, &dedtej[0], 1, &outarray[j][0], 1, &outarray[j][0], 1);
1452 }
1453}
1454
1455void MMFSystem::UpwindMaxwellFlux1D(
1456 Array<OneD, Array<OneD, NekDouble>> &physfield,
1457 Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1458 Array<OneD, Array<OneD, NekDouble>> &numfluxBwd)
1459{
1460 int i;
1461 int nTraceNumPoints = GetTraceTotPoints();
1462 int nvar = 2;
1463
1464 // get temporary arrays
1465 Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1466 Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1467
1468 for (i = 0; i < nvar; ++i)
1469 {
1470 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1471 Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1472 }
1473
1474 // get the physical values at the trace from the dependent variables
1475 for (i = 0; i < nvar; ++i)
1476 {
1477 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1478 }
1479
1480 // E = 0 at the boundaries
1481 CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PEC");
1482
1483 // d H / d n = 0 at the boundaries
1484 CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC");
1485
1486 Array<OneD, NekDouble> dE(nTraceNumPoints);
1487 Array<OneD, NekDouble> dH(nTraceNumPoints);
1488
1489 Vmath::Vsub(nTraceNumPoints, &Fwd[0][0], 1, &Bwd[0][0], 1, &dE[0], 1);
1490 Vmath::Vsub(nTraceNumPoints, &Fwd[1][0], 1, &Bwd[1][0], 1, &dH[0], 1);
1491
1492 NekDouble nx, AverZ, AverY, AverZH, AverYE;
1493 for (i = 0; i < nTraceNumPoints; ++i)
1494 {
1495 nx = m_traceNormals[0][i];
1496 AverZ = m_ZimFwd[0][i] + m_ZimBwd[0][i];
1497 AverY = m_YimFwd[0][i] + m_YimBwd[0][i];
1498 AverZH = m_ZimFwd[0][i] * Fwd[1][i] + m_ZimBwd[0][i] * Bwd[1][i];
1499 AverYE = m_YimFwd[0][i] * Fwd[0][i] + m_YimBwd[0][i] * Bwd[0][i];
1500
1501 numfluxFwd[0][i] = nx / AverZ * (AverZH - nx * dE[i]);
1502 numfluxFwd[1][i] = nx / AverY * (AverYE - nx * dH[i]);
1503
1504 numfluxBwd[0][i] = nx / AverZ * (AverZH - nx * dE[i]);
1505 numfluxBwd[1][i] = nx / AverY * (AverYE - nx * dH[i]);
1506 }
1507}
1508
1509void MMFSystem::LaxFriedrichMaxwellFlux1D(
1510 Array<OneD, Array<OneD, NekDouble>> &physfield,
1511 Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1512 Array<OneD, Array<OneD, NekDouble>> &numfluxBwd)
1513{
1514 int i;
1515 int nTraceNumPoints = GetTraceTotPoints();
1516 int nvar = 2;
1517
1518 // get temporary arrays
1519 Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1520 Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1521
1522 for (i = 0; i < nvar; ++i)
1523 {
1524 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1525 Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1526 }
1527
1528 // get the physical values at the trace from the dependent variables
1529 for (i = 0; i < nvar; ++i)
1530 {
1531 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1532 }
1533
1534 // E = 0 at the boundaries
1535 CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PEC");
1536
1537 // d H / d n = 0 at the boundaries
1538 CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC");
1539
1540 Array<OneD, NekDouble> dE(nTraceNumPoints);
1541 Array<OneD, NekDouble> dH(nTraceNumPoints);
1542
1543 Vmath::Vsub(nTraceNumPoints, &Fwd[0][0], 1, &Bwd[0][0], 1, &dE[0], 1);
1544 Vmath::Vsub(nTraceNumPoints, &Fwd[1][0], 1, &Bwd[1][0], 1, &dH[0], 1);
1545
1546 NekDouble nx;
1547 for (i = 0; i < nTraceNumPoints; ++i)
1548 {
1549 nx = m_traceNormals[0][i];
1550 numfluxFwd[0][i] =
1551 0.5 * nx * ((Fwd[1][i] + Bwd[1][i]) - nx * m_YimFwd[0][i] * dE[i]);
1552 numfluxFwd[1][i] =
1553 0.5 * nx * ((Fwd[0][i] + Bwd[0][i]) - nx * m_ZimFwd[0][i] * dH[i]);
1554
1555 numfluxBwd[0][i] =
1556 0.5 * nx * ((Fwd[1][i] + Bwd[1][i]) - nx * m_YimFwd[0][i] * dE[i]);
1557 numfluxBwd[1][i] =
1558 0.5 * nx * ((Fwd[0][i] + Bwd[0][i]) - nx * m_ZimFwd[0][i] * dH[i]);
1559 }
1560}
1561
1562void MMFSystem::AverageMaxwellFlux1D(
1563 Array<OneD, Array<OneD, NekDouble>> &physfield,
1564 Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1565 Array<OneD, Array<OneD, NekDouble>> &numfluxBwd)
1566{
1567 int i;
1568 int nTraceNumPoints = GetTraceTotPoints();
1569 int nvar = 2;
1570
1571 // get temporary arrays
1572 Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1573 Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1574
1575 for (i = 0; i < nvar; ++i)
1576 {
1577 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1578 Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1579 }
1580
1581 // get the physical values at the trace from the dependent variables
1582 for (i = 0; i < nvar; ++i)
1583 {
1584 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1585 }
1586
1587 // E = 0 at the boundaries
1588 CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PEC");
1589
1590 // d H / d n = 0 at the boundaries
1591 CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC");
1592
1593 for (i = 0; i < nTraceNumPoints; ++i)
1594 {
1595 numfluxFwd[0][i] = 0.5 * m_traceNormals[0][i] * (Fwd[1][i] + Bwd[1][i]);
1596 numfluxFwd[1][i] = 0.5 * m_traceNormals[0][i] * (Fwd[0][i] + Bwd[0][i]);
1597
1598 numfluxBwd[0][i] = 0.5 * m_traceNormals[0][i] * (Fwd[1][i] + Bwd[1][i]);
1599 numfluxBwd[1][i] = 0.5 * m_traceNormals[0][i] * (Fwd[0][i] + Bwd[0][i]);
1600 }
1601}
1602
1603void MMFSystem::GetMaxwellFluxVector(
1604 const int var, const Array<OneD, const Array<OneD, NekDouble>> &physfield,
1605 Array<OneD, Array<OneD, NekDouble>> &flux)
1606{
1607 switch (m_TestMaxwellType)
1608 {
1609 case eMaxwell1D:
1610 case eScatField1D:
1611 {
1612 GetMaxwellFlux1D(var, physfield, flux);
1613 }
1614 break;
1615
1616 case eTestMaxwell2DPEC:
1617 case eTestMaxwell2DPECAVGFLUX:
1618 case eTestMaxwell2DPMC:
1619 case eScatField2D:
1620 case eTotField2D:
1621 case eMaxwellSphere:
1622 case eELF2DSurface:
1623 {
1624 GetMaxwellFlux2D(var, physfield, flux);
1625 }
1626 break;
1627
1628 default:
1629 break;
1630 }
1631}
1632
1633void MMFSystem::GetMaxwellFlux1D(
1634 const int var, const Array<OneD, const Array<OneD, NekDouble>> &physfield,
1635 Array<OneD, Array<OneD, NekDouble>> &flux)
1636{
1637 int nq = m_fields[0]->GetTotPoints();
1638
1639 switch (var)
1640 {
1641 case 0:
1642 {
1643 // H in flux 0
1644 Vmath::Vcopy(nq, physfield[1], 1, flux[0], 1);
1645
1646 // E in flux 1
1647 Vmath::Zero(nq, flux[1], 1);
1648 }
1649 break;
1650
1651 case 1:
1652 {
1653 // E in flux 0
1654 Vmath::Vcopy(nq, physfield[0], 1, flux[0], 1);
1655
1656 // H in flux 1
1657 Vmath::Zero(nq, flux[1], 1);
1658 }
1659 break;
1660 //----------------------------------------------------
1661
1662 default:
1663 break;
1664 }
1665}
1666
1667void MMFSystem::GetMaxwellFlux2D(
1668 const int var, const Array<OneD, const Array<OneD, NekDouble>> &physfield,
1669 Array<OneD, Array<OneD, NekDouble>> &flux)
1670{
1671 int nq = m_fields[0]->GetTotPoints();
1672
1673 NekDouble sign = 1.0;
1674 switch (m_PolType)
1675 {
1676 // TransMagnetic
1677 case 0:
1678 {
1679 sign = -1.0;
1680 }
1681 break;
1682
1683 // TransElectric
1684 case 1:
1685 {
1686 sign = 1.0;
1687 }
1688 break;
1689
1690 default:
1691 break;
1692 }
1693
1694 switch (var)
1695 {
1696 case 0:
1697 {
1698 // -Ez in flux 1
1699 Vmath::Smul(nq, sign, physfield[2], 1, flux[0], 1);
1700 Vmath::Zero(nq, flux[1], 1);
1701 }
1702 break;
1703
1704 case 1:
1705 {
1706 // Ez in flux 0
1707 Vmath::Zero(nq, flux[0], 1);
1708 Vmath::Smul(nq, -sign, physfield[2], 1, flux[1], 1);
1709 }
1710 break;
1711
1712 case 2:
1713 {
1714 Vmath::Smul(nq, sign, physfield[0], 1, flux[0], 1);
1715 Vmath::Smul(nq, -sign, physfield[1], 1, flux[1], 1);
1716 }
1717 break;
1718
1719 default:
1720 ASSERTL0(false, "GetFluxVector2D: illegal vector index");
1721 }
1722}
1723
1724void MMFSystem::NumericalMaxwellFlux(
1725 Array<OneD, Array<OneD, NekDouble>> &physfield,
1726 Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1727 Array<OneD, Array<OneD, NekDouble>> &numfluxBwd, const NekDouble time)
1728{
1729
1730 switch (m_TestMaxwellType)
1731 {
1732 case eMaxwell1D:
1733 case eScatField1D:
1734 {
1735 switch (m_upwindType)
1736 {
1737 case eAverage:
1738 {
1739 AverageMaxwellFlux1D(physfield, numfluxFwd, numfluxBwd);
1740 }
1741 break;
1742
1743 case eLaxFriedrich:
1744 {
1745 LaxFriedrichMaxwellFlux1D(physfield, numfluxFwd,
1746 numfluxBwd);
1747 }
1748 break;
1749
1750 case eUpwind:
1751 {
1752 UpwindMaxwellFlux1D(physfield, numfluxFwd, numfluxBwd);
1753 }
1754 break;
1755
1756 default:
1757 {
1758 ASSERTL0(false,
1759 "populate switch statement for upwind flux");
1760 }
1761 break; // upwindType
1762 }
1763 }
1764 break; // eMaxwell1D
1765
1766 case eTestMaxwell2DPEC:
1767 case eTestMaxwell2DPECAVGFLUX:
1768 case eTestMaxwell2DPMC:
1769 case eScatField2D:
1770 case eTotField2D:
1771 case eMaxwellSphere:
1772 case eELF2DSurface:
1773 {
1774 switch (m_PolType)
1775 {
1776 case eTransMagnetic:
1777 {
1778 NumericalMaxwellFluxTM(physfield, numfluxFwd, numfluxBwd,
1779 time);
1780 }
1781 break;
1782
1783 case eTransElectric:
1784 {
1785 NumericalMaxwellFluxTE(physfield, numfluxFwd, numfluxBwd,
1786 time);
1787 }
1788 break;
1789
1790 default:
1791 break;
1792 }
1793 }
1794 break; // eMaxwell2D
1795
1796 default:
1797 break;
1798 } // m_TestMaxwellType
1799}
1800
1801void MMFSystem::NumericalMaxwellFluxTM(
1802 Array<OneD, Array<OneD, NekDouble>> &physfield,
1803 Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1804 Array<OneD, Array<OneD, NekDouble>> &numfluxBwd, const NekDouble time)
1805{
1806 int nq = m_fields[0]->GetNpoints();
1807 int nTraceNumPoints = GetTraceTotPoints();
1808 int nvar = physfield.size();
1809
1810 // get temporary arrays
1811 Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1812 Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1813 for (int i = 0; i < nvar; ++i)
1814 {
1815 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1816 Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1817
1818 // get the physical values at the trace
1819 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1820 }
1821
1822 // E^|| = 0 at the PEC boundaries vs. H^|| = 0 at PMC boundaries
1823 Array<OneD, NekDouble> IncField(nq, 0.0);
1824 Array<OneD, NekDouble> IncFieldBwd(nTraceNumPoints, 0.0);
1825 Array<OneD, Array<OneD, NekDouble>> IncFieldFwd(m_spacedim);
1826 for (int i = 0; i < m_spacedim; ++i)
1827 {
1828 IncFieldFwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1829
1830 IncField = GetIncidentField(i, time);
1831 m_fields[0]->GetFwdBwdTracePhys(IncField, IncFieldFwd[i], IncFieldBwd);
1832
1833 Vmath::Svtvp(nTraceNumPoints, 2.0, &IncFieldFwd[i][0], 1, &Fwd[i][0], 1,
1834 &IncFieldFwd[i][0], 1);
1835 }
1836
1837 // Total Field Formulation
1838 CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0, "PEC_Forces");
1839 CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC_Forces");
1840 CopyBoundaryTrace(IncFieldFwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2,
1841 "PEC_Forces");
1842
1843 CopyBoundaryTrace(IncFieldFwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0,
1844 "PMC_Forces");
1845 CopyBoundaryTrace(IncFieldFwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1,
1846 "PMC_Forces");
1847 CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2, "PMC_Forces");
1848
1849 CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0, "PEC");
1850 CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC");
1851 CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2, "PEC");
1852
1853 CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PMC");
1854 CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1, "PMC");
1855 CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2, "PMC");
1856
1857 Array<OneD, NekDouble> e1Fwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1858 Array<OneD, NekDouble> e1Bwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1859 Array<OneD, NekDouble> e2Fwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1860 Array<OneD, NekDouble> e2Bwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1861 Array<OneD, NekDouble> e3Fwd_cdot_dEe3(nTraceNumPoints, 0.0);
1862 Array<OneD, NekDouble> e3Bwd_cdot_dEe3(nTraceNumPoints, 0.0);
1863
1864 // Compute numfluxFwd[dir] = (eFwd^[dir] \cdot n \times e^3) * (YimFwd *
1865 // EFwd^3 + YimBwd * EBwd^3 )
1866 ComputeNtimesFz(0, Fwd, Bwd, m_YimFwd[0], m_YimBwd[0], numfluxFwd[0],
1867 numfluxBwd[0]);
1868 ComputeNtimesFz(1, Fwd, Bwd, m_YimFwd[1], m_YimBwd[1], numfluxFwd[1],
1869 numfluxBwd[1]);
1870
1871 // Compute numfluxFwd[2] = eFwd^3 \cdot ( n1e1 \times ( ZimFwd HFwd + ZimBwd
1872 // HBwd ) ) / 2 {{Z_i}}
1873 ComputeNtimesF12(Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0], m_ZimFwd[1],
1874 m_ZimBwd[1], numfluxFwd[2], numfluxBwd[2]);
1875
1876 // Compute e1Fwd_cdot_ncrossdE = eFwd[dir] \cdot \alpha n \times n \times
1877 // [H] / 2 {{YimFwd}}
1878 ComputeNtimestimesdFz(0, Fwd, Bwd, m_YimFwd[0], m_YimBwd[0],
1879 e1Fwd_cdot_ncrossdH, e1Bwd_cdot_ncrossdH);
1880 ComputeNtimestimesdFz(1, Fwd, Bwd, m_YimFwd[1], m_YimBwd[1],
1881 e2Fwd_cdot_ncrossdH, e2Bwd_cdot_ncrossdH);
1882
1883 // Compute \alpha [E3] * ( 1/2{{Zim1}} + 1/2{{Zim2}} )
1884 ComputeNtimestimesdF12(Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0], m_ZimFwd[1],
1885 m_ZimBwd[1], e3Fwd_cdot_dEe3, e3Bwd_cdot_dEe3);
1886
1887 Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxFwd[0], 1, e1Fwd_cdot_ncrossdH,
1888 1, numfluxFwd[0], 1);
1889 Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxFwd[1], 1, e2Fwd_cdot_ncrossdH,
1890 1, numfluxFwd[1], 1);
1891 Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxFwd[2], 1, e3Fwd_cdot_dEe3, 1,
1892 numfluxFwd[2], 1);
1893
1894 Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxBwd[0], 1, e1Bwd_cdot_ncrossdH,
1895 1, numfluxBwd[0], 1);
1896 Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxBwd[1], 1, e2Bwd_cdot_ncrossdH,
1897 1, numfluxBwd[1], 1);
1898 Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxBwd[2], 1, e3Bwd_cdot_dEe3, 1,
1899 numfluxBwd[2], 1);
1900}
1901
1902void MMFSystem::NumericalMaxwellFluxTE(
1903 Array<OneD, Array<OneD, NekDouble>> &physfield,
1904 Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1905 Array<OneD, Array<OneD, NekDouble>> &numfluxBwd, const NekDouble time)
1906
1907{
1908 int nq = m_fields[0]->GetNpoints();
1909 int nTraceNumPoints = GetTraceTotPoints();
1910 int nvar = physfield.size();
1911
1912 // Get temporary arrays
1913 Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1914 Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1915 for (int i = 0; i < nvar; ++i)
1916 {
1917 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1918 Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1919
1920 // get the physical values at the trace
1921 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1922 }
1923
1924 // E = 0 at the PEC boundaries:
1925 Array<OneD, NekDouble> IncField(nq, 0.0);
1926 Array<OneD, NekDouble> IncFieldBwd(nTraceNumPoints, 0.0);
1927 Array<OneD, Array<OneD, NekDouble>> IncFieldFwd(m_spacedim);
1928 for (int i = 0; i < m_spacedim; ++i)
1929 {
1930 IncFieldFwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1931
1932 IncField = GetIncidentField(i, time);
1933 m_fields[0]->GetFwdBwdTracePhys(IncField, IncFieldFwd[i], IncFieldBwd);
1934
1935 Vmath::Svtvp(nTraceNumPoints, 2.0, &IncFieldFwd[i][0], 1, &Fwd[i][0], 1,
1936 &IncFieldFwd[i][0], 1);
1937 }
1938
1939 CopyBoundaryTrace(IncFieldFwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0,
1940 "PEC_Forces");
1941 CopyBoundaryTrace(IncFieldFwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1,
1942 "PEC_Forces");
1943 CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2, "PEC_Forces");
1944
1945 CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0, "PMC_Forces");
1946 CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PMC_Forces");
1947 CopyBoundaryTrace(IncFieldFwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2,
1948 "PMC_Forces");
1949
1950 CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PEC");
1951 CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1, "PEC");
1952 CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2, "PEC");
1953
1954 CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0, "PMC");
1955 CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PMC");
1956 CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2, "PMC");
1957
1958 Array<OneD, NekDouble> e1Fwd_cdot_ncrossdE(nTraceNumPoints);
1959 Array<OneD, NekDouble> e1Bwd_cdot_ncrossdE(nTraceNumPoints);
1960 Array<OneD, NekDouble> e2Fwd_cdot_ncrossdE(nTraceNumPoints);
1961 Array<OneD, NekDouble> e2Bwd_cdot_ncrossdE(nTraceNumPoints);
1962 Array<OneD, NekDouble> e3Fwd_cdot_dHe3(nTraceNumPoints);
1963 Array<OneD, NekDouble> e3Bwd_cdot_dHe3(nTraceNumPoints);
1964
1965 // Compute numfluxFwd[dir] = (eFwd^[dir] \cdot n \times e^3) * (ZimFwd *
1966 // HFwd^3 + ZimBwd * HBwd^3 )
1967 // Compute numfluxBwd[dir] = (eBwd^[dir] \cdot n \times e^3) * (ZimFwd *
1968 // HFwd^3 + ZimBwd * HBwd^3 )
1969 ComputeNtimesFz(0, Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0], numfluxFwd[0],
1970 numfluxBwd[0]);
1971 ComputeNtimesFz(1, Fwd, Bwd, m_ZimFwd[1], m_ZimBwd[1], numfluxFwd[1],
1972 numfluxBwd[1]);
1973
1974 // Compute numfluxFwd[2] = eFwd^3 \cdot ( n1e1 \times ( imFwd EFwd + imBwd
1975 // EBwd ) ) / 2 {{Y_i}}
1976 ComputeNtimesF12(Fwd, Bwd, m_YimFwd[0], m_YimBwd[0], m_YimFwd[1],
1977 m_YimBwd[1], numfluxFwd[2], numfluxBwd[2]);
1978
1979 // Compute e1Fwd_cdot_ncrossdE = eFwd[dir] \cdot \alpha n \times n \times
1980 // [E] / 2 {{ZimFwd}}
1981 ComputeNtimestimesdFz(0, Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0],
1982 e1Fwd_cdot_ncrossdE, e1Bwd_cdot_ncrossdE);
1983 ComputeNtimestimesdFz(1, Fwd, Bwd, m_ZimFwd[1], m_ZimBwd[1],
1984 e2Fwd_cdot_ncrossdE, e2Bwd_cdot_ncrossdE);
1985
1986 // Compute - \alpha [H3] * ( 1/2{{Yim1}} + 1/2{{Yim2}} )
1987 ComputeNtimestimesdF12(Fwd, Bwd, m_YimFwd[0], m_YimBwd[0], m_YimFwd[1],
1988 m_YimBwd[1], e3Fwd_cdot_dHe3, e3Bwd_cdot_dHe3);
1989
1990 Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxFwd[0], 1, e1Fwd_cdot_ncrossdE, 1,
1991 numfluxFwd[0], 1);
1992 Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxFwd[1], 1, e2Fwd_cdot_ncrossdE, 1,
1993 numfluxFwd[1], 1);
1994 Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxFwd[2], 1, e3Fwd_cdot_dHe3, 1,
1995 numfluxFwd[2], 1);
1996
1997 Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxBwd[0], 1, e1Bwd_cdot_ncrossdE, 1,
1998 numfluxBwd[0], 1);
1999 Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxBwd[1], 1, e2Bwd_cdot_ncrossdE, 1,
2000 numfluxBwd[1], 1);
2001 Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxBwd[2], 1, e3Bwd_cdot_dHe3, 1,
2002 numfluxBwd[2], 1);
2003}
2004
2005Array<OneD, NekDouble> MMFSystem::GetIncidentField(const int var,
2006 const NekDouble time)
2007{
2008 int nq = m_fields[0]->GetNpoints();
2009
2010 Array<OneD, NekDouble> x0(nq);
2011 Array<OneD, NekDouble> x1(nq);
2012 Array<OneD, NekDouble> x2(nq);
2013
2014 m_fields[0]->GetCoords(x0, x1, x2);
2015
2016 // GetSmoothFactor such that wave propages from the left to the object.
2017 // a = 0.1, ta = 1, f = 1.0./(1.0 + exp( -0.5.*(time-ta)/a ));
2018 Array<OneD, NekDouble> SmoothFactor(nq, 1.0);
2019 switch (m_SmoothFactor)
2020 {
2021 case 0:
2022 {
2023 for (int i = 0; i < nq; i++)
2024 {
2025 SmoothFactor[i] = 1.0 / (1.0 + exp(-1.0 * (time - 1.0) / 0.1));
2026 }
2027 }
2028 break;
2029
2030 case 1:
2031 {
2032 NekDouble xmin = Vmath::Vmin(nq, x0, 1);
2033 NekDouble xp;
2034 for (int i = 0; i < nq; i++)
2035 {
2036 xp = x0[i] - xmin - time - m_SFinit;
2037 if (xp > 0.0)
2038 {
2039 SmoothFactor[i] =
2040 2.0 / (1.0 + exp(0.5 * (sqrt(xp * xp) - 0.1)));
2041 }
2042
2043 else
2044 {
2045 SmoothFactor[i] = 1.0;
2046 }
2047 }
2048 }
2049 break;
2050
2051 default:
2052 break;
2053 }
2054
2055 // Generate a factor for smoothly increasing wave
2056 Array<OneD, NekDouble> F1(nq);
2057 Array<OneD, NekDouble> F2(nq);
2058 Array<OneD, NekDouble> F3(nq);
2059
2060 Array<OneD, NekDouble> dF1dt(nq);
2061 Array<OneD, NekDouble> dF2dt(nq);
2062 Array<OneD, NekDouble> dF3dt(nq);
2063
2064 Array<OneD, NekDouble> F1int(nq);
2065 Array<OneD, NekDouble> F2int(nq);
2066 Array<OneD, NekDouble> F3int(nq);
2067
2068 Array<OneD, NekDouble> outarray(nq);
2069
2070 switch (m_IncType)
2071 {
2072 case ePlaneWave:
2073 {
2074 NekDouble cs, sn;
2075 NekDouble e1y, e2y;
2076 switch (m_PolType)
2077 {
2078 case eTransMagnetic:
2079 {
2080 // H = 0 \hat{x} - e^{ikx} \hat{y}
2081 // E = e^{ikx} \hat{z}
2082 // outarray1 = Hr1inc
2083 // outarray2 = Hr2inc
2084 // outarray3 = Ezrinc
2085 for (int i = 0; i < nq; i++)
2086 {
2087 e1y = m_movingframes[0][nq + i];
2088 e2y = m_movingframes[1][nq + i];
2089
2090 cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2091 sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2092
2093 F1[i] = -1.0 * cs * e1y;
2094 F2[i] = -1.0 * cs * e2y;
2095 F3[i] = cs;
2096
2097 dF1dt[i] = -1.0 * m_Incfreq * sn * e1y;
2098 dF2dt[i] = -1.0 * m_Incfreq * sn * e2y;
2099 dF3dt[i] = 1.0 * m_Incfreq * sn;
2100
2101 F1int[i] = (1.0 / m_Incfreq) * sn * e1y;
2102 F2int[i] = (1.0 / m_Incfreq) * sn * e2y;
2103 F3int[i] = (-1.0 / m_Incfreq) * sn;
2104 }
2105 }
2106 break;
2107
2108 case eTransElectric:
2109 {
2110 // E = 0 \hat{x} + e^{ikx} \hat{y}
2111 // H = e^{ikx} \hat{z}
2112 // outarray1 = Er1inc
2113 // outarray2 = Er2inc
2114 // outarray3 = Hzrinc
2115 // outarray1 = Ei1inc
2116 // outarray2 = Ei2inc
2117 // outarray3 = Hziinc
2118 for (int i = 0; i < nq; i++)
2119 {
2120 e1y = m_movingframes[0][nq + i];
2121 e2y = m_movingframes[1][nq + i];
2122
2123 cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2124 sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2125
2126 F1[i] = cs * e1y;
2127 F2[i] = cs * e2y;
2128 F3[i] = cs;
2129
2130 dF1dt[i] = m_Incfreq * sn * e1y;
2131 dF2dt[i] = m_Incfreq * sn * e2y;
2132 dF3dt[i] = m_Incfreq * sn;
2133
2134 F1int[i] = (-1.0 / m_Incfreq) * sn * e1y;
2135 F2int[i] = (-1.0 / m_Incfreq) * sn * e2y;
2136 F3int[i] = (-1.0 / m_Incfreq) * sn;
2137 }
2138 }
2139 break;
2140
2141 default:
2142 break;
2143 }
2144 }
2145 break;
2146
2147 case ePlaneWaveImag:
2148 {
2149 NekDouble cs, sn;
2150 NekDouble e1y, e2y;
2151 switch (m_PolType)
2152 {
2153 case eTransMagnetic:
2154 {
2155 // H = 0 \hat{x} - e^{ikx} \hat{y}
2156 // E = e^{ikx} \hat{z}
2157 // outarray1 = Hr1inc
2158 // outarray2 = Hr2inc
2159 // outarray3 = Ezrinc
2160 for (int i = 0; i < nq; i++)
2161 {
2162 e1y = m_movingframes[0][nq + i];
2163 e2y = m_movingframes[1][nq + i];
2164
2165 cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2166 sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2167
2168 F1[i] = -1.0 * sn * e1y;
2169 F2[i] = -1.0 * sn * e2y;
2170 F3[i] = sn;
2171
2172 dF1dt[i] = m_Incfreq * cs * e1y;
2173 dF2dt[i] = m_Incfreq * cs * e2y;
2174 dF3dt[i] = -1.0 * m_Incfreq * cs;
2175
2176 F1int[i] = (-1.0 / m_Incfreq) * cs * e1y;
2177 F2int[i] = (-1.0 / m_Incfreq) * cs * e2y;
2178 F3int[i] = (1.0 / m_Incfreq) * cs;
2179 }
2180 }
2181 break;
2182
2183 case eTransElectric:
2184 {
2185 // E = 0 \hat{x} + e^{ikx} \hat{y}
2186 // H = e^{ikx} \hat{z}
2187 // outarray1 = Er1inc
2188 // outarray2 = Er2inc
2189 // outarray3 = Hzrinc
2190 // outarray1 = Ei1inc
2191 // outarray2 = Ei2inc
2192 // outarray3 = Hziinc
2193 for (int i = 0; i < nq; i++)
2194 {
2195 e1y = m_movingframes[0][nq + i];
2196 e2y = m_movingframes[1][nq + i];
2197
2198 cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2199 sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2200
2201 F1[i] = sn * e1y;
2202 F2[i] = sn * e2y;
2203 F3[i] = sn;
2204
2205 dF1dt[i] = -1.0 * m_Incfreq * cs * e1y;
2206 dF2dt[i] = -1.0 * m_Incfreq * cs * e2y;
2207 dF3dt[i] = -1.0 * m_Incfreq * cs;
2208
2209 F1int[i] = (1.0 / m_Incfreq) * cs * e1y;
2210 F2int[i] = (1.0 / m_Incfreq) * cs * e2y;
2211 F3int[i] = (1.0 / m_Incfreq) * cs;
2212 }
2213 }
2214 break;
2215
2216 default:
2217 break;
2218 }
2219 }
2220 break;
2221
2222 default:
2223 break;
2224 }
2225
2226 switch (var)
2227 {
2228 case 0:
2229 {
2230 outarray = F1;
2231 }
2232 break;
2233
2234 case 1:
2235 {
2236 outarray = F2;
2237 }
2238 break;
2239
2240 case 2:
2241 {
2242 outarray = F3;
2243 }
2244 break;
2245
2246 case 10:
2247 {
2248 outarray = dF1dt;
2249 }
2250 break;
2251
2252 case 11:
2253 {
2254 outarray = dF2dt;
2255 }
2256 break;
2257
2258 case 12:
2259 {
2260 outarray = dF3dt;
2261 }
2262 break;
2263
2264 case 20:
2265 {
2266 outarray = F1int;
2267 }
2268 break;
2269
2270 case 21:
2271 {
2272 outarray = F2int;
2273 }
2274 break;
2275
2276 case 22:
2277 {
2278 outarray = F3int;
2279 }
2280 break;
2281
2282 default:
2283 {
2284 Vmath::Zero(nq, outarray, 1);
2285 }
2286 break;
2287 }
2288
2289 return outarray;
2290}
2291
2292NekDouble MMFSystem::AvgInt(const Array<OneD, const NekDouble> &inarray)
2293{
2294 int nq = m_fields[0]->GetNpoints();
2295 Array<OneD, NekDouble> Ones(nq, 1.0);
2296
2297 if (inarray.size() != nq)
2298 {
2299 ASSERTL0(false, "AvgInt Error: Vector size is not correct");
2300 }
2301
2302 NekDouble jac = m_fields[0]->Integral(Ones);
2303
2304 return (m_fields[0]->Integral(inarray)) / jac;
2305}
2306
2307NekDouble MMFSystem::AvgAbsInt(const Array<OneD, const NekDouble> &inarray)
2308{
2309 int nq = m_fields[0]->GetNpoints();
2310 Array<OneD, NekDouble> Ones(nq, 1.0);
2311 Array<OneD, NekDouble> tmp(nq);
2312
2313 if (inarray.size() != nq)
2314 {
2315 ASSERTL0(false, "AvgAbsInt Error: Vector size is not correct");
2316 }
2317
2318 NekDouble jac = m_fields[0]->Integral(Ones);
2319
2320 Vmath::Vabs(nq, inarray, 1, tmp, 1);
2321 return (m_fields[0]->Integral(tmp)) / jac;
2322}
2323
2324NekDouble MMFSystem::AbsIntegral(const Array<OneD, const NekDouble> &inarray)
2325{
2326 int nq = m_fields[0]->GetNpoints();
2327 Array<OneD, NekDouble> tmp(nq);
2328
2329 if (inarray.size() != nq)
2330 {
2331 ASSERTL0(false, "AbsIntegral Error: Vector size is not correct");
2332 }
2333
2334 Vmath::Vabs(nq, inarray, 1, tmp, 1);
2335 return m_fields[0]->Integral(tmp);
2336}
2337
2338NekDouble MMFSystem::RootMeanSquare(const Array<OneD, const NekDouble> &inarray)
2339{
2340 int Ntot = inarray.size();
2341
2342 NekDouble reval = 0.0;
2343 for (int i = 0; i < Ntot; ++i)
2344 {
2345 reval = reval + inarray[i] * inarray[i];
2346 }
2347 reval = sqrt(reval / Ntot);
2348
2349 return reval;
2350}
2351
2352NekDouble MMFSystem::VectorAvgMagnitude(
2353 const Array<OneD, const Array<OneD, NekDouble>> &inarray)
2354{
2355 int nq = inarray[0].size();
2356
2357 Array<OneD, NekDouble> tmp(nq, 0.0);
2358 for (int k = 0; k < m_spacedim; k++)
2359 {
2360 Vmath::Vvtvp(nq, &inarray[k][0], 1, &inarray[k][0], 1, &tmp[0], 1,
2361 &tmp[0], 1);
2362 }
2363 Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
2364
2365 return RootMeanSquare(tmp);
2366}
2367
2368void MMFSystem::BubbleSort(Array<OneD, NekDouble> &refarray,
2369 Array<OneD, NekDouble> &sortarray)
2370{
2371 int nq = refarray.size();
2372
2373 bool swapped = true;
2374 int j = 0;
2375 NekDouble tmp;
2376
2377 while (swapped)
2378 {
2379 swapped = false;
2380 j++;
2381 for (int i = 0; i < nq - j; i++)
2382 {
2383 if (refarray[i] > refarray[i + 1])
2384 {
2385 tmp = refarray[i];
2386 refarray[i] = refarray[i + 1];
2387 refarray[i + 1] = tmp;
2388
2389 tmp = sortarray[i];
2390 sortarray[i] = sortarray[i + 1];
2391 sortarray[i + 1] = tmp;
2392
2393 swapped = true;
2394 }
2395 }
2396 }
2397}
2398
2399void MMFSystem::GramSchumitz(
2400 const Array<OneD, const Array<OneD, NekDouble>> &v1,
2401 const Array<OneD, const Array<OneD, NekDouble>> &v2,
2402 Array<OneD, Array<OneD, NekDouble>> &outarray, bool KeepTheMagnitude)
2403{
2404
2405 int nq = v1[0].size();
2406 Array<OneD, NekDouble> tmp(nq, 0.0);
2407 Array<OneD, NekDouble> mag(nq, 0.0);
2408
2409 for (int i = 0; i < m_spacedim; ++i)
2410 {
2411 // u2 = v2 - < u1 , v2 > ( u1 / < u1, u1 > )
2412 Vmath::Vvtvp(nq, &v1[i][0], 1, &v2[i][0], 1, &tmp[0], 1, &tmp[0], 1);
2413 Vmath::Vvtvp(nq, &v1[i][0], 1, &v1[i][0], 1, &mag[0], 1, &mag[0], 1);
2414 }
2415 Vmath::Vdiv(nq, &tmp[0], 1, &mag[0], 1, &tmp[0], 1);
2416 Vmath::Neg(nq, &tmp[0], 1);
2417
2418 // outarray = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
2419
2420 // u2 = v2 - < u1 , v2 > ( u1 / < u1, u1 > )
2421 for (int i = 0; i < m_spacedim; ++i)
2422 {
2423 outarray[i] = Array<OneD, NekDouble>(nq, 0.0);
2424 Vmath::Vvtvp(nq, &tmp[0], 1, &v1[i][0], 1, &v2[i][0], 1,
2425 &outarray[i][0], 1);
2426 }
2427
2428 if (KeepTheMagnitude)
2429 {
2430 Array<OneD, NekDouble> magorig(nq, 0.0);
2431 Array<OneD, NekDouble> magnew(nq, 0.0);
2432
2433 for (int i = 0; i < m_spacedim; ++i)
2434 {
2435 Vmath::Vmul(nq, &v2[0][0], 1, &v2[0][0], 1, &magorig[0], 1);
2436 Vmath::Vvtvp(nq, &v2[1][0], 1, &v2[1][0], 1, &magorig[0], 1,
2437 &magorig[0], 1);
2438 Vmath::Vvtvp(nq, &v2[2][0], 1, &v2[2][0], 1, &magorig[0], 1,
2439 &magorig[0], 1);
2440
2441 Vmath::Vmul(nq, &outarray[0][0], 1, &outarray[0][0], 1, &magnew[0],
2442 1);
2443 Vmath::Vvtvp(nq, &outarray[1][0], 1, &outarray[1][0], 1, &magnew[0],
2444 1, &magnew[0], 1);
2445 Vmath::Vvtvp(nq, &outarray[2][0], 1, &outarray[2][0], 1, &magnew[0],
2446 1, &magnew[0], 1);
2447 }
2448
2449 for (int i = 0; i < m_spacedim; ++i)
2450 {
2451 for (int j = 0; j < nq; ++j)
2452 {
2453 if (fabs(magnew[j]) > 0.000000001)
2454 {
2455 outarray[i][j] =
2456 outarray[i][j] * sqrt(magorig[j] / magnew[j]);
2457 }
2458 }
2459 }
2460 }
2461}
2462
2463void MMFSystem::v_GenerateSummary(SummaryList &s)
2464{
2465 int nq = m_fields[0]->GetNpoints();
2466 UnsteadySystem::v_GenerateSummary(s);
2467
2468 AddSummaryItem(s, "Total grids", nq);
2469 AddSummaryItem(s, "Shape Dimension", m_shapedim);
2470 AddSummaryItem(s, "Surface", SurfaceTypeMap[m_surfaceType]);
2471 if (m_surfaceType == eSphere)
2472 {
2473 NekDouble MeshError;
2474
2475 Array<OneD, NekDouble> x(nq);
2476 Array<OneD, NekDouble> y(nq);
2477 Array<OneD, NekDouble> z(nq);
2478 Array<OneD, NekDouble> rad(nq, 0.0);
2479
2480 m_fields[0]->GetCoords(x, y, z);
2481
2482 Vmath::Vvtvp(nq, x, 1, x, 1, rad, 1, rad, 1);
2483 Vmath::Vvtvp(nq, y, 1, y, 1, rad, 1, rad, 1);
2484 Vmath::Vvtvp(nq, z, 1, z, 1, rad, 1, rad, 1);
2485 Vmath::Vsqrt(nq, rad, 1, rad, 1);
2486
2487 Vmath::Sadd(nq, -1.0, rad, 1, rad, 1);
2488 Vmath::Vabs(nq, rad, 1, rad, 1);
2489
2490 MeshError = m_fields[0]->Integral(rad);
2491 SolverUtils::AddSummaryItem(s, "Mesh Error", MeshError);
2492 }
2493
2494 AddSummaryItem(s, "MMFdir", SpatialDomains::GeomMMFMap[m_MMFdir]);
2495 AddSummaryItem(s, "SmoothFactor", m_SmoothFactor);
2496 AddSummaryItem(s, "SFinit", m_SFinit);
2497
2498 if (fabs(m_alpha - 1.0) > 0.001)
2499 {
2500 AddSummaryItem(s, "Alpha", m_alpha);
2501 }
2502
2503 if (m_Incfreq > 0.0)
2504 {
2505 AddSummaryItem(s, "Incfreq", m_Incfreq);
2506 }
2507
2508 if (m_MMFdir == 4)
2509 {
2510 AddSummaryItem(s, "MMFCircAxisX", m_MMFfactors[0]);
2511 AddSummaryItem(s, "MMFCircAxisY", m_MMFfactors[1]);
2512 AddSummaryItem(s, "MMFCircCentreX", m_MMFfactors[2]);
2513 AddSummaryItem(s, "MMFCircCentreY", m_MMFfactors[3]);
2514 }
2515}
2516} // namespace Nektar::SolverUtils
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:47
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:396
Nektar::SolverUtils::EquationSystem::GetTraceNpoints
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Definition: EquationSystem.h:744
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:739
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:359
Nektar::SolverUtils::EquationSystem::GetExpSize
SOLVER_UTILS_EXPORT int GetExpSize()
Definition: EquationSystem.h:749
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:352
Nektar::SolverUtils::EquationSystem::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
Definition: EquationSystem.h:421
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:764
Nektar::SolverUtils::MMFSystem::LaxFriedrichMaxwellFlux1D
SOLVER_UTILS_EXPORT void LaxFriedrichMaxwellFlux1D(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxFwd, Array< OneD, Array< OneD, NekDouble > > &numfluxBwd)
Definition: MMFSystem.cpp:1509
Nektar::SolverUtils::MMFSystem::CheckMovingFrames
void CheckMovingFrames(const Array< OneD, const Array< OneD, NekDouble > > &movingframes)
Definition: MMFSystem.cpp:230
Nektar::SolverUtils::MMFSystem::ComputeNtimesF12
SOLVER_UTILS_EXPORT void ComputeNtimesF12(const Array< OneD, Array< OneD, NekDouble > > &Fwd, const Array< OneD, Array< OneD, NekDouble > > &Bwd, const Array< OneD, const NekDouble > &im1Fwd, const Array< OneD, const NekDouble > &im1Bwd, const Array< OneD, const NekDouble > &im2Fwd, const Array< OneD, const NekDouble > &im2Bwd, Array< OneD, NekDouble > &outarrayFwd, Array< OneD, NekDouble > &outarrayBwd)
Definition: MMFSystem.cpp:945
Nektar::SolverUtils::MMFSystem::VectorDotProd
SOLVER_UTILS_EXPORT void VectorDotProd(const Array< OneD, const Array< OneD, NekDouble > > &v1, const Array< OneD, const Array< OneD, NekDouble > > &v2, Array< OneD, NekDouble > &v3)
Definition: MMFSystem.cpp:673
Nektar::SolverUtils::MMFSystem::ComputeCurl
SOLVER_UTILS_EXPORT void ComputeCurl(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSystem.cpp:731
Nektar::SolverUtils::MMFSystem::m_dedxi_cdot_e
Array< OneD, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > > m_dedxi_cdot_e
Definition: MMFSystem.h:217
Nektar::SolverUtils::MMFSystem::m_ZimFwd
Array< OneD, Array< OneD, NekDouble > > m_ZimFwd
Definition: MMFSystem.h:204
Nektar::SolverUtils::MMFSystem::BubbleSort
SOLVER_UTILS_EXPORT void BubbleSort(Array< OneD, NekDouble > &refarray, Array< OneD, NekDouble > &sortarray)
Definition: MMFSystem.cpp:2368
Nektar::SolverUtils::MMFSystem::m_ntimes_ntimesMFBwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimes_ntimesMFBwd
Definition: MMFSystem.h:202
Nektar::SolverUtils::MMFSystem::GetMaxwellFluxVector
SOLVER_UTILS_EXPORT void GetMaxwellFluxVector(const int var, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
Definition: MMFSystem.cpp:1603
Nektar::SolverUtils::MMFSystem::m_DivMF
Array< OneD, Array< OneD, NekDouble > > m_DivMF
Definition: MMFSystem.h:192
Nektar::SolverUtils::MMFSystem::ComputeNtimesFz
SOLVER_UTILS_EXPORT void ComputeNtimesFz(const int dir, const Array< OneD, Array< OneD, NekDouble > > &Fwd, const Array< OneD, Array< OneD, NekDouble > > &Bwd, const Array< OneD, const NekDouble > &imFwd, const Array< OneD, const NekDouble > &imBwd, Array< OneD, NekDouble > &outarrayFwd, Array< OneD, NekDouble > &outarrayBwd)
Definition: MMFSystem.cpp:909
Nektar::SolverUtils::MMFSystem::NumericalMaxwellFluxTM
SOLVER_UTILS_EXPORT void NumericalMaxwellFluxTM(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxFwd, Array< OneD, Array< OneD, NekDouble > > &numfluxBwd, const NekDouble time)
Definition: MMFSystem.cpp:1801
Nektar::SolverUtils::MMFSystem::m_nperpcdotMFFwd
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFFwd
Definition: MMFSystem.h:189
Nektar::SolverUtils::MMFSystem::GetMaxwellFlux1D
SOLVER_UTILS_EXPORT void GetMaxwellFlux1D(const int var, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
Definition: MMFSystem.cpp:1633
Nektar::SolverUtils::MMFSystem::v_GenerateSummary
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Virtual function for generating summary information.
Definition: MMFSystem.cpp:2463
Nektar::SolverUtils::MMFSystem::AbsIntegral
SOLVER_UTILS_EXPORT NekDouble AbsIntegral(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2324
Nektar::SolverUtils::MMFSystem::DeriveCrossProductMF
SOLVER_UTILS_EXPORT void DeriveCrossProductMF(Array< OneD, Array< OneD, NekDouble > > &CrossProductMF)
Definition: MMFSystem.cpp:568
Nektar::SolverUtils::MMFSystem::GetMaxwellFlux2D
SOLVER_UTILS_EXPORT void GetMaxwellFlux2D(const int var, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
Definition: MMFSystem.cpp:1667
Nektar::SolverUtils::MMFSystem::m_ntimes_ntimesMFFwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimes_ntimesMFFwd
Definition: MMFSystem.h:201
Nektar::SolverUtils::MMFSystem::VectorCrossProd
SOLVER_UTILS_EXPORT void VectorCrossProd(const Array< OneD, const Array< OneD, NekDouble > > &v1, const Array< OneD, const Array< OneD, NekDouble > > &v2, Array< OneD, Array< OneD, NekDouble > > &v3)
Definition: MMFSystem.cpp:696
Nektar::SolverUtils::MMFSystem::NumericalMaxwellFlux
SOLVER_UTILS_EXPORT void NumericalMaxwellFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxFwd, Array< OneD, Array< OneD, NekDouble > > &numfluxBwd, const NekDouble time=0.0)
Definition: MMFSystem.cpp:1724
Nektar::SolverUtils::MMFSystem::AvgInt
SOLVER_UTILS_EXPORT NekDouble AvgInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2292
Nektar::SolverUtils::MMFSystem::AdddedtMaxwell
SOLVER_UTILS_EXPORT void AdddedtMaxwell(const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSystem.cpp:1364
Nektar::SolverUtils::MMFSystem::ComputeZimYim
SOLVER_UTILS_EXPORT void ComputeZimYim(Array< OneD, Array< OneD, NekDouble > > &epsvec, Array< OneD, Array< OneD, NekDouble > > &muvec)
Definition: MMFSystem.cpp:1102
Nektar::SolverUtils::MMFSystem::MMFSystem
SOLVER_UTILS_EXPORT MMFSystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Definition: MMFSystem.cpp:41
Nektar::SolverUtils::MMFSystem::AverageMaxwellFlux1D
SOLVER_UTILS_EXPORT void AverageMaxwellFlux1D(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxFwd, Array< OneD, Array< OneD, NekDouble > > &numfluxBwd)
Definition: MMFSystem.cpp:1562
Nektar::SolverUtils::MMFSystem::m_nperpcdotMFBwd
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFBwd
Definition: MMFSystem.h:190
Nektar::SolverUtils::MMFSystem::m_YimBwd
Array< OneD, Array< OneD, NekDouble > > m_YimBwd
Definition: MMFSystem.h:207
Nektar::SolverUtils::MMFSystem::~MMFSystem
SOLVER_UTILS_EXPORT ~MMFSystem() override
Definition: MMFSystem.cpp:47
Nektar::SolverUtils::MMFSystem::VectorAvgMagnitude
SOLVER_UTILS_EXPORT NekDouble VectorAvgMagnitude(const Array< OneD, const Array< OneD, NekDouble > > &inarray)
Definition: MMFSystem.cpp:2352
Nektar::SolverUtils::MMFSystem::m_muvec
Array< OneD, Array< OneD, NekDouble > > m_muvec
Definition: MMFSystem.h:210
Nektar::SolverUtils::MMFSystem::m_MMFdir
SpatialDomains::GeomMMF m_MMFdir
Definition: MMFSystem.h:219
Nektar::SolverUtils::MMFSystem::AvgAbsInt
SOLVER_UTILS_EXPORT NekDouble AvgAbsInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2307
Nektar::SolverUtils::MMFSystem::m_movingframes
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:183
Nektar::SolverUtils::MMFSystem::m_MMFfactors
Array< OneD, NekDouble > m_MMFfactors
Definition: MMFSystem.h:157
Nektar::SolverUtils::MMFSystem::m_ntimesMFBwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFBwd
Definition: MMFSystem.h:200
Nektar::SolverUtils::MMFSystem::m_MFtraceFwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceFwd
Definition: MMFSystem.h:196
Nektar::SolverUtils::MMFSystem::m_MFtraceBwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceBwd
Definition: MMFSystem.h:197
Nektar::SolverUtils::MMFSystem::GetIncidentField
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > GetIncidentField(const int var, const NekDouble time)
Definition: MMFSystem.cpp:2005
Nektar::SolverUtils::MMFSystem::SetUpMovingFrames
void SetUpMovingFrames(const Array< OneD, const Array< OneD, NekDouble > > &Anisotropy, const int TangentXelem)
Definition: MMFSystem.cpp:118
Nektar::SolverUtils::MMFSystem::m_ntimesMFFwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFFwd
Definition: MMFSystem.h:199
Nektar::SolverUtils::MMFSystem::ComputeDivCurlMF
SOLVER_UTILS_EXPORT void ComputeDivCurlMF()
Definition: MMFSystem.cpp:442
Nektar::SolverUtils::MMFSystem::GramSchumitz
SOLVER_UTILS_EXPORT void GramSchumitz(const Array< OneD, const Array< OneD, NekDouble > > &v1, const Array< OneD, const Array< OneD, NekDouble > > &v2, Array< OneD, Array< OneD, NekDouble > > &outarray, bool KeepTheMagnitude=true)
Definition: MMFSystem.cpp:2399
Nektar::SolverUtils::MMFSystem::m_ZimBwd
Array< OneD, Array< OneD, NekDouble > > m_ZimBwd
Definition: MMFSystem.h:205
Nektar::SolverUtils::MMFSystem::CartesianToMovingframes
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble > > &uvec, unsigned int field)
Definition: MMFSystem.cpp:771
Nektar::SolverUtils::MMFSystem::ComputeNtimestimesdF12
SOLVER_UTILS_EXPORT void ComputeNtimestimesdF12(const Array< OneD, Array< OneD, NekDouble > > &Fwd, const Array< OneD, Array< OneD, NekDouble > > &Bwd, const Array< OneD, const NekDouble > &im1Fwd, const Array< OneD, const NekDouble > &im1Bwd, const Array< OneD, const NekDouble > &im2Fwd, const Array< OneD, const NekDouble > &im2Bwd, Array< OneD, NekDouble > &outarrayFwd, Array< OneD, NekDouble > &outarrayBwd)
Definition: MMFSystem.cpp:1074
Nektar::SolverUtils::MMFSystem::CopyBoundaryTrace
SOLVER_UTILS_EXPORT void CopyBoundaryTrace(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const BoundaryCopyType BDCopyType, const int var=0, const std::string btype="NoUserDefined")
Definition: MMFSystem.cpp:835
Nektar::SolverUtils::MMFSystem::m_ncdotMFBwd
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
Definition: MMFSystem.h:187
Nektar::SolverUtils::MMFSystem::ComputeMFtrace
SOLVER_UTILS_EXPORT void ComputeMFtrace()
Definition: MMFSystem.cpp:519
Nektar::SolverUtils::MMFSystem::m_YimFwd
Array< OneD, Array< OneD, NekDouble > > m_YimFwd
Definition: MMFSystem.h:206
Nektar::SolverUtils::MMFSystem::m_epsvec
Array< OneD, Array< OneD, NekDouble > > m_epsvec
Definition: MMFSystem.h:209
Nektar::SolverUtils::MMFSystem::UpwindMaxwellFlux1D
SOLVER_UTILS_EXPORT void UpwindMaxwellFlux1D(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxFwd, Array< OneD, Array< OneD, NekDouble > > &numfluxBwd)
Definition: MMFSystem.cpp:1455
Nektar::SolverUtils::MMFSystem::ComputencdotMF
SOLVER_UTILS_EXPORT void ComputencdotMF()
Definition: MMFSystem.cpp:317
Nektar::SolverUtils::MMFSystem::CartesianToSpherical
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
Definition: MMFSystem.cpp:792
Nektar::SolverUtils::MMFSystem::ComputeNtimestimesdFz
SOLVER_UTILS_EXPORT void ComputeNtimestimesdFz(const int dir, const Array< OneD, Array< OneD, NekDouble > > &Fwd, const Array< OneD, Array< OneD, NekDouble > > &Bwd, const Array< OneD, const NekDouble > &imFwd, const Array< OneD, const NekDouble > &imBwd, Array< OneD, NekDouble > &outarrayFwd, Array< OneD, NekDouble > &outarrayBwd)
Definition: MMFSystem.cpp:1008
Nektar::SolverUtils::MMFSystem::NumericalMaxwellFluxTE
SOLVER_UTILS_EXPORT void NumericalMaxwellFluxTE(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxFwd, Array< OneD, Array< OneD, NekDouble > > &numfluxBwd, const NekDouble time)
Definition: MMFSystem.cpp:1902
Nektar::SolverUtils::MMFSystem::ComputeNtimesMF
SOLVER_UTILS_EXPORT void ComputeNtimesMF()
Definition: MMFSystem.cpp:615
Nektar::SolverUtils::MMFSystem::RootMeanSquare
SOLVER_UTILS_EXPORT NekDouble RootMeanSquare(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2338
Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
Definition: MMFSystem.h:186
Nektar::SolverUtils::MMFSystem::MMFInitObject
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble > > &Anisotropy, const int TangentXelem=-1)
Definition: MMFSystem.cpp:51
Nektar::SolverUtils::MMFSystem::m_TestMaxwellType
TestMaxwellType m_TestMaxwellType
Definition: MMFSystem.h:153
Nektar::SolverUtils::MMFSystem::m_CurlMF
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
Definition: MMFSystem.h:193
Nektar::SolverUtils::MMFSystem::Computedemdxicdote
SOLVER_UTILS_EXPORT void Computedemdxicdote()
Definition: MMFSystem.cpp:1274
Nektar::SolverUtils::UnsteadySystem
Base class for unsteady solvers.
Definition: UnsteadySystem.h:46
Nektar::SolverUtils::UnsteadySystem::v_GenerateSummary
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
Definition: UnsteadySystem.cpp:606
Nektar::LibUtilities::rad
static NekDouble rad(NekDouble x, NekDouble y)
Definition: Interpreter/Interpreter.cpp:76
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:113
Nektar::LibUtilities::EquationSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:125
Nektar::SolverUtils::UpwindTypeMap
const char *const UpwindTypeMap[]
Definition: MMFSystem.h:86
Nektar::SolverUtils::SummaryList
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
Nektar::SolverUtils::SIZE_UpwindType
@ SIZE_UpwindType
Length of enum list.
Definition: MMFSystem.h:83
Nektar::SolverUtils::eAverage
@ eAverage
averaged (or centred) flux
Definition: MMFSystem.h:77
Nektar::SolverUtils::eLaxFriedrich
@ eLaxFriedrich
Lax-Friedrich flux.
Definition: MMFSystem.h:78
Nektar::SolverUtils::AddSummaryItem
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:47
Nektar::SolverUtils::SurfaceTypeMap
const char *const SurfaceTypeMap[]
Definition: MMFSystem.h:54
Nektar::SpatialDomains::GeomMMFMap
const char *const GeomMMFMap[]
Session file names associated with tangent principle directions.
Definition: SpatialDomains.hpp:87
Nektar::SpatialDomains::eLOCAL
@ eLOCAL
No Principal direction.
Definition: SpatialDomains.hpp:80
Nektar::SpatialDomains::eTangentIrregular
@ eTangentIrregular
Circular around the centre of domain.
Definition: SpatialDomains.hpp:78
Nektar::SpatialDomains::eTangentX
@ eTangentX
X coordinate direction.
Definition: SpatialDomains.hpp:73
Nektar::SpatialDomains::eTangentCircular
@ eTangentCircular
Circular around the centre of domain.
Definition: SpatialDomains.hpp:77
Nektar::SpatialDomains::eTangentNonconvex
@ eTangentNonconvex
Circular around the centre of domain.
Definition: SpatialDomains.hpp:79
Nektar::SpatialDomains::eTangentXY
@ eTangentXY
XY direction.
Definition: SpatialDomains.hpp:75
Nektar::SpatialDomains::eTangentZ
@ eTangentZ
Z coordinate direction.
Definition: SpatialDomains.hpp:76
Nektar::SpatialDomains::eTangentY
@ eTangentY
Y coordinate direction.
Definition: SpatialDomains.hpp:74
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
Nektar::UnitTests::z
std::vector< double > z(NPUPPER)
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvp (scalar times vector plus vector): z = alpha*x + y.
Definition: Vmath.hpp:396
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.hpp:352
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
Vmath::Vmin
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.hpp:725
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.hpp:366
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.hpp:608
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
Vmath::Vvtvm
void Vvtvm(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvm (vector times vector minus vector): z = w*x - y
Definition: Vmath.hpp:381
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
Vmath::Imax
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.hpp:623
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.hpp:126
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.hpp:194
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.hpp:644
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.hpp:220
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294
Generated on Wed Jan 31 2024 09:35:56 for Nektar++ by doxygen 1.9.4