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