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