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