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