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