Nektar++
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 0, auxArray1 = u1_hat + phys_offset,
1016 auxArray2 = m_tmp1[i][j] + phys_offset);
1017
1018 fields[i]->GetExp(n)->StdPhysDeriv(
1019 1, auxArray1 = u2_hat + phys_offset,
1020 auxArray2 = m_tmp2[i][j] + phys_offset);
1021 }
1022
1023 Vmath::Vadd(nSolutionPts, &m_tmp1[i][j][0], 1,
1024 &m_tmp2[i][j][0], 1, &m_DU1[i][j][0], 1);
1025
1026 // Divide by the metric jacobian
1027 Vmath::Vdiv(nSolutionPts, &m_DU1[i][j][0], 1, &m_jac[0], 1,
1028 &m_DU1[i][j][0], 1);
1029
1030 // Computing the standard first-order correction
1031 // derivatives
1032 DerCFlux_2D(nConvectiveFields, j, fields, inarray[i],
1033 m_IF1[i][j], m_DFC1[i][j]);
1034 }
1035
1036 // Multiplying derivatives by B matrix to get auxiliary
1037 // variables
1038 for (j = 0; j < nSolutionPts; ++j)
1039 {
1040 // std(q1)
1041 m_BD1[i][0][j] = (m_gmat[0][j] * m_gmat[0][j] +
1042 m_gmat[2][j] * m_gmat[2][j]) *
1043 m_DFC1[i][0][j] +
1044 (m_gmat[1][j] * m_gmat[0][j] +
1045 m_gmat[3][j] * m_gmat[2][j]) *
1046 m_DFC1[i][1][j];
1047
1048 // std(q2)
1049 m_BD1[i][1][j] = (m_gmat[1][j] * m_gmat[0][j] +
1050 m_gmat[3][j] * m_gmat[2][j]) *
1051 m_DFC1[i][0][j] +
1052 (m_gmat[1][j] * m_gmat[1][j] +
1053 m_gmat[3][j] * m_gmat[3][j]) *
1054 m_DFC1[i][1][j];
1055 }
1056
1057 // Multiplying derivatives by A^(-1) to get back
1058 // into the physical space
1059 for (j = 0; j < nSolutionPts; j++)
1060 {
1061 // q1 = A11^(-1)*std(q1) + A12^(-1)*std(q2)
1062 m_DFC1[i][0][j] = m_gmat[3][j] * m_BD1[i][0][j] -
1063 m_gmat[2][j] * m_BD1[i][1][j];
1064
1065 // q2 = A21^(-1)*std(q1) + A22^(-1)*std(q2)
1066 m_DFC1[i][1][j] = m_gmat[0][j] * m_BD1[i][1][j] -
1067 m_gmat[1][j] * m_BD1[i][0][j];
1068 }
1069
1070 // Computing the physical first-order derivatives
1071 for (j = 0; j < nDim; ++j)
1072 {
1073 Vmath::Vadd(nSolutionPts, &m_DU1[i][j][0], 1,
1074 &m_DFC1[i][j][0], 1, &m_D1[i][j][0], 1);
1075 }
1076 }
1077
1078 // Computing interface numerical fluxes for m_D1
1079 // in physical space
1080 NumFluxforVector(fields, inarray, m_D1, m_IF2);
1081
1082 // Computing the standard second-order discontinuous
1083 // derivatives
1084 for (i = 0; i < nConvectiveFields; ++i)
1085 {
1086 Array<OneD, NekDouble> f_hat(nSolutionPts, 0.0);
1087 Array<OneD, NekDouble> g_hat(nSolutionPts, 0.0);
1088
1089 for (j = 0; j < nSolutionPts; j++)
1090 {
1091 f_hat[j] = (m_D1[i][0][j] * m_gmat[0][j] +
1092 m_D1[i][1][j] * m_gmat[2][j]) *
1093 m_jac[j];
1094
1095 g_hat[j] = (m_D1[i][0][j] * m_gmat[1][j] +
1096 m_D1[i][1][j] * m_gmat[3][j]) *
1097 m_jac[j];
1098 }
1099
1100 for (n = 0; n < nElements; n++)
1101 {
1102 phys_offset = fields[0]->GetPhys_Offset(n);
1103
1104 fields[0]->GetExp(n)->StdPhysDeriv(
1105 0, auxArray1 = f_hat + phys_offset,
1106 auxArray2 = m_DD1[i][0] + phys_offset);
1107
1108 fields[0]->GetExp(n)->StdPhysDeriv(
1109 1, auxArray1 = g_hat + phys_offset,
1110 auxArray2 = m_DD1[i][1] + phys_offset);
1111 }
1112
1113 // Divergence of the standard discontinuous flux
1114 Vmath::Vadd(nSolutionPts, &m_DD1[i][0][0], 1, &m_DD1[i][1][0],
1115 1, &m_divFD[i][0], 1);
1116
1117 // Divergence of the standard correction flux
1118 if (Basis[0]->GetPointsType() ==
1119 LibUtilities::eGaussGaussLegendre &&
1120 Basis[1]->GetPointsType() ==
1121 LibUtilities::eGaussGaussLegendre)
1122 {
1123 DivCFlux_2D_Gauss(nConvectiveFields, fields, f_hat, g_hat,
1124 m_IF2[i], m_divFC[i]);
1125 }
1126 else
1127 {
1128 DivCFlux_2D(nConvectiveFields, fields, m_D1[i][0],
1129 m_D1[i][1], m_IF2[i], m_divFC[i]);
1130 }
1131
1132 // Divergence of the standard final flux
1133 Vmath::Vadd(nSolutionPts, &m_divFD[i][0], 1, &m_divFC[i][0], 1,
1134 &outarray[i][0], 1);
1135
1136 // Dividing by the jacobian to get back into
1137 // physical space
1138 Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1, &m_jac[0], 1,
1139 &outarray[i][0], 1);
1140
1141 // Primitive Dealiasing 2D
1142 if (!(Basis[0]->Collocation()))
1143 {
1144 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1145 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1146 }
1147 } // Close loop on nConvectiveFields
1148 break;
1149 }
1150 // 3D problems
1151 case 3:
1152 {
1153 ASSERTL0(false, "3D FRDG case not implemented yet");
1154 break;
1155 }
1156 }
1157}
1158
1159/**
1160 * @brief Builds the numerical flux for the 1st order derivatives
1161 *
1162 */
1163void DiffusionLFR::NumFluxforScalar(
1164 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1165 const Array<OneD, Array<OneD, NekDouble>> &ufield,
1166 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &uflux)
1167{
1168 int i, j;
1169 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1170 int nvariables = fields.size();
1171 int nDim = fields[0]->GetCoordim(0);
1172
1173 Array<OneD, NekDouble> Fwd(nTracePts);
1174 Array<OneD, NekDouble> Bwd(nTracePts);
1175 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
1176 Array<OneD, NekDouble> fluxtemp(nTracePts, 0.0);
1177
1178 // Get the normal velocity Vn
1179 for (i = 0; i < nDim; ++i)
1180 {
1181 Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1, Vn, 1, Vn, 1);
1182 }
1183
1184 // Get the sign of (v \cdot n), v = an arbitrary vector
1185 // Evaluate upwind flux:
1186 // uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
1187 for (j = 0; j < nDim; ++j)
1188 {
1189 for (i = 0; i < nvariables; ++i)
1190 {
1191 // Compute Fwd and Bwd value of ufield of i direction
1192 fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1193
1194 // if Vn >= 0, flux = uFwd, i.e.,
1195 // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
1196 // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
1197
1198 // else if Vn < 0, flux = uBwd, i.e.,
1199 // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
1200 // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
1201
1202 fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, fluxtemp);
1203
1204 // Imposing weak boundary condition with flux
1205 // if Vn >= 0, uflux = uBwd at Neumann, i.e.,
1206 // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
1207 // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
1208
1209 // if Vn >= 0, uflux = uFwd at Neumann, i.e.,
1210 // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
1211 // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
1212
1213 if (fields[0]->GetBndCondExpansions().size())
1214 {
1215 WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
1216 }
1217
1218 // if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
1219 // i.e,
1220 // edge::eForward, uFwd \‍(\tan_{\xi}^Fwd \cdot \vec{n})
1221 // edge::eBackward, uFwd \‍(\tan_{\xi}^Bwd \cdot \vec{n})
1222
1223 // else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}),
1224 // i.e,
1225 // edge::eForward, uBwd \‍(\tan_{\xi}^Fwd \cdot \vec{n})
1226 // edge::eBackward, uBwd \‍(\tan_{\xi}^Bwd \cdot \vec{n})
1227
1228 Vmath::Vcopy(nTracePts, fluxtemp, 1, uflux[i][j], 1);
1229 }
1230 }
1231}
1232
1233/**
1234 * @brief Imposes appropriate bcs for the 1st order derivatives
1235 *
1236 */
1237void DiffusionLFR::WeakPenaltyforScalar(
1238 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields, const int var,
1239 const Array<OneD, const NekDouble> &ufield,
1240 Array<OneD, NekDouble> &penaltyflux)
1241{
1242 // Variable initialisation
1243 int i, e, id1, id2;
1244 int nBndEdgePts, nBndEdges;
1245 int cnt = 0;
1246 int nBndRegions = fields[var]->GetBndCondExpansions().size();
1247 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1248
1249 // Extract physical values of the fields at the trace space
1250 Array<OneD, NekDouble> uplus(nTracePts);
1251 fields[var]->ExtractTracePhys(ufield, uplus);
1252
1253 // Impose boundary conditions
1254 for (i = 0; i < nBndRegions; ++i)
1255 {
1256 // Number of boundary edges related to bcRegion 'i'
1257 nBndEdges = fields[var]->GetBndCondExpansions()[i]->GetExpSize();
1258
1259 // Weakly impose boundary conditions by modifying flux values
1260 for (e = 0; e < nBndEdges; ++e)
1261 {
1262 // Number of points on the expansion
1263 nBndEdgePts = fields[var]
1264 ->GetBndCondExpansions()[i]
1265 ->GetExp(e)
1266 ->GetTotPoints();
1267
1268 // Offset of the boundary expansion
1269 id1 = fields[var]->GetBndCondExpansions()[i]->GetPhys_Offset(e);
1270
1271 // Offset of the trace space related to boundary expansion
1272 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1273 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1274
1275 // Dirichlet bcs ==> uflux = gD
1276 if (fields[var]
1277 ->GetBndConditions()[i]
1278 ->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1279 {
1280 Vmath::Vcopy(
1281 nBndEdgePts,
1282 &(fields[var]->GetBndCondExpansions()[i]->GetPhys())[id1],
1283 1, &penaltyflux[id2], 1);
1284 }
1285 // Neumann bcs ==> uflux = u+
1286 else if ((fields[var]->GetBndConditions()[i])
1287 ->GetBoundaryConditionType() ==
1288 SpatialDomains::eNeumann)
1289 {
1290 Vmath::Vcopy(nBndEdgePts, &uplus[id2], 1, &penaltyflux[id2], 1);
1291 }
1292 }
1293 }
1294}
1295
1296/**
1297 * @brief Build the numerical flux for the 2nd order derivatives
1298 *
1299 */
1300void DiffusionLFR::NumFluxforVector(
1301 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1302 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &ufield,
1303 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &qfield,
1304 Array<OneD, Array<OneD, NekDouble>> &qflux)
1305{
1306 int i, j;
1307 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1308 int nvariables = fields.size();
1309 int nDim = fields[0]->GetCoordim(0);
1310
1311 NekDouble C11 = 0.0;
1312 Array<OneD, NekDouble> Fwd(nTracePts);
1313 Array<OneD, NekDouble> Bwd(nTracePts);
1314 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
1315
1316 Array<OneD, NekDouble> qFwd(nTracePts);
1317 Array<OneD, NekDouble> qBwd(nTracePts);
1318 Array<OneD, NekDouble> qfluxtemp(nTracePts, 0.0);
1319 Array<OneD, NekDouble> uterm(nTracePts);
1320
1321 // Get the normal velocity Vn
1322 for (i = 0; i < nDim; ++i)
1323 {
1324 Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1, Vn, 1, Vn, 1);
1325 }
1326
1327 // Evaulate upwind flux:
1328 // qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
1329 for (i = 0; i < nvariables; ++i)
1330 {
1331 qflux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1332 for (j = 0; j < nDim; ++j)
1333 {
1334 // Compute Fwd and Bwd value of ufield of jth direction
1335 fields[i]->GetFwdBwdTracePhys(qfield[i][j], qFwd, qBwd);
1336
1337 // if Vn >= 0, flux = uFwd, i.e.,
1338 // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick
1339 // qflux = qBwd = q+
1340 // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick
1341 // qflux = qBwd = q-
1342
1343 // else if Vn < 0, flux = uBwd, i.e.,
1344 // edge::eForward, if V*n<0 <=> V*n_F<0, pick
1345 // qflux = qFwd = q-
1346 // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
1347 // qflux = qFwd = q+
1348
1349 fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1350
1351 Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
1352 qfluxtemp, 1);
1353
1354 // Imposing weak boundary condition with flux
1355 if (fields[0]->GetBndCondExpansions().size())
1356 {
1357 WeakPenaltyforVector(fields, i, j, qfield[i][j], qfluxtemp,
1358 C11);
1359 }
1360
1361 // q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
1362 // n_xi = n_x * tan_xi_x + n_y * tan_xi_y + n_z * tan_xi_z
1363 // n_xi = n_x * tan_eta_x + n_y * tan_eta_y + n_z*tan_eta_z
1364 Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1, qflux[i], 1);
1365 }
1366 }
1367}
1368
1369/**
1370 * @brief Imposes appropriate bcs for the 2nd order derivatives
1371 *
1372 */
1373void DiffusionLFR::WeakPenaltyforVector(
1374 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields, const int var,
1375 const int dir, const Array<OneD, const NekDouble> &qfield,
1376 Array<OneD, NekDouble> &penaltyflux, [[maybe_unused]] NekDouble C11)
1377{
1378 int i, e, id1, id2;
1379 int nBndEdges, nBndEdgePts;
1380 int nBndRegions = fields[var]->GetBndCondExpansions().size();
1381 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1382
1383 Array<OneD, NekDouble> uterm(nTracePts);
1384 Array<OneD, NekDouble> qtemp(nTracePts);
1385 int cnt = 0;
1386
1387 fields[var]->ExtractTracePhys(qfield, qtemp);
1388
1389 for (i = 0; i < nBndRegions; ++i)
1390 {
1391 nBndEdges = fields[var]->GetBndCondExpansions()[i]->GetExpSize();
1392
1393 // Weakly impose boundary conditions by modifying flux values
1394 for (e = 0; e < nBndEdges; ++e)
1395 {
1396 nBndEdgePts = fields[var]
1397 ->GetBndCondExpansions()[i]
1398 ->GetExp(e)
1399 ->GetTotPoints();
1400
1401 id1 = fields[var]->GetBndCondExpansions()[i]->GetPhys_Offset(e);
1402
1403 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1404 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1405
1406 // For Dirichlet boundary condition:
1407 // qflux = q+ - C_11 (u+ - g_D) (nx, ny)
1408 if (fields[var]
1409 ->GetBndConditions()[i]
1410 ->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1411 {
1412 Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1413 &qtemp[id2], 1, &penaltyflux[id2], 1);
1414 }
1415 // For Neumann boundary condition: qflux = g_N
1416 else if ((fields[var]->GetBndConditions()[i])
1417 ->GetBoundaryConditionType() ==
1418 SpatialDomains::eNeumann)
1419 {
1420 Vmath::Vmul(
1421 nBndEdgePts, &m_traceNormals[dir][id2], 1,
1422 &(fields[var]->GetBndCondExpansions()[i]->GetPhys())[id1],
1423 1, &penaltyflux[id2], 1);
1424 }
1425 }
1426 }
1427}
1428
1429/**
1430 * @brief Compute the derivative of the corrective flux for 1D problems.
1431 *
1432 * @param nConvectiveFields Number of fields.
1433 * @param fields Pointer to fields.
1434 * @param flux Volumetric flux in the physical space.
1435 * @param iFlux Numerical interface flux in physical space.
1436 * @param derCFlux Derivative of the corrective flux.
1437 *
1438 */
1439void DiffusionLFR::DerCFlux_1D(
1440 [[maybe_unused]] const int nConvectiveFields,
1441 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1442 const Array<OneD, const NekDouble> &flux,
1443 const Array<OneD, const NekDouble> &iFlux, Array<OneD, NekDouble> &derCFlux)
1444{
1445 int n;
1446 int nLocalSolutionPts, phys_offset, t_offset;
1447
1448 Array<OneD, NekDouble> auxArray1, auxArray2;
1449 Array<TwoD, const NekDouble> gmat;
1450 Array<OneD, const NekDouble> jac;
1451
1452 LibUtilities::PointsKeyVector ptsKeys;
1453 LibUtilities::BasisSharedPtr Basis;
1454 Basis = fields[0]->GetExp(0)->GetBasis(0);
1455
1456 int nElements = fields[0]->GetExpSize();
1457 int nSolutionPts = fields[0]->GetTotPoints();
1458
1459 vector<bool> leftAdjacentTraces =
1460 std::static_pointer_cast<MultiRegions::DisContField>(fields[0])
1461 ->GetLeftAdjacentTraces();
1462
1463 // Arrays to store the derivatives of the correction flux
1464 Array<OneD, NekDouble> DCL(nSolutionPts / nElements, 0.0);
1465 Array<OneD, NekDouble> DCR(nSolutionPts / nElements, 0.0);
1466
1467 // Arrays to store the intercell numerical flux jumps
1468 Array<OneD, NekDouble> JumpL(nElements);
1469 Array<OneD, NekDouble> JumpR(nElements);
1470
1471 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1472 fields[0]->GetTraceMap()->GetElmtToTrace();
1473
1474 for (n = 0; n < nElements; ++n)
1475 {
1476 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1477 phys_offset = fields[0]->GetPhys_Offset(n);
1478
1479 Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1480 Array<OneD, NekDouble> tmpFluxVertex(1, 0.0);
1481 Vmath::Vcopy(nLocalSolutionPts, &flux[phys_offset], 1, &tmparrayX1[0],
1482 1);
1483
1484 fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0], tmparrayX1,
1485 tmpFluxVertex);
1486
1487 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1488 elmtToTrace[n][0]->GetElmtId());
1489
1490 if (leftAdjacentTraces[2 * n])
1491 {
1492 JumpL[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1493 }
1494 else
1495 {
1496 JumpL[n] = iFlux[t_offset] - tmpFluxVertex[0];
1497 }
1498
1499 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1500 elmtToTrace[n][1]->GetElmtId());
1501
1502 fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1], tmparrayX1,
1503 tmpFluxVertex);
1504 if (leftAdjacentTraces[2 * n + 1])
1505 {
1506 JumpR[n] = iFlux[t_offset] - tmpFluxVertex[0];
1507 }
1508 else
1509 {
1510 JumpR[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1511 }
1512 }
1513
1514 for (n = 0; n < nElements; ++n)
1515 {
1516 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1517 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1518 phys_offset = fields[0]->GetPhys_Offset(n);
1519 jac = fields[0]
1520 ->GetExp(n)
1521 ->as<LocalRegions::Expansion1D>()
1522 ->GetGeom1D()
1523 ->GetMetricInfo()
1524 ->GetJac(ptsKeys);
1525
1526 JumpL[n] = JumpL[n] * jac[0];
1527 JumpR[n] = JumpR[n] * jac[0];
1528
1529 // Left jump multiplied by left derivative of C function
1530 Vmath::Smul(nLocalSolutionPts, JumpL[n], &m_dGL_xi1[n][0], 1, &DCL[0],
1531 1);
1532
1533 // Right jump multiplied by right derivative of C function
1534 Vmath::Smul(nLocalSolutionPts, JumpR[n], &m_dGR_xi1[n][0], 1, &DCR[0],
1535 1);
1536
1537 // Assembling derivative of the correction flux
1538 Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1539 &derCFlux[phys_offset], 1);
1540 }
1541}
1542
1543/**
1544 * @brief Compute the derivative of the corrective flux wrt a given
1545 * coordinate for 2D problems.
1546 *
1547 * There could be a bug for deformed elements since first derivatives
1548 * are not exactly the same results of DiffusionLDG as expected
1549 *
1550 * @param nConvectiveFields Number of fields.
1551 * @param fields Pointer to fields.
1552 * @param flux Volumetric flux in the physical space.
1553 * @param numericalFlux Numerical interface flux in physical space.
1554 * @param derCFlux Derivative of the corrective flux.
1555 *
1556 * \todo: Switch on shapes eventually here.
1557 */
1558void DiffusionLFR::DerCFlux_2D(
1559 [[maybe_unused]] const int nConvectiveFields, const int direction,
1560 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1561 const Array<OneD, const NekDouble> &flux,
1562 const Array<OneD, NekDouble> &iFlux, Array<OneD, NekDouble> &derCFlux)
1563{
1564 int n, e, i, j, cnt;
1565
1566 Array<OneD, const NekDouble> jac;
1567
1568 int nElements = fields[0]->GetExpSize();
1569 int trace_offset, phys_offset;
1570 int nLocalSolutionPts;
1571 int nquad0, nquad1;
1572 int nEdgePts;
1573
1574 LibUtilities::PointsKeyVector ptsKeys;
1575 Array<OneD, NekDouble> auxArray1, auxArray2;
1576 Array<OneD, LibUtilities::BasisSharedPtr> base;
1577 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1578 fields[0]->GetTraceMap()->GetElmtToTrace();
1579
1580 // Loop on the elements
1581 for (n = 0; n < nElements; ++n)
1582 {
1583 // Offset of the element on the global vector
1584 phys_offset = fields[0]->GetPhys_Offset(n);
1585 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1586 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1587
1588 jac = fields[0]
1589 ->GetExp(n)
1590 ->as<LocalRegions::Expansion2D>()
1591 ->GetGeom2D()
1592 ->GetMetricInfo()
1593 ->GetJac(ptsKeys);
1594
1595 base = fields[0]->GetExp(n)->GetBase();
1596 nquad0 = base[0]->GetNumPoints();
1597 nquad1 = base[1]->GetNumPoints();
1598
1599 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1600 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1601 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1602 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1603
1604 // Loop on the edges
1605 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1606 {
1607 // Number of edge points of edge 'e'
1608 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1609
1610 // Array for storing volumetric fluxes on each edge
1611 Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1612
1613 // Offset of the trace space correspondent to edge 'e'
1614 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1615 elmtToTrace[n][e]->GetElmtId());
1616
1617 // Extract the edge values of the volumetric fluxes
1618 // on edge 'e' and order them accordingly to the order
1619 // of the trace space
1620 fields[0]->GetExp(n)->GetTracePhysVals(
1621 e, elmtToTrace[n][e], flux + phys_offset, auxArray1 = tmparray);
1622
1623 // Compute the fluxJumps per each edge 'e' and each
1624 // flux point
1625 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1626 Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1, &tmparray[0], 1,
1627 &fluxJumps[0], 1);
1628
1629 // Check the ordering of the fluxJumps and reverse
1630 // it in case of backward definition of edge 'e'
1631 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1632 StdRegions::eBackwards)
1633 {
1634 Vmath::Reverse(nEdgePts, &fluxJumps[0], 1, &fluxJumps[0], 1);
1635 }
1636
1637 // Deformed elements
1638 if (fields[0]
1639 ->GetExp(n)
1640 ->as<LocalRegions::Expansion2D>()
1641 ->GetGeom2D()
1642 ->GetMetricInfo()
1643 ->GetGtype() == SpatialDomains::eDeformed)
1644 {
1645 // Extract the Jacobians along edge 'e'
1646 Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1647
1648 fields[0]->GetExp(n)->GetTracePhysVals(
1649 e, elmtToTrace[n][e], jac, auxArray1 = jacEdge);
1650
1651 // Check the ordering of the fluxJumps and reverse
1652 // it in case of backward definition of edge 'e'
1653 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1654 StdRegions::eBackwards)
1655 {
1656 Vmath::Reverse(nEdgePts, &jacEdge[0], 1, &jacEdge[0], 1);
1657 }
1658
1659 // Multiply the fluxJumps by the edge 'e' Jacobians
1660 // to bring the fluxJumps into the standard space
1661 for (j = 0; j < nEdgePts; j++)
1662 {
1663 fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1664 }
1665 }
1666 // Non-deformed elements
1667 else
1668 {
1669 // Multiply the fluxJumps by the edge 'e' Jacobians
1670 // to bring the fluxJumps into the standard space
1671 Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1, fluxJumps, 1);
1672 }
1673
1674 // Multiply jumps by derivatives of the correction functions
1675 // All the quntities at this level should be defined into
1676 // the standard space
1677 switch (e)
1678 {
1679 case 0:
1680 for (i = 0; i < nquad0; ++i)
1681 {
1682 for (j = 0; j < nquad1; ++j)
1683 {
1684 cnt = i + j * nquad0;
1685 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1686 }
1687 }
1688 break;
1689 case 1:
1690 for (i = 0; i < nquad1; ++i)
1691 {
1692 for (j = 0; j < nquad0; ++j)
1693 {
1694 cnt = (nquad0)*i + j;
1695 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1696 }
1697 }
1698 break;
1699 case 2:
1700 for (i = 0; i < nquad0; ++i)
1701 {
1702 for (j = 0; j < nquad1; ++j)
1703 {
1704 cnt = j * nquad0 + i;
1705 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1706 }
1707 }
1708 break;
1709 case 3:
1710 for (i = 0; i < nquad1; ++i)
1711 {
1712 for (j = 0; j < nquad0; ++j)
1713 {
1714 cnt = j + i * nquad0;
1715 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1716 }
1717 }
1718 break;
1719
1720 default:
1721 ASSERTL0(false, "edge value (< 3) is out of range");
1722 break;
1723 }
1724 }
1725
1726 // Summing the component of the flux for calculating the
1727 // derivatives per each direction
1728 if (direction == 0)
1729 {
1730 Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1, &divCFluxE3[0], 1,
1731 &derCFlux[phys_offset], 1);
1732 }
1733 else if (direction == 1)
1734 {
1735 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE2[0], 1,
1736 &derCFlux[phys_offset], 1);
1737 }
1738 }
1739}
1740
1741/**
1742 * @brief Compute the divergence of the corrective flux for 2D problems.
1743 *
1744 * @param nConvectiveFields Number of fields.
1745 * @param fields Pointer to fields.
1746 * @param fluxX1 X1 - volumetric flux in physical space.
1747 * @param fluxX2 X2 - volumetric flux in physical space.
1748 * @param numericalFlux Interface flux in physical space.
1749 * @param divCFlux Divergence of the corrective flux for
1750 * 2D problems.
1751 *
1752 * \todo: Switch on shapes eventually here.
1753 */
1754void DiffusionLFR::DivCFlux_2D(
1755 [[maybe_unused]] const int nConvectiveFields,
1756 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1757 const Array<OneD, const NekDouble> &fluxX1,
1758 const Array<OneD, const NekDouble> &fluxX2,
1759 const Array<OneD, const NekDouble> &numericalFlux,
1760 Array<OneD, NekDouble> &divCFlux)
1761{
1762 int n, e, i, j, cnt;
1763
1764 int nElements = fields[0]->GetExpSize();
1765
1766 int nLocalSolutionPts;
1767 int nEdgePts;
1768 int trace_offset;
1769 int phys_offset;
1770 int nquad0;
1771 int nquad1;
1772
1773 Array<OneD, NekDouble> auxArray1, auxArray2;
1774 Array<OneD, LibUtilities::BasisSharedPtr> base;
1775
1776 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1777 fields[0]->GetTraceMap()->GetElmtToTrace();
1778
1779 // Loop on the elements
1780 for (n = 0; n < nElements; ++n)
1781 {
1782 // Offset of the element on the global vector
1783 phys_offset = fields[0]->GetPhys_Offset(n);
1784 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1785
1786 base = fields[0]->GetExp(n)->GetBase();
1787 nquad0 = base[0]->GetNumPoints();
1788 nquad1 = base[1]->GetNumPoints();
1789
1790 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1791 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1792 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1793 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1794
1795 // Loop on the edges
1796 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1797 {
1798 // Number of edge points of edge e
1799 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1800
1801 Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1802 Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1803 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1804 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1805 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1806
1807 // Offset of the trace space correspondent to edge e
1808 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1809 elmtToTrace[n][e]->GetElmtId());
1810
1811 // Get the normals of edge e
1812 const Array<OneD, const Array<OneD, NekDouble>> &normals =
1813 fields[0]->GetExp(n)->GetTraceNormal(e);
1814
1815 // Extract the edge values of flux-x on edge e and order
1816 // them accordingly to the order of the trace space
1817 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1818 fluxX1 + phys_offset,
1819 auxArray1 = tmparrayX1);
1820
1821 // Extract the edge values of flux-y on edge e and order
1822 // them accordingly to the order of the trace space
1823 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1824 fluxX2 + phys_offset,
1825 auxArray1 = tmparrayX2);
1826
1827 // Multiply the edge components of the flux by the normal
1828 for (i = 0; i < nEdgePts; ++i)
1829 {
1830 fluxN[i] = tmparrayX1[i] * m_traceNormals[0][trace_offset + i] +
1831 tmparrayX2[i] * m_traceNormals[1][trace_offset + i];
1832 }
1833
1834 // Subtract to the Riemann flux the discontinuous flux
1835 Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1, &fluxN[0], 1,
1836 &fluxJumps[0], 1);
1837
1838 // Check the ordering of the jump vectors
1839 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1840 StdRegions::eBackwards)
1841 {
1842 Vmath::Reverse(nEdgePts, auxArray1 = fluxJumps, 1,
1843 auxArray2 = fluxJumps, 1);
1844 }
1845
1846 for (i = 0; i < nEdgePts; ++i)
1847 {
1848 if (m_traceNormals[0][trace_offset + i] != normals[0][i] ||
1849 m_traceNormals[1][trace_offset + i] != normals[1][i])
1850 {
1851 fluxJumps[i] = -fluxJumps[i];
1852 }
1853 }
1854
1855 // Multiply jumps by derivatives of the correction functions
1856 switch (e)
1857 {
1858 case 0:
1859 for (i = 0; i < nquad0; ++i)
1860 {
1861 // Multiply fluxJumps by Q factors
1862 fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1863
1864 for (j = 0; j < nquad1; ++j)
1865 {
1866 cnt = i + j * nquad0;
1867 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1868 }
1869 }
1870 break;
1871 case 1:
1872 for (i = 0; i < nquad1; ++i)
1873 {
1874 // Multiply fluxJumps by Q factors
1875 fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1876
1877 for (j = 0; j < nquad0; ++j)
1878 {
1879 cnt = (nquad0)*i + j;
1880 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1881 }
1882 }
1883 break;
1884 case 2:
1885 for (i = 0; i < nquad0; ++i)
1886 {
1887 // Multiply fluxJumps by Q factors
1888 fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1889
1890 for (j = 0; j < nquad1; ++j)
1891 {
1892 cnt = j * nquad0 + i;
1893 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1894 }
1895 }
1896 break;
1897 case 3:
1898 for (i = 0; i < nquad1; ++i)
1899 {
1900 // Multiply fluxJumps by Q factors
1901 fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1902 for (j = 0; j < nquad0; ++j)
1903 {
1904 cnt = j + i * nquad0;
1905 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1906 }
1907 }
1908 break;
1909
1910 default:
1911 ASSERTL0(false, "edge value (< 3) is out of range");
1912 break;
1913 }
1914 }
1915
1916 // Sum each edge contribution
1917 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
1918 &divCFlux[phys_offset], 1);
1919
1920 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1921 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1922
1923 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1924 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1925 }
1926}
1927
1928/**
1929 * @brief Compute the divergence of the corrective flux for 2D problems
1930 * where POINTSTYPE="GaussGaussLegendre"
1931 *
1932 * @param nConvectiveFields Number of fields.
1933 * @param fields Pointer to fields.
1934 * @param fluxX1 X1-volumetric flux in physical space.
1935 * @param fluxX2 X2-volumetric flux in physical space.
1936 * @param numericalFlux Interface flux in physical space.
1937 * @param divCFlux Divergence of the corrective flux.
1938 *
1939 * \todo: Switch on shapes eventually here.
1940 */
1941
1942void DiffusionLFR::DivCFlux_2D_Gauss(
1943 [[maybe_unused]] const int nConvectiveFields,
1944 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1945 const Array<OneD, const NekDouble> &fluxX1,
1946 const Array<OneD, const NekDouble> &fluxX2,
1947 const Array<OneD, const NekDouble> &numericalFlux,
1948 Array<OneD, NekDouble> &divCFlux)
1949{
1950 int n, e, i, j, cnt;
1951
1952 int nElements = fields[0]->GetExpSize();
1953 int nLocalSolutionPts;
1954 int nEdgePts;
1955 int trace_offset;
1956 int phys_offset;
1957 int nquad0;
1958 int nquad1;
1959
1960 Array<OneD, NekDouble> auxArray1, auxArray2;
1961 Array<OneD, LibUtilities::BasisSharedPtr> base;
1962
1963 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1964 fields[0]->GetTraceMap()->GetElmtToTrace();
1965
1966 // Loop on the elements
1967 for (n = 0; n < nElements; ++n)
1968 {
1969 // Offset of the element on the global vector
1970 phys_offset = fields[0]->GetPhys_Offset(n);
1971 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1972
1973 base = fields[0]->GetExp(n)->GetBase();
1974 nquad0 = base[0]->GetNumPoints();
1975 nquad1 = base[1]->GetNumPoints();
1976
1977 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1978 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1979 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1980 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1981
1982 // Loop on the edges
1983 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1984 {
1985 // Number of edge points of edge e
1986 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1987
1988 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1989 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1990 Array<OneD, NekDouble> fluxN_R(nEdgePts, 0.0);
1991 Array<OneD, NekDouble> fluxN_D(nEdgePts, 0.0);
1992 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1993
1994 // Offset of the trace space correspondent to edge e
1995 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1996 elmtToTrace[n][e]->GetElmtId());
1997
1998 // Get the normals of edge e
1999 const Array<OneD, const Array<OneD, NekDouble>> &normals =
2000 fields[0]->GetExp(n)->GetTraceNormal(e);
2001
2002 // Extract the trasformed normal flux at each edge
2003 switch (e)
2004 {
2005 case 0:
2006 // Extract the edge values of transformed flux-y on
2007 // edge e and order them accordingly to the order of
2008 // the trace space
2009 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2010 fluxX2 + phys_offset,
2011 auxArray1 = fluxN_D);
2012
2013 Vmath::Neg(nEdgePts, fluxN_D, 1);
2014
2015 // Extract the physical Rieamann flux at each edge
2016 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2017 &fluxN[0], 1);
2018
2019 // Check the ordering of vectors
2020 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2021 StdRegions::eBackwards)
2022 {
2023 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2024 auxArray2 = fluxN, 1);
2025
2026 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2027 auxArray2 = fluxN_D, 1);
2028 }
2029
2030 // Transform Riemann Fluxes in the standard element
2031 for (i = 0; i < nquad0; ++i)
2032 {
2033 // Multiply Riemann Flux by Q factors
2034 fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
2035 }
2036
2037 for (i = 0; i < nEdgePts; ++i)
2038 {
2039 if (m_traceNormals[0][trace_offset + i] !=
2040 normals[0][i] ||
2041 m_traceNormals[1][trace_offset + i] !=
2042 normals[1][i])
2043 {
2044 fluxN_R[i] = -fluxN_R[i];
2045 }
2046 }
2047
2048 // Subtract to the Riemann flux the discontinuous
2049 // flux
2050 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2051 &fluxJumps[0], 1);
2052
2053 // Multiplicate the flux jumps for the correction
2054 // function
2055 for (i = 0; i < nquad0; ++i)
2056 {
2057 for (j = 0; j < nquad1; ++j)
2058 {
2059 cnt = i + j * nquad0;
2060 divCFluxE0[cnt] = -fluxJumps[i] * m_dGL_xi2[n][j];
2061 }
2062 }
2063 break;
2064 case 1:
2065 // Extract the edge values of transformed flux-x on
2066 // edge e and order them accordingly to the order of
2067 // the trace space
2068 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2069 fluxX1 + phys_offset,
2070 auxArray1 = fluxN_D);
2071
2072 // Extract the physical Rieamann flux at each edge
2073 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2074 &fluxN[0], 1);
2075
2076 // Check the ordering of vectors
2077 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2078 StdRegions::eBackwards)
2079 {
2080 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2081 auxArray2 = fluxN, 1);
2082
2083 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2084 auxArray2 = fluxN_D, 1);
2085 }
2086
2087 // Transform Riemann Fluxes in the standard element
2088 for (i = 0; i < nquad1; ++i)
2089 {
2090 // Multiply Riemann Flux by Q factors
2091 fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
2092 }
2093
2094 for (i = 0; i < nEdgePts; ++i)
2095 {
2096 if (m_traceNormals[0][trace_offset + i] !=
2097 normals[0][i] ||
2098 m_traceNormals[1][trace_offset + i] !=
2099 normals[1][i])
2100 {
2101 fluxN_R[i] = -fluxN_R[i];
2102 }
2103 }
2104
2105 // Subtract to the Riemann flux the discontinuous
2106 // flux
2107 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2108 &fluxJumps[0], 1);
2109
2110 // Multiplicate the flux jumps for the correction
2111 // function
2112 for (i = 0; i < nquad1; ++i)
2113 {
2114 for (j = 0; j < nquad0; ++j)
2115 {
2116 cnt = (nquad0)*i + j;
2117 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
2118 }
2119 }
2120 break;
2121 case 2:
2122
2123 // Extract the edge values of transformed flux-y on
2124 // edge e and order them accordingly to the order of
2125 // the trace space
2126
2127 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2128 fluxX2 + phys_offset,
2129 auxArray1 = fluxN_D);
2130
2131 // Extract the physical Rieamann flux at each edge
2132 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2133 &fluxN[0], 1);
2134
2135 // Check the ordering of vectors
2136 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2137 StdRegions::eBackwards)
2138 {
2139 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2140 auxArray2 = fluxN, 1);
2141
2142 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2143 auxArray2 = fluxN_D, 1);
2144 }
2145
2146 // Transform Riemann Fluxes in the standard element
2147 for (i = 0; i < nquad0; ++i)
2148 {
2149 // Multiply Riemann Flux by Q factors
2150 fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
2151 }
2152
2153 for (i = 0; i < nEdgePts; ++i)
2154 {
2155 if (m_traceNormals[0][trace_offset + i] !=
2156 normals[0][i] ||
2157 m_traceNormals[1][trace_offset + i] !=
2158 normals[1][i])
2159 {
2160 fluxN_R[i] = -fluxN_R[i];
2161 }
2162 }
2163
2164 // Subtract to the Riemann flux the discontinuous
2165 // flux
2166
2167 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2168 &fluxJumps[0], 1);
2169
2170 // Multiplicate the flux jumps for the correction
2171 // function
2172 for (i = 0; i < nquad0; ++i)
2173 {
2174 for (j = 0; j < nquad1; ++j)
2175 {
2176 cnt = j * nquad0 + i;
2177 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
2178 }
2179 }
2180 break;
2181 case 3:
2182 // Extract the edge values of transformed flux-x on
2183 // edge e and order them accordingly to the order of
2184 // the trace space
2185
2186 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2187 fluxX1 + phys_offset,
2188 auxArray1 = fluxN_D);
2189 Vmath::Neg(nEdgePts, fluxN_D, 1);
2190
2191 // Extract the physical Rieamann flux at each edge
2192 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2193 &fluxN[0], 1);
2194
2195 // Check the ordering of vectors
2196 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2197 StdRegions::eBackwards)
2198 {
2199 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2200 auxArray2 = fluxN, 1);
2201
2202 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2203 auxArray2 = fluxN_D, 1);
2204 }
2205
2206 // Transform Riemann Fluxes in the standard element
2207 for (i = 0; i < nquad1; ++i)
2208 {
2209 // Multiply Riemann Flux by Q factors
2210 fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
2211 }
2212
2213 for (i = 0; i < nEdgePts; ++i)
2214 {
2215 if (m_traceNormals[0][trace_offset + i] !=
2216 normals[0][i] ||
2217 m_traceNormals[1][trace_offset + i] !=
2218 normals[1][i])
2219 {
2220 fluxN_R[i] = -fluxN_R[i];
2221 }
2222 }
2223
2224 // Subtract to the Riemann flux the discontinuous
2225 // flux
2226
2227 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2228 &fluxJumps[0], 1);
2229
2230 // Multiplicate the flux jumps for the correction
2231 // function
2232 for (i = 0; i < nquad1; ++i)
2233 {
2234 for (j = 0; j < nquad0; ++j)
2235 {
2236 cnt = j + i * nquad0;
2237 divCFluxE3[cnt] = -fluxJumps[i] * m_dGL_xi1[n][j];
2238 }
2239 }
2240 break;
2241 default:
2242 ASSERTL0(false, "edge value (< 3) is out of range");
2243 break;
2244 }
2245 }
2246
2247 // Sum each edge contribution
2248 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
2249 &divCFlux[phys_offset], 1);
2250
2251 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2252 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2253
2254 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2255 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2256 }
2257}
2258
2259} // namespace Nektar::SolverUtils
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
WARNINGL0
#define WARNINGL0(condition, msg)
Definition: ErrorUtil.hpp:215
Nektar::LibUtilities::Basis
Represents a basis of a given type.
Definition: Basis.h:198
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: BasicUtils/NekFactory.hpp:197
Nektar::LocalRegions::Expansion::GetMetricInfo
const SpatialDomains::GeomFactorsSharedPtr & GetMetricInfo() const
Definition: Expansion.cpp:246
Nektar::SolverUtils::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:1439
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:1558
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:1237
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:1373
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:1754
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:1300
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:1163
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:1942
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::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Polylib::jacobd
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition: Polylib.cpp:1378
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
Vmath::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 Jun 15 2025 07:52:32 for Nektar++ by doxygen 1.9.4