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