Nektar++
Loading...
Searching...
No Matches
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 auxArray1 = u1_hat + phys_offset,
1069 auxArray2 = m_tmp1[i][j] + phys_offset);
1070
1071 fields[i]->GetExp(n)->StdPhysDeriv(
1072 auxArray1 = u2_hat + phys_offset,
1073 NullNekDouble1DArray,
1074 auxArray2 = m_tmp2[i][j] + phys_offset);
1075 }
1076
1077 Vmath::Vadd(nSolutionPts, &m_tmp1[i][j][0], 1,
1078 &m_tmp2[i][j][0], 1, &m_DU1[i][j][0], 1);
1079
1080 // Divide by the metric jacobian
1081 Vmath::Vdiv(nSolutionPts, &m_DU1[i][j][0], 1, &m_jac[0], 1,
1082 &m_DU1[i][j][0], 1);
1083
1084 // Computing the standard first-order correction
1085 // derivatives
1086 DerCFlux_2D(nConvectiveFields, j, fields, inarray[i],
1087 m_IF1[i][j], m_DFC1[i][j]);
1088 }
1089
1090 // Multiplying derivatives by B matrix to get auxiliary
1091 // variables
1092 for (j = 0; j < nSolutionPts; ++j)
1093 {
1094 // std(q1)
1095 m_BD1[i][0][j] = (m_gmat[0][j] * m_gmat[0][j] +
1096 m_gmat[2][j] * m_gmat[2][j]) *
1097 m_DFC1[i][0][j] +
1098 (m_gmat[1][j] * m_gmat[0][j] +
1099 m_gmat[3][j] * m_gmat[2][j]) *
1100 m_DFC1[i][1][j];
1101
1102 // std(q2)
1103 m_BD1[i][1][j] = (m_gmat[1][j] * m_gmat[0][j] +
1104 m_gmat[3][j] * m_gmat[2][j]) *
1105 m_DFC1[i][0][j] +
1106 (m_gmat[1][j] * m_gmat[1][j] +
1107 m_gmat[3][j] * m_gmat[3][j]) *
1108 m_DFC1[i][1][j];
1109 }
1110
1111 // Multiplying derivatives by A^(-1) to get back
1112 // into the physical space
1113 for (j = 0; j < nSolutionPts; j++)
1114 {
1115 // q1 = A11^(-1)*std(q1) + A12^(-1)*std(q2)
1116 m_DFC1[i][0][j] = m_gmat[3][j] * m_BD1[i][0][j] -
1117 m_gmat[2][j] * m_BD1[i][1][j];
1118
1119 // q2 = A21^(-1)*std(q1) + A22^(-1)*std(q2)
1120 m_DFC1[i][1][j] = m_gmat[0][j] * m_BD1[i][1][j] -
1121 m_gmat[1][j] * m_BD1[i][0][j];
1122 }
1123
1124 // Computing the physical first-order derivatives
1125 for (j = 0; j < nDim; ++j)
1126 {
1127 Vmath::Vadd(nSolutionPts, &m_DU1[i][j][0], 1,
1128 &m_DFC1[i][j][0], 1, &m_D1[j][i][0], 1);
1129 }
1130 } // Close loop on nScalars
1131
1132 // For 3D Homogeneous 1D only take derivatives in 3rd direction
1133 if (m_diffDim == 1)
1134 {
1135 for (i = 0; i < nScalars; ++i)
1136 {
1137 m_D1[2][i] = m_homoDerivs[i];
1138 }
1139 }
1140
1141 // Computing the viscous tensor
1142 m_fluxVectorNS(inarray, m_D1, m_viscTensor);
1143
1144 // Compute u from q_{\eta} and q_{\xi}
1145 // Obtain numerical fluxes
1146 NumericalFluxO2(fields, inarray, m_viscTensor, m_viscFlux);
1147
1148 // Computing the standard second-order discontinuous
1149 // derivatives
1150 for (i = 0; i < nConvectiveFields; ++i)
1151 {
1152 // Temporary vectors
1153 Array<OneD, NekDouble> f_hat(nSolutionPts, 0.0);
1154 Array<OneD, NekDouble> g_hat(nSolutionPts, 0.0);
1155
1156 for (j = 0; j < nSolutionPts; j++)
1157 {
1158 f_hat[j] = (m_viscTensor[0][i][j] * m_gmat[0][j] +
1159 m_viscTensor[1][i][j] * m_gmat[2][j]) *
1160 m_jac[j];
1161
1162 g_hat[j] = (m_viscTensor[0][i][j] * m_gmat[1][j] +
1163 m_viscTensor[1][i][j] * m_gmat[3][j]) *
1164 m_jac[j];
1165 }
1166
1167 for (n = 0; n < nElements; n++)
1168 {
1169 phys_offset = fields[0]->GetPhys_Offset(n);
1170
1171 fields[0]->GetExp(n)->StdPhysDeriv(
1172 auxArray1 = f_hat + phys_offset,
1173 auxArray2 = m_DD1[i][0] + phys_offset);
1174
1175 fields[0]->GetExp(n)->StdPhysDeriv(
1176 auxArray1 = g_hat + phys_offset, NullNekDouble1DArray,
1177 auxArray2 = m_DD1[i][1] + phys_offset);
1178 }
1179
1180 // Divergence of the standard discontinuous flux
1181 Vmath::Vadd(nSolutionPts, &m_DD1[i][0][0], 1, &m_DD1[i][1][0],
1182 1, &m_divFD[i][0], 1);
1183
1184 // Divergence of the standard correction flux
1185 if (Basis[0]->GetPointsType() ==
1186 LibUtilities::eGaussGaussLegendre &&
1187 Basis[1]->GetPointsType() ==
1188 LibUtilities::eGaussGaussLegendre)
1189 {
1190 DivCFlux_2D_Gauss(nConvectiveFields, fields, f_hat, g_hat,
1191 m_viscFlux[i], m_divFC[i]);
1192 }
1193 else
1194 {
1195 DivCFlux_2D(nConvectiveFields, fields, m_viscTensor[0][i],
1196 m_viscTensor[1][i], m_viscFlux[i], m_divFC[i]);
1197 }
1198
1199 // Divergence of the standard final flux
1200 Vmath::Vadd(nSolutionPts, &m_divFD[i][0], 1, &m_divFC[i][0], 1,
1201 &outarray[i][0], 1);
1202
1203 // Dividing by the jacobian to get back into
1204 // physical space
1205 Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1, &m_jac[0], 1,
1206 &outarray[i][0], 1);
1207
1208 // Primitive Dealiasing 2D
1209 if (!(Basis[0]->Collocation()))
1210 {
1211 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1212 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1213 }
1214 }
1215 break;
1216 }
1217 // 3D problems
1218 case 3:
1219 {
1220 ASSERTL0(false, "3D FRDG case not implemented yet");
1221 break;
1222 }
1223 }
1224}
1225
1226/**
1227 * @brief Builds the numerical flux for the 1st order derivatives
1228 *
1229 */
1230void DiffusionLFRNS::NumericalFluxO1(
1231 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1232 const Array<OneD, Array<OneD, NekDouble>> &inarray,
1233 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &numericalFluxO1)
1234{
1235 int i, j;
1236 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1237 int nScalars = inarray.size();
1238 int nDim = fields[0]->GetCoordim(0);
1239
1240 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
1241 Array<OneD, NekDouble> fluxtemp(nTracePts, 0.0);
1242
1243 // Get the normal velocity Vn
1244 for (i = 0; i < nDim; ++i)
1245 {
1246 fields[0]->ExtractTracePhys(inarray[i], m_traceVel[i]);
1247 Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, m_traceVel[i], 1, Vn, 1,
1248 Vn, 1);
1249 }
1250
1251 // Store forwards/backwards space along trace space
1252 Array<OneD, Array<OneD, NekDouble>> Fwd(nScalars);
1253 Array<OneD, Array<OneD, NekDouble>> Bwd(nScalars);
1254 Array<OneD, Array<OneD, NekDouble>> numflux(nScalars);
1255
1256 for (i = 0; i < nScalars; ++i)
1257 {
1258 Fwd[i] = Array<OneD, NekDouble>(nTracePts);
1259 Bwd[i] = Array<OneD, NekDouble>(nTracePts);
1260 numflux[i] = Array<OneD, NekDouble>(nTracePts);
1261 fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
1262 fields[0]->GetTrace()->Upwind(Vn, Fwd[i], Bwd[i], numflux[i]);
1263 }
1264
1265 // Modify the values in case of boundary interfaces
1266 if (fields[0]->GetBndCondExpansions().size())
1267 {
1268 WeakPenaltyO1(fields, inarray, numflux);
1269 }
1270
1271 // Splitting the numerical flux into the dimensions
1272 for (j = 0; j < nDim; ++j)
1273 {
1274 for (i = 0; i < nScalars; ++i)
1275 {
1276 Vmath::Vcopy(nTracePts, numflux[i], 1, numericalFluxO1[i][j], 1);
1277 }
1278 }
1279}
1280
1281/**
1282 * @brief Imposes appropriate bcs for the 1st order derivatives
1283 *
1284 */
1285void DiffusionLFRNS::WeakPenaltyO1(
1286 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1287 const Array<OneD, Array<OneD, NekDouble>> &inarray,
1288 Array<OneD, Array<OneD, NekDouble>> &penaltyfluxO1)
1289{
1290 int cnt;
1291 int i, j, e;
1292 int id1, id2;
1293
1294 int nBndEdgePts, nBndEdges, nBndRegions;
1295
1296 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1297 int nScalars = inarray.size();
1298
1299 Array<OneD, NekDouble> tmp1(nTracePts, 0.0);
1300 Array<OneD, NekDouble> tmp2(nTracePts, 0.0);
1301 Array<OneD, NekDouble> Tw(nTracePts, m_Twall);
1302
1303 Array<OneD, Array<OneD, NekDouble>> scalarVariables(nScalars);
1304 Array<OneD, Array<OneD, NekDouble>> uplus(nScalars);
1305
1306 // Extract internal values of the scalar variables for Neumann bcs
1307 for (i = 0; i < nScalars; ++i)
1308 {
1309 scalarVariables[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1310
1311 uplus[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1312 fields[i]->ExtractTracePhys(inarray[i], uplus[i]);
1313 }
1314
1315 // Compute boundary conditions for velocities
1316 for (i = 0; i < nScalars - 1; ++i)
1317 {
1318 // Note that cnt has to loop on nBndRegions and nBndEdges
1319 // and has to be reset to zero per each equation
1320 cnt = 0;
1321 nBndRegions = fields[i + 1]->GetBndCondExpansions().size();
1322 for (j = 0; j < nBndRegions; ++j)
1323 {
1324 if (fields[i + 1]
1325 ->GetBndConditions()[j]
1326 ->GetBoundaryConditionType() == SpatialDomains::ePeriodic)
1327 {
1328 continue;
1329 }
1330
1331 nBndEdges = fields[i + 1]->GetBndCondExpansions()[j]->GetExpSize();
1332 for (e = 0; e < nBndEdges; ++e)
1333 {
1334 nBndEdgePts = fields[i + 1]
1335 ->GetBndCondExpansions()[j]
1336 ->GetExp(e)
1337 ->GetTotPoints();
1338
1339 id1 =
1340 fields[i + 1]->GetBndCondExpansions()[j]->GetPhys_Offset(e);
1341
1342 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1343 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
1344 cnt++));
1345
1346 // Reinforcing bcs for velocity in case of Wall bcs
1347 if (boost::iequals(
1348 fields[i]->GetBndConditions()[j]->GetUserDefined(),
1349 "WallViscous") ||
1350 boost::iequals(
1351 fields[i]->GetBndConditions()[j]->GetUserDefined(),
1352 "WallAdiabatic"))
1353 {
1354 Vmath::Zero(nBndEdgePts, &scalarVariables[i][id2], 1);
1355 }
1356
1357 // Imposing velocity bcs if not Wall
1358 else if (fields[i]
1359 ->GetBndConditions()[j]
1360 ->GetBoundaryConditionType() ==
1361 SpatialDomains::eDirichlet)
1362 {
1363 Vmath::Vdiv(nBndEdgePts,
1364 &(fields[i + 1]
1365 ->GetBndCondExpansions()[j]
1366 ->UpdatePhys())[id1],
1367 1,
1368 &(fields[0]
1369 ->GetBndCondExpansions()[j]
1370 ->UpdatePhys())[id1],
1371 1, &scalarVariables[i][id2], 1);
1372 }
1373
1374 // For Dirichlet boundary condition: uflux = u_bcs
1375 if (fields[i]
1376 ->GetBndConditions()[j]
1377 ->GetBoundaryConditionType() ==
1378 SpatialDomains::eDirichlet)
1379 {
1380 Vmath::Vcopy(nBndEdgePts, &scalarVariables[i][id2], 1,
1381 &penaltyfluxO1[i][id2], 1);
1382 }
1383
1384 // For Neumann boundary condition: uflux = u_+
1385 else if ((fields[i]->GetBndConditions()[j])
1386 ->GetBoundaryConditionType() ==
1387 SpatialDomains::eNeumann)
1388 {
1389 Vmath::Vcopy(nBndEdgePts, &uplus[i][id2], 1,
1390 &penaltyfluxO1[i][id2], 1);
1391 }
1392
1393 // Building kinetic energy to be used for T bcs
1394 Vmath::Vmul(nBndEdgePts, &scalarVariables[i][id2], 1,
1395 &scalarVariables[i][id2], 1, &tmp1[id2], 1);
1396
1397 Vmath::Smul(nBndEdgePts, 0.5, &tmp1[id2], 1, &tmp1[id2], 1);
1398
1399 Vmath::Vadd(nBndEdgePts, &tmp2[id2], 1, &tmp1[id2], 1,
1400 &tmp2[id2], 1);
1401 }
1402 }
1403 }
1404
1405 // Compute boundary conditions for temperature
1406 cnt = 0;
1407 nBndRegions = fields[nScalars]->GetBndCondExpansions().size();
1408 for (j = 0; j < nBndRegions; ++j)
1409 {
1410 nBndEdges = fields[nScalars]->GetBndCondExpansions()[j]->GetExpSize();
1411
1412 if (fields[nScalars]
1413 ->GetBndConditions()[j]
1414 ->GetBoundaryConditionType() == SpatialDomains::ePeriodic)
1415 {
1416 continue;
1417 }
1418
1419 for (e = 0; e < nBndEdges; ++e)
1420 {
1421 nBndEdgePts = fields[nScalars]
1422 ->GetBndCondExpansions()[j]
1423 ->GetExp(e)
1424 ->GetTotPoints();
1425
1426 id1 =
1427 fields[nScalars]->GetBndCondExpansions()[j]->GetPhys_Offset(e);
1428
1429 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1430 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1431
1432 // Imposing Temperature Twall at the wall
1433 if (boost::iequals(
1434 fields[i]->GetBndConditions()[j]->GetUserDefined(),
1435 "WallViscous"))
1436 {
1437 Vmath::Vcopy(nBndEdgePts, &Tw[0], 1,
1438 &scalarVariables[nScalars - 1][id2], 1);
1439 }
1440 // Imposing Temperature through condition on the Energy
1441 // for no wall boundaries (e.g. farfield)
1442 else if (fields[i]
1443 ->GetBndConditions()[j]
1444 ->GetBoundaryConditionType() ==
1445 SpatialDomains::eDirichlet)
1446 {
1447 // Divide E by rho
1448 Vmath::Vdiv(
1449 nBndEdgePts,
1450 &(fields[nScalars]
1451 ->GetBndCondExpansions()[j]
1452 ->GetPhys())[id1],
1453 1, &(fields[0]->GetBndCondExpansions()[j]->GetPhys())[id1],
1454 1, &scalarVariables[nScalars - 1][id2], 1);
1455
1456 // Subtract kinetic energy to E/rho
1457 Vmath::Vsub(nBndEdgePts, &scalarVariables[nScalars - 1][id2], 1,
1458 &tmp2[id2], 1, &scalarVariables[nScalars - 1][id2],
1459 1);
1460
1461 // Multiply by constant factor (gamma-1)/R
1462 Vmath::Smul(nBndEdgePts, (m_gamma - 1) / m_gasConstant,
1463 &scalarVariables[nScalars - 1][id2], 1,
1464 &scalarVariables[nScalars - 1][id2], 1);
1465 }
1466
1467 // For Dirichlet boundary condition: uflux = u_bcs
1468 if (fields[nScalars]
1469 ->GetBndConditions()[j]
1470 ->GetBoundaryConditionType() ==
1471 SpatialDomains::eDirichlet &&
1472 !boost::iequals(
1473 fields[nScalars]->GetBndConditions()[j]->GetUserDefined(),
1474 "WallAdiabatic"))
1475 {
1476 Vmath::Vcopy(nBndEdgePts, &scalarVariables[nScalars - 1][id2],
1477 1, &penaltyfluxO1[nScalars - 1][id2], 1);
1478 }
1479
1480 // For Neumann boundary condition: uflux = u_+
1481 else if (((fields[nScalars]->GetBndConditions()[j])
1482 ->GetBoundaryConditionType() ==
1483 SpatialDomains::eNeumann) ||
1484 boost::iequals(fields[nScalars]
1485 ->GetBndConditions()[j]
1486 ->GetUserDefined(),
1487 "WallAdiabatic"))
1488 {
1489 Vmath::Vcopy(nBndEdgePts, &uplus[nScalars - 1][id2], 1,
1490 &penaltyfluxO1[nScalars - 1][id2], 1);
1491 }
1492 }
1493 }
1494}
1495
1496/**
1497 * @brief Build the numerical flux for the 2nd order derivatives
1498 *
1499 */
1500void DiffusionLFRNS::NumericalFluxO2(
1501 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1502 const Array<OneD, Array<OneD, NekDouble>> &ufield,
1503 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &qfield,
1504 Array<OneD, Array<OneD, NekDouble>> &qflux)
1505{
1506 int i, j;
1507 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1508 int nVariables = fields.size();
1509 int nDim = fields[0]->GetCoordim(0);
1510
1511 Array<OneD, NekDouble> Fwd(nTracePts);
1512 Array<OneD, NekDouble> Bwd(nTracePts);
1513 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
1514
1515 Array<OneD, NekDouble> qFwd(nTracePts);
1516 Array<OneD, NekDouble> qBwd(nTracePts);
1517 Array<OneD, NekDouble> qfluxtemp(nTracePts, 0.0);
1518
1519 // Get the normal velocity Vn
1520 for (i = 0; i < nDim; ++i)
1521 {
1522 fields[0]->ExtractTracePhys(ufield[i], m_traceVel[i]);
1523 Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, m_traceVel[i], 1, Vn, 1,
1524 Vn, 1);
1525 }
1526
1527 // Evaulate Riemann flux
1528 // qflux = \hat{q} \cdot u = q \cdot n
1529 // Notice: i = 1 (first row of the viscous tensor is zero)
1530 for (i = 1; i < nVariables; ++i)
1531 {
1532 qflux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1533 for (j = 0; j < nDim; ++j)
1534 {
1535 // Compute qFwd and qBwd value of qfield in position 'ji'
1536 fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
1537
1538 // Get Riemann flux of qflux --> LDG implies upwind
1539 fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1540
1541 // Multiply the Riemann flux by the trace normals
1542 Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
1543 qfluxtemp, 1);
1544
1545 // Impose weak boundary condition with flux
1546 if (fields[0]->GetBndCondExpansions().size())
1547 {
1548 WeakPenaltyO2(fields, i, j, qfield[j][i], qfluxtemp);
1549 }
1550
1551 // Store the final flux into qflux
1552 Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1, qflux[i], 1);
1553 }
1554 }
1555}
1556
1557/**
1558 * @brief Imposes appropriate bcs for the 2nd order derivatives
1559 *
1560 */
1561void DiffusionLFRNS::WeakPenaltyO2(
1562 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields, const int var,
1563 const int dir, const Array<OneD, const NekDouble> &qfield,
1564 Array<OneD, NekDouble> &penaltyflux)
1565{
1566 int cnt = 0;
1567 int nBndEdges, nBndEdgePts;
1568 int i, e;
1569 int id2;
1570
1571 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1572 int nBndRegions = fields[var]->GetBndCondExpansions().size();
1573
1574 Array<OneD, NekDouble> uterm(nTracePts);
1575 Array<OneD, NekDouble> qtemp(nTracePts);
1576
1577 // Extract the physical values of the solution at the boundaries
1578 fields[var]->ExtractTracePhys(qfield, qtemp);
1579
1580 // Loop on the boundary regions to apply appropriate bcs
1581 for (i = 0; i < nBndRegions; ++i)
1582 {
1583 // Number of boundary regions related to region 'i'
1584 nBndEdges = fields[var]->GetBndCondExpansions()[i]->GetExpSize();
1585
1586 if (fields[var]->GetBndConditions()[i]->GetBoundaryConditionType() ==
1587 SpatialDomains::ePeriodic)
1588 {
1589 continue;
1590 }
1591
1592 // Weakly impose bcs by modifying flux values
1593 for (e = 0; e < nBndEdges; ++e)
1594 {
1595 nBndEdgePts = fields[var]
1596 ->GetBndCondExpansions()[i]
1597 ->GetExp(e)
1598 ->GetTotPoints();
1599
1600 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1601 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1602
1603 // In case of Dirichlet bcs:
1604 // uflux = gD
1605 if (fields[var]
1606 ->GetBndConditions()[i]
1607 ->GetBoundaryConditionType() ==
1608 SpatialDomains::eDirichlet &&
1609 !boost::iequals(
1610 fields[var]->GetBndConditions()[i]->GetUserDefined(),
1611 "WallAdiabatic"))
1612 {
1613 Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1614 &qtemp[id2], 1, &penaltyflux[id2], 1);
1615 }
1616 // 3.4) In case of Neumann bcs:
1617 // uflux = u+
1618 else if ((fields[var]->GetBndConditions()[i])
1619 ->GetBoundaryConditionType() ==
1620 SpatialDomains::eNeumann)
1621 {
1622 ASSERTL0(false, "Neumann bcs not implemented for LFRNS");
1623 }
1624 else if (boost::iequals(
1625 fields[var]->GetBndConditions()[i]->GetUserDefined(),
1626 "WallAdiabatic"))
1627 {
1628 if ((var == m_spaceDim + 1))
1629 {
1630 Vmath::Zero(nBndEdgePts, &penaltyflux[id2], 1);
1631 }
1632 else
1633 {
1634
1635 Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1636 &qtemp[id2], 1, &penaltyflux[id2], 1);
1637 }
1638 }
1639 }
1640 }
1641}
1642
1643/**
1644 * @brief Compute the derivative of the corrective flux for 1D problems.
1645 *
1646 * @param nConvectiveFields Number of fields.
1647 * @param fields Pointer to fields.
1648 * @param flux Volumetric flux in the physical space.
1649 * @param iFlux Numerical interface flux in physical space.
1650 * @param derCFlux Derivative of the corrective flux.
1651 *
1652 */
1653void DiffusionLFRNS::DerCFlux_1D(
1654 [[maybe_unused]] const int nConvectiveFields,
1655 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1656 const Array<OneD, const NekDouble> &flux,
1657 const Array<OneD, const NekDouble> &iFlux, Array<OneD, NekDouble> &derCFlux)
1658{
1659 int n;
1660 int nLocalSolutionPts, phys_offset;
1661
1662 Array<OneD, NekDouble> auxArray1, auxArray2;
1663 Array<TwoD, const NekDouble> gmat;
1664 Array<OneD, const NekDouble> jac;
1665
1666 LibUtilities::PointsKeyVector ptsKeys;
1667 LibUtilities::BasisSharedPtr Basis;
1668 Basis = fields[0]->GetExp(0)->GetBasis(0);
1669
1670 int nElements = fields[0]->GetExpSize();
1671 int nPts = fields[0]->GetTotPoints();
1672
1673 // Arrays to store the derivatives of the correction flux
1674 Array<OneD, NekDouble> DCL(nPts / nElements, 0.0);
1675 Array<OneD, NekDouble> DCR(nPts / nElements, 0.0);
1676
1677 // Arrays to store the intercell numerical flux jumps
1678 Array<OneD, NekDouble> JumpL(nElements);
1679 Array<OneD, NekDouble> JumpR(nElements);
1680
1681 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1682 fields[0]->GetTraceMap()->GetElmtToTrace();
1683
1684 for (n = 0; n < nElements; ++n)
1685 {
1686 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1687 phys_offset = fields[0]->GetPhys_Offset(n);
1688
1689 Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1690 Array<OneD, NekDouble> tmpFluxVertex(1, 0.0);
1691 Vmath::Vcopy(nLocalSolutionPts, &flux[phys_offset], 1, &tmparrayX1[0],
1692 1);
1693
1694 fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0], tmparrayX1,
1695 tmpFluxVertex);
1696 JumpL[n] = iFlux[n] - tmpFluxVertex[0];
1697
1698 fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1], tmparrayX1,
1699 tmpFluxVertex);
1700 JumpR[n] = iFlux[n + 1] - tmpFluxVertex[0];
1701 }
1702
1703 for (n = 0; n < nElements; ++n)
1704 {
1705 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1706 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1707 phys_offset = fields[0]->GetPhys_Offset(n);
1708 jac = fields[0]
1709 ->GetExp(n)
1710 ->as<LocalRegions::Expansion1D>()
1711 ->GetGeom1D()
1712 ->GetMetricInfo()
1713 ->GetJac(ptsKeys);
1714
1715 JumpL[n] = JumpL[n] * jac[0];
1716 JumpR[n] = JumpR[n] * jac[0];
1717
1718 // Left jump multiplied by left derivative of C function
1719 Vmath::Smul(nLocalSolutionPts, JumpL[n], &m_dGL_xi1[n][0], 1, &DCL[0],
1720 1);
1721
1722 // Right jump multiplied by right derivative of C function
1723 Vmath::Smul(nLocalSolutionPts, JumpR[n], &m_dGR_xi1[n][0], 1, &DCR[0],
1724 1);
1725
1726 // Assembling derivative of the correction flux
1727 Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1728 &derCFlux[phys_offset], 1);
1729 }
1730}
1731
1732/**
1733 * @brief Compute the derivative of the corrective flux wrt a given
1734 * coordinate for 2D problems.
1735 *
1736 * There could be a bug for deformed elements since first derivatives
1737 * are not exactly the same results of DiffusionLDG as expected
1738 *
1739 * @param nConvectiveFields Number of fields.
1740 * @param fields Pointer to fields.
1741 * @param flux Volumetric flux in the physical space.
1742 * @param numericalFlux Numerical interface flux in physical space.
1743 * @param derCFlux Derivative of the corrective flux.
1744 *
1745 * \todo: Switch on shapes eventually here.
1746 */
1747void DiffusionLFRNS::DerCFlux_2D(
1748 [[maybe_unused]] const int nConvectiveFields, const int direction,
1749 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1750 const Array<OneD, const NekDouble> &flux,
1751 const Array<OneD, NekDouble> &iFlux, Array<OneD, NekDouble> &derCFlux)
1752{
1753 int n, e, i, j, cnt;
1754
1755 Array<OneD, const NekDouble> jac;
1756
1757 int nElements = fields[0]->GetExpSize();
1758 int trace_offset, phys_offset;
1759 int nLocalSolutionPts;
1760 int nquad0, nquad1;
1761 int nEdgePts;
1762
1763 LibUtilities::PointsKeyVector ptsKeys;
1764 Array<OneD, NekDouble> auxArray1, auxArray2;
1765 Array<OneD, LibUtilities::BasisSharedPtr> base;
1766 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1767 fields[0]->GetTraceMap()->GetElmtToTrace();
1768
1769 // Loop on the elements
1770 for (n = 0; n < nElements; ++n)
1771 {
1772 // Offset of the element on the global vector
1773 phys_offset = fields[0]->GetPhys_Offset(n);
1774 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1775 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1776
1777 jac = fields[0]
1778 ->GetExp(n)
1779 ->as<LocalRegions::Expansion2D>()
1780 ->GetGeom2D()
1781 ->GetMetricInfo()
1782 ->GetJac(ptsKeys);
1783
1784 base = fields[0]->GetExp(n)->GetBase();
1785 nquad0 = base[0]->GetNumPoints();
1786 nquad1 = base[1]->GetNumPoints();
1787
1788 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1789 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1790 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1791 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1792
1793 // Loop on the edges
1794 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1795 {
1796 // Number of edge points of edge 'e'
1797 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1798
1799 // Array for storing volumetric fluxes on each edge
1800 Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1801
1802 // Offset of the trace space correspondent to edge 'e'
1803 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1804 elmtToTrace[n][e]->GetElmtId());
1805
1806 // Extract the edge values of the volumetric fluxes
1807 // on edge 'e' and order them accordingly to the order
1808 // of the trace space
1809 fields[0]->GetExp(n)->GetTracePhysVals(
1810 e, elmtToTrace[n][e], flux + phys_offset, auxArray1 = tmparray);
1811
1812 // Compute the fluxJumps per each edge 'e' and each
1813 // flux point
1814 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1815 Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1, &tmparray[0], 1,
1816 &fluxJumps[0], 1);
1817
1818 // Check the ordering of the fluxJumps and reverse
1819 // it in case of backward definition of edge 'e'
1820 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1821 StdRegions::eBackwards)
1822 {
1823 Vmath::Reverse(nEdgePts, &fluxJumps[0], 1, &fluxJumps[0], 1);
1824 }
1825
1826 // Deformed elements
1827 if (fields[0]
1828 ->GetExp(n)
1829 ->as<LocalRegions::Expansion2D>()
1830 ->GetGeom2D()
1831 ->GetMetricInfo()
1832 ->GetGtype() == SpatialDomains::eDeformed)
1833 {
1834 // Extract the Jacobians along edge 'e'
1835 Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1836
1837 fields[0]->GetExp(n)->GetTracePhysVals(
1838 e, elmtToTrace[n][e], jac, auxArray1 = jacEdge);
1839
1840 // Check the ordering of the fluxJumps and reverse
1841 // it in case of backward definition of edge 'e'
1842 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1843 StdRegions::eBackwards)
1844 {
1845 Vmath::Reverse(nEdgePts, &jacEdge[0], 1, &jacEdge[0], 1);
1846 }
1847
1848 // Multiply the fluxJumps by the edge 'e' Jacobians
1849 // to bring the fluxJumps into the standard space
1850 for (j = 0; j < nEdgePts; j++)
1851 {
1852 fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1853 }
1854 }
1855 // Non-deformed elements
1856 else
1857 {
1858 // Multiply the fluxJumps by the edge 'e' Jacobians
1859 // to bring the fluxJumps into the standard space
1860 Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1, fluxJumps, 1);
1861 }
1862
1863 // Multiply jumps by derivatives of the correction functions
1864 // All the quntities at this level should be defined into
1865 // the standard space
1866 switch (e)
1867 {
1868 case 0:
1869 for (i = 0; i < nquad0; ++i)
1870 {
1871 for (j = 0; j < nquad1; ++j)
1872 {
1873 cnt = i + j * nquad0;
1874 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1875 }
1876 }
1877 break;
1878 case 1:
1879 for (i = 0; i < nquad1; ++i)
1880 {
1881 for (j = 0; j < nquad0; ++j)
1882 {
1883 cnt = (nquad0)*i + j;
1884 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1885 }
1886 }
1887 break;
1888 case 2:
1889 for (i = 0; i < nquad0; ++i)
1890 {
1891 for (j = 0; j < nquad1; ++j)
1892 {
1893 cnt = j * nquad0 + i;
1894 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1895 }
1896 }
1897 break;
1898 case 3:
1899 for (i = 0; i < nquad1; ++i)
1900 {
1901 for (j = 0; j < nquad0; ++j)
1902 {
1903 cnt = j + i * nquad0;
1904 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1905 }
1906 }
1907 break;
1908
1909 default:
1910 ASSERTL0(false, "edge value (< 3) is out of range");
1911 break;
1912 }
1913 }
1914
1915 // Summing the component of the flux for calculating the
1916 // derivatives per each direction
1917 if (direction == 0)
1918 {
1919 Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1, &divCFluxE3[0], 1,
1920 &derCFlux[phys_offset], 1);
1921 }
1922 else if (direction == 1)
1923 {
1924 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE2[0], 1,
1925 &derCFlux[phys_offset], 1);
1926 }
1927 }
1928}
1929
1930/**
1931 * @brief Compute the divergence of the corrective flux for 2D problems.
1932 *
1933 * @param nConvectiveFields Number of fields.
1934 * @param fields Pointer to fields.
1935 * @param fluxX1 X1 - volumetric flux in physical space.
1936 * @param fluxX2 X2 - volumetric flux in physical space.
1937 * @param numericalFlux Interface flux in physical space.
1938 * @param divCFlux Divergence of the corrective flux for
1939 * 2D problems.
1940 *
1941 * \todo: Switch on shapes eventually here.
1942 */
1943void DiffusionLFRNS::DivCFlux_2D(
1944 [[maybe_unused]] const int nConvectiveFields,
1945 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1946 const Array<OneD, const NekDouble> &fluxX1,
1947 const Array<OneD, const NekDouble> &fluxX2,
1948 const Array<OneD, const NekDouble> &numericalFlux,
1949 Array<OneD, NekDouble> &divCFlux)
1950{
1951 int n, e, i, j, cnt;
1952
1953 int nElements = fields[0]->GetExpSize();
1954
1955 int nLocalSolutionPts;
1956 int nEdgePts;
1957 int trace_offset;
1958 int phys_offset;
1959 int nquad0;
1960 int nquad1;
1961
1962 Array<OneD, NekDouble> auxArray1, auxArray2;
1963 Array<OneD, LibUtilities::BasisSharedPtr> base;
1964
1965 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1966 fields[0]->GetTraceMap()->GetElmtToTrace();
1967
1968 // Loop on the elements
1969 for (n = 0; n < nElements; ++n)
1970 {
1971 // Offset of the element on the global vector
1972 phys_offset = fields[0]->GetPhys_Offset(n);
1973 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1974
1975 base = fields[0]->GetExp(n)->GetBase();
1976 nquad0 = base[0]->GetNumPoints();
1977 nquad1 = base[1]->GetNumPoints();
1978
1979 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1980 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1981 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1982 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1983
1984 // Loop on the edges
1985 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1986 {
1987 // Number of edge points of edge e
1988 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1989
1990 Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1991 Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1992 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1993 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1994 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1995
1996 // Offset of the trace space correspondent to edge e
1997 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1998 elmtToTrace[n][e]->GetElmtId());
1999
2000 // Get the normals of edge e
2001 const Array<OneD, const Array<OneD, NekDouble>> &normals =
2002 fields[0]->GetExp(n)->GetTraceNormal(e);
2003
2004 // Extract the edge values of flux-x on edge e and order
2005 // them accordingly to the order of the trace space
2006 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2007 fluxX1 + phys_offset,
2008 auxArray1 = tmparrayX1);
2009
2010 // Extract the edge values of flux-y on edge e and order
2011 // them accordingly to the order of the trace space
2012 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2013 fluxX2 + phys_offset,
2014 auxArray1 = tmparrayX2);
2015
2016 // Multiply the edge components of the flux by the normal
2017 for (i = 0; i < nEdgePts; ++i)
2018 {
2019 fluxN[i] = tmparrayX1[i] * m_traceNormals[0][trace_offset + i] +
2020 tmparrayX2[i] * m_traceNormals[1][trace_offset + i];
2021 }
2022
2023 // Subtract to the Riemann flux the discontinuous flux
2024 Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1, &fluxN[0], 1,
2025 &fluxJumps[0], 1);
2026
2027 // Check the ordering of the jump vectors
2028 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2029 StdRegions::eBackwards)
2030 {
2031 Vmath::Reverse(nEdgePts, auxArray1 = fluxJumps, 1,
2032 auxArray2 = fluxJumps, 1);
2033 }
2034
2035 for (i = 0; i < nEdgePts; ++i)
2036 {
2037 if (m_traceNormals[0][trace_offset + i] != normals[0][i] ||
2038 m_traceNormals[1][trace_offset + i] != normals[1][i])
2039 {
2040 fluxJumps[i] = -fluxJumps[i];
2041 }
2042 }
2043
2044 // Multiply jumps by derivatives of the correction functions
2045 switch (e)
2046 {
2047 case 0:
2048 for (i = 0; i < nquad0; ++i)
2049 {
2050 // Multiply fluxJumps by Q factors
2051 fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
2052
2053 for (j = 0; j < nquad1; ++j)
2054 {
2055 cnt = i + j * nquad0;
2056 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
2057 }
2058 }
2059 break;
2060 case 1:
2061 for (i = 0; i < nquad1; ++i)
2062 {
2063 // Multiply fluxJumps by Q factors
2064 fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
2065
2066 for (j = 0; j < nquad0; ++j)
2067 {
2068 cnt = (nquad0)*i + j;
2069 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
2070 }
2071 }
2072 break;
2073 case 2:
2074 for (i = 0; i < nquad0; ++i)
2075 {
2076 // Multiply fluxJumps by Q factors
2077 fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
2078
2079 for (j = 0; j < nquad1; ++j)
2080 {
2081 cnt = j * nquad0 + i;
2082 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
2083 }
2084 }
2085 break;
2086 case 3:
2087 for (i = 0; i < nquad1; ++i)
2088 {
2089 // Multiply fluxJumps by Q factors
2090 fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
2091 for (j = 0; j < nquad0; ++j)
2092 {
2093 cnt = j + i * nquad0;
2094 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
2095 }
2096 }
2097 break;
2098
2099 default:
2100 ASSERTL0(false, "edge value (< 3) is out of range");
2101 break;
2102 }
2103 }
2104
2105 // Sum each edge contribution
2106 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
2107 &divCFlux[phys_offset], 1);
2108
2109 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2110 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2111
2112 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2113 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2114 }
2115}
2116
2117/**
2118 * @brief Compute the divergence of the corrective flux for 2D problems
2119 * where POINTSTYPE="GaussGaussLegendre"
2120 *
2121 * @param nConvectiveFields Number of fields.
2122 * @param fields Pointer to fields.
2123 * @param fluxX1 X1-volumetric flux in physical space.
2124 * @param fluxX2 X2-volumetric flux in physical space.
2125 * @param numericalFlux Interface flux in physical space.
2126 * @param divCFlux Divergence of the corrective flux.
2127 *
2128 * \todo: Switch on shapes eventually here.
2129 */
2130
2131void DiffusionLFRNS::DivCFlux_2D_Gauss(
2132 [[maybe_unused]] const int nConvectiveFields,
2133 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
2134 const Array<OneD, const NekDouble> &fluxX1,
2135 const Array<OneD, const NekDouble> &fluxX2,
2136 const Array<OneD, const NekDouble> &numericalFlux,
2137 Array<OneD, NekDouble> &divCFlux)
2138{
2139 int n, e, i, j, cnt;
2140
2141 int nElements = fields[0]->GetExpSize();
2142 int nLocalSolutionPts;
2143 int nEdgePts;
2144 int trace_offset;
2145 int phys_offset;
2146 int nquad0;
2147 int nquad1;
2148
2149 Array<OneD, NekDouble> auxArray1, auxArray2;
2150 Array<OneD, LibUtilities::BasisSharedPtr> base;
2151
2152 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
2153 fields[0]->GetTraceMap()->GetElmtToTrace();
2154
2155 // Loop on the elements
2156 for (n = 0; n < nElements; ++n)
2157 {
2158 // Offset of the element on the global vector
2159 phys_offset = fields[0]->GetPhys_Offset(n);
2160 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
2161
2162 base = fields[0]->GetExp(n)->GetBase();
2163 nquad0 = base[0]->GetNumPoints();
2164 nquad1 = base[1]->GetNumPoints();
2165
2166 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
2167 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
2168 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
2169 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
2170
2171 // Loop on the edges
2172 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
2173 {
2174 // Number of edge points of edge e
2175 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
2176
2177 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
2178 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
2179 Array<OneD, NekDouble> fluxN_R(nEdgePts, 0.0);
2180 Array<OneD, NekDouble> fluxN_D(nEdgePts, 0.0);
2181 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
2182
2183 // Offset of the trace space correspondent to edge e
2184 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
2185 elmtToTrace[n][e]->GetElmtId());
2186
2187 // Get the normals of edge e
2188 const Array<OneD, const Array<OneD, NekDouble>> &normals =
2189 fields[0]->GetExp(n)->GetTraceNormal(e);
2190
2191 // Extract the trasformed normal flux at each edge
2192 switch (e)
2193 {
2194 case 0:
2195 // Extract the edge values of transformed flux-y on
2196 // edge e and order them accordingly to the order of
2197 // the trace space
2198 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2199 fluxX2 + phys_offset,
2200 auxArray1 = fluxN_D);
2201
2202 Vmath::Neg(nEdgePts, fluxN_D, 1);
2203
2204 // Extract the physical Rieamann flux at each edge
2205 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2206 &fluxN[0], 1);
2207
2208 // Check the ordering of vectors
2209 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2210 StdRegions::eBackwards)
2211 {
2212 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2213 auxArray2 = fluxN, 1);
2214
2215 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2216 auxArray2 = fluxN_D, 1);
2217 }
2218
2219 // Transform Riemann Fluxes in the standard element
2220 for (i = 0; i < nquad0; ++i)
2221 {
2222 // Multiply Riemann Flux by Q factors
2223 fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
2224 }
2225
2226 for (i = 0; i < nEdgePts; ++i)
2227 {
2228 if (m_traceNormals[0][trace_offset + i] !=
2229 normals[0][i] ||
2230 m_traceNormals[1][trace_offset + i] !=
2231 normals[1][i])
2232 {
2233 fluxN_R[i] = -fluxN_R[i];
2234 }
2235 }
2236
2237 // Subtract to the Riemann flux the discontinuous
2238 // flux
2239 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2240 &fluxJumps[0], 1);
2241
2242 // Multiplicate the flux jumps for the correction
2243 // function
2244 for (i = 0; i < nquad0; ++i)
2245 {
2246 for (j = 0; j < nquad1; ++j)
2247 {
2248 cnt = i + j * nquad0;
2249 divCFluxE0[cnt] = -fluxJumps[i] * m_dGL_xi2[n][j];
2250 }
2251 }
2252 break;
2253 case 1:
2254 // Extract the edge values of transformed flux-x on
2255 // edge e and order them accordingly to the order of
2256 // the trace space
2257 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2258 fluxX1 + phys_offset,
2259 auxArray1 = fluxN_D);
2260
2261 // Extract the physical Rieamann flux at each edge
2262 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2263 &fluxN[0], 1);
2264
2265 // Check the ordering of vectors
2266 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2267 StdRegions::eBackwards)
2268 {
2269 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2270 auxArray2 = fluxN, 1);
2271
2272 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2273 auxArray2 = fluxN_D, 1);
2274 }
2275
2276 // Transform Riemann Fluxes in the standard element
2277 for (i = 0; i < nquad1; ++i)
2278 {
2279 // Multiply Riemann Flux by Q factors
2280 fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
2281 }
2282
2283 for (i = 0; i < nEdgePts; ++i)
2284 {
2285 if (m_traceNormals[0][trace_offset + i] !=
2286 normals[0][i] ||
2287 m_traceNormals[1][trace_offset + i] !=
2288 normals[1][i])
2289 {
2290 fluxN_R[i] = -fluxN_R[i];
2291 }
2292 }
2293
2294 // Subtract to the Riemann flux the discontinuous
2295 // flux
2296 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2297 &fluxJumps[0], 1);
2298
2299 // Multiplicate the flux jumps for the correction
2300 // function
2301 for (i = 0; i < nquad1; ++i)
2302 {
2303 for (j = 0; j < nquad0; ++j)
2304 {
2305 cnt = (nquad0)*i + j;
2306 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
2307 }
2308 }
2309 break;
2310 case 2:
2311
2312 // Extract the edge values of transformed flux-y on
2313 // edge e and order them accordingly to the order of
2314 // the trace space
2315
2316 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2317 fluxX2 + phys_offset,
2318 auxArray1 = fluxN_D);
2319
2320 // Extract the physical Rieamann flux at each edge
2321 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2322 &fluxN[0], 1);
2323
2324 // Check the ordering of vectors
2325 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2326 StdRegions::eBackwards)
2327 {
2328 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2329 auxArray2 = fluxN, 1);
2330
2331 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2332 auxArray2 = fluxN_D, 1);
2333 }
2334
2335 // Transform Riemann Fluxes in the standard element
2336 for (i = 0; i < nquad0; ++i)
2337 {
2338 // Multiply Riemann Flux by Q factors
2339 fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
2340 }
2341
2342 for (i = 0; i < nEdgePts; ++i)
2343 {
2344 if (m_traceNormals[0][trace_offset + i] !=
2345 normals[0][i] ||
2346 m_traceNormals[1][trace_offset + i] !=
2347 normals[1][i])
2348 {
2349 fluxN_R[i] = -fluxN_R[i];
2350 }
2351 }
2352
2353 // Subtract to the Riemann flux the discontinuous
2354 // flux
2355
2356 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2357 &fluxJumps[0], 1);
2358
2359 // Multiplicate the flux jumps for the correction
2360 // function
2361 for (i = 0; i < nquad0; ++i)
2362 {
2363 for (j = 0; j < nquad1; ++j)
2364 {
2365 cnt = j * nquad0 + i;
2366 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
2367 }
2368 }
2369 break;
2370 case 3:
2371 // Extract the edge values of transformed flux-x on
2372 // edge e and order them accordingly to the order of
2373 // the trace space
2374
2375 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2376 fluxX1 + phys_offset,
2377 auxArray1 = fluxN_D);
2378 Vmath::Neg(nEdgePts, fluxN_D, 1);
2379
2380 // Extract the physical Rieamann flux at each edge
2381 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2382 &fluxN[0], 1);
2383
2384 // Check the ordering of vectors
2385 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2386 StdRegions::eBackwards)
2387 {
2388 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2389 auxArray2 = fluxN, 1);
2390
2391 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2392 auxArray2 = fluxN_D, 1);
2393 }
2394
2395 // Transform Riemann Fluxes in the standard element
2396 for (i = 0; i < nquad1; ++i)
2397 {
2398 // Multiply Riemann Flux by Q factors
2399 fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
2400 }
2401
2402 for (i = 0; i < nEdgePts; ++i)
2403 {
2404 if (m_traceNormals[0][trace_offset + i] !=
2405 normals[0][i] ||
2406 m_traceNormals[1][trace_offset + i] !=
2407 normals[1][i])
2408 {
2409 fluxN_R[i] = -fluxN_R[i];
2410 }
2411 }
2412
2413 // Subtract to the Riemann flux the discontinuous
2414 // flux
2415
2416 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2417 &fluxJumps[0], 1);
2418
2419 // Multiplicate the flux jumps for the correction
2420 // function
2421 for (i = 0; i < nquad1; ++i)
2422 {
2423 for (j = 0; j < nquad0; ++j)
2424 {
2425 cnt = j + i * nquad0;
2426 divCFluxE3[cnt] = -fluxJumps[i] * m_dGL_xi1[n][j];
2427 }
2428 }
2429 break;
2430 default:
2431 ASSERTL0(false, "edge value (< 3) is out of range");
2432 break;
2433 }
2434 }
2435
2436 // Sum each edge contribution
2437 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
2438 &divCFlux[phys_offset], 1);
2439
2440 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2441 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2442
2443 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2444 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2445 }
2446}
2447
2448} // namespace Nektar::SolverUtils
ASSERTL0
#define ASSERTL0(condition, msg)
Definition ErrorUtil.hpp:216
WARNINGL0
#define WARNINGL0(condition, msg)
Definition ErrorUtil.hpp:223
Nektar::LibUtilities::Basis
Represents a basis of a given type.
Definition Basis.h:198
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition BasicUtils/NekFactory.hpp:196
Nektar::LocalRegions::Expansion::GetMetricInfo
const SpatialDomains::GeomFactorsSharedPtr & GetMetricInfo() const
Definition Expansion.cpp:246
Nektar::SolverUtils::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:1653
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:1230
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:1285
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:1943
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:1747
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:2131
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:1500
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:1561
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::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition BasicUtils/SharedArray.hpp:928
Polylib::jacobd
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition Polylib.cpp:1378
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition Vmath.hpp:72
Vmath::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 Sun Nov 9 2025 20:16:59 for Nektar++ by doxygen 1.9.8