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