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