Nektar++
AdvectionFR.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: AdvectionFR.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: FR advection class.
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <LibUtilities/Foundations/Basis.h>
36#include <LibUtilities/Foundations/ManagerAccess.h>
37#include <LibUtilities/Foundations/Points.h>
38#include <LibUtilities/Polylib/Polylib.h>
39#include <LocalRegions/Expansion1D.h>
40#include <LocalRegions/Expansion2D.h>
41#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
42#include <MultiRegions/DisContField.h>
43#include <SolverUtils/Advection/AdvectionFR.h>
44#include <StdRegions/StdSegExp.h>
45
46#include <iomanip>
47#include <iostream>
48
49using namespace std;
50
51namespace Nektar::SolverUtils
52{
53std::string AdvectionFR::type[] = {
54 GetAdvectionFactory().RegisterCreatorFunction("FRDG", AdvectionFR::create),
55 GetAdvectionFactory().RegisterCreatorFunction("FRSD", AdvectionFR::create),
56 GetAdvectionFactory().RegisterCreatorFunction("FRHU", AdvectionFR::create),
57 GetAdvectionFactory().RegisterCreatorFunction("FRcmin",
58 AdvectionFR::create),
59 GetAdvectionFactory().RegisterCreatorFunction("FRcinf",
60 AdvectionFR::create)};
61
62/**
63 * @brief AdvectionFR uses the Flux Reconstruction (FR) approach to
64 * compute the advection term. The implementation is only for segments,
65 * quadrilaterals and hexahedra at the moment.
66 *
67 * \todo Extension to triangles, tetrahedra and other shapes.
68 * (Long term objective)
69 */
70AdvectionFR::AdvectionFR(std::string advType) : m_advType(advType)
71{
72}
73
74/**
75 * @brief Initiliase AdvectionFR objects and store them before starting
76 * the time-stepping.
77 *
78 * This routine calls the functions #SetupMetrics and
79 * #SetupCFunctions to initialise the objects needed by AdvectionFR.
80 *
81 * @param pSession Pointer to session reader.
82 * @param pFields Pointer to fields.
83 */
84void AdvectionFR::v_InitObject(
85 LibUtilities::SessionReaderSharedPtr pSession,
86 Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
87{
88 WARNINGL0(false, "FR is deprecated, use DG instead");
89
90 Advection::v_InitObject(pSession, pFields);
91 SetupMetrics(pSession, pFields);
92 SetupCFunctions(pSession, pFields);
93}
94
95/**
96 * @brief Setup the metric terms to compute the contravariant
97 * fluxes. (i.e. this special metric terms transform the fluxes
98 * at the interfaces of each element from the physical space to
99 * the standard space).
100 *
101 * This routine calls the function #GetEdgeQFactors to compute and
102 * store the metric factors following an anticlockwise conventions
103 * along the edges/faces of each element. Note: for 1D problem
104 * the transformation is not needed.
105 *
106 * @param pSession Pointer to session reader.
107 * @param pFields Pointer to fields.
108 *
109 * \todo Add the metric terms for 3D Hexahedra.
110 */
111void AdvectionFR::SetupMetrics(
112 [[maybe_unused]] LibUtilities::SessionReaderSharedPtr pSession,
113 Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
114{
115 int i, n;
116 int nquad0, nquad1;
117 int phys_offset;
118 int nLocalSolutionPts;
119 int nElements = pFields[0]->GetExpSize();
120 int nDimensions = pFields[0]->GetCoordim(0);
121 int nSolutionPts = pFields[0]->GetTotPoints();
122 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
123
124 Array<TwoD, const NekDouble> gmat;
125 Array<OneD, const NekDouble> jac;
126
127 m_jac = Array<OneD, NekDouble>(nSolutionPts);
128
129 Array<OneD, NekDouble> auxArray1;
130 Array<OneD, LibUtilities::BasisSharedPtr> base;
131 LibUtilities::PointsKeyVector ptsKeys;
132
133 switch (nDimensions)
134 {
135 case 1:
136 {
137 for (n = 0; n < nElements; ++n)
138 {
139 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
140 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
141 phys_offset = pFields[0]->GetPhys_Offset(n);
142 jac = pFields[0]
143 ->GetExp(n)
144 ->as<LocalRegions::Expansion1D>()
145 ->GetGeom1D()
146 ->GetMetricInfo()
147 ->GetJac(ptsKeys);
148 for (i = 0; i < nLocalSolutionPts; ++i)
149 {
150 m_jac[i + phys_offset] = jac[0];
151 }
152 }
153 break;
154 }
155 case 2:
156 {
157 m_gmat = Array<OneD, Array<OneD, NekDouble>>(4);
158 m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
159 m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
160 m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
161 m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
162
163 m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble>>(nElements);
164 m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
165 m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
166 m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
167
168 LibUtilities::PointsKeyVector ptsKeys;
169
170 for (n = 0; n < nElements; ++n)
171 {
172 base = pFields[0]->GetExp(n)->GetBase();
173 nquad0 = base[0]->GetNumPoints();
174 nquad1 = base[1]->GetNumPoints();
175
176 m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
177 m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
178 m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
179 m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
180
181 // Extract the Q factors at each edge point
182 pFields[0]->GetExp(n)->GetTraceQFactors(0, auxArray1 =
183 m_Q2D_e0[n]);
184 pFields[0]->GetExp(n)->GetTraceQFactors(1, auxArray1 =
185 m_Q2D_e1[n]);
186 pFields[0]->GetExp(n)->GetTraceQFactors(2, auxArray1 =
187 m_Q2D_e2[n]);
188 pFields[0]->GetExp(n)->GetTraceQFactors(3, auxArray1 =
189 m_Q2D_e3[n]);
190
191 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
192 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
193 phys_offset = pFields[0]->GetPhys_Offset(n);
194
195 jac = pFields[0]
196 ->GetExp(n)
197 ->as<LocalRegions::Expansion2D>()
198 ->GetGeom2D()
199 ->GetMetricInfo()
200 ->GetJac(ptsKeys);
201 gmat = pFields[0]
202 ->GetExp(n)
203 ->as<LocalRegions::Expansion2D>()
204 ->GetGeom2D()
205 ->GetMetricInfo()
206 ->GetDerivFactors(ptsKeys);
207
208 if (pFields[0]
209 ->GetExp(n)
210 ->as<LocalRegions::Expansion2D>()
211 ->GetGeom2D()
212 ->GetMetricInfo()
213 ->GetGtype() == SpatialDomains::eDeformed)
214 {
215 for (i = 0; i < nLocalSolutionPts; ++i)
216 {
217 m_jac[i + phys_offset] = jac[i];
218 m_gmat[0][i + phys_offset] = gmat[0][i];
219 m_gmat[1][i + phys_offset] = gmat[1][i];
220 m_gmat[2][i + phys_offset] = gmat[2][i];
221 m_gmat[3][i + phys_offset] = gmat[3][i];
222 }
223 }
224 else
225 {
226 for (i = 0; i < nLocalSolutionPts; ++i)
227 {
228 m_jac[i + phys_offset] = jac[0];
229 m_gmat[0][i + phys_offset] = gmat[0][0];
230 m_gmat[1][i + phys_offset] = gmat[1][0];
231 m_gmat[2][i + phys_offset] = gmat[2][0];
232 m_gmat[3][i + phys_offset] = gmat[3][0];
233 }
234 }
235 }
236
237 m_traceNormals = Array<OneD, Array<OneD, NekDouble>>(nDimensions);
238 for (i = 0; i < nDimensions; ++i)
239 {
240 m_traceNormals[i] = Array<OneD, NekDouble>(nTracePts);
241 }
242 pFields[0]->GetTrace()->GetNormals(m_traceNormals);
243
244 break;
245 }
246 case 3:
247 {
248 ASSERTL0(false, "3DFR Metric terms not implemented yet");
249 break;
250 }
251 default:
252 {
253 ASSERTL0(false, "Expansion dimension not recognised");
254 break;
255 }
256 }
257}
258
259/**
260 * @brief Setup the derivatives of the correction functions. For more
261 * details see J Sci Comput (2011) 47: 50–72
262 *
263 * This routine calls 3 different bases:
264 * #eDG_DG_Left - #eDG_DG_Left which recovers a nodal DG scheme,
265 * #eDG_SD_Left - #eDG_SD_Left which recovers the SD scheme,
266 * #eDG_HU_Left - #eDG_HU_Left which recovers the Huynh scheme.
267 * The values of the derivatives of the correction function are then
268 * stored into global variables and reused into the functions
269 * #DivCFlux_1D, #DivCFlux_2D, #DivCFlux_3D to compute the
270 * the divergence of the correction flux for 1D, 2D or 3D problems.
271 *
272 * @param pSession Pointer to session reader.
273 * @param pFields Pointer to fields.
274 */
275void AdvectionFR::SetupCFunctions(
276 [[maybe_unused]] LibUtilities::SessionReaderSharedPtr pSession,
277 Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
278{
279 int i, n;
280 NekDouble c0 = 0.0;
281 NekDouble c1 = 0.0;
282 NekDouble c2 = 0.0;
283 int nquad0, nquad1, nquad2;
284 int nmodes0, nmodes1, nmodes2;
285 Array<OneD, LibUtilities::BasisSharedPtr> base;
286
287 int nElements = pFields[0]->GetExpSize();
288 int nDimensions = pFields[0]->GetCoordim(0);
289
290 switch (nDimensions)
291 {
292 case 1:
293 {
294 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
295 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
296
297 for (n = 0; n < nElements; ++n)
298 {
299 base = pFields[0]->GetExp(n)->GetBase();
300 nquad0 = base[0]->GetNumPoints();
301 nmodes0 = base[0]->GetNumModes();
302 Array<OneD, const NekDouble> z0;
303 Array<OneD, const NekDouble> w0;
304
305 base[0]->GetZW(z0, w0);
306
307 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
308 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
309
310 // Auxiliary vectors to build up the auxiliary
311 // derivatives of the Legendre polynomials
312 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
313 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
314 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
315
316 // Degree of the correction functions
317 int p0 = nmodes0 - 1;
318
319 // Function sign to compute left correction function
320 NekDouble sign0 = pow(-1.0, p0);
321
322 // Factorial factor to build the scheme
323 NekDouble ap0 =
324 std::tgamma(2 * p0 + 1) /
325 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
326
327 // Scalar parameter which recovers the FR schemes
328 if (m_advType == "FRDG")
329 {
330 c0 = 0.0;
331 }
332 else if (m_advType == "FRSD")
333 {
334 c0 = 2.0 * p0 /
335 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
336 (ap0 * std::tgamma(p0 + 1)) *
337 (ap0 * std::tgamma(p0 + 1)));
338 }
339 else if (m_advType == "FRHU")
340 {
341 c0 = 2.0 * (p0 + 1.0) /
342 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
343 (ap0 * std::tgamma(p0 + 1)));
344 }
345 else if (m_advType == "FRcmin")
346 {
347 c0 =
348 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
349 (ap0 * std::tgamma(p0 + 1)));
350 }
351 else if (m_advType == "FRcinf")
352 {
353 c0 = 10000000000000000.0;
354 }
355
356 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
357 (ap0 * std::tgamma(p0 + 1)) *
358 (ap0 * std::tgamma(p0 + 1));
359
360 NekDouble overeta0 = 1.0 / (1.0 + etap0);
361
362 // Derivative of the Legendre polynomials
363 // dLp = derivative of the Legendre polynomial of order p
364 // dLpp = derivative of the Legendre polynomial of order p+1
365 // dLpm = derivative of the Legendre polynomial of order p-1
366 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
367 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
368 0.0);
369 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
370 0.0);
371
372 // Building the DG_c_Left
373 for (i = 0; i < nquad0; ++i)
374 {
375 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
376 m_dGL_xi1[n][i] += dLpp0[i];
377 m_dGL_xi1[n][i] *= overeta0;
378 m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
379 m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
380 }
381
382 // Building the DG_c_Right
383 for (i = 0; i < nquad0; ++i)
384 {
385 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
386 m_dGR_xi1[n][i] += dLpp0[i];
387 m_dGR_xi1[n][i] *= overeta0;
388 m_dGR_xi1[n][i] += dLp0[i];
389 m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
390 }
391 }
392 break;
393 }
394 case 2:
395 {
396 LibUtilities::BasisSharedPtr BasisFR_Left0;
397 LibUtilities::BasisSharedPtr BasisFR_Right0;
398 LibUtilities::BasisSharedPtr BasisFR_Left1;
399 LibUtilities::BasisSharedPtr BasisFR_Right1;
400
401 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
402 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
403 m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
404 m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
405
406 for (n = 0; n < nElements; ++n)
407 {
408 base = pFields[0]->GetExp(n)->GetBase();
409 nquad0 = base[0]->GetNumPoints();
410 nquad1 = base[1]->GetNumPoints();
411 nmodes0 = base[0]->GetNumModes();
412 nmodes1 = base[1]->GetNumModes();
413
414 Array<OneD, const NekDouble> z0;
415 Array<OneD, const NekDouble> w0;
416 Array<OneD, const NekDouble> z1;
417 Array<OneD, const NekDouble> w1;
418
419 base[0]->GetZW(z0, w0);
420 base[1]->GetZW(z1, w1);
421
422 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
423 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
424 m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
425 m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
426
427 // Auxiliary vectors to build up the auxiliary
428 // derivatives of the Legendre polynomials
429 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
430 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
431 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
432 Array<OneD, NekDouble> dLp1(nquad1, 0.0);
433 Array<OneD, NekDouble> dLpp1(nquad1, 0.0);
434 Array<OneD, NekDouble> dLpm1(nquad1, 0.0);
435
436 // Degree of the correction functions
437 int p0 = nmodes0 - 1;
438 int p1 = nmodes1 - 1;
439
440 // Function sign to compute left correction function
441 NekDouble sign0 = pow(-1.0, p0);
442 NekDouble sign1 = pow(-1.0, p1);
443
444 // Factorial factor to build the scheme
445 NekDouble ap0 =
446 std::tgamma(2 * p0 + 1) /
447 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
448
449 NekDouble ap1 =
450 std::tgamma(2 * p1 + 1) /
451 (pow(2.0, p1) * std::tgamma(p1 + 1) * std::tgamma(p1 + 1));
452
453 // Scalar parameter which recovers the FR schemes
454 if (m_advType == "FRDG")
455 {
456 c0 = 0.0;
457 c1 = 0.0;
458 }
459 else if (m_advType == "FRSD")
460 {
461 c0 = 2.0 * p0 /
462 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
463 (ap0 * std::tgamma(p0 + 1)) *
464 (ap0 * std::tgamma(p0 + 1)));
465
466 c1 = 2.0 * p1 /
467 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
468 (ap1 * std::tgamma(p1 + 1)) *
469 (ap1 * std::tgamma(p1 + 1)));
470 }
471 else if (m_advType == "FRHU")
472 {
473 c0 = 2.0 * (p0 + 1.0) /
474 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
475 (ap0 * std::tgamma(p0 + 1)));
476
477 c1 = 2.0 * (p1 + 1.0) /
478 ((2.0 * p1 + 1.0) * p1 * (ap1 * std::tgamma(p1 + 1)) *
479 (ap1 * std::tgamma(p1 + 1)));
480 }
481 else if (m_advType == "FRcmin")
482 {
483 c0 =
484 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
485 (ap0 * std::tgamma(p0 + 1)));
486
487 c1 =
488 -2.0 / ((2.0 * p1 + 1.0) * (ap1 * std::tgamma(p1 + 1)) *
489 (ap1 * std::tgamma(p1 + 1)));
490 }
491 else if (m_advType == "FRcinf")
492 {
493 c0 = 10000000000000000.0;
494 c1 = 10000000000000000.0;
495 }
496
497 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
498 (ap0 * std::tgamma(p0 + 1)) *
499 (ap0 * std::tgamma(p0 + 1));
500
501 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
502 (ap1 * std::tgamma(p1 + 1)) *
503 (ap1 * std::tgamma(p1 + 1));
504
505 NekDouble overeta0 = 1.0 / (1.0 + etap0);
506 NekDouble overeta1 = 1.0 / (1.0 + etap1);
507
508 // Derivative of the Legendre polynomials
509 // dLp = derivative of the Legendre polynomial of order p
510 // dLpp = derivative of the Legendre polynomial of order p+1
511 // dLpm = derivative of the Legendre polynomial of order p-1
512 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
513 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
514 0.0);
515 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
516 0.0);
517
518 Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
519 Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1 + 1, 0.0,
520 0.0);
521 Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1 - 1, 0.0,
522 0.0);
523
524 // Building the DG_c_Left
525 for (i = 0; i < nquad0; ++i)
526 {
527 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
528 m_dGL_xi1[n][i] += dLpp0[i];
529 m_dGL_xi1[n][i] *= overeta0;
530 m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
531 m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
532 }
533
534 // Building the DG_c_Left
535 for (i = 0; i < nquad1; ++i)
536 {
537 m_dGL_xi2[n][i] = etap1 * dLpm1[i];
538 m_dGL_xi2[n][i] += dLpp1[i];
539 m_dGL_xi2[n][i] *= overeta1;
540 m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
541 m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
542 }
543
544 // Building the DG_c_Right
545 for (i = 0; i < nquad0; ++i)
546 {
547 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
548 m_dGR_xi1[n][i] += dLpp0[i];
549 m_dGR_xi1[n][i] *= overeta0;
550 m_dGR_xi1[n][i] += dLp0[i];
551 m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
552 }
553
554 // Building the DG_c_Right
555 for (i = 0; i < nquad1; ++i)
556 {
557 m_dGR_xi2[n][i] = etap1 * dLpm1[i];
558 m_dGR_xi2[n][i] += dLpp1[i];
559 m_dGR_xi2[n][i] *= overeta1;
560 m_dGR_xi2[n][i] += dLp1[i];
561 m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
562 }
563 }
564 break;
565 }
566 case 3:
567 {
568 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
569 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
570 m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
571 m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
572 m_dGL_xi3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
573 m_dGR_xi3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
574
575 for (n = 0; n < nElements; ++n)
576 {
577 base = pFields[0]->GetExp(n)->GetBase();
578 nquad0 = base[0]->GetNumPoints();
579 nquad1 = base[1]->GetNumPoints();
580 nquad2 = base[2]->GetNumPoints();
581 nmodes0 = base[0]->GetNumModes();
582 nmodes1 = base[1]->GetNumModes();
583 nmodes2 = base[2]->GetNumModes();
584
585 Array<OneD, const NekDouble> z0;
586 Array<OneD, const NekDouble> w0;
587 Array<OneD, const NekDouble> z1;
588 Array<OneD, const NekDouble> w1;
589 Array<OneD, const NekDouble> z2;
590 Array<OneD, const NekDouble> w2;
591
592 base[0]->GetZW(z0, w0);
593 base[1]->GetZW(z1, w1);
594 base[1]->GetZW(z2, w2);
595
596 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
597 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
598 m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
599 m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
600 m_dGL_xi3[n] = Array<OneD, NekDouble>(nquad2);
601 m_dGR_xi3[n] = Array<OneD, NekDouble>(nquad2);
602
603 // Auxiliary vectors to build up the auxiliary
604 // derivatives of the Legendre polynomials
605 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
606 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
607 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
608 Array<OneD, NekDouble> dLp1(nquad1, 0.0);
609 Array<OneD, NekDouble> dLpp1(nquad1, 0.0);
610 Array<OneD, NekDouble> dLpm1(nquad1, 0.0);
611 Array<OneD, NekDouble> dLp2(nquad2, 0.0);
612 Array<OneD, NekDouble> dLpp2(nquad2, 0.0);
613 Array<OneD, NekDouble> dLpm2(nquad2, 0.0);
614
615 // Degree of the correction functions
616 int p0 = nmodes0 - 1;
617 int p1 = nmodes1 - 1;
618 int p2 = nmodes2 - 1;
619
620 // Function sign to compute left correction function
621 NekDouble sign0 = pow(-1.0, p0);
622 NekDouble sign1 = pow(-1.0, p1);
623
624 // Factorial factor to build the scheme
625 NekDouble ap0 =
626 std::tgamma(2 * p0 + 1) /
627 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
628
629 // Factorial factor to build the scheme
630 NekDouble ap1 =
631 std::tgamma(2 * p1 + 1) /
632 (pow(2.0, p1) * std::tgamma(p1 + 1) * std::tgamma(p1 + 1));
633
634 // Factorial factor to build the scheme
635 NekDouble ap2 =
636 std::tgamma(2 * p2 + 1) /
637 (pow(2.0, p2) * std::tgamma(p2 + 1) * std::tgamma(p2 + 1));
638
639 // Scalar parameter which recovers the FR schemes
640 if (m_advType == "FRDG")
641 {
642 c0 = 0.0;
643 c1 = 0.0;
644 c2 = 0.0;
645 }
646 else if (m_advType == "FRSD")
647 {
648 c0 = 2.0 * p0 /
649 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
650 (ap0 * std::tgamma(p0 + 1)) *
651 (ap0 * std::tgamma(p0 + 1)));
652
653 c1 = 2.0 * p1 /
654 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
655 (ap1 * std::tgamma(p1 + 1)) *
656 (ap1 * std::tgamma(p1 + 1)));
657
658 c2 = 2.0 * p2 /
659 ((2.0 * p2 + 1.0) * (p2 + 1.0) *
660 (ap2 * std::tgamma(p2 + 1)) *
661 (ap2 * std::tgamma(p2 + 1)));
662 }
663 else if (m_advType == "FRHU")
664 {
665 c0 = 2.0 * (p0 + 1.0) /
666 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
667 (ap0 * std::tgamma(p0 + 1)));
668
669 c1 = 2.0 * (p1 + 1.0) /
670 ((2.0 * p1 + 1.0) * p1 * (ap1 * std::tgamma(p1 + 1)) *
671 (ap1 * std::tgamma(p1 + 1)));
672
673 c2 = 2.0 * (p2 + 1.0) /
674 ((2.0 * p2 + 1.0) * p2 * (ap2 * std::tgamma(p2 + 1)) *
675 (ap2 * std::tgamma(p2 + 1)));
676 }
677 else if (m_advType == "FRcmin")
678 {
679 c0 =
680 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
681 (ap0 * std::tgamma(p0 + 1)));
682
683 c1 =
684 -2.0 / ((2.0 * p1 + 1.0) * (ap1 * std::tgamma(p1 + 1)) *
685 (ap1 * std::tgamma(p1 + 1)));
686
687 c2 =
688 -2.0 / ((2.0 * p2 + 1.0) * (ap2 * std::tgamma(p2 + 1)) *
689 (ap2 * std::tgamma(p2 + 1)));
690 }
691 else if (m_advType == "FRcinf")
692 {
693 c0 = 10000000000000000.0;
694 c1 = 10000000000000000.0;
695 c2 = 10000000000000000.0;
696 }
697
698 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
699 (ap0 * std::tgamma(p0 + 1)) *
700 (ap0 * std::tgamma(p0 + 1));
701
702 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
703 (ap1 * std::tgamma(p1 + 1)) *
704 (ap1 * std::tgamma(p1 + 1));
705
706 NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0) *
707 (ap2 * std::tgamma(p2 + 1)) *
708 (ap2 * std::tgamma(p2 + 1));
709
710 NekDouble overeta0 = 1.0 / (1.0 + etap0);
711 NekDouble overeta1 = 1.0 / (1.0 + etap1);
712 NekDouble overeta2 = 1.0 / (1.0 + etap2);
713
714 // Derivative of the Legendre polynomials
715 // dLp = derivative of the Legendre polynomial of order p
716 // dLpp = derivative of the Legendre polynomial of order p+1
717 // dLpm = derivative of the Legendre polynomial of order p-1
718 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
719 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
720 0.0);
721 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
722 0.0);
723
724 Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
725 Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1 + 1, 0.0,
726 0.0);
727 Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1 - 1, 0.0,
728 0.0);
729
730 Polylib::jacobd(nquad2, z2.data(), &(dLp2[0]), p2, 0.0, 0.0);
731 Polylib::jacobd(nquad2, z2.data(), &(dLpp2[0]), p2 + 1, 0.0,
732 0.0);
733 Polylib::jacobd(nquad2, z2.data(), &(dLpm2[0]), p2 - 1, 0.0,
734 0.0);
735
736 // Building the DG_c_Left
737 for (i = 0; i < nquad0; ++i)
738 {
739 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
740 m_dGL_xi1[n][i] += dLpp0[i];
741 m_dGL_xi1[n][i] *= overeta0;
742 m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
743 m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
744 }
745
746 // Building the DG_c_Left
747 for (i = 0; i < nquad1; ++i)
748 {
749 m_dGL_xi2[n][i] = etap1 * dLpm1[i];
750 m_dGL_xi2[n][i] += dLpp1[i];
751 m_dGL_xi2[n][i] *= overeta1;
752 m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
753 m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
754 }
755
756 // Building the DG_c_Left
757 for (i = 0; i < nquad2; ++i)
758 {
759 m_dGL_xi3[n][i] = etap2 * dLpm2[i];
760 m_dGL_xi3[n][i] += dLpp2[i];
761 m_dGL_xi3[n][i] *= overeta2;
762 m_dGL_xi3[n][i] = dLp2[i] - m_dGL_xi3[n][i];
763 m_dGL_xi3[n][i] = 0.5 * sign1 * m_dGL_xi3[n][i];
764 }
765
766 // Building the DG_c_Right
767 for (i = 0; i < nquad0; ++i)
768 {
769 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
770 m_dGR_xi1[n][i] += dLpp0[i];
771 m_dGR_xi1[n][i] *= overeta0;
772 m_dGR_xi1[n][i] += dLp0[i];
773 m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
774 }
775
776 // Building the DG_c_Right
777 for (i = 0; i < nquad1; ++i)
778 {
779 m_dGR_xi2[n][i] = etap1 * dLpm1[i];
780 m_dGR_xi2[n][i] += dLpp1[i];
781 m_dGR_xi2[n][i] *= overeta1;
782 m_dGR_xi2[n][i] += dLp1[i];
783 m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
784 }
785
786 // Building the DG_c_Right
787 for (i = 0; i < nquad2; ++i)
788 {
789 m_dGR_xi3[n][i] = etap2 * dLpm2[i];
790 m_dGR_xi3[n][i] += dLpp2[i];
791 m_dGR_xi3[n][i] *= overeta2;
792 m_dGR_xi3[n][i] += dLp2[i];
793 m_dGR_xi3[n][i] = 0.5 * m_dGR_xi3[n][i];
794 }
795 }
796 break;
797 }
798 default:
799 {
800 ASSERTL0(false, "Expansion dimension not recognised");
801 break;
802 }
803 }
804}
805
806/**
807 * @brief Compute the advection term at each time-step using the Flux
808 * Reconstruction approach (FR).
809 *
810 * @param nConvectiveFields Number of fields.
811 * @param fields Pointer to fields.
812 * @param advVel Advection velocities.
813 * @param inarray Solution at the previous time-step.
814 * @param outarray Advection term to be passed at the
815 * time integration class.
816 *
817 */
818void AdvectionFR::v_Advect(
819 const int nConvectiveFields,
820 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
821 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &advVel,
822 const Array<OneD, Array<OneD, NekDouble>> &inarray,
823 Array<OneD, Array<OneD, NekDouble>> &outarray,
824 [[maybe_unused]] const NekDouble &time,
825 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pFwd,
826 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pBwd)
827{
828 int i, j, n;
829 int phys_offset;
830
831 Array<OneD, NekDouble> auxArray1, auxArray2;
832
833 Array<OneD, LibUtilities::BasisSharedPtr> Basis;
834 Basis = fields[0]->GetExp(0)->GetBase();
835
836 int nElements = fields[0]->GetExpSize();
837 int nDimensions = fields[0]->GetCoordim(0);
838 int nSolutionPts = fields[0]->GetTotPoints();
839 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
840 int nCoeffs = fields[0]->GetNcoeffs();
841
842 // Storage for flux vector.
843 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> fluxvector(
844 nConvectiveFields);
845 // Outarray for Galerkin projection in case of primitive dealising
846 Array<OneD, Array<OneD, NekDouble>> outarrayCoeff(nConvectiveFields);
847 // Store forwards/backwards space along trace space
848 Array<OneD, Array<OneD, NekDouble>> Fwd(nConvectiveFields);
849 Array<OneD, Array<OneD, NekDouble>> Bwd(nConvectiveFields);
850 Array<OneD, Array<OneD, NekDouble>> numflux(nConvectiveFields);
851
852 // Set up storage for flux vector.
853 for (i = 0; i < nConvectiveFields; ++i)
854 {
855 fluxvector[i] = Array<OneD, Array<OneD, NekDouble>>(m_spaceDim);
856
857 for (j = 0; j < m_spaceDim; ++j)
858 {
859 fluxvector[i][j] = Array<OneD, NekDouble>(nSolutionPts);
860 }
861 }
862
863 for (i = 0; i < nConvectiveFields; ++i)
864 {
865 outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
866 Fwd[i] = Array<OneD, NekDouble>(nTracePts);
867 Bwd[i] = Array<OneD, NekDouble>(nTracePts);
868 numflux[i] = Array<OneD, NekDouble>(nTracePts);
869 fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
870 }
871
872 // Computing the interface flux at each trace point
873 m_riemann->Solve(m_spaceDim, Fwd, Bwd, numflux);
874
875 // Divergence of the flux (computing the RHS)
876 switch (nDimensions)
877 {
878 // 1D problems
879 case 1:
880 {
881 Array<OneD, NekDouble> DfluxvectorX1(nSolutionPts);
882 Array<OneD, NekDouble> divFC(nSolutionPts);
883
884 // Get the advection flux vector (solver specific)
885 m_fluxVector(inarray, fluxvector);
886
887 // Get the discontinuous flux divergence
888 for (i = 0; i < nConvectiveFields; ++i)
889 {
890 for (n = 0; n < nElements; n++)
891 {
892 phys_offset = fields[0]->GetPhys_Offset(n);
893
894 fields[i]->GetExp(n)->PhysDeriv(
895 0, fluxvector[i][0] + phys_offset,
896 auxArray2 = DfluxvectorX1 + phys_offset);
897 }
898
899 // Get the correction flux divergence
900 DivCFlux_1D(nConvectiveFields, fields, fluxvector[i][0],
901 numflux[i], divFC);
902
903 // Back to the physical space using global operations
904 Vmath::Vdiv(nSolutionPts, &divFC[0], 1, &m_jac[0], 1,
905 &outarray[i][0], 1);
906
907 // Adding the total divergence to outarray (RHS)
908 Vmath::Vadd(nSolutionPts, &outarray[i][0], 1, &DfluxvectorX1[0],
909 1, &outarray[i][0], 1);
910
911 // Primitive Dealiasing 1D
912 if (!(Basis[0]->Collocation()))
913 {
914 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
915 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
916 }
917 }
918 break;
919 }
920 // 2D problems
921 case 2:
922 {
923 Array<OneD, NekDouble> DfluxvectorX1(nSolutionPts);
924 Array<OneD, NekDouble> DfluxvectorX2(nSolutionPts);
925 Array<OneD, NekDouble> divFD(nSolutionPts);
926 Array<OneD, NekDouble> divFC(nSolutionPts);
927
928 // Get the advection flux vector (solver specific)
929 m_fluxVector(inarray, fluxvector);
930
931 for (i = 0; i < nConvectiveFields; ++i)
932 {
933 // Temporary vectors
934 Array<OneD, NekDouble> f_hat(nSolutionPts);
935 Array<OneD, NekDouble> g_hat(nSolutionPts);
936
937 Vmath::Vvtvvtp(nSolutionPts, &m_gmat[0][0], 1,
938 &fluxvector[i][0][0], 1, &m_gmat[2][0], 1,
939 &fluxvector[i][1][0], 1, &f_hat[0], 1);
940
941 Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &f_hat[0], 1, &f_hat[0],
942 1);
943
944 Vmath::Vvtvvtp(nSolutionPts, &m_gmat[1][0], 1,
945 &fluxvector[i][0][0], 1, &m_gmat[3][0], 1,
946 &fluxvector[i][1][0], 1, &g_hat[0], 1);
947
948 Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &g_hat[0], 1, &g_hat[0],
949 1);
950
951 // Get the discontinuous flux derivatives
952 for (n = 0; n < nElements; n++)
953 {
954 phys_offset = fields[0]->GetPhys_Offset(n);
955 fields[0]->GetExp(n)->StdPhysDeriv(
956 0, auxArray1 = f_hat + phys_offset,
957 auxArray2 = DfluxvectorX1 + phys_offset);
958 fields[0]->GetExp(n)->StdPhysDeriv(
959 1, auxArray1 = g_hat + phys_offset,
960 auxArray2 = DfluxvectorX2 + phys_offset);
961 }
962
963 // Divergence of the discontinuous flux
964 Vmath::Vadd(nSolutionPts, DfluxvectorX1, 1, DfluxvectorX2, 1,
965 divFD, 1);
966
967 // Divergence of the correction flux
968 if (Basis[0]->GetPointsType() ==
969 LibUtilities::eGaussGaussLegendre &&
970 Basis[1]->GetPointsType() ==
971 LibUtilities::eGaussGaussLegendre)
972 {
973 DivCFlux_2D_Gauss(nConvectiveFields, fields, f_hat, g_hat,
974 numflux[i], divFC);
975 }
976 else
977 {
978 DivCFlux_2D(nConvectiveFields, fields, fluxvector[i][0],
979 fluxvector[i][1], numflux[i], divFC);
980 }
981
982 // Divergence of the final flux
983 Vmath::Vadd(nSolutionPts, divFD, 1, divFC, 1, outarray[i], 1);
984
985 // Back to the physical space using a global operation
986 Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1, &m_jac[0], 1,
987 &outarray[i][0], 1);
988
989 // Primitive Dealiasing 2D
990 if (!(Basis[0]->Collocation()))
991 {
992 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
993 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
994 }
995 } // close nConvectiveFields loop
996 break;
997 }
998 // 3D problems
999 case 3:
1000 {
1001 ASSERTL0(false, "3D FRDG case not implemented yet");
1002 break;
1003 }
1004 }
1005}
1006
1007/**
1008 * @brief Compute the divergence of the corrective flux for 1D problems.
1009 *
1010 * @param nConvectiveFields Number of fields.
1011 * @param fields Pointer to fields.
1012 * @param fluxX1 X1-volumetric flux in physical space.
1013 * @param numericalFlux Interface flux in physical space.
1014 * @param divCFlux Divergence of the corrective flux.
1015 *
1016 */
1017void AdvectionFR::DivCFlux_1D(
1018 [[maybe_unused]] const int nConvectiveFields,
1019 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1020 const Array<OneD, const NekDouble> &fluxX1,
1021 const Array<OneD, const NekDouble> &numericalFlux,
1022 Array<OneD, NekDouble> &divCFlux)
1023{
1024 int n;
1025 int nLocalSolutionPts, phys_offset, t_offset;
1026
1027 LibUtilities::BasisSharedPtr Basis;
1028 Basis = fields[0]->GetExp(0)->GetBasis(0);
1029
1030 int nElements = fields[0]->GetExpSize();
1031 int nSolutionPts = fields[0]->GetTotPoints();
1032
1033 vector<bool> leftAdjacentTraces =
1034 std::static_pointer_cast<MultiRegions::DisContField>(fields[0])
1035 ->GetLeftAdjacentTraces();
1036
1037 // Arrays to store the derivatives of the correction flux
1038 Array<OneD, NekDouble> DCL(nSolutionPts / nElements, 0.0);
1039 Array<OneD, NekDouble> DCR(nSolutionPts / nElements, 0.0);
1040
1041 // Arrays to store the intercell numerical flux jumps
1042 Array<OneD, NekDouble> JumpL(nElements);
1043 Array<OneD, NekDouble> JumpR(nElements);
1044
1045 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1046 fields[0]->GetTraceMap()->GetElmtToTrace();
1047
1048 for (n = 0; n < nElements; ++n)
1049 {
1050 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1051 phys_offset = fields[0]->GetPhys_Offset(n);
1052
1053 Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1054 Array<OneD, NekDouble> tmpFluxVertex(1, 0.0);
1055 Vmath::Vcopy(nLocalSolutionPts, &fluxX1[phys_offset], 1, &tmparrayX1[0],
1056 1);
1057
1058 fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0], tmparrayX1,
1059 tmpFluxVertex);
1060
1061 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1062 elmtToTrace[n][0]->GetElmtId());
1063
1064 if (leftAdjacentTraces[2 * n])
1065 {
1066 JumpL[n] = -numericalFlux[t_offset] - tmpFluxVertex[0];
1067 }
1068 else
1069 {
1070 JumpL[n] = numericalFlux[t_offset] - tmpFluxVertex[0];
1071 }
1072
1073 fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1], tmparrayX1,
1074 tmpFluxVertex);
1075
1076 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1077 elmtToTrace[n][1]->GetElmtId());
1078
1079 if (leftAdjacentTraces[2 * n + 1])
1080 {
1081 JumpR[n] = numericalFlux[t_offset] - tmpFluxVertex[0];
1082 }
1083 else
1084 {
1085 JumpR[n] = -numericalFlux[t_offset] - tmpFluxVertex[0];
1086 }
1087 }
1088
1089 for (n = 0; n < nElements; ++n)
1090 {
1091 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1092 phys_offset = fields[0]->GetPhys_Offset(n);
1093
1094 // Left jump multiplied by left derivative of C function
1095 Vmath::Smul(nLocalSolutionPts, JumpL[n], &m_dGL_xi1[n][0], 1, &DCL[0],
1096 1);
1097
1098 // Right jump multiplied by right derivative of C function
1099 Vmath::Smul(nLocalSolutionPts, JumpR[n], &m_dGR_xi1[n][0], 1, &DCR[0],
1100 1);
1101
1102 // Assembling divergence of the correction flux
1103 Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1104 &divCFlux[phys_offset], 1);
1105 }
1106}
1107
1108/**
1109 * @brief Compute the divergence of the corrective flux for 2D problems.
1110 *
1111 * @param nConvectiveFields Number of fields.
1112 * @param fields Pointer to fields.
1113 * @param fluxX1 X1-volumetric flux in physical space.
1114 * @param fluxX2 X2-volumetric flux in physical space.
1115 * @param numericalFlux Interface flux in physical space.
1116 * @param divCFlux Divergence of the corrective flux.
1117 *
1118 * \todo: Switch on shapes eventually here.
1119 */
1120void AdvectionFR::DivCFlux_2D(
1121 [[maybe_unused]] const int nConvectiveFields,
1122 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1123 const Array<OneD, const NekDouble> &fluxX1,
1124 const Array<OneD, const NekDouble> &fluxX2,
1125 const Array<OneD, const NekDouble> &numericalFlux,
1126 Array<OneD, NekDouble> &divCFlux)
1127{
1128 int n, e, i, j, cnt;
1129
1130 int nElements = fields[0]->GetExpSize();
1131
1132 int nLocalSolutionPts, nEdgePts;
1133 int trace_offset, phys_offset;
1134 int nquad0, nquad1;
1135
1136 Array<OneD, NekDouble> auxArray1;
1137 Array<OneD, LibUtilities::BasisSharedPtr> base;
1138
1139 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1140 fields[0]->GetTraceMap()->GetElmtToTrace();
1141
1142 // Loop on the elements
1143 for (n = 0; n < nElements; ++n)
1144 {
1145 // Offset of the element on the global vector
1146 phys_offset = fields[0]->GetPhys_Offset(n);
1147 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1148
1149 base = fields[0]->GetExp(n)->GetBase();
1150 nquad0 = base[0]->GetNumPoints();
1151 nquad1 = base[1]->GetNumPoints();
1152
1153 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1154 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1155 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1156 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1157
1158 // Loop on the edges
1159 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1160 {
1161 // Number of edge points of edge e
1162 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1163
1164 Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1165 Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1166 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1167 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1168 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1169
1170 // Offset of the trace space correspondent to edge e
1171 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1172 elmtToTrace[n][e]->GetElmtId());
1173
1174 // Get the normals of edge e
1175 const Array<OneD, const Array<OneD, NekDouble>> &normals =
1176 fields[0]->GetExp(n)->GetTraceNormal(e);
1177
1178 // Extract the edge values of flux-x on edge e and order
1179 // them accordingly to the order of the trace space
1180 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1181 fluxX1 + phys_offset,
1182 auxArray1 = tmparrayX1);
1183
1184 // Extract the edge values of flux-y on edge e and order
1185 // them accordingly to the order of the trace space
1186 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1187 fluxX2 + phys_offset,
1188 auxArray1 = tmparrayX2);
1189
1190 // Multiply the edge components of the flux by the normal
1191 Vmath::Vvtvvtp(nEdgePts, &tmparrayX1[0], 1,
1192 &m_traceNormals[0][trace_offset], 1, &tmparrayX2[0],
1193 1, &m_traceNormals[1][trace_offset], 1, &fluxN[0],
1194 1);
1195
1196 // Subtract to the Riemann flux the discontinuous flux
1197 Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1, &fluxN[0], 1,
1198 &fluxJumps[0], 1);
1199
1200 // Check the ordering of the jump vectors
1201 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1202 StdRegions::eBackwards)
1203 {
1204 Vmath::Reverse(nEdgePts, &fluxJumps[0], 1, &fluxJumps[0], 1);
1205 }
1206
1207 for (i = 0; i < nEdgePts; ++i)
1208 {
1209 if (m_traceNormals[0][trace_offset + i] != normals[0][i] ||
1210 m_traceNormals[1][trace_offset + i] != normals[1][i])
1211 {
1212 fluxJumps[i] = -fluxJumps[i];
1213 }
1214 }
1215
1216 // Multiply jumps by derivatives of the correction functions
1217 switch (e)
1218 {
1219 case 0:
1220 for (i = 0; i < nquad0; ++i)
1221 {
1222 // Multiply fluxJumps by Q factors
1223 fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1224
1225 for (j = 0; j < nquad1; ++j)
1226 {
1227 cnt = i + j * nquad0;
1228 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1229 }
1230 }
1231 break;
1232 case 1:
1233 for (i = 0; i < nquad1; ++i)
1234 {
1235 // Multiply fluxJumps by Q factors
1236 fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1237
1238 for (j = 0; j < nquad0; ++j)
1239 {
1240 cnt = (nquad0)*i + j;
1241 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1242 }
1243 }
1244 break;
1245 case 2:
1246 for (i = 0; i < nquad0; ++i)
1247 {
1248 // Multiply fluxJumps by Q factors
1249 fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1250
1251 for (j = 0; j < nquad1; ++j)
1252 {
1253 cnt = j * nquad0 + i;
1254 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1255 }
1256 }
1257 break;
1258 case 3:
1259 for (i = 0; i < nquad1; ++i)
1260 {
1261 // Multiply fluxJumps by Q factors
1262 fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1263 for (j = 0; j < nquad0; ++j)
1264 {
1265 cnt = j + i * nquad0;
1266 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1267 }
1268 }
1269 break;
1270 default:
1271 ASSERTL0(false, "edge value (< 3) is out of range");
1272 break;
1273 }
1274 }
1275
1276 // Sum each edge contribution
1277 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
1278 &divCFlux[phys_offset], 1);
1279
1280 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1281 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1282
1283 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1284 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1285 }
1286}
1287
1288/**
1289 * @brief Compute the divergence of the corrective flux for 2D problems
1290 * where POINTSTYPE="GaussGaussLegendre"
1291 *
1292 * @param nConvectiveFields Number of fields.
1293 * @param fields Pointer to fields.
1294 * @param fluxX1 X1-volumetric flux in physical space.
1295 * @param fluxX2 X2-volumetric flux in physical space.
1296 * @param numericalFlux Interface flux in physical space.
1297 * @param divCFlux Divergence of the corrective flux.
1298 *
1299 * \todo: Switch on shapes eventually here.
1300 */
1301
1302void AdvectionFR::DivCFlux_2D_Gauss(
1303 [[maybe_unused]] const int nConvectiveFields,
1304 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1305 const Array<OneD, const NekDouble> &fluxX1,
1306 const Array<OneD, const NekDouble> &fluxX2,
1307 const Array<OneD, const NekDouble> &numericalFlux,
1308 Array<OneD, NekDouble> &divCFlux)
1309{
1310 int n, e, i, j, cnt;
1311
1312 int nElements = fields[0]->GetExpSize();
1313 int nLocalSolutionPts;
1314 int nEdgePts;
1315 int trace_offset;
1316 int phys_offset;
1317 int nquad0;
1318 int nquad1;
1319
1320 Array<OneD, NekDouble> auxArray1, auxArray2;
1321 Array<OneD, LibUtilities::BasisSharedPtr> base;
1322
1323 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1324 fields[0]->GetTraceMap()->GetElmtToTrace();
1325
1326 // Loop on the elements
1327 for (n = 0; n < nElements; ++n)
1328 {
1329 // Offset of the element on the global vector
1330 phys_offset = fields[0]->GetPhys_Offset(n);
1331 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1332
1333 base = fields[0]->GetExp(n)->GetBase();
1334 nquad0 = base[0]->GetNumPoints();
1335 nquad1 = base[1]->GetNumPoints();
1336
1337 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1338 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1339 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1340 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1341
1342 // Loop on the edges
1343 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1344 {
1345 // Number of edge points of edge e
1346 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1347
1348 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1349 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1350 Array<OneD, NekDouble> fluxN_R(nEdgePts, 0.0);
1351 Array<OneD, NekDouble> fluxN_D(nEdgePts, 0.0);
1352 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1353
1354 // Offset of the trace space correspondent to edge e
1355 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1356 elmtToTrace[n][e]->GetElmtId());
1357
1358 // Get the normals of edge e
1359 const Array<OneD, const Array<OneD, NekDouble>> &normals =
1360 fields[0]->GetExp(n)->GetTraceNormal(e);
1361
1362 // Extract the trasformed normal flux at each edge
1363 switch (e)
1364 {
1365 case 0:
1366 // Extract the edge values of transformed flux-y on
1367 // edge e and order them accordingly to the order of
1368 // the trace space
1369 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1370 fluxX2 + phys_offset,
1371 auxArray1 = fluxN_D);
1372
1373 Vmath::Neg(nEdgePts, fluxN_D, 1);
1374
1375 // Extract the physical Rieamann flux at each edge
1376 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
1377 &fluxN[0], 1);
1378
1379 // Check the ordering of vectors
1380 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1381 StdRegions::eBackwards)
1382 {
1383 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
1384 auxArray2 = fluxN, 1);
1385
1386 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
1387 auxArray2 = fluxN_D, 1);
1388 }
1389
1390 // Transform Riemann Fluxes in the standard element
1391 for (i = 0; i < nquad0; ++i)
1392 {
1393 // Multiply Riemann Flux by Q factors
1394 fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
1395 }
1396
1397 for (i = 0; i < nEdgePts; ++i)
1398 {
1399 if (m_traceNormals[0][trace_offset + i] !=
1400 normals[0][i] ||
1401 m_traceNormals[1][trace_offset + i] !=
1402 normals[1][i])
1403 {
1404 fluxN_R[i] = -fluxN_R[i];
1405 }
1406 }
1407
1408 // Subtract to the Riemann flux the discontinuous
1409 // flux
1410 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
1411 &fluxJumps[0], 1);
1412
1413 // Multiplicate the flux jumps for the correction
1414 // function
1415 for (i = 0; i < nquad0; ++i)
1416 {
1417 for (j = 0; j < nquad1; ++j)
1418 {
1419 cnt = i + j * nquad0;
1420 divCFluxE0[cnt] = -fluxJumps[i] * m_dGL_xi2[n][j];
1421 }
1422 }
1423 break;
1424 case 1:
1425 // Extract the edge values of transformed flux-x on
1426 // edge e and order them accordingly to the order of
1427 // the trace space
1428 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1429 fluxX1 + phys_offset,
1430 auxArray1 = fluxN_D);
1431
1432 // Extract the physical Rieamann flux at each edge
1433 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
1434 &fluxN[0], 1);
1435
1436 // Check the ordering of vectors
1437 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1438 StdRegions::eBackwards)
1439 {
1440 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
1441 auxArray2 = fluxN, 1);
1442
1443 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
1444 auxArray2 = fluxN_D, 1);
1445 }
1446
1447 // Transform Riemann Fluxes in the standard element
1448 for (i = 0; i < nquad1; ++i)
1449 {
1450 // Multiply Riemann Flux by Q factors
1451 fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
1452 }
1453
1454 for (i = 0; i < nEdgePts; ++i)
1455 {
1456 if (m_traceNormals[0][trace_offset + i] !=
1457 normals[0][i] ||
1458 m_traceNormals[1][trace_offset + i] !=
1459 normals[1][i])
1460 {
1461 fluxN_R[i] = -fluxN_R[i];
1462 }
1463 }
1464
1465 // Subtract to the Riemann flux the discontinuous
1466 // flux
1467 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
1468 &fluxJumps[0], 1);
1469
1470 // Multiplicate the flux jumps for the correction
1471 // function
1472 for (i = 0; i < nquad1; ++i)
1473 {
1474 for (j = 0; j < nquad0; ++j)
1475 {
1476 cnt = (nquad0)*i + j;
1477 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1478 }
1479 }
1480 break;
1481 case 2:
1482
1483 // Extract the edge values of transformed flux-y on
1484 // edge e and order them accordingly to the order of
1485 // the trace space
1486
1487 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1488 fluxX2 + phys_offset,
1489 auxArray1 = fluxN_D);
1490
1491 // Extract the physical Rieamann flux at each edge
1492 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
1493 &fluxN[0], 1);
1494
1495 // Check the ordering of vectors
1496 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1497 StdRegions::eBackwards)
1498 {
1499 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
1500 auxArray2 = fluxN, 1);
1501
1502 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
1503 auxArray2 = fluxN_D, 1);
1504 }
1505
1506 // Transform Riemann Fluxes in the standard element
1507 for (i = 0; i < nquad0; ++i)
1508 {
1509 // Multiply Riemann Flux by Q factors
1510 fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
1511 }
1512
1513 for (i = 0; i < nEdgePts; ++i)
1514 {
1515 if (m_traceNormals[0][trace_offset + i] !=
1516 normals[0][i] ||
1517 m_traceNormals[1][trace_offset + i] !=
1518 normals[1][i])
1519 {
1520 fluxN_R[i] = -fluxN_R[i];
1521 }
1522 }
1523
1524 // Subtract to the Riemann flux the discontinuous
1525 // flux
1526
1527 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
1528 &fluxJumps[0], 1);
1529
1530 // Multiplicate the flux jumps for the correction
1531 // function
1532 for (i = 0; i < nquad0; ++i)
1533 {
1534 for (j = 0; j < nquad1; ++j)
1535 {
1536 cnt = j * nquad0 + i;
1537 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1538 }
1539 }
1540 break;
1541 case 3:
1542 // Extract the edge values of transformed flux-x on
1543 // edge e and order them accordingly to the order of
1544 // the trace space
1545
1546 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1547 fluxX1 + phys_offset,
1548 auxArray1 = fluxN_D);
1549 Vmath::Neg(nEdgePts, fluxN_D, 1);
1550
1551 // Extract the physical Rieamann flux at each edge
1552 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
1553 &fluxN[0], 1);
1554
1555 // Check the ordering of vectors
1556 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1557 StdRegions::eBackwards)
1558 {
1559 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
1560 auxArray2 = fluxN, 1);
1561
1562 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
1563 auxArray2 = fluxN_D, 1);
1564 }
1565
1566 // Transform Riemann Fluxes in the standard element
1567 for (i = 0; i < nquad1; ++i)
1568 {
1569 // Multiply Riemann Flux by Q factors
1570 fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
1571 }
1572
1573 for (i = 0; i < nEdgePts; ++i)
1574 {
1575 if (m_traceNormals[0][trace_offset + i] !=
1576 normals[0][i] ||
1577 m_traceNormals[1][trace_offset + i] !=
1578 normals[1][i])
1579 {
1580 fluxN_R[i] = -fluxN_R[i];
1581 }
1582 }
1583
1584 // Subtract to the Riemann flux the discontinuous
1585 // flux
1586
1587 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
1588 &fluxJumps[0], 1);
1589
1590 // Multiplicate the flux jumps for the correction
1591 // function
1592 for (i = 0; i < nquad1; ++i)
1593 {
1594 for (j = 0; j < nquad0; ++j)
1595 {
1596 cnt = j + i * nquad0;
1597 divCFluxE3[cnt] = -fluxJumps[i] * m_dGL_xi1[n][j];
1598 }
1599 }
1600 break;
1601 default:
1602 ASSERTL0(false, "edge value (< 3) is out of range");
1603 break;
1604 }
1605 }
1606
1607 // Sum each edge contribution
1608 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
1609 &divCFlux[phys_offset], 1);
1610
1611 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1612 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1613
1614 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1615 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1616 }
1617}
1618
1619/**
1620 * @brief Compute the divergence of the corrective flux for 3D problems.
1621 *
1622 * @param nConvectiveFields Number of fields.
1623 * @param fields Pointer to fields.
1624 * @param fluxX1 X1-volumetric flux in physical space.
1625 * @param fluxX2 X2-volumetric flux in physical space.
1626 * @param fluxX3 X3-volumetric flux in physical space.
1627 * @param numericalFlux Interface flux in physical space.
1628 * @param divCFlux Divergence of the corrective flux.
1629 *
1630 * \todo: To be implemented. Switch on shapes eventually here.
1631 */
1632void AdvectionFR::DivCFlux_3D(
1633 [[maybe_unused]] const int nConvectiveFields,
1634 [[maybe_unused]] const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1635 [[maybe_unused]] const Array<OneD, const NekDouble> &fluxX1,
1636 [[maybe_unused]] const Array<OneD, const NekDouble> &fluxX2,
1637 [[maybe_unused]] const Array<OneD, const NekDouble> &fluxX3,
1638 [[maybe_unused]] const Array<OneD, const NekDouble> &numericalFlux,
1639 [[maybe_unused]] Array<OneD, NekDouble> &divCFlux)
1640{
1641}
1642
1643} // namespace Nektar::SolverUtils
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
WARNINGL0
#define WARNINGL0(condition, msg)
Definition: ErrorUtil.hpp:215
Nektar::LibUtilities::Basis
Represents a basis of a given type.
Definition: Basis.h:198
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: BasicUtils/NekFactory.hpp:197
Nektar::LocalRegions::Expansion::GetMetricInfo
const SpatialDomains::GeomFactorsSharedPtr & GetMetricInfo() const
Definition: Expansion.cpp:246
Nektar::SolverUtils::AdvectionFR::create
static AdvectionSharedPtr create(std::string advType)
Definition: AdvectionFR.h:47
Nektar::SolverUtils::AdvectionFR::v_Advect
void v_Advect(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray) override
Compute the advection term at each time-step using the Flux Reconstruction approach (FR).
Definition: AdvectionFR.cpp:818
Nektar::SolverUtils::AdvectionFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: AdvectionFR.h:85
Nektar::SolverUtils::AdvectionFR::m_dGR_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition: AdvectionFR.h:90
Nektar::SolverUtils::AdvectionFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: AdvectionFR.h:86
Nektar::SolverUtils::AdvectionFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: AdvectionFR.h:88
Nektar::SolverUtils::AdvectionFR::SetupMetrics
void SetupMetrics(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform...
Definition: AdvectionFR.cpp:111
Nektar::SolverUtils::AdvectionFR::m_dGL_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition: AdvectionFR.h:89
Nektar::SolverUtils::AdvectionFR::m_jac
Array< OneD, NekDouble > m_jac
Definition: AdvectionFR.h:77
Nektar::SolverUtils::AdvectionFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: AdvectionFR.h:80
Nektar::SolverUtils::AdvectionFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: AdvectionFR.h:83
Nektar::SolverUtils::AdvectionFR::AdvectionFR
AdvectionFR(std::string advType)
AdvectionFR uses the Flux Reconstruction (FR) approach to compute the advection term....
Definition: AdvectionFR.cpp:70
Nektar::SolverUtils::AdvectionFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: AdvectionFR.h:81
Nektar::SolverUtils::AdvectionFR::v_InitObject
void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
Initiliase AdvectionFR objects and store them before starting the time-stepping.
Definition: AdvectionFR.cpp:84
Nektar::SolverUtils::AdvectionFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: AdvectionFR.h:87
Nektar::SolverUtils::AdvectionFR::DivCFlux_2D_Gauss
void DivCFlux_2D_Gauss(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre".
Definition: AdvectionFR.cpp:1302
Nektar::SolverUtils::AdvectionFR::SetupCFunctions
void SetupCFunctions(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–72...
Definition: AdvectionFR.cpp:275
Nektar::SolverUtils::AdvectionFR::DivCFlux_2D
void DivCFlux_2D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems.
Definition: AdvectionFR.cpp:1120
Nektar::SolverUtils::AdvectionFR::DivCFlux_3D
void DivCFlux_3D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &fluxX3, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 3D problems.
Definition: AdvectionFR.cpp:1632
Nektar::SolverUtils::AdvectionFR::type
static std::string type[]
Definition: AdvectionFR.h:52
Nektar::SolverUtils::AdvectionFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: AdvectionFR.h:75
Nektar::SolverUtils::AdvectionFR::DivCFlux_1D
void DivCFlux_1D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 1D problems.
Definition: AdvectionFR.cpp:1017
Nektar::SolverUtils::AdvectionFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: AdvectionFR.h:82
Nektar::SolverUtils::AdvectionFR::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: AdvectionFR.h:78
Nektar::SolverUtils::Advection::m_spaceDim
int m_spaceDim
Storage for space dimension. Used for homogeneous extension.
Definition: Advection.h:179
Nektar::SolverUtils::Advection::m_fluxVector
AdvectionFluxVecCB m_fluxVector
Callback function to the flux vector (set when advection is in conservative form).
Definition: Advection.h:175
Nektar::SolverUtils::Advection::v_InitObject
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Definition: Advection.cpp:264
Nektar::SolverUtils::Advection::m_riemann
RiemannSolverSharedPtr m_riemann
Riemann solver for DG-type schemes.
Definition: Advection.h:177
Nektar::LibUtilities::BasisSharedPtr
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:69
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:119
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:231
Nektar::LibUtilities::eGaussGaussLegendre
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:46
Nektar::SolverUtils::GetAdvectionFactory
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43
Nektar::SpatialDomains::eDeformed
@ eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:54
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Polylib::jacobd
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition: Polylib.cpp:1378
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.hpp:126
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:844
Vmath::Vvtvvtp
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.hpp:439
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.hpp:220
std
STL namespace.
Generated on Sun Jun 15 2025 07:52:32 for Nektar++ by doxygen 1.9.4