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