Nektar++
Loading...
Searching...
No Matches
DiffusionLFR.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: DiffusionLFR.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: LFR diffusion class.
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <LibUtilities/Polylib/Polylib.h>
36#include <LocalRegions/Expansion1D.h>
37#include <LocalRegions/Expansion2D.h>
38#include <MultiRegions/DisContField.h>
39#include <SolverUtils/Diffusion/DiffusionLFR.h>
40#include <iomanip>
41#include <iostream>
42
43using namespace std;
44
45namespace Nektar::SolverUtils
46{
47std::string DiffusionLFR::type[] = {
48 GetDiffusionFactory().RegisterCreatorFunction("LFRDG", DiffusionLFR::create,
49 "LFRDG"),
50 GetDiffusionFactory().RegisterCreatorFunction("LFRSD", DiffusionLFR::create,
51 "LRFSD"),
52 GetDiffusionFactory().RegisterCreatorFunction("LFRHU", DiffusionLFR::create,
53 "LFRHU"),
54 GetDiffusionFactory().RegisterCreatorFunction(
55 "LFRcmin", DiffusionLFR::create, "LFRcmin"),
56 GetDiffusionFactory().RegisterCreatorFunction(
57 "LFRcinf", DiffusionLFR::create, "LFRcinf")};
58
59/**
60 * @brief DiffusionLFR uses the Flux Reconstruction (FR) approach to
61 * compute the diffusion term. The implementation is only for segments,
62 * quadrilaterals and hexahedra at the moment.
63 *
64 * \todo Extension to triangles, tetrahedra and other shapes.
65 * (Long term objective)
66 */
67DiffusionLFR::DiffusionLFR(std::string diffType) : m_diffType(diffType)
68{
69}
70
71/**
72 * @brief Initiliase DiffusionLFR objects and store them before starting
73 * the time-stepping.
74 *
75 * This routine calls the functions #SetupMetrics and
76 * #SetupCFunctions to initialise the objects needed by DiffusionLFR.
77 *
78 * @param pSession Pointer to session reader.
79 * @param pFields Pointer to fields.
80 */
81void DiffusionLFR::v_InitObject(
82 LibUtilities::SessionReaderSharedPtr pSession,
83 Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
84{
85 WARNINGL0(false, "LFR is deprecated, use LDG instead");
86
87 m_session = pSession;
88 SetupMetrics(pSession, pFields);
89 SetupCFunctions(pSession, pFields);
90
91 // Initialising arrays
92 int i, j;
93 int nConvectiveFields = pFields.size();
94 int nDim = pFields[0]->GetCoordim(0);
95 int nSolutionPts = pFields[0]->GetTotPoints();
96 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
97
98 m_IF1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
99 m_DU1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
100 m_DFC1 =
101 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
102 m_BD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
103 m_D1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
104 m_DD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
105 m_IF2 = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
106 m_DFC2 =
107 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
108 m_divFD = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
109 m_divFC = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
110
111 m_tmp1 =
112 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
113 m_tmp2 =
114 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
115
116 for (i = 0; i < nConvectiveFields; ++i)
117 {
118 m_IF1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
119 m_DU1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
120 m_DFC1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
121 m_BD1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
122 m_D1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
123 m_DD1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
124 m_IF2[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
125 m_DFC2[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
126 m_divFD[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
127 m_divFC[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
128
129 m_tmp1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
130 m_tmp2[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
131
132 for (j = 0; j < nDim; ++j)
133 {
134 m_IF1[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
135 m_DU1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
136 m_DFC1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
137 m_BD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
138 m_D1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
139 m_DD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
140 m_DFC2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
141
142 m_tmp1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
143 m_tmp2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
144 }
145 }
146}
147
148/**
149 * @brief Setup the metric terms to compute the contravariant
150 * fluxes. (i.e. this special metric terms transform the fluxes
151 * at the interfaces of each element from the physical space to
152 * the standard space).
153 *
154 * This routine calls the function #GetTraceQFactors to compute and
155 * store the metric factors following an anticlockwise conventions
156 * along the edges/faces of each element. Note: for 1D problem
157 * the transformation is not needed.
158 *
159 * @param pSession Pointer to session reader.
160 * @param pFields Pointer to fields.
161 *
162 * \todo Add the metric terms for 3D Hexahedra.
163 */
164void DiffusionLFR::SetupMetrics(
165 [[maybe_unused]] LibUtilities::SessionReaderSharedPtr pSession,
166 Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
167{
168 int i, n;
169 int nquad0, nquad1;
170 int phys_offset;
171 int nLocalSolutionPts;
172 int nElements = pFields[0]->GetExpSize();
173 int nDimensions = pFields[0]->GetCoordim(0);
174 int nSolutionPts = pFields[0]->GetTotPoints();
175 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
176
177 m_traceNormals = Array<OneD, Array<OneD, NekDouble>>(nDimensions);
178 for (i = 0; i < nDimensions; ++i)
179 {
180 m_traceNormals[i] = Array<OneD, NekDouble>(nTracePts);
181 }
182 pFields[0]->GetTrace()->GetNormals(m_traceNormals);
183
184 Array<TwoD, const NekDouble> gmat;
185 Array<OneD, const NekDouble> jac;
186
187 m_jac = Array<OneD, NekDouble>(nSolutionPts);
188
189 Array<OneD, NekDouble> auxArray1;
190 Array<OneD, LibUtilities::BasisSharedPtr> base;
191 LibUtilities::PointsKeyVector ptsKeys;
192
193 switch (nDimensions)
194 {
195 case 1:
196 {
197 for (n = 0; n < nElements; ++n)
198 {
199 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
200 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
201 phys_offset = pFields[0]->GetPhys_Offset(n);
202 jac = pFields[0]
203 ->GetExp(n)
204 ->as<LocalRegions::Expansion1D>()
205 ->GetGeom1D()
206 ->GetMetricInfo()
207 ->GetJac(ptsKeys);
208 for (i = 0; i < nLocalSolutionPts; ++i)
209 {
210 m_jac[i + phys_offset] = jac[0];
211 }
212 }
213 break;
214 }
215 case 2:
216 {
217 m_gmat = Array<OneD, Array<OneD, NekDouble>>(4);
218 m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
219 m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
220 m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
221 m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
222
223 m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble>>(nElements);
224 m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
225 m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
226 m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
227
228 for (n = 0; n < nElements; ++n)
229 {
230 base = pFields[0]->GetExp(n)->GetBase();
231 nquad0 = base[0]->GetNumPoints();
232 nquad1 = base[1]->GetNumPoints();
233
234 m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
235 m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
236 m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
237 m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
238
239 // Extract the Q factors at each edge point
240 pFields[0]->GetExp(n)->GetTraceQFactors(0, auxArray1 =
241 m_Q2D_e0[n]);
242 pFields[0]->GetExp(n)->GetTraceQFactors(1, auxArray1 =
243 m_Q2D_e1[n]);
244 pFields[0]->GetExp(n)->GetTraceQFactors(2, auxArray1 =
245 m_Q2D_e2[n]);
246 pFields[0]->GetExp(n)->GetTraceQFactors(3, auxArray1 =
247 m_Q2D_e3[n]);
248
249 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
250 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
251 phys_offset = pFields[0]->GetPhys_Offset(n);
252
253 jac = pFields[0]
254 ->GetExp(n)
255 ->as<LocalRegions::Expansion2D>()
256 ->GetGeom2D()
257 ->GetMetricInfo()
258 ->GetJac(ptsKeys);
259 gmat = pFields[0]
260 ->GetExp(n)
261 ->as<LocalRegions::Expansion2D>()
262 ->GetGeom2D()
263 ->GetMetricInfo()
264 ->GetDerivFactors(ptsKeys);
265
266 if (pFields[0]
267 ->GetExp(n)
268 ->as<LocalRegions::Expansion2D>()
269 ->GetGeom2D()
270 ->GetMetricInfo()
271 ->GetGtype() == SpatialDomains::eDeformed)
272 {
273 for (i = 0; i < nLocalSolutionPts; ++i)
274 {
275 m_jac[i + phys_offset] = jac[i];
276 m_gmat[0][i + phys_offset] = gmat[0][i];
277 m_gmat[1][i + phys_offset] = gmat[1][i];
278 m_gmat[2][i + phys_offset] = gmat[2][i];
279 m_gmat[3][i + phys_offset] = gmat[3][i];
280 }
281 }
282 else
283 {
284 for (i = 0; i < nLocalSolutionPts; ++i)
285 {
286 m_jac[i + phys_offset] = jac[0];
287 m_gmat[0][i + phys_offset] = gmat[0][0];
288 m_gmat[1][i + phys_offset] = gmat[1][0];
289 m_gmat[2][i + phys_offset] = gmat[2][0];
290 m_gmat[3][i + phys_offset] = gmat[3][0];
291 }
292 }
293 }
294 break;
295 }
296 case 3:
297 {
298 ASSERTL0(false, "3DFR Metric terms not implemented yet");
299 break;
300 }
301 default:
302 {
303 ASSERTL0(false, "Expansion dimension not recognised");
304 break;
305 }
306 }
307}
308
309/**
310 * @brief Setup the derivatives of the correction functions. For more
311 * details see J Sci Comput (2011) 47: 50–72
312 *
313 * This routine calls 3 different bases:
314 * #eDG_DG_Left - #eDG_DG_Left which recovers a nodal DG scheme,
315 * #eDG_SD_Left - #eDG_SD_Left which recovers the SD scheme,
316 * #eDG_HU_Left - #eDG_HU_Left which recovers the Huynh scheme.
317 *
318 * The values of the derivatives of the correction function are then
319 * stored into global variables and reused into the Functions
320 * #DivCFlux_1D, #DivCFlux_2D, #DivCFlux_3D to compute the
321 * the divergence of the correction flux for 1D, 2D or 3D problems.
322 *
323 * @param pSession Pointer to session reader.
324 * @param pFields Pointer to fields.
325 */
326void DiffusionLFR::SetupCFunctions(
327 [[maybe_unused]] LibUtilities::SessionReaderSharedPtr pSession,
328 Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
329{
330 int i, n;
331 NekDouble c0 = 0.0;
332 NekDouble c1 = 0.0;
333 NekDouble c2 = 0.0;
334 int nquad0, nquad1, nquad2;
335 int nmodes0, nmodes1, nmodes2;
336 Array<OneD, LibUtilities::BasisSharedPtr> base;
337
338 int nElements = pFields[0]->GetExpSize();
339 int nDim = pFields[0]->GetCoordim(0);
340
341 switch (nDim)
342 {
343 case 1:
344 {
345 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
346 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
347
348 for (n = 0; n < nElements; ++n)
349 {
350 base = pFields[0]->GetExp(n)->GetBase();
351 nquad0 = base[0]->GetNumPoints();
352 nmodes0 = base[0]->GetNumModes();
353 Array<OneD, const NekDouble> z0;
354 Array<OneD, const NekDouble> w0;
355
356 base[0]->GetZW(z0, w0);
357
358 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
359 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
360
361 // Auxiliary vectors to build up the auxiliary
362 // derivatives of the Legendre polynomials
363 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
364 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
365 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
366
367 // Degree of the correction functions
368 int p0 = nmodes0 - 1;
369
370 // Function sign to compute left correction function
371 NekDouble sign0 = pow(-1.0, p0);
372
373 // Factorial factor to build the scheme
374 NekDouble ap0 =
375 std::tgamma(2 * p0 + 1) /
376 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
377
378 // Scalar parameter which recovers the FR schemes
379 if (m_diffType == "LFRDG")
380 {
381 c0 = 0.0;
382 }
383 else if (m_diffType == "LFRSD")
384 {
385 c0 = 2.0 * p0 /
386 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
387 (ap0 * std::tgamma(p0 + 1)) *
388 (ap0 * std::tgamma(p0 + 1)));
389 }
390 else if (m_diffType == "LFRHU")
391 {
392 c0 = 2.0 * (p0 + 1.0) /
393 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
394 (ap0 * std::tgamma(p0 + 1)));
395 }
396 else if (m_diffType == "LFRcmin")
397 {
398 c0 =
399 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
400 (ap0 * std::tgamma(p0 + 1)));
401 }
402 else if (m_diffType == "LFRinf")
403 {
404 c0 = 10000000000000000.0;
405 }
406
407 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
408 (ap0 * std::tgamma(p0 + 1)) *
409 (ap0 * std::tgamma(p0 + 1));
410
411 NekDouble overeta0 = 1.0 / (1.0 + etap0);
412
413 // Derivative of the Legendre polynomials
414 // dLp = derivative of the Legendre polynomial of order p
415 // dLpp = derivative of the Legendre polynomial of order p+1
416 // dLpm = derivative of the Legendre polynomial of order p-1
417 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
418 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
419 0.0);
420 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
421 0.0);
422
423 // Building the DG_c_Left
424 for (i = 0; i < nquad0; ++i)
425 {
426 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
427 m_dGL_xi1[n][i] += dLpp0[i];
428 m_dGL_xi1[n][i] *= overeta0;
429 m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
430 m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
431 }
432
433 // Building the DG_c_Right
434 for (i = 0; i < nquad0; ++i)
435 {
436 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
437 m_dGR_xi1[n][i] += dLpp0[i];
438 m_dGR_xi1[n][i] *= overeta0;
439 m_dGR_xi1[n][i] += dLp0[i];
440 m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
441 }
442 }
443 break;
444 }
445 case 2:
446 {
447 LibUtilities::BasisSharedPtr BasisFR_Left0;
448 LibUtilities::BasisSharedPtr BasisFR_Right0;
449 LibUtilities::BasisSharedPtr BasisFR_Left1;
450 LibUtilities::BasisSharedPtr BasisFR_Right1;
451
452 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
453 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
454 m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
455 m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
456
457 for (n = 0; n < nElements; ++n)
458 {
459 base = pFields[0]->GetExp(n)->GetBase();
460 nquad0 = base[0]->GetNumPoints();
461 nquad1 = base[1]->GetNumPoints();
462 nmodes0 = base[0]->GetNumModes();
463 nmodes1 = base[1]->GetNumModes();
464
465 Array<OneD, const NekDouble> z0;
466 Array<OneD, const NekDouble> w0;
467 Array<OneD, const NekDouble> z1;
468 Array<OneD, const NekDouble> w1;
469
470 base[0]->GetZW(z0, w0);
471 base[1]->GetZW(z1, w1);
472
473 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
474 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
475 m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
476 m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
477
478 // Auxiliary vectors to build up the auxiliary
479 // derivatives of the Legendre polynomials
480 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
481 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
482 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
483 Array<OneD, NekDouble> dLp1(nquad1, 0.0);
484 Array<OneD, NekDouble> dLpp1(nquad1, 0.0);
485 Array<OneD, NekDouble> dLpm1(nquad1, 0.0);
486
487 // Degree of the correction functions
488 int p0 = nmodes0 - 1;
489 int p1 = nmodes1 - 1;
490
491 // Function sign to compute left correction function
492 NekDouble sign0 = pow(-1.0, p0);
493 NekDouble sign1 = pow(-1.0, p1);
494
495 // Factorial factor to build the scheme
496 NekDouble ap0 =
497 std::tgamma(2 * p0 + 1) /
498 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
499
500 NekDouble ap1 =
501 std::tgamma(2 * p1 + 1) /
502 (pow(2.0, p1) * std::tgamma(p1 + 1) * std::tgamma(p1 + 1));
503
504 // Scalar parameter which recovers the FR schemes
505 if (m_diffType == "LFRDG")
506 {
507 c0 = 0.0;
508 c1 = 0.0;
509 }
510 else if (m_diffType == "LFRSD")
511 {
512 c0 = 2.0 * p0 /
513 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
514 (ap0 * std::tgamma(p0 + 1)) *
515 (ap0 * std::tgamma(p0 + 1)));
516
517 c1 = 2.0 * p1 /
518 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
519 (ap1 * std::tgamma(p1 + 1)) *
520 (ap1 * std::tgamma(p1 + 1)));
521 }
522 else if (m_diffType == "LFRHU")
523 {
524 c0 = 2.0 * (p0 + 1.0) /
525 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
526 (ap0 * std::tgamma(p0 + 1)));
527
528 c1 = 2.0 * (p1 + 1.0) /
529 ((2.0 * p1 + 1.0) * p1 * (ap1 * std::tgamma(p1 + 1)) *
530 (ap1 * std::tgamma(p1 + 1)));
531 }
532 else if (m_diffType == "LFRcmin")
533 {
534 c0 =
535 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
536 (ap0 * std::tgamma(p0 + 1)));
537
538 c1 =
539 -2.0 / ((2.0 * p1 + 1.0) * (ap1 * std::tgamma(p1 + 1)) *
540 (ap1 * std::tgamma(p1 + 1)));
541 }
542 else if (m_diffType == "LFRinf")
543 {
544 c0 = 10000000000000000.0;
545 c1 = 10000000000000000.0;
546 }
547
548 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
549 (ap0 * std::tgamma(p0 + 1)) *
550 (ap0 * std::tgamma(p0 + 1));
551
552 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
553 (ap1 * std::tgamma(p1 + 1)) *
554 (ap1 * std::tgamma(p1 + 1));
555
556 NekDouble overeta0 = 1.0 / (1.0 + etap0);
557 NekDouble overeta1 = 1.0 / (1.0 + etap1);
558
559 // Derivative of the Legendre polynomials
560 // dLp = derivative of the Legendre polynomial of order p
561 // dLpp = derivative of the Legendre polynomial of order p+1
562 // dLpm = derivative of the Legendre polynomial of order p-1
563 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
564 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
565 0.0);
566 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
567 0.0);
568
569 Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
570 Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1 + 1, 0.0,
571 0.0);
572 Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1 - 1, 0.0,
573 0.0);
574
575 // Building the DG_c_Left
576 for (i = 0; i < nquad0; ++i)
577 {
578 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
579 m_dGL_xi1[n][i] += dLpp0[i];
580 m_dGL_xi1[n][i] *= overeta0;
581 m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
582 m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
583 }
584
585 // Building the DG_c_Left
586 for (i = 0; i < nquad1; ++i)
587 {
588 m_dGL_xi2[n][i] = etap1 * dLpm1[i];
589 m_dGL_xi2[n][i] += dLpp1[i];
590 m_dGL_xi2[n][i] *= overeta1;
591 m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
592 m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
593 }
594
595 // Building the DG_c_Right
596 for (i = 0; i < nquad0; ++i)
597 {
598 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
599 m_dGR_xi1[n][i] += dLpp0[i];
600 m_dGR_xi1[n][i] *= overeta0;
601 m_dGR_xi1[n][i] += dLp0[i];
602 m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
603 }
604
605 // Building the DG_c_Right
606 for (i = 0; i < nquad1; ++i)
607 {
608 m_dGR_xi2[n][i] = etap1 * dLpm1[i];
609 m_dGR_xi2[n][i] += dLpp1[i];
610 m_dGR_xi2[n][i] *= overeta1;
611 m_dGR_xi2[n][i] += dLp1[i];
612 m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
613 }
614 }
615 break;
616 }
617 case 3:
618 {
619 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
620 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
621 m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
622 m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
623 m_dGL_xi3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
624 m_dGR_xi3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
625
626 for (n = 0; n < nElements; ++n)
627 {
628 base = pFields[0]->GetExp(n)->GetBase();
629 nquad0 = base[0]->GetNumPoints();
630 nquad1 = base[1]->GetNumPoints();
631 nquad2 = base[2]->GetNumPoints();
632 nmodes0 = base[0]->GetNumModes();
633 nmodes1 = base[1]->GetNumModes();
634 nmodes2 = base[2]->GetNumModes();
635
636 Array<OneD, const NekDouble> z0;
637 Array<OneD, const NekDouble> w0;
638 Array<OneD, const NekDouble> z1;
639 Array<OneD, const NekDouble> w1;
640 Array<OneD, const NekDouble> z2;
641 Array<OneD, const NekDouble> w2;
642
643 base[0]->GetZW(z0, w0);
644 base[1]->GetZW(z1, w1);
645 base[1]->GetZW(z2, w2);
646
647 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
648 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
649 m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
650 m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
651 m_dGL_xi3[n] = Array<OneD, NekDouble>(nquad2);
652 m_dGR_xi3[n] = Array<OneD, NekDouble>(nquad2);
653
654 // Auxiliary vectors to build up the auxiliary
655 // derivatives of the Legendre polynomials
656 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
657 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
658 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
659 Array<OneD, NekDouble> dLp1(nquad1, 0.0);
660 Array<OneD, NekDouble> dLpp1(nquad1, 0.0);
661 Array<OneD, NekDouble> dLpm1(nquad1, 0.0);
662 Array<OneD, NekDouble> dLp2(nquad2, 0.0);
663 Array<OneD, NekDouble> dLpp2(nquad2, 0.0);
664 Array<OneD, NekDouble> dLpm2(nquad2, 0.0);
665
666 // Degree of the correction functions
667 int p0 = nmodes0 - 1;
668 int p1 = nmodes1 - 1;
669 int p2 = nmodes2 - 1;
670
671 // Function sign to compute left correction function
672 NekDouble sign0 = pow(-1.0, p0);
673 NekDouble sign1 = pow(-1.0, p1);
674 NekDouble sign2 = pow(-1.0, p2);
675
676 // Factorial factor to build the scheme
677 NekDouble ap0 =
678 std::tgamma(2 * p0 + 1) /
679 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
680
681 // Factorial factor to build the scheme
682 NekDouble ap1 =
683 std::tgamma(2 * p1 + 1) /
684 (pow(2.0, p1) * std::tgamma(p1 + 1) * std::tgamma(p1 + 1));
685
686 // Factorial factor to build the scheme
687 NekDouble ap2 =
688 std::tgamma(2 * p2 + 1) /
689 (pow(2.0, p2) * std::tgamma(p2 + 1) * std::tgamma(p2 + 1));
690
691 // Scalar parameter which recovers the FR schemes
692 if (m_diffType == "LFRDG")
693 {
694 c0 = 0.0;
695 c1 = 0.0;
696 c2 = 0.0;
697 }
698 else if (m_diffType == "LFRSD")
699 {
700 c0 = 2.0 * p0 /
701 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
702 (ap0 * std::tgamma(p0 + 1)) *
703 (ap0 * std::tgamma(p0 + 1)));
704
705 c1 = 2.0 * p1 /
706 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
707 (ap1 * std::tgamma(p1 + 1)) *
708 (ap1 * std::tgamma(p1 + 1)));
709
710 c2 = 2.0 * p2 /
711 ((2.0 * p2 + 1.0) * (p2 + 1.0) *
712 (ap2 * std::tgamma(p2 + 1)) *
713 (ap2 * std::tgamma(p2 + 1)));
714 }
715 else if (m_diffType == "LFRHU")
716 {
717 c0 = 2.0 * (p0 + 1.0) /
718 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
719 (ap0 * std::tgamma(p0 + 1)));
720
721 c1 = 2.0 * (p1 + 1.0) /
722 ((2.0 * p1 + 1.0) * p1 * (ap1 * std::tgamma(p1 + 1)) *
723 (ap1 * std::tgamma(p1 + 1)));
724
725 c2 = 2.0 * (p2 + 1.0) /
726 ((2.0 * p2 + 1.0) * p2 * (ap2 * std::tgamma(p2 + 1)) *
727 (ap2 * std::tgamma(p2 + 1)));
728 }
729 else if (m_diffType == "LFRcmin")
730 {
731 c0 =
732 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
733 (ap0 * std::tgamma(p0 + 1)));
734
735 c1 =
736 -2.0 / ((2.0 * p1 + 1.0) * (ap1 * std::tgamma(p1 + 1)) *
737 (ap1 * std::tgamma(p1 + 1)));
738
739 c2 =
740 -2.0 / ((2.0 * p2 + 1.0) * (ap2 * std::tgamma(p2 + 1)) *
741 (ap2 * std::tgamma(p2 + 1)));
742 }
743 else if (m_diffType == "LFRinf")
744 {
745 c0 = 10000000000000000.0;
746 c1 = 10000000000000000.0;
747 c2 = 10000000000000000.0;
748 }
749
750 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
751 (ap0 * std::tgamma(p0 + 1)) *
752 (ap0 * std::tgamma(p0 + 1));
753
754 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
755 (ap1 * std::tgamma(p1 + 1)) *
756 (ap1 * std::tgamma(p1 + 1));
757
758 NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0) *
759 (ap2 * std::tgamma(p2 + 1)) *
760 (ap2 * std::tgamma(p2 + 1));
761
762 NekDouble overeta0 = 1.0 / (1.0 + etap0);
763 NekDouble overeta1 = 1.0 / (1.0 + etap1);
764 NekDouble overeta2 = 1.0 / (1.0 + etap2);
765
766 // Derivative of the Legendre polynomials
767 // dLp = derivative of the Legendre polynomial of order p
768 // dLpp = derivative of the Legendre polynomial of order p+1
769 // dLpm = derivative of the Legendre polynomial of order p-1
770 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
771 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
772 0.0);
773 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
774 0.0);
775
776 Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
777 Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1 + 1, 0.0,
778 0.0);
779 Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1 - 1, 0.0,
780 0.0);
781
782 Polylib::jacobd(nquad2, z2.data(), &(dLp2[0]), p2, 0.0, 0.0);
783 Polylib::jacobd(nquad2, z2.data(), &(dLpp2[0]), p2 + 1, 0.0,
784 0.0);
785 Polylib::jacobd(nquad2, z2.data(), &(dLpm2[0]), p2 - 1, 0.0,
786 0.0);
787
788 // Building the DG_c_Left
789 for (i = 0; i < nquad0; ++i)
790 {
791 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
792 m_dGL_xi1[n][i] += dLpp0[i];
793 m_dGL_xi1[n][i] *= overeta0;
794 m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
795 m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
796 }
797
798 // Building the DG_c_Left
799 for (i = 0; i < nquad1; ++i)
800 {
801 m_dGL_xi2[n][i] = etap1 * dLpm1[i];
802 m_dGL_xi2[n][i] += dLpp1[i];
803 m_dGL_xi2[n][i] *= overeta1;
804 m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
805 m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
806 }
807
808 // Building the DG_c_Left
809 for (i = 0; i < nquad2; ++i)
810 {
811 m_dGL_xi3[n][i] = etap2 * dLpm2[i];
812 m_dGL_xi3[n][i] += dLpp2[i];
813 m_dGL_xi3[n][i] *= overeta2;
814 m_dGL_xi3[n][i] = dLp2[i] - m_dGL_xi3[n][i];
815 m_dGL_xi3[n][i] = 0.5 * sign2 * m_dGL_xi3[n][i];
816 }
817
818 // Building the DG_c_Right
819 for (i = 0; i < nquad0; ++i)
820 {
821 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
822 m_dGR_xi1[n][i] += dLpp0[i];
823 m_dGR_xi1[n][i] *= overeta0;
824 m_dGR_xi1[n][i] += dLp0[i];
825 m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
826 }
827
828 // Building the DG_c_Right
829 for (i = 0; i < nquad1; ++i)
830 {
831 m_dGR_xi2[n][i] = etap1 * dLpm1[i];
832 m_dGR_xi2[n][i] += dLpp1[i];
833 m_dGR_xi2[n][i] *= overeta1;
834 m_dGR_xi2[n][i] += dLp1[i];
835 m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
836 }
837
838 // Building the DG_c_Right
839 for (i = 0; i < nquad2; ++i)
840 {
841 m_dGR_xi3[n][i] = etap2 * dLpm2[i];
842 m_dGR_xi3[n][i] += dLpp2[i];
843 m_dGR_xi3[n][i] *= overeta2;
844 m_dGR_xi3[n][i] += dLp2[i];
845 m_dGR_xi3[n][i] = 0.5 * m_dGR_xi3[n][i];
846 }
847 }
848 break;
849 }
850 default:
851 {
852 ASSERTL0(false, "Expansion dimension not recognised");
853 break;
854 }
855 }
856}
857
858/**
859 * @brief Calculate FR Diffusion for the linear problems
860 * using an LDG interface flux.
861 *
862 */
863void DiffusionLFR::v_Diffuse(
864 const std::size_t nConvectiveFields,
865 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
866 const Array<OneD, Array<OneD, NekDouble>> &inarray,
867 Array<OneD, Array<OneD, NekDouble>> &outarray,
868 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pFwd,
869 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pBwd)
870{
871 int i, j, n;
872 int phys_offset;
873 Array<OneD, NekDouble> auxArray1, auxArray2;
874
875 Array<OneD, LibUtilities::BasisSharedPtr> Basis;
876 Basis = fields[0]->GetExp(0)->GetBase();
877
878 int nElements = fields[0]->GetExpSize();
879 int nDim = fields[0]->GetCoordim(0);
880 int nSolutionPts = fields[0]->GetTotPoints();
881 int nCoeffs = fields[0]->GetNcoeffs();
882
883 Array<OneD, Array<OneD, NekDouble>> outarrayCoeff(nConvectiveFields);
884 for (i = 0; i < nConvectiveFields; ++i)
885 {
886 outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
887 }
888
889 // Compute interface numerical fluxes for inarray in physical space
890 NumFluxforScalar(fields, inarray, m_IF1);
891
892 switch (nDim)
893 {
894 // 1D problems
895 case 1:
896 {
897 for (i = 0; i < nConvectiveFields; ++i)
898 {
899 // Computing the physical first-order discountinuous
900 // derivative
901 for (n = 0; n < nElements; n++)
902 {
903 phys_offset = fields[0]->GetPhys_Offset(n);
904
905 fields[i]->GetExp(n)->PhysDeriv(
906 0, auxArray1 = inarray[i] + phys_offset,
907 auxArray2 = m_DU1[i][0] + phys_offset);
908 }
909
910 // Computing the standard first-order correction
911 // derivative
912 DerCFlux_1D(nConvectiveFields, fields, inarray[i], m_IF1[i][0],
913 m_DFC1[i][0]);
914
915 // Back to the physical space using global operations
916 Vmath::Vdiv(nSolutionPts, &m_DFC1[i][0][0], 1, &m_jac[0], 1,
917 &m_DFC1[i][0][0], 1);
918 Vmath::Vdiv(nSolutionPts, &m_DFC1[i][0][0], 1, &m_jac[0], 1,
919 &m_DFC1[i][0][0], 1);
920
921 // Computing total first order derivatives
922 Vmath::Vadd(nSolutionPts, &m_DFC1[i][0][0], 1, &m_DU1[i][0][0],
923 1, &m_D1[i][0][0], 1);
924
925 Vmath::Vcopy(nSolutionPts, &m_D1[i][0][0], 1, &m_tmp1[i][0][0],
926 1);
927 }
928
929 // Computing interface numerical fluxes for m_D1
930 // in physical space
931 NumFluxforVector(fields, inarray, m_tmp1, m_IF2);
932
933 for (i = 0; i < nConvectiveFields; ++i)
934 {
935 // Computing the physical second-order discountinuous
936 // derivative
937 for (n = 0; n < nElements; n++)
938 {
939 phys_offset = fields[0]->GetPhys_Offset(n);
940
941 fields[i]->GetExp(n)->PhysDeriv(
942 0, auxArray1 = m_D1[i][0] + phys_offset,
943 auxArray2 = m_DD1[i][0] + phys_offset);
944 }
945
946 // Computing the standard second-order correction
947 // derivative
948 DerCFlux_1D(nConvectiveFields, fields, m_D1[i][0], m_IF2[i],
949 m_DFC2[i][0]);
950
951 // Back to the physical space using global operations
952 Vmath::Vdiv(nSolutionPts, &m_DFC2[i][0][0], 1, &m_jac[0], 1,
953 &m_DFC2[i][0][0], 1);
954 Vmath::Vdiv(nSolutionPts, &m_DFC2[i][0][0], 1, &m_jac[0], 1,
955 &m_DFC2[i][0][0], 1);
956
957 // Adding the total divergence to outarray (RHS)
958 Vmath::Vadd(nSolutionPts, &m_DFC2[i][0][0], 1, &m_DD1[i][0][0],
959 1, &outarray[i][0], 1);
960
961 // Primitive Dealiasing 1D
962 if (!(Basis[0]->Collocation()))
963 {
964 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
965 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
966 }
967 }
968 break;
969 }
970 // 2D problems
971 case 2:
972 {
973 for (i = 0; i < nConvectiveFields; ++i)
974 {
975 for (j = 0; j < nDim; ++j)
976 {
977 // Temporary vectors
978 Array<OneD, NekDouble> u1_hat(nSolutionPts, 0.0);
979 Array<OneD, NekDouble> u2_hat(nSolutionPts, 0.0);
980
981 if (j == 0)
982 {
983 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
984 &m_gmat[0][0], 1, &u1_hat[0], 1);
985
986 Vmath::Vmul(nSolutionPts, &u1_hat[0], 1, &m_jac[0], 1,
987 &u1_hat[0], 1);
988
989 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
990 &m_gmat[1][0], 1, &u2_hat[0], 1);
991
992 Vmath::Vmul(nSolutionPts, &u2_hat[0], 1, &m_jac[0], 1,
993 &u2_hat[0], 1);
994 }
995 else if (j == 1)
996 {
997 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
998 &m_gmat[2][0], 1, &u1_hat[0], 1);
999
1000 Vmath::Vmul(nSolutionPts, &u1_hat[0], 1, &m_jac[0], 1,
1001 &u1_hat[0], 1);
1002
1003 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
1004 &m_gmat[3][0], 1, &u2_hat[0], 1);
1005
1006 Vmath::Vmul(nSolutionPts, &u2_hat[0], 1, &m_jac[0], 1,
1007 &u2_hat[0], 1);
1008 }
1009
1010 for (n = 0; n < nElements; n++)
1011 {
1012 phys_offset = fields[0]->GetPhys_Offset(n);
1013
1014 fields[i]->GetExp(n)->StdPhysDeriv(
1015 auxArray1 = u1_hat + phys_offset,
1016 auxArray2 = m_tmp1[i][j] + phys_offset);
1017
1018 fields[i]->GetExp(n)->StdPhysDeriv(
1019 auxArray1 = u2_hat + phys_offset,
1020 NullNekDouble1DArray,
1021 auxArray2 = m_tmp2[i][j] + phys_offset);
1022 }
1023
1024 Vmath::Vadd(nSolutionPts, &m_tmp1[i][j][0], 1,
1025 &m_tmp2[i][j][0], 1, &m_DU1[i][j][0], 1);
1026
1027 // Divide by the metric jacobian
1028 Vmath::Vdiv(nSolutionPts, &m_DU1[i][j][0], 1, &m_jac[0], 1,
1029 &m_DU1[i][j][0], 1);
1030
1031 // Computing the standard first-order correction
1032 // derivatives
1033 DerCFlux_2D(nConvectiveFields, j, fields, inarray[i],
1034 m_IF1[i][j], m_DFC1[i][j]);
1035 }
1036
1037 // Multiplying derivatives by B matrix to get auxiliary
1038 // variables
1039 for (j = 0; j < nSolutionPts; ++j)
1040 {
1041 // std(q1)
1042 m_BD1[i][0][j] = (m_gmat[0][j] * m_gmat[0][j] +
1043 m_gmat[2][j] * m_gmat[2][j]) *
1044 m_DFC1[i][0][j] +
1045 (m_gmat[1][j] * m_gmat[0][j] +
1046 m_gmat[3][j] * m_gmat[2][j]) *
1047 m_DFC1[i][1][j];
1048
1049 // std(q2)
1050 m_BD1[i][1][j] = (m_gmat[1][j] * m_gmat[0][j] +
1051 m_gmat[3][j] * m_gmat[2][j]) *
1052 m_DFC1[i][0][j] +
1053 (m_gmat[1][j] * m_gmat[1][j] +
1054 m_gmat[3][j] * m_gmat[3][j]) *
1055 m_DFC1[i][1][j];
1056 }
1057
1058 // Multiplying derivatives by A^(-1) to get back
1059 // into the physical space
1060 for (j = 0; j < nSolutionPts; j++)
1061 {
1062 // q1 = A11^(-1)*std(q1) + A12^(-1)*std(q2)
1063 m_DFC1[i][0][j] = m_gmat[3][j] * m_BD1[i][0][j] -
1064 m_gmat[2][j] * m_BD1[i][1][j];
1065
1066 // q2 = A21^(-1)*std(q1) + A22^(-1)*std(q2)
1067 m_DFC1[i][1][j] = m_gmat[0][j] * m_BD1[i][1][j] -
1068 m_gmat[1][j] * m_BD1[i][0][j];
1069 }
1070
1071 // Computing the physical first-order derivatives
1072 for (j = 0; j < nDim; ++j)
1073 {
1074 Vmath::Vadd(nSolutionPts, &m_DU1[i][j][0], 1,
1075 &m_DFC1[i][j][0], 1, &m_D1[i][j][0], 1);
1076 }
1077 }
1078
1079 // Computing interface numerical fluxes for m_D1
1080 // in physical space
1081 NumFluxforVector(fields, inarray, m_D1, m_IF2);
1082
1083 // Computing the standard second-order discontinuous
1084 // derivatives
1085 for (i = 0; i < nConvectiveFields; ++i)
1086 {
1087 Array<OneD, NekDouble> f_hat(nSolutionPts, 0.0);
1088 Array<OneD, NekDouble> g_hat(nSolutionPts, 0.0);
1089
1090 for (j = 0; j < nSolutionPts; j++)
1091 {
1092 f_hat[j] = (m_D1[i][0][j] * m_gmat[0][j] +
1093 m_D1[i][1][j] * m_gmat[2][j]) *
1094 m_jac[j];
1095
1096 g_hat[j] = (m_D1[i][0][j] * m_gmat[1][j] +
1097 m_D1[i][1][j] * m_gmat[3][j]) *
1098 m_jac[j];
1099 }
1100
1101 for (n = 0; n < nElements; n++)
1102 {
1103 phys_offset = fields[0]->GetPhys_Offset(n);
1104
1105 fields[0]->GetExp(n)->StdPhysDeriv(
1106 auxArray1 = f_hat + phys_offset,
1107 auxArray2 = m_DD1[i][0] + phys_offset);
1108
1109 fields[0]->GetExp(n)->StdPhysDeriv(
1110 auxArray1 = g_hat + phys_offset, NullNekDouble1DArray,
1111 auxArray2 = m_DD1[i][1] + phys_offset);
1112 }
1113
1114 // Divergence of the standard discontinuous flux
1115 Vmath::Vadd(nSolutionPts, &m_DD1[i][0][0], 1, &m_DD1[i][1][0],
1116 1, &m_divFD[i][0], 1);
1117
1118 // Divergence of the standard correction flux
1119 if (Basis[0]->GetPointsType() ==
1120 LibUtilities::eGaussGaussLegendre &&
1121 Basis[1]->GetPointsType() ==
1122 LibUtilities::eGaussGaussLegendre)
1123 {
1124 DivCFlux_2D_Gauss(nConvectiveFields, fields, f_hat, g_hat,
1125 m_IF2[i], m_divFC[i]);
1126 }
1127 else
1128 {
1129 DivCFlux_2D(nConvectiveFields, fields, m_D1[i][0],
1130 m_D1[i][1], m_IF2[i], m_divFC[i]);
1131 }
1132
1133 // Divergence of the standard final flux
1134 Vmath::Vadd(nSolutionPts, &m_divFD[i][0], 1, &m_divFC[i][0], 1,
1135 &outarray[i][0], 1);
1136
1137 // Dividing by the jacobian to get back into
1138 // physical space
1139 Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1, &m_jac[0], 1,
1140 &outarray[i][0], 1);
1141
1142 // Primitive Dealiasing 2D
1143 if (!(Basis[0]->Collocation()))
1144 {
1145 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1146 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1147 }
1148 } // Close loop on nConvectiveFields
1149 break;
1150 }
1151 // 3D problems
1152 case 3:
1153 {
1154 ASSERTL0(false, "3D FRDG case not implemented yet");
1155 break;
1156 }
1157 }
1158}
1159
1160/**
1161 * @brief Builds the numerical flux for the 1st order derivatives
1162 *
1163 */
1164void DiffusionLFR::NumFluxforScalar(
1165 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1166 const Array<OneD, Array<OneD, NekDouble>> &ufield,
1167 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &uflux)
1168{
1169 int i, j;
1170 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1171 int nvariables = fields.size();
1172 int nDim = fields[0]->GetCoordim(0);
1173
1174 Array<OneD, NekDouble> Fwd(nTracePts);
1175 Array<OneD, NekDouble> Bwd(nTracePts);
1176 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
1177 Array<OneD, NekDouble> fluxtemp(nTracePts, 0.0);
1178
1179 // Get the normal velocity Vn
1180 for (i = 0; i < nDim; ++i)
1181 {
1182 Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1, Vn, 1, Vn, 1);
1183 }
1184
1185 // Get the sign of (v \cdot n), v = an arbitrary vector
1186 // Evaluate upwind flux:
1187 // uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
1188 for (j = 0; j < nDim; ++j)
1189 {
1190 for (i = 0; i < nvariables; ++i)
1191 {
1192 // Compute Fwd and Bwd value of ufield of i direction
1193 fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1194
1195 // if Vn >= 0, flux = uFwd, i.e.,
1196 // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
1197 // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
1198
1199 // else if Vn < 0, flux = uBwd, i.e.,
1200 // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
1201 // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
1202
1203 fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, fluxtemp);
1204
1205 // Imposing weak boundary condition with flux
1206 // if Vn >= 0, uflux = uBwd at Neumann, i.e.,
1207 // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
1208 // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
1209
1210 // if Vn >= 0, uflux = uFwd at Neumann, i.e.,
1211 // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
1212 // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
1213
1214 if (fields[0]->GetBndCondExpansions().size())
1215 {
1216 WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
1217 }
1218
1219 // if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
1220 // i.e,
1221 // edge::eForward, uFwd \‍(\tan_{\xi}^Fwd \cdot \vec{n})
1222 // edge::eBackward, uFwd \‍(\tan_{\xi}^Bwd \cdot \vec{n})
1223
1224 // else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}),
1225 // i.e,
1226 // edge::eForward, uBwd \‍(\tan_{\xi}^Fwd \cdot \vec{n})
1227 // edge::eBackward, uBwd \‍(\tan_{\xi}^Bwd \cdot \vec{n})
1228
1229 Vmath::Vcopy(nTracePts, fluxtemp, 1, uflux[i][j], 1);
1230 }
1231 }
1232}
1233
1234/**
1235 * @brief Imposes appropriate bcs for the 1st order derivatives
1236 *
1237 */
1238void DiffusionLFR::WeakPenaltyforScalar(
1239 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields, const int var,
1240 const Array<OneD, const NekDouble> &ufield,
1241 Array<OneD, NekDouble> &penaltyflux)
1242{
1243 // Variable initialisation
1244 int i, e, id1, id2;
1245 int nBndEdgePts, nBndEdges;
1246 int cnt = 0;
1247 int nBndRegions = fields[var]->GetBndCondExpansions().size();
1248 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1249
1250 // Extract physical values of the fields at the trace space
1251 Array<OneD, NekDouble> uplus(nTracePts);
1252 fields[var]->ExtractTracePhys(ufield, uplus);
1253
1254 // Impose boundary conditions
1255 for (i = 0; i < nBndRegions; ++i)
1256 {
1257 // Number of boundary edges related to bcRegion 'i'
1258 nBndEdges = fields[var]->GetBndCondExpansions()[i]->GetExpSize();
1259
1260 // Weakly impose boundary conditions by modifying flux values
1261 for (e = 0; e < nBndEdges; ++e)
1262 {
1263 // Number of points on the expansion
1264 nBndEdgePts = fields[var]
1265 ->GetBndCondExpansions()[i]
1266 ->GetExp(e)
1267 ->GetTotPoints();
1268
1269 // Offset of the boundary expansion
1270 id1 = fields[var]->GetBndCondExpansions()[i]->GetPhys_Offset(e);
1271
1272 // Offset of the trace space related to boundary expansion
1273 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1274 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1275
1276 // Dirichlet bcs ==> uflux = gD
1277 if (fields[var]
1278 ->GetBndConditions()[i]
1279 ->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1280 {
1281 Vmath::Vcopy(
1282 nBndEdgePts,
1283 &(fields[var]->GetBndCondExpansions()[i]->GetPhys())[id1],
1284 1, &penaltyflux[id2], 1);
1285 }
1286 // Neumann bcs ==> uflux = u+
1287 else if ((fields[var]->GetBndConditions()[i])
1288 ->GetBoundaryConditionType() ==
1289 SpatialDomains::eNeumann)
1290 {
1291 Vmath::Vcopy(nBndEdgePts, &uplus[id2], 1, &penaltyflux[id2], 1);
1292 }
1293 }
1294 }
1295}
1296
1297/**
1298 * @brief Build the numerical flux for the 2nd order derivatives
1299 *
1300 */
1301void DiffusionLFR::NumFluxforVector(
1302 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1303 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &ufield,
1304 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &qfield,
1305 Array<OneD, Array<OneD, NekDouble>> &qflux)
1306{
1307 int i, j;
1308 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1309 int nvariables = fields.size();
1310 int nDim = fields[0]->GetCoordim(0);
1311
1312 NekDouble C11 = 0.0;
1313 Array<OneD, NekDouble> Fwd(nTracePts);
1314 Array<OneD, NekDouble> Bwd(nTracePts);
1315 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
1316
1317 Array<OneD, NekDouble> qFwd(nTracePts);
1318 Array<OneD, NekDouble> qBwd(nTracePts);
1319 Array<OneD, NekDouble> qfluxtemp(nTracePts, 0.0);
1320 Array<OneD, NekDouble> uterm(nTracePts);
1321
1322 // Get the normal velocity Vn
1323 for (i = 0; i < nDim; ++i)
1324 {
1325 Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1, Vn, 1, Vn, 1);
1326 }
1327
1328 // Evaulate upwind flux:
1329 // qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
1330 for (i = 0; i < nvariables; ++i)
1331 {
1332 qflux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1333 for (j = 0; j < nDim; ++j)
1334 {
1335 // Compute Fwd and Bwd value of ufield of jth direction
1336 fields[i]->GetFwdBwdTracePhys(qfield[i][j], qFwd, qBwd);
1337
1338 // if Vn >= 0, flux = uFwd, i.e.,
1339 // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick
1340 // qflux = qBwd = q+
1341 // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick
1342 // qflux = qBwd = q-
1343
1344 // else if Vn < 0, flux = uBwd, i.e.,
1345 // edge::eForward, if V*n<0 <=> V*n_F<0, pick
1346 // qflux = qFwd = q-
1347 // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
1348 // qflux = qFwd = q+
1349
1350 fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1351
1352 Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
1353 qfluxtemp, 1);
1354
1355 // Imposing weak boundary condition with flux
1356 if (fields[0]->GetBndCondExpansions().size())
1357 {
1358 WeakPenaltyforVector(fields, i, j, qfield[i][j], qfluxtemp,
1359 C11);
1360 }
1361
1362 // q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
1363 // n_xi = n_x * tan_xi_x + n_y * tan_xi_y + n_z * tan_xi_z
1364 // n_xi = n_x * tan_eta_x + n_y * tan_eta_y + n_z*tan_eta_z
1365 Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1, qflux[i], 1);
1366 }
1367 }
1368}
1369
1370/**
1371 * @brief Imposes appropriate bcs for the 2nd order derivatives
1372 *
1373 */
1374void DiffusionLFR::WeakPenaltyforVector(
1375 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields, const int var,
1376 const int dir, const Array<OneD, const NekDouble> &qfield,
1377 Array<OneD, NekDouble> &penaltyflux, [[maybe_unused]] NekDouble C11)
1378{
1379 int i, e, id1, id2;
1380 int nBndEdges, nBndEdgePts;
1381 int nBndRegions = fields[var]->GetBndCondExpansions().size();
1382 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1383
1384 Array<OneD, NekDouble> uterm(nTracePts);
1385 Array<OneD, NekDouble> qtemp(nTracePts);
1386 int cnt = 0;
1387
1388 fields[var]->ExtractTracePhys(qfield, qtemp);
1389
1390 for (i = 0; i < nBndRegions; ++i)
1391 {
1392 nBndEdges = fields[var]->GetBndCondExpansions()[i]->GetExpSize();
1393
1394 // Weakly impose boundary conditions by modifying flux values
1395 for (e = 0; e < nBndEdges; ++e)
1396 {
1397 nBndEdgePts = fields[var]
1398 ->GetBndCondExpansions()[i]
1399 ->GetExp(e)
1400 ->GetTotPoints();
1401
1402 id1 = fields[var]->GetBndCondExpansions()[i]->GetPhys_Offset(e);
1403
1404 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1405 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1406
1407 // For Dirichlet boundary condition:
1408 // qflux = q+ - C_11 (u+ - g_D) (nx, ny)
1409 if (fields[var]
1410 ->GetBndConditions()[i]
1411 ->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1412 {
1413 Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1414 &qtemp[id2], 1, &penaltyflux[id2], 1);
1415 }
1416 // For Neumann boundary condition: qflux = g_N
1417 else if ((fields[var]->GetBndConditions()[i])
1418 ->GetBoundaryConditionType() ==
1419 SpatialDomains::eNeumann)
1420 {
1421 Vmath::Vmul(
1422 nBndEdgePts, &m_traceNormals[dir][id2], 1,
1423 &(fields[var]->GetBndCondExpansions()[i]->GetPhys())[id1],
1424 1, &penaltyflux[id2], 1);
1425 }
1426 }
1427 }
1428}
1429
1430/**
1431 * @brief Compute the derivative of the corrective flux for 1D problems.
1432 *
1433 * @param nConvectiveFields Number of fields.
1434 * @param fields Pointer to fields.
1435 * @param flux Volumetric flux in the physical space.
1436 * @param iFlux Numerical interface flux in physical space.
1437 * @param derCFlux Derivative of the corrective flux.
1438 *
1439 */
1440void DiffusionLFR::DerCFlux_1D(
1441 [[maybe_unused]] const int nConvectiveFields,
1442 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1443 const Array<OneD, const NekDouble> &flux,
1444 const Array<OneD, const NekDouble> &iFlux, Array<OneD, NekDouble> &derCFlux)
1445{
1446 int n;
1447 int nLocalSolutionPts, phys_offset, t_offset;
1448
1449 Array<OneD, NekDouble> auxArray1, auxArray2;
1450 Array<TwoD, const NekDouble> gmat;
1451 Array<OneD, const NekDouble> jac;
1452
1453 LibUtilities::PointsKeyVector ptsKeys;
1454 LibUtilities::BasisSharedPtr Basis;
1455 Basis = fields[0]->GetExp(0)->GetBasis(0);
1456
1457 int nElements = fields[0]->GetExpSize();
1458 int nSolutionPts = fields[0]->GetTotPoints();
1459
1460 vector<bool> leftAdjacentTraces =
1461 std::static_pointer_cast<MultiRegions::DisContField>(fields[0])
1462 ->GetLeftAdjacentTraces();
1463
1464 // Arrays to store the derivatives of the correction flux
1465 Array<OneD, NekDouble> DCL(nSolutionPts / nElements, 0.0);
1466 Array<OneD, NekDouble> DCR(nSolutionPts / nElements, 0.0);
1467
1468 // Arrays to store the intercell numerical flux jumps
1469 Array<OneD, NekDouble> JumpL(nElements);
1470 Array<OneD, NekDouble> JumpR(nElements);
1471
1472 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1473 fields[0]->GetTraceMap()->GetElmtToTrace();
1474
1475 for (n = 0; n < nElements; ++n)
1476 {
1477 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1478 phys_offset = fields[0]->GetPhys_Offset(n);
1479
1480 Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1481 Array<OneD, NekDouble> tmpFluxVertex(1, 0.0);
1482 Vmath::Vcopy(nLocalSolutionPts, &flux[phys_offset], 1, &tmparrayX1[0],
1483 1);
1484
1485 fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0], tmparrayX1,
1486 tmpFluxVertex);
1487
1488 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1489 elmtToTrace[n][0]->GetElmtId());
1490
1491 if (leftAdjacentTraces[2 * n])
1492 {
1493 JumpL[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1494 }
1495 else
1496 {
1497 JumpL[n] = iFlux[t_offset] - tmpFluxVertex[0];
1498 }
1499
1500 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1501 elmtToTrace[n][1]->GetElmtId());
1502
1503 fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1], tmparrayX1,
1504 tmpFluxVertex);
1505 if (leftAdjacentTraces[2 * n + 1])
1506 {
1507 JumpR[n] = iFlux[t_offset] - tmpFluxVertex[0];
1508 }
1509 else
1510 {
1511 JumpR[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1512 }
1513 }
1514
1515 for (n = 0; n < nElements; ++n)
1516 {
1517 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1518 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1519 phys_offset = fields[0]->GetPhys_Offset(n);
1520 jac = fields[0]
1521 ->GetExp(n)
1522 ->as<LocalRegions::Expansion1D>()
1523 ->GetGeom1D()
1524 ->GetMetricInfo()
1525 ->GetJac(ptsKeys);
1526
1527 JumpL[n] = JumpL[n] * jac[0];
1528 JumpR[n] = JumpR[n] * jac[0];
1529
1530 // Left jump multiplied by left derivative of C function
1531 Vmath::Smul(nLocalSolutionPts, JumpL[n], &m_dGL_xi1[n][0], 1, &DCL[0],
1532 1);
1533
1534 // Right jump multiplied by right derivative of C function
1535 Vmath::Smul(nLocalSolutionPts, JumpR[n], &m_dGR_xi1[n][0], 1, &DCR[0],
1536 1);
1537
1538 // Assembling derivative of the correction flux
1539 Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1540 &derCFlux[phys_offset], 1);
1541 }
1542}
1543
1544/**
1545 * @brief Compute the derivative of the corrective flux wrt a given
1546 * coordinate for 2D problems.
1547 *
1548 * There could be a bug for deformed elements since first derivatives
1549 * are not exactly the same results of DiffusionLDG as expected
1550 *
1551 * @param nConvectiveFields Number of fields.
1552 * @param fields Pointer to fields.
1553 * @param flux Volumetric flux in the physical space.
1554 * @param numericalFlux Numerical interface flux in physical space.
1555 * @param derCFlux Derivative of the corrective flux.
1556 *
1557 * \todo: Switch on shapes eventually here.
1558 */
1559void DiffusionLFR::DerCFlux_2D(
1560 [[maybe_unused]] const int nConvectiveFields, const int direction,
1561 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1562 const Array<OneD, const NekDouble> &flux,
1563 const Array<OneD, NekDouble> &iFlux, Array<OneD, NekDouble> &derCFlux)
1564{
1565 int n, e, i, j, cnt;
1566
1567 Array<OneD, const NekDouble> jac;
1568
1569 int nElements = fields[0]->GetExpSize();
1570 int trace_offset, phys_offset;
1571 int nLocalSolutionPts;
1572 int nquad0, nquad1;
1573 int nEdgePts;
1574
1575 LibUtilities::PointsKeyVector ptsKeys;
1576 Array<OneD, NekDouble> auxArray1, auxArray2;
1577 Array<OneD, LibUtilities::BasisSharedPtr> base;
1578 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1579 fields[0]->GetTraceMap()->GetElmtToTrace();
1580
1581 // Loop on the elements
1582 for (n = 0; n < nElements; ++n)
1583 {
1584 // Offset of the element on the global vector
1585 phys_offset = fields[0]->GetPhys_Offset(n);
1586 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1587 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1588
1589 jac = fields[0]
1590 ->GetExp(n)
1591 ->as<LocalRegions::Expansion2D>()
1592 ->GetGeom2D()
1593 ->GetMetricInfo()
1594 ->GetJac(ptsKeys);
1595
1596 base = fields[0]->GetExp(n)->GetBase();
1597 nquad0 = base[0]->GetNumPoints();
1598 nquad1 = base[1]->GetNumPoints();
1599
1600 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1601 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1602 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1603 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1604
1605 // Loop on the edges
1606 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1607 {
1608 // Number of edge points of edge 'e'
1609 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1610
1611 // Array for storing volumetric fluxes on each edge
1612 Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1613
1614 // Offset of the trace space correspondent to edge 'e'
1615 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1616 elmtToTrace[n][e]->GetElmtId());
1617
1618 // Extract the edge values of the volumetric fluxes
1619 // on edge 'e' and order them accordingly to the order
1620 // of the trace space
1621 fields[0]->GetExp(n)->GetTracePhysVals(
1622 e, elmtToTrace[n][e], flux + phys_offset, auxArray1 = tmparray);
1623
1624 // Compute the fluxJumps per each edge 'e' and each
1625 // flux point
1626 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1627 Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1, &tmparray[0], 1,
1628 &fluxJumps[0], 1);
1629
1630 // Check the ordering of the fluxJumps and reverse
1631 // it in case of backward definition of edge 'e'
1632 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1633 StdRegions::eBackwards)
1634 {
1635 Vmath::Reverse(nEdgePts, &fluxJumps[0], 1, &fluxJumps[0], 1);
1636 }
1637
1638 // Deformed elements
1639 if (fields[0]
1640 ->GetExp(n)
1641 ->as<LocalRegions::Expansion2D>()
1642 ->GetGeom2D()
1643 ->GetMetricInfo()
1644 ->GetGtype() == SpatialDomains::eDeformed)
1645 {
1646 // Extract the Jacobians along edge 'e'
1647 Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1648
1649 fields[0]->GetExp(n)->GetTracePhysVals(
1650 e, elmtToTrace[n][e], jac, auxArray1 = jacEdge);
1651
1652 // Check the ordering of the fluxJumps and reverse
1653 // it in case of backward definition of edge 'e'
1654 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1655 StdRegions::eBackwards)
1656 {
1657 Vmath::Reverse(nEdgePts, &jacEdge[0], 1, &jacEdge[0], 1);
1658 }
1659
1660 // Multiply the fluxJumps by the edge 'e' Jacobians
1661 // to bring the fluxJumps into the standard space
1662 for (j = 0; j < nEdgePts; j++)
1663 {
1664 fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1665 }
1666 }
1667 // Non-deformed elements
1668 else
1669 {
1670 // Multiply the fluxJumps by the edge 'e' Jacobians
1671 // to bring the fluxJumps into the standard space
1672 Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1, fluxJumps, 1);
1673 }
1674
1675 // Multiply jumps by derivatives of the correction functions
1676 // All the quntities at this level should be defined into
1677 // the standard space
1678 switch (e)
1679 {
1680 case 0:
1681 for (i = 0; i < nquad0; ++i)
1682 {
1683 for (j = 0; j < nquad1; ++j)
1684 {
1685 cnt = i + j * nquad0;
1686 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1687 }
1688 }
1689 break;
1690 case 1:
1691 for (i = 0; i < nquad1; ++i)
1692 {
1693 for (j = 0; j < nquad0; ++j)
1694 {
1695 cnt = (nquad0)*i + j;
1696 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1697 }
1698 }
1699 break;
1700 case 2:
1701 for (i = 0; i < nquad0; ++i)
1702 {
1703 for (j = 0; j < nquad1; ++j)
1704 {
1705 cnt = j * nquad0 + i;
1706 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1707 }
1708 }
1709 break;
1710 case 3:
1711 for (i = 0; i < nquad1; ++i)
1712 {
1713 for (j = 0; j < nquad0; ++j)
1714 {
1715 cnt = j + i * nquad0;
1716 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1717 }
1718 }
1719 break;
1720
1721 default:
1722 ASSERTL0(false, "edge value (< 3) is out of range");
1723 break;
1724 }
1725 }
1726
1727 // Summing the component of the flux for calculating the
1728 // derivatives per each direction
1729 if (direction == 0)
1730 {
1731 Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1, &divCFluxE3[0], 1,
1732 &derCFlux[phys_offset], 1);
1733 }
1734 else if (direction == 1)
1735 {
1736 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE2[0], 1,
1737 &derCFlux[phys_offset], 1);
1738 }
1739 }
1740}
1741
1742/**
1743 * @brief Compute the divergence of the corrective flux for 2D problems.
1744 *
1745 * @param nConvectiveFields Number of fields.
1746 * @param fields Pointer to fields.
1747 * @param fluxX1 X1 - volumetric flux in physical space.
1748 * @param fluxX2 X2 - volumetric flux in physical space.
1749 * @param numericalFlux Interface flux in physical space.
1750 * @param divCFlux Divergence of the corrective flux for
1751 * 2D problems.
1752 *
1753 * \todo: Switch on shapes eventually here.
1754 */
1755void DiffusionLFR::DivCFlux_2D(
1756 [[maybe_unused]] const int nConvectiveFields,
1757 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1758 const Array<OneD, const NekDouble> &fluxX1,
1759 const Array<OneD, const NekDouble> &fluxX2,
1760 const Array<OneD, const NekDouble> &numericalFlux,
1761 Array<OneD, NekDouble> &divCFlux)
1762{
1763 int n, e, i, j, cnt;
1764
1765 int nElements = fields[0]->GetExpSize();
1766
1767 int nLocalSolutionPts;
1768 int nEdgePts;
1769 int trace_offset;
1770 int phys_offset;
1771 int nquad0;
1772 int nquad1;
1773
1774 Array<OneD, NekDouble> auxArray1, auxArray2;
1775 Array<OneD, LibUtilities::BasisSharedPtr> base;
1776
1777 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1778 fields[0]->GetTraceMap()->GetElmtToTrace();
1779
1780 // Loop on the elements
1781 for (n = 0; n < nElements; ++n)
1782 {
1783 // Offset of the element on the global vector
1784 phys_offset = fields[0]->GetPhys_Offset(n);
1785 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1786
1787 base = fields[0]->GetExp(n)->GetBase();
1788 nquad0 = base[0]->GetNumPoints();
1789 nquad1 = base[1]->GetNumPoints();
1790
1791 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1792 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1793 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1794 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1795
1796 // Loop on the edges
1797 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1798 {
1799 // Number of edge points of edge e
1800 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1801
1802 Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1803 Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1804 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1805 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1806 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1807
1808 // Offset of the trace space correspondent to edge e
1809 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1810 elmtToTrace[n][e]->GetElmtId());
1811
1812 // Get the normals of edge e
1813 const Array<OneD, const Array<OneD, NekDouble>> &normals =
1814 fields[0]->GetExp(n)->GetTraceNormal(e);
1815
1816 // Extract the edge values of flux-x on edge e and order
1817 // them accordingly to the order of the trace space
1818 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1819 fluxX1 + phys_offset,
1820 auxArray1 = tmparrayX1);
1821
1822 // Extract the edge values of flux-y on edge e and order
1823 // them accordingly to the order of the trace space
1824 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1825 fluxX2 + phys_offset,
1826 auxArray1 = tmparrayX2);
1827
1828 // Multiply the edge components of the flux by the normal
1829 for (i = 0; i < nEdgePts; ++i)
1830 {
1831 fluxN[i] = tmparrayX1[i] * m_traceNormals[0][trace_offset + i] +
1832 tmparrayX2[i] * m_traceNormals[1][trace_offset + i];
1833 }
1834
1835 // Subtract to the Riemann flux the discontinuous flux
1836 Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1, &fluxN[0], 1,
1837 &fluxJumps[0], 1);
1838
1839 // Check the ordering of the jump vectors
1840 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1841 StdRegions::eBackwards)
1842 {
1843 Vmath::Reverse(nEdgePts, auxArray1 = fluxJumps, 1,
1844 auxArray2 = fluxJumps, 1);
1845 }
1846
1847 for (i = 0; i < nEdgePts; ++i)
1848 {
1849 if (m_traceNormals[0][trace_offset + i] != normals[0][i] ||
1850 m_traceNormals[1][trace_offset + i] != normals[1][i])
1851 {
1852 fluxJumps[i] = -fluxJumps[i];
1853 }
1854 }
1855
1856 // Multiply jumps by derivatives of the correction functions
1857 switch (e)
1858 {
1859 case 0:
1860 for (i = 0; i < nquad0; ++i)
1861 {
1862 // Multiply fluxJumps by Q factors
1863 fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1864
1865 for (j = 0; j < nquad1; ++j)
1866 {
1867 cnt = i + j * nquad0;
1868 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1869 }
1870 }
1871 break;
1872 case 1:
1873 for (i = 0; i < nquad1; ++i)
1874 {
1875 // Multiply fluxJumps by Q factors
1876 fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1877
1878 for (j = 0; j < nquad0; ++j)
1879 {
1880 cnt = (nquad0)*i + j;
1881 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1882 }
1883 }
1884 break;
1885 case 2:
1886 for (i = 0; i < nquad0; ++i)
1887 {
1888 // Multiply fluxJumps by Q factors
1889 fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1890
1891 for (j = 0; j < nquad1; ++j)
1892 {
1893 cnt = j * nquad0 + i;
1894 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1895 }
1896 }
1897 break;
1898 case 3:
1899 for (i = 0; i < nquad1; ++i)
1900 {
1901 // Multiply fluxJumps by Q factors
1902 fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1903 for (j = 0; j < nquad0; ++j)
1904 {
1905 cnt = j + i * nquad0;
1906 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1907 }
1908 }
1909 break;
1910
1911 default:
1912 ASSERTL0(false, "edge value (< 3) is out of range");
1913 break;
1914 }
1915 }
1916
1917 // Sum each edge contribution
1918 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
1919 &divCFlux[phys_offset], 1);
1920
1921 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1922 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1923
1924 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1925 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1926 }
1927}
1928
1929/**
1930 * @brief Compute the divergence of the corrective flux for 2D problems
1931 * where POINTSTYPE="GaussGaussLegendre"
1932 *
1933 * @param nConvectiveFields Number of fields.
1934 * @param fields Pointer to fields.
1935 * @param fluxX1 X1-volumetric flux in physical space.
1936 * @param fluxX2 X2-volumetric flux in physical space.
1937 * @param numericalFlux Interface flux in physical space.
1938 * @param divCFlux Divergence of the corrective flux.
1939 *
1940 * \todo: Switch on shapes eventually here.
1941 */
1942
1943void DiffusionLFR::DivCFlux_2D_Gauss(
1944 [[maybe_unused]] const int nConvectiveFields,
1945 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1946 const Array<OneD, const NekDouble> &fluxX1,
1947 const Array<OneD, const NekDouble> &fluxX2,
1948 const Array<OneD, const NekDouble> &numericalFlux,
1949 Array<OneD, NekDouble> &divCFlux)
1950{
1951 int n, e, i, j, cnt;
1952
1953 int nElements = fields[0]->GetExpSize();
1954 int nLocalSolutionPts;
1955 int nEdgePts;
1956 int trace_offset;
1957 int phys_offset;
1958 int nquad0;
1959 int nquad1;
1960
1961 Array<OneD, NekDouble> auxArray1, auxArray2;
1962 Array<OneD, LibUtilities::BasisSharedPtr> base;
1963
1964 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1965 fields[0]->GetTraceMap()->GetElmtToTrace();
1966
1967 // Loop on the elements
1968 for (n = 0; n < nElements; ++n)
1969 {
1970 // Offset of the element on the global vector
1971 phys_offset = fields[0]->GetPhys_Offset(n);
1972 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1973
1974 base = fields[0]->GetExp(n)->GetBase();
1975 nquad0 = base[0]->GetNumPoints();
1976 nquad1 = base[1]->GetNumPoints();
1977
1978 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1979 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1980 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1981 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1982
1983 // Loop on the edges
1984 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1985 {
1986 // Number of edge points of edge e
1987 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1988
1989 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1990 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1991 Array<OneD, NekDouble> fluxN_R(nEdgePts, 0.0);
1992 Array<OneD, NekDouble> fluxN_D(nEdgePts, 0.0);
1993 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1994
1995 // Offset of the trace space correspondent to edge e
1996 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1997 elmtToTrace[n][e]->GetElmtId());
1998
1999 // Get the normals of edge e
2000 const Array<OneD, const Array<OneD, NekDouble>> &normals =
2001 fields[0]->GetExp(n)->GetTraceNormal(e);
2002
2003 // Extract the trasformed normal flux at each edge
2004 switch (e)
2005 {
2006 case 0:
2007 // Extract the edge values of transformed flux-y on
2008 // edge e and order them accordingly to the order of
2009 // the trace space
2010 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2011 fluxX2 + phys_offset,
2012 auxArray1 = fluxN_D);
2013
2014 Vmath::Neg(nEdgePts, fluxN_D, 1);
2015
2016 // Extract the physical Rieamann flux at each edge
2017 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2018 &fluxN[0], 1);
2019
2020 // Check the ordering of vectors
2021 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2022 StdRegions::eBackwards)
2023 {
2024 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2025 auxArray2 = fluxN, 1);
2026
2027 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2028 auxArray2 = fluxN_D, 1);
2029 }
2030
2031 // Transform Riemann Fluxes in the standard element
2032 for (i = 0; i < nquad0; ++i)
2033 {
2034 // Multiply Riemann Flux by Q factors
2035 fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
2036 }
2037
2038 for (i = 0; i < nEdgePts; ++i)
2039 {
2040 if (m_traceNormals[0][trace_offset + i] !=
2041 normals[0][i] ||
2042 m_traceNormals[1][trace_offset + i] !=
2043 normals[1][i])
2044 {
2045 fluxN_R[i] = -fluxN_R[i];
2046 }
2047 }
2048
2049 // Subtract to the Riemann flux the discontinuous
2050 // flux
2051 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2052 &fluxJumps[0], 1);
2053
2054 // Multiplicate the flux jumps for the correction
2055 // function
2056 for (i = 0; i < nquad0; ++i)
2057 {
2058 for (j = 0; j < nquad1; ++j)
2059 {
2060 cnt = i + j * nquad0;
2061 divCFluxE0[cnt] = -fluxJumps[i] * m_dGL_xi2[n][j];
2062 }
2063 }
2064 break;
2065 case 1:
2066 // Extract the edge values of transformed flux-x on
2067 // edge e and order them accordingly to the order of
2068 // the trace space
2069 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2070 fluxX1 + phys_offset,
2071 auxArray1 = fluxN_D);
2072
2073 // Extract the physical Rieamann flux at each edge
2074 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2075 &fluxN[0], 1);
2076
2077 // Check the ordering of vectors
2078 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2079 StdRegions::eBackwards)
2080 {
2081 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2082 auxArray2 = fluxN, 1);
2083
2084 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2085 auxArray2 = fluxN_D, 1);
2086 }
2087
2088 // Transform Riemann Fluxes in the standard element
2089 for (i = 0; i < nquad1; ++i)
2090 {
2091 // Multiply Riemann Flux by Q factors
2092 fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
2093 }
2094
2095 for (i = 0; i < nEdgePts; ++i)
2096 {
2097 if (m_traceNormals[0][trace_offset + i] !=
2098 normals[0][i] ||
2099 m_traceNormals[1][trace_offset + i] !=
2100 normals[1][i])
2101 {
2102 fluxN_R[i] = -fluxN_R[i];
2103 }
2104 }
2105
2106 // Subtract to the Riemann flux the discontinuous
2107 // flux
2108 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2109 &fluxJumps[0], 1);
2110
2111 // Multiplicate the flux jumps for the correction
2112 // function
2113 for (i = 0; i < nquad1; ++i)
2114 {
2115 for (j = 0; j < nquad0; ++j)
2116 {
2117 cnt = (nquad0)*i + j;
2118 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
2119 }
2120 }
2121 break;
2122 case 2:
2123
2124 // Extract the edge values of transformed flux-y on
2125 // edge e and order them accordingly to the order of
2126 // the trace space
2127
2128 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2129 fluxX2 + phys_offset,
2130 auxArray1 = fluxN_D);
2131
2132 // Extract the physical Rieamann flux at each edge
2133 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2134 &fluxN[0], 1);
2135
2136 // Check the ordering of vectors
2137 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2138 StdRegions::eBackwards)
2139 {
2140 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2141 auxArray2 = fluxN, 1);
2142
2143 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2144 auxArray2 = fluxN_D, 1);
2145 }
2146
2147 // Transform Riemann Fluxes in the standard element
2148 for (i = 0; i < nquad0; ++i)
2149 {
2150 // Multiply Riemann Flux by Q factors
2151 fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
2152 }
2153
2154 for (i = 0; i < nEdgePts; ++i)
2155 {
2156 if (m_traceNormals[0][trace_offset + i] !=
2157 normals[0][i] ||
2158 m_traceNormals[1][trace_offset + i] !=
2159 normals[1][i])
2160 {
2161 fluxN_R[i] = -fluxN_R[i];
2162 }
2163 }
2164
2165 // Subtract to the Riemann flux the discontinuous
2166 // flux
2167
2168 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2169 &fluxJumps[0], 1);
2170
2171 // Multiplicate the flux jumps for the correction
2172 // function
2173 for (i = 0; i < nquad0; ++i)
2174 {
2175 for (j = 0; j < nquad1; ++j)
2176 {
2177 cnt = j * nquad0 + i;
2178 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
2179 }
2180 }
2181 break;
2182 case 3:
2183 // Extract the edge values of transformed flux-x on
2184 // edge e and order them accordingly to the order of
2185 // the trace space
2186
2187 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2188 fluxX1 + phys_offset,
2189 auxArray1 = fluxN_D);
2190 Vmath::Neg(nEdgePts, fluxN_D, 1);
2191
2192 // Extract the physical Rieamann flux at each edge
2193 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2194 &fluxN[0], 1);
2195
2196 // Check the ordering of vectors
2197 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2198 StdRegions::eBackwards)
2199 {
2200 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2201 auxArray2 = fluxN, 1);
2202
2203 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2204 auxArray2 = fluxN_D, 1);
2205 }
2206
2207 // Transform Riemann Fluxes in the standard element
2208 for (i = 0; i < nquad1; ++i)
2209 {
2210 // Multiply Riemann Flux by Q factors
2211 fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
2212 }
2213
2214 for (i = 0; i < nEdgePts; ++i)
2215 {
2216 if (m_traceNormals[0][trace_offset + i] !=
2217 normals[0][i] ||
2218 m_traceNormals[1][trace_offset + i] !=
2219 normals[1][i])
2220 {
2221 fluxN_R[i] = -fluxN_R[i];
2222 }
2223 }
2224
2225 // Subtract to the Riemann flux the discontinuous
2226 // flux
2227
2228 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2229 &fluxJumps[0], 1);
2230
2231 // Multiplicate the flux jumps for the correction
2232 // function
2233 for (i = 0; i < nquad1; ++i)
2234 {
2235 for (j = 0; j < nquad0; ++j)
2236 {
2237 cnt = j + i * nquad0;
2238 divCFluxE3[cnt] = -fluxJumps[i] * m_dGL_xi1[n][j];
2239 }
2240 }
2241 break;
2242 default:
2243 ASSERTL0(false, "edge value (< 3) is out of range");
2244 break;
2245 }
2246 }
2247
2248 // Sum each edge contribution
2249 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
2250 &divCFlux[phys_offset], 1);
2251
2252 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2253 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2254
2255 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2256 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2257 }
2258}
2259
2260} // namespace Nektar::SolverUtils
ASSERTL0
#define ASSERTL0(condition, msg)
Definition ErrorUtil.hpp:216
WARNINGL0
#define WARNINGL0(condition, msg)
Definition ErrorUtil.hpp:223
Nektar::LibUtilities::Basis
Represents a basis of a given type.
Definition Basis.h:198
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition BasicUtils/NekFactory.hpp:196
Nektar::LocalRegions::Expansion::GetMetricInfo
const SpatialDomains::GeomFactorsSharedPtr & GetMetricInfo() const
Definition Expansion.cpp:246
Nektar::SolverUtils::DiffusionLFR::m_tmp2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Definition DiffusionLFR.h:101
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition DiffusionLFR.h:82
Nektar::SolverUtils::DiffusionLFR::v_Diffuse
void v_Diffuse(const size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd) override
Calculate FR Diffusion for the linear problems using an LDG interface flux.
Definition DiffusionLFR.cpp:863
Nektar::SolverUtils::DiffusionLFR::m_IF1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
Definition DiffusionLFR.h:89
Nektar::SolverUtils::DiffusionLFR::m_jac
Array< OneD, NekDouble > m_jac
Definition DiffusionLFR.h:72
Nektar::SolverUtils::DiffusionLFR::DerCFlux_1D
void DerCFlux_1D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, const NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux for 1D problems.
Definition DiffusionLFR.cpp:1440
Nektar::SolverUtils::DiffusionLFR::m_DFC1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Definition DiffusionLFR.h:91
Nektar::SolverUtils::DiffusionLFR::m_DD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
Definition DiffusionLFR.h:94
Nektar::SolverUtils::DiffusionLFR::DiffusionLFR
DiffusionLFR(std::string diffType)
DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion term....
Definition DiffusionLFR.cpp:67
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition DiffusionLFR.h:81
Nektar::SolverUtils::DiffusionLFR::m_divFC
Array< OneD, Array< OneD, NekDouble > > m_divFC
Definition DiffusionLFR.h:98
Nektar::SolverUtils::DiffusionLFR::DerCFlux_2D
void DerCFlux_2D(const int nConvectiveFields, const int direction, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux wrt a given coordinate for 2D problems.
Definition DiffusionLFR.cpp:1559
Nektar::SolverUtils::DiffusionLFR::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 DiffusionLFR.cpp:164
Nektar::SolverUtils::DiffusionLFR::m_session
LibUtilities::SessionReaderSharedPtr m_session
Definition DiffusionLFR.h:70
Nektar::SolverUtils::DiffusionLFR::WeakPenaltyforScalar
void WeakPenaltyforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const Array< OneD, const NekDouble > &ufield, Array< OneD, NekDouble > &penaltyflux)
Imposes appropriate bcs for the 1st order derivatives.
Definition DiffusionLFR.cpp:1238
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition DiffusionLFR.h:83
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition DiffusionLFR.h:80
Nektar::SolverUtils::DiffusionLFR::WeakPenaltyforVector
void WeakPenaltyforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux, NekDouble C11)
Imposes appropriate bcs for the 2nd order derivatives.
Definition DiffusionLFR.cpp:1374
Nektar::SolverUtils::DiffusionLFR::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 DiffusionLFR.cpp:326
Nektar::SolverUtils::DiffusionLFR::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 DiffusionLFR.cpp:1755
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition DiffusionLFR.h:85
Nektar::SolverUtils::DiffusionLFR::NumFluxforVector
void NumFluxforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
Build the numerical flux for the 2nd order derivatives.
Definition DiffusionLFR.cpp:1301
Nektar::SolverUtils::DiffusionLFR::create
static DiffusionSharedPtr create(std::string diffType)
Definition DiffusionLFR.h:45
Nektar::SolverUtils::DiffusionLFR::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition DiffusionLFR.h:73
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition DiffusionLFR.h:69
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition DiffusionLFR.h:75
Nektar::SolverUtils::DiffusionLFR::m_DU1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Definition DiffusionLFR.h:90
Nektar::SolverUtils::DiffusionLFR::m_DFC2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Definition DiffusionLFR.h:96
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition DiffusionLFR.h:78
Nektar::SolverUtils::DiffusionLFR::m_tmp1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Definition DiffusionLFR.h:100
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition DiffusionLFR.h:77
Nektar::SolverUtils::DiffusionLFR::m_D1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
Definition DiffusionLFR.h:93
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition DiffusionLFR.h:84
Nektar::SolverUtils::DiffusionLFR::NumFluxforScalar
void NumFluxforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
Builds the numerical flux for the 1st order derivatives.
Definition DiffusionLFR.cpp:1164
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition DiffusionLFR.h:76
Nektar::SolverUtils::DiffusionLFR::v_InitObject
void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
Initiliase DiffusionLFR objects and store them before starting the time-stepping.
Definition DiffusionLFR.cpp:81
Nektar::SolverUtils::DiffusionLFR::m_IF2
Array< OneD, Array< OneD, NekDouble > > m_IF2
Definition DiffusionLFR.h:95
Nektar::SolverUtils::DiffusionLFR::m_divFD
Array< OneD, Array< OneD, NekDouble > > m_divFD
Definition DiffusionLFR.h:97
Nektar::SolverUtils::DiffusionLFR::m_BD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
Definition DiffusionLFR.h:92
Nektar::SolverUtils::DiffusionLFR::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 DiffusionLFR.cpp:1943
Nektar::LibUtilities::BasisSharedPtr
std::shared_ptr< Basis > BasisSharedPtr
Definition FoundationsFwd.hpp:69
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition SessionReader.h:119
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition Points.h:231
Nektar::LibUtilities::eGaussGaussLegendre
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition PointsType.h:46
Nektar::SolverUtils::GetDiffusionFactory
DiffusionFactory & GetDiffusionFactory()
Definition Diffusion.cpp:39
Nektar::SpatialDomains::eDeformed
@ eDeformed
Geometry is curved or has non-constant factors.
Definition SpatialDomains.hpp:54
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition BasicUtils/SharedArray.hpp:928
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::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvp (scalar times vector plus vector): z = alpha*x + y.
Definition Vmath.hpp:396
Vmath::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::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition Vmath.hpp:825
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition Vmath.hpp:220
std
STL namespace.
Generated on Sun Nov 9 2025 20:16:59 for Nektar++ by doxygen 1.9.8