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