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