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.size();
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)->GetTraceQFactors(
300  0, auxArray1 = m_Q2D_e0[n]);
301  pFields[0]->GetExp(n)->GetTraceQFactors(
302  1, auxArray1 = m_Q2D_e1[n]);
303  pFields[0]->GetExp(n)->GetTraceQFactors(
304  2, auxArray1 = m_Q2D_e2[n]);
305  pFields[0]->GetExp(n)->GetTraceQFactors(
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.size();
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.size();
1251  int nDim = fields[0]->GetCoordim(0);
1252 
1253  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1254  Array<OneD, NekDouble > fluxtemp(nTracePts, 0.0);
1255 
1256  // Get the normal velocity Vn
1257  for (i = 0; i < nDim; ++i)
1258  {
1259  fields[0]->ExtractTracePhys(inarray[i], m_traceVel[i]);
1260  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1,
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().size())
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.size();
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().size();
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  GetBndCondIDToGlobalTraceID(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().size();
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  GetBndCondIDToGlobalTraceID(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.size();
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().size())
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().size();
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  GetBndCondIDToGlobalTraceID(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  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1711  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1712 
1713  for (n = 0; n < nElements; ++n)
1714  {
1715  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1716  phys_offset = fields[0]->GetPhys_Offset(n);
1717 
1718  Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1719  Array<OneD, NekDouble> tmpFluxVertex(1,0.0);
1720  Vmath::Vcopy(nLocalSolutionPts,
1721  &flux[phys_offset], 1,
1722  &tmparrayX1[0], 1);
1723 
1724  fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0],
1725  tmparrayX1,
1726  tmpFluxVertex);
1727  JumpL[n] = iFlux[n] - tmpFluxVertex[0];
1728 
1729  fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1],
1730  tmparrayX1,
1731  tmpFluxVertex);
1732  JumpR[n] = iFlux[n+1] - tmpFluxVertex[0];
1733  }
1734 
1735  for (n = 0; n < nElements; ++n)
1736  {
1737  ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1738  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1739  phys_offset = fields[0]->GetPhys_Offset(n);
1740  jac = fields[0]->GetExp(n)
1741  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
1742  ->GetMetricInfo()->GetJac(ptsKeys);
1743 
1744  JumpL[n] = JumpL[n] * jac[0];
1745  JumpR[n] = JumpR[n] * jac[0];
1746 
1747  // Left jump multiplied by left derivative of C function
1748  Vmath::Smul(nLocalSolutionPts, JumpL[n],
1749  &m_dGL_xi1[n][0], 1, &DCL[0], 1);
1750 
1751  // Right jump multiplied by right derivative of C function
1752  Vmath::Smul(nLocalSolutionPts, JumpR[n],
1753  &m_dGR_xi1[n][0], 1, &DCR[0], 1);
1754 
1755  // Assembling derivative of the correction flux
1756  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1757  &derCFlux[phys_offset], 1);
1758  }
1759  }
1760 
1761  /**
1762  * @brief Compute the derivative of the corrective flux wrt a given
1763  * coordinate for 2D problems.
1764  *
1765  * There could be a bug for deformed elements since first derivatives
1766  * are not exactly the same results of DiffusionLDG as expected
1767  *
1768  * @param nConvectiveFields Number of fields.
1769  * @param fields Pointer to fields.
1770  * @param flux Volumetric flux in the physical space.
1771  * @param numericalFlux Numerical interface flux in physical space.
1772  * @param derCFlux Derivative of the corrective flux.
1773  *
1774  * \todo: Switch on shapes eventually here.
1775  */
1776  void DiffusionLFRNS::v_DerCFlux_2D(
1777  const int nConvectiveFields,
1778  const int direction,
1779  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1780  const Array<OneD, const NekDouble> &flux,
1781  const Array<OneD, NekDouble> &iFlux,
1782  Array<OneD, NekDouble> &derCFlux)
1783  {
1784  boost::ignore_unused(nConvectiveFields);
1785 
1786  int n, e, i, j, cnt;
1787 
1788  Array<OneD, const NekDouble> jac;
1789 
1790  int nElements = fields[0]->GetExpSize();
1791  int trace_offset, phys_offset;
1792  int nLocalSolutionPts;
1793  int nquad0, nquad1;
1794  int nEdgePts;
1795 
1796  LibUtilities::PointsKeyVector ptsKeys;
1797  Array<OneD, NekDouble> auxArray1, auxArray2;
1798  Array<OneD, LibUtilities::BasisSharedPtr> base;
1799  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1800  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1801 
1802  // Loop on the elements
1803  for (n = 0; n < nElements; ++n)
1804  {
1805  // Offset of the element on the global vector
1806  phys_offset = fields[0]->GetPhys_Offset(n);
1807  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1808  ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1809 
1810  jac = fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1811  ->GetGeom2D()->GetMetricInfo()->GetJac(ptsKeys);
1812 
1813  base = fields[0]->GetExp(n)->GetBase();
1814  nquad0 = base[0]->GetNumPoints();
1815  nquad1 = base[1]->GetNumPoints();
1816 
1817  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1818  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1819  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1820  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1821 
1822  // Loop on the edges
1823  for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1824  {
1825  // Number of edge points of edge 'e'
1826  nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1827 
1828  // Array for storing volumetric fluxes on each edge
1829  Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1830 
1831  // Offset of the trace space correspondent to edge 'e'
1832  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1833  elmtToTrace[n][e]->GetElmtId());
1834 
1835  // Get the normals of edge 'e'
1836  //const Array<OneD, const Array<OneD, NekDouble> > &normals =
1837  //fields[0]->GetExp(n)->GetTraceNormal(e);
1838 
1839  // Extract the edge values of the volumetric fluxes
1840  // on edge 'e' and order them accordingly to the order
1841  // of the trace space
1842  fields[0]->GetExp(n)->GetTracePhysVals(
1843  e, elmtToTrace[n][e],
1844  flux + phys_offset,
1845  auxArray1 = tmparray);
1846 
1847  // Splitting inarray into the 'j' direction
1848  /*Vmath::Vmul(nEdgePts,
1849  &m_traceNormals[direction][trace_offset], 1,
1850  &tmparray[0], 1, &tmparray[0], 1);*/
1851 
1852  // Compute the fluxJumps per each edge 'e' and each
1853  // flux point
1854  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1855  Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1,
1856  &tmparray[0], 1, &fluxJumps[0], 1);
1857 
1858  // Check the ordering of the fluxJumps and reverse
1859  // it in case of backward definition of edge 'e'
1860  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1861  StdRegions::eBackwards)
1862  {
1863  Vmath::Reverse(nEdgePts,
1864  &fluxJumps[0], 1,
1865  &fluxJumps[0], 1);
1866  }
1867 
1868  // Deformed elements
1869  if (fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1870  ->GetGeom2D()->GetMetricInfo()->GetGtype()
1871  == SpatialDomains::eDeformed)
1872  {
1873  // Extract the Jacobians along edge 'e'
1874  Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1875  fields[0]->GetExp(n)->GetTracePhysVals(
1876  e, elmtToTrace[n][e],
1877  jac, auxArray1 = jacEdge);
1878 
1879  // Check the ordering of the fluxJumps and reverse
1880  // it in case of backward definition of edge 'e'
1881  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1882  StdRegions::eBackwards)
1883  {
1884  Vmath::Reverse(nEdgePts, &jacEdge[0], 1,
1885  &jacEdge[0], 1);
1886  }
1887 
1888  // Multiply the fluxJumps by the edge 'e' Jacobians
1889  // to bring the fluxJumps into the standard space
1890  for (j = 0; j < nEdgePts; j++)
1891  {
1892  fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1893  }
1894  }
1895  // Non-deformed elements
1896  else
1897  {
1898  // Multiply the fluxJumps by the edge 'e' Jacobians
1899  // to bring the fluxJumps into the standard space
1900  Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1,
1901  fluxJumps, 1);
1902  }
1903 
1904  // Multiply jumps by derivatives of the correction functions
1905  // All the quntities at this level should be defined into
1906  // the standard space
1907  switch (e)
1908  {
1909  case 0:
1910  for (i = 0; i < nquad0; ++i)
1911  {
1912  // Multiply fluxJumps by Q factors
1913  //fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1914 
1915  for (j = 0; j < nquad1; ++j)
1916  {
1917  cnt = i + j*nquad0;
1918  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1919  }
1920  }
1921  break;
1922  case 1:
1923  for (i = 0; i < nquad1; ++i)
1924  {
1925  // Multiply fluxJumps by Q factors
1926  //fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1927 
1928  for (j = 0; j < nquad0; ++j)
1929  {
1930  cnt = (nquad0)*i + j;
1931  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1932  }
1933  }
1934  break;
1935  case 2:
1936  for (i = 0; i < nquad0; ++i)
1937  {
1938  // Multiply fluxJumps by Q factors
1939  //fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1940 
1941  for (j = 0; j < nquad1; ++j)
1942  {
1943  cnt = j*nquad0 + i;
1944  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1945  }
1946  }
1947  break;
1948  case 3:
1949  for (i = 0; i < nquad1; ++i)
1950  {
1951  // Multiply fluxJumps by Q factors
1952  //fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1953  for (j = 0; j < nquad0; ++j)
1954  {
1955  cnt = j + i*nquad0;
1956  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1957  }
1958  }
1959  break;
1960 
1961  default:
1962  ASSERTL0(false, "edge value (< 3) is out of range");
1963  break;
1964  }
1965  }
1966 
1967 
1968  // Sum all the edge contributions since I am passing the
1969  // component of the flux x and y already. So I should not
1970  // need to sum E0-E2 to get the derivative wrt xi2 and E1-E3
1971  // to get the derivative wrt xi1
1972 
1973  if (direction == 0)
1974  {
1975  Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1,
1976  &divCFluxE3[0], 1, &derCFlux[phys_offset], 1);
1977  }
1978  else if (direction == 1)
1979  {
1980  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1981  &divCFluxE2[0], 1, &derCFlux[phys_offset], 1);
1982  }
1983  }
1984  }
1985 
1986  /**
1987  * @brief Compute the divergence of the corrective flux for 2D problems.
1988  *
1989  * @param nConvectiveFields Number of fields.
1990  * @param fields Pointer to fields.
1991  * @param fluxX1 X1 - volumetric flux in physical space.
1992  * @param fluxX2 X2 - volumetric flux in physical space.
1993  * @param numericalFlux Interface flux in physical space.
1994  * @param divCFlux Divergence of the corrective flux for
1995  * 2D problems.
1996  *
1997  * \todo: Switch on shapes eventually here.
1998  */
1999  void DiffusionLFRNS::v_DivCFlux_2D(
2000  const int nConvectiveFields,
2001  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
2002  const Array<OneD, const NekDouble> &fluxX1,
2003  const Array<OneD, const NekDouble> &fluxX2,
2004  const Array<OneD, const NekDouble> &numericalFlux,
2005  Array<OneD, NekDouble> &divCFlux)
2006  {
2007  boost::ignore_unused(nConvectiveFields);
2008 
2009  int n, e, i, j, cnt;
2010 
2011  int nElements = fields[0]->GetExpSize();
2012  int nLocalSolutionPts;
2013  int nEdgePts;
2014  int trace_offset;
2015  int phys_offset;
2016  int nquad0;
2017  int nquad1;
2018 
2019  Array<OneD, NekDouble> auxArray1, auxArray2;
2020  Array<OneD, LibUtilities::BasisSharedPtr> base;
2021 
2022  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
2023  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
2024 
2025  // Loop on the elements
2026  for(n = 0; n < nElements; ++n)
2027  {
2028  // Offset of the element on the global vector
2029  phys_offset = fields[0]->GetPhys_Offset(n);
2030  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
2031 
2032  base = fields[0]->GetExp(n)->GetBase();
2033  nquad0 = base[0]->GetNumPoints();
2034  nquad1 = base[1]->GetNumPoints();
2035 
2036  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
2037  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
2038  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
2039  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
2040 
2041  // Loop on the edges
2042  for(e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
2043  {
2044  // Number of edge points of edge e
2045  nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
2046 
2047  Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
2048  Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
2049  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
2050  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
2051  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
2052 
2053  // Offset of the trace space correspondent to edge e
2054  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
2055  elmtToTrace[n][e]->GetElmtId());
2056 
2057  // Get the normals of edge e
2058  const Array<OneD, const Array<OneD, NekDouble> > &normals =
2059  fields[0]->GetExp(n)->GetTraceNormal(e);
2060 
2061  // Extract the edge values of flux-x on edge e and order
2062  // them accordingly to the order of the trace space
2063  fields[0]->GetExp(n)->GetTracePhysVals(
2064  e, elmtToTrace[n][e],
2065  fluxX1 + phys_offset,
2066  auxArray1 = tmparrayX1);
2067 
2068  // Extract the edge values of flux-y on edge e and order
2069  // them accordingly to the order of the trace space
2070  fields[0]->GetExp(n)->GetTracePhysVals(
2071  e, elmtToTrace[n][e],
2072  fluxX2 + phys_offset,
2073  auxArray1 = tmparrayX2);
2074 
2075  // Multiply the edge components of the flux by the normal
2076  for (i = 0; i < nEdgePts; ++i)
2077  {
2078  fluxN[i] =
2079  tmparrayX1[i]*m_traceNormals[0][trace_offset+i] +
2080  tmparrayX2[i]*m_traceNormals[1][trace_offset+i];
2081  }
2082 
2083  // Subtract to the Riemann flux the discontinuous flux
2084  Vmath::Vsub(nEdgePts,
2085  &numericalFlux[trace_offset], 1,
2086  &fluxN[0], 1, &fluxJumps[0], 1);
2087 
2088  // Check the ordering of the jump vectors
2089  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2090  StdRegions::eBackwards)
2091  {
2092  Vmath::Reverse(nEdgePts,
2093  auxArray1 = fluxJumps, 1,
2094  auxArray2 = fluxJumps, 1);
2095  }
2096 
2097  for (i = 0; i < nEdgePts; ++i)
2098  {
2099  if (m_traceNormals[0][trace_offset+i] != normals[0][i]
2100  || m_traceNormals[1][trace_offset+i] != normals[1][i])
2101  {
2102  fluxJumps[i] = -fluxJumps[i];
2103  }
2104  }
2105 
2106  // Multiply jumps by derivatives of the correction functions
2107  switch (e)
2108  {
2109  case 0:
2110  for (i = 0; i < nquad0; ++i)
2111  {
2112  // Multiply fluxJumps by Q factors
2113  fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
2114 
2115  for (j = 0; j < nquad1; ++j)
2116  {
2117  cnt = i + j*nquad0;
2118  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
2119  }
2120  }
2121  break;
2122  case 1:
2123  for (i = 0; i < nquad1; ++i)
2124  {
2125  // Multiply fluxJumps by Q factors
2126  fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
2127 
2128  for (j = 0; j < nquad0; ++j)
2129  {
2130  cnt = (nquad0)*i + j;
2131  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
2132  }
2133  }
2134  break;
2135  case 2:
2136  for (i = 0; i < nquad0; ++i)
2137  {
2138  // Multiply fluxJumps by Q factors
2139  fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
2140 
2141  for (j = 0; j < nquad1; ++j)
2142  {
2143  cnt = j*nquad0 + i;
2144  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
2145  }
2146  }
2147  break;
2148  case 3:
2149  for (i = 0; i < nquad1; ++i)
2150  {
2151  // Multiply fluxJumps by Q factors
2152  fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
2153  for (j = 0; j < nquad0; ++j)
2154  {
2155  cnt = j + i*nquad0;
2156  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
2157 
2158  }
2159  }
2160  break;
2161 
2162  default:
2163  ASSERTL0(false,"edge value (< 3) is out of range");
2164  break;
2165  }
2166  }
2167 
2168  // Sum each edge contribution
2169  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
2170  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
2171 
2172  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2173  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2174 
2175  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2176  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2177  }
2178  }
2179 
2180  /**
2181  * @brief Compute the divergence of the corrective flux for 2D problems
2182  * where POINTSTYPE="GaussGaussLegendre"
2183  *
2184  * @param nConvectiveFields Number of fields.
2185  * @param fields Pointer to fields.
2186  * @param fluxX1 X1-volumetric flux in physical space.
2187  * @param fluxX2 X2-volumetric flux in physical space.
2188  * @param numericalFlux Interface flux in physical space.
2189  * @param divCFlux Divergence of the corrective flux.
2190  *
2191  * \todo: Switch on shapes eventually here.
2192  */
2193 
2194  void DiffusionLFRNS::v_DivCFlux_2D_Gauss(
2195  const int nConvectiveFields,
2196  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
2197  const Array<OneD, const NekDouble> &fluxX1,
2198  const Array<OneD, const NekDouble> &fluxX2,
2199  const Array<OneD, const NekDouble> &numericalFlux,
2200  Array<OneD, NekDouble> &divCFlux)
2201  {
2202  boost::ignore_unused(nConvectiveFields);
2203 
2204  int n, e, i, j, cnt;
2205 
2206  int nElements = fields[0]->GetExpSize();
2207  int nLocalSolutionPts;
2208  int nEdgePts;
2209  int trace_offset;
2210  int phys_offset;
2211  int nquad0;
2212  int nquad1;
2213 
2214  Array<OneD, NekDouble> auxArray1, auxArray2;
2215  Array<OneD, LibUtilities::BasisSharedPtr> base;
2216 
2217  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
2218  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
2219 
2220  // Loop on the elements
2221  for(n = 0; n < nElements; ++n)
2222  {
2223  // Offset of the element on the global vector
2224  phys_offset = fields[0]->GetPhys_Offset(n);
2225  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
2226 
2227  base = fields[0]->GetExp(n)->GetBase();
2228  nquad0 = base[0]->GetNumPoints();
2229  nquad1 = base[1]->GetNumPoints();
2230 
2231  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
2232  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
2233  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
2234  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
2235 
2236  // Loop on the edges
2237  for(e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
2238  {
2239  // Number of edge points of edge e
2240  nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
2241 
2242  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
2243  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
2244  Array<OneD, NekDouble> fluxN_R (nEdgePts, 0.0);
2245  Array<OneD, NekDouble> fluxN_D (nEdgePts, 0.0);
2246  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
2247 
2248  // Offset of the trace space correspondent to edge e
2249  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
2250  elmtToTrace[n][e]->GetElmtId());
2251 
2252  // Get the normals of edge e
2253  const Array<OneD, const Array<OneD, NekDouble> > &normals =
2254  fields[0]->GetExp(n)->GetTraceNormal(e);
2255 
2256  // Extract the trasformed normal flux at each edge
2257  switch (e)
2258  {
2259  case 0:
2260  // Extract the edge values of transformed flux-y on
2261  // edge e and order them accordingly to the order of
2262  // the trace space
2263  fields[0]->GetExp(n)->GetTracePhysVals(
2264  e, elmtToTrace[n][e],
2265  fluxX2 + phys_offset,
2266  auxArray1 = fluxN_D);
2267 
2268  Vmath::Neg (nEdgePts, fluxN_D, 1);
2269 
2270  // Extract the physical Rieamann flux at each edge
2271  Vmath::Vcopy(nEdgePts,
2272  &numericalFlux[trace_offset], 1,
2273  &fluxN[0], 1);
2274 
2275  // Check the ordering of vectors
2276  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2277  StdRegions::eBackwards)
2278  {
2279  Vmath::Reverse(nEdgePts,
2280  auxArray1 = fluxN, 1,
2281  auxArray2 = fluxN, 1);
2282 
2283  Vmath::Reverse(nEdgePts,
2284  auxArray1 = fluxN_D, 1,
2285  auxArray2 = fluxN_D, 1);
2286  }
2287 
2288  // Transform Riemann Fluxes in the standard element
2289  for (i = 0; i < nquad0; ++i)
2290  {
2291  // Multiply Riemann Flux by Q factors
2292  fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
2293  }
2294 
2295  for (i = 0; i < nEdgePts; ++i)
2296  {
2297  if (m_traceNormals[0][trace_offset+i]
2298  != normals[0][i] ||
2299  m_traceNormals[1][trace_offset+i]
2300  != normals[1][i])
2301  {
2302  fluxN_R[i] = -fluxN_R[i];
2303  }
2304  }
2305 
2306  // Subtract to the Riemann flux the discontinuous
2307  // flux
2308  Vmath::Vsub(nEdgePts,
2309  &fluxN_R[0], 1,
2310  &fluxN_D[0], 1, &fluxJumps[0], 1);
2311 
2312  // Multiplicate the flux jumps for the correction
2313  // function
2314  for (i = 0; i < nquad0; ++i)
2315  {
2316  for (j = 0; j < nquad1; ++j)
2317  {
2318  cnt = i + j*nquad0;
2319  divCFluxE0[cnt] =
2320  -fluxJumps[i] * m_dGL_xi2[n][j];
2321  }
2322  }
2323  break;
2324  case 1:
2325  // Extract the edge values of transformed flux-x on
2326  // edge e and order them accordingly to the order of
2327  // the trace space
2328  fields[0]->GetExp(n)->GetTracePhysVals(
2329  e, elmtToTrace[n][e],
2330  fluxX1 + phys_offset,
2331  auxArray1 = fluxN_D);
2332 
2333  // Extract the physical Rieamann flux at each edge
2334  Vmath::Vcopy(nEdgePts,
2335  &numericalFlux[trace_offset], 1,
2336  &fluxN[0], 1);
2337 
2338  // Check the ordering of vectors
2339  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2340  StdRegions::eBackwards)
2341  {
2342  Vmath::Reverse(nEdgePts,
2343  auxArray1 = fluxN, 1,
2344  auxArray2 = fluxN, 1);
2345 
2346  Vmath::Reverse(nEdgePts,
2347  auxArray1 = fluxN_D, 1,
2348  auxArray2 = fluxN_D, 1);
2349  }
2350 
2351  // Transform Riemann Fluxes in the standard element
2352  for (i = 0; i < nquad1; ++i)
2353  {
2354  // Multiply Riemann Flux by Q factors
2355  fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
2356  }
2357 
2358  for (i = 0; i < nEdgePts; ++i)
2359  {
2360  if (m_traceNormals[0][trace_offset+i]
2361  != normals[0][i] ||
2362  m_traceNormals[1][trace_offset+i]
2363  != normals[1][i])
2364  {
2365  fluxN_R[i] = -fluxN_R[i];
2366  }
2367  }
2368 
2369  // Subtract to the Riemann flux the discontinuous
2370  // flux
2371  Vmath::Vsub(nEdgePts,
2372  &fluxN_R[0], 1,
2373  &fluxN_D[0], 1, &fluxJumps[0], 1);
2374 
2375  // Multiplicate the flux jumps for the correction
2376  // function
2377  for (i = 0; i < nquad1; ++i)
2378  {
2379  for (j = 0; j < nquad0; ++j)
2380  {
2381  cnt = (nquad0)*i + j;
2382  divCFluxE1[cnt] =
2383  fluxJumps[i] * m_dGR_xi1[n][j];
2384  }
2385  }
2386  break;
2387  case 2:
2388 
2389  // Extract the edge values of transformed flux-y on
2390  // edge e and order them accordingly to the order of
2391  // the trace space
2392 
2393  fields[0]->GetExp(n)->GetTracePhysVals(
2394  e, elmtToTrace[n][e],
2395  fluxX2 + phys_offset,
2396  auxArray1 = fluxN_D);
2397 
2398  // Extract the physical Rieamann flux at each edge
2399  Vmath::Vcopy(nEdgePts,
2400  &numericalFlux[trace_offset], 1,
2401  &fluxN[0], 1);
2402 
2403  // Check the ordering of vectors
2404  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2405  StdRegions::eBackwards)
2406  {
2407  Vmath::Reverse(nEdgePts,
2408  auxArray1 = fluxN, 1,
2409  auxArray2 = fluxN, 1);
2410 
2411  Vmath::Reverse(nEdgePts,
2412  auxArray1 = fluxN_D, 1,
2413  auxArray2 = fluxN_D, 1);
2414  }
2415 
2416  // Transform Riemann Fluxes in the standard element
2417  for (i = 0; i < nquad0; ++i)
2418  {
2419  // Multiply Riemann Flux by Q factors
2420  fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
2421  }
2422 
2423  for (i = 0; i < nEdgePts; ++i)
2424  {
2425  if (m_traceNormals[0][trace_offset+i]
2426  != normals[0][i] ||
2427  m_traceNormals[1][trace_offset+i]
2428  != normals[1][i])
2429  {
2430  fluxN_R[i] = -fluxN_R[i];
2431  }
2432  }
2433 
2434  // Subtract to the Riemann flux the discontinuous
2435  // flux
2436 
2437  Vmath::Vsub(nEdgePts,
2438  &fluxN_R[0], 1,
2439  &fluxN_D[0], 1, &fluxJumps[0], 1);
2440 
2441  // Multiplicate the flux jumps for the correction
2442  // function
2443  for (i = 0; i < nquad0; ++i)
2444  {
2445  for (j = 0; j < nquad1; ++j)
2446  {
2447  cnt = j*nquad0 + i;
2448  divCFluxE2[cnt] =
2449  fluxJumps[i] * m_dGR_xi2[n][j];
2450  }
2451  }
2452  break;
2453  case 3:
2454  // Extract the edge values of transformed flux-x on
2455  // edge e and order them accordingly to the order of
2456  // the trace space
2457 
2458  fields[0]->GetExp(n)->GetTracePhysVals(
2459  e, elmtToTrace[n][e],
2460  fluxX1 + phys_offset,
2461  auxArray1 = fluxN_D);
2462  Vmath::Neg (nEdgePts, fluxN_D, 1);
2463 
2464  // Extract the physical Rieamann flux at each edge
2465  Vmath::Vcopy(nEdgePts,
2466  &numericalFlux[trace_offset], 1,
2467  &fluxN[0], 1);
2468 
2469  // Check the ordering of vectors
2470  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2471  StdRegions::eBackwards)
2472  {
2473  Vmath::Reverse(nEdgePts,
2474  auxArray1 = fluxN, 1,
2475  auxArray2 = fluxN, 1);
2476 
2477  Vmath::Reverse(nEdgePts,
2478  auxArray1 = fluxN_D, 1,
2479  auxArray2 = fluxN_D, 1);
2480  }
2481 
2482  // Transform Riemann Fluxes in the standard element
2483  for (i = 0; i < nquad1; ++i)
2484  {
2485  // Multiply Riemann Flux by Q factors
2486  fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
2487  }
2488 
2489  for (i = 0; i < nEdgePts; ++i)
2490  {
2491  if (m_traceNormals[0][trace_offset+i]
2492  != normals[0][i] ||
2493  m_traceNormals[1][trace_offset+i]
2494  != normals[1][i])
2495  {
2496  fluxN_R[i] = -fluxN_R[i];
2497  }
2498  }
2499 
2500  // Subtract to the Riemann flux the discontinuous
2501  // flux
2502 
2503  Vmath::Vsub(nEdgePts,
2504  &fluxN_R[0], 1,
2505  &fluxN_D[0], 1, &fluxJumps[0], 1);
2506 
2507  // Multiplicate the flux jumps for the correction
2508  // function
2509  for (i = 0; i < nquad1; ++i)
2510  {
2511  for (j = 0; j < nquad0; ++j)
2512  {
2513  cnt = j + i*nquad0;
2514  divCFluxE3[cnt] =
2515  -fluxJumps[i] * m_dGL_xi1[n][j];
2516  }
2517  }
2518  break;
2519  default:
2520  ASSERTL0(false,"edge value (< 3) is out of range");
2521  break;
2522  }
2523  }
2524 
2525 
2526  // Sum each edge contribution
2527  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
2528  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
2529 
2530  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2531  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2532 
2533  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2534  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2535  }
2536  }
2537 
2538  }//end of namespace SolverUtils
2539 }//end of namespace Nektar
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::LibUtilities::Basis
Represents a basis of a given type.
Definition: Basis.h:212
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:200
Nektar::LocalRegions::Expansion::GetMetricInfo
const SpatialDomains::GeomFactorsSharedPtr & GetMetricInfo() const
Definition: Expansion.cpp:251
Nektar::SolverUtils::Diffusion::m_fluxVectorNS
DiffusionFluxVecCBNS m_fluxVectorNS
Definition: Diffusion.h:467
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_traceVel
Array< OneD, Array< OneD, NekDouble > > m_traceVel
Definition: DiffusionLFRNS.h:74
Nektar::SolverUtils::DiffusionLFRNS::v_DivCFlux_2D_Gauss
virtual void v_DivCFlux_2D_Gauss(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre".
Definition: DiffusionLFRNS.cpp:2194
Nektar::SolverUtils::DiffusionLFRNS::m_divFD
Array< OneD, Array< OneD, NekDouble > > m_divFD
Definition: DiffusionLFRNS.h:96
Nektar::SolverUtils::DiffusionLFRNS::m_jac
Array< OneD, NekDouble > m_jac
Definition: DiffusionLFRNS.h:54
Nektar::SolverUtils::DiffusionLFRNS::v_DerCFlux_2D
virtual void v_DerCFlux_2D(const int nConvectiveFields, const int direction, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux wrt a given coordinate for 2D problems.
Definition: DiffusionLFRNS.cpp:1776
Nektar::SolverUtils::DiffusionLFRNS::m_viscFlux
Array< OneD, Array< OneD, NekDouble > > m_viscFlux
Definition: DiffusionLFRNS.h:94
Nektar::SolverUtils::DiffusionLFRNS::m_tmp2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Definition: DiffusionLFRNS.h:100
Nektar::SolverUtils::DiffusionLFRNS::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFRNS.h:58
Nektar::SolverUtils::DiffusionLFRNS::m_DFC2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Definition: DiffusionLFRNS.h:95
Nektar::SolverUtils::DiffusionLFRNS::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFRNS.h:65
Nektar::SolverUtils::DiffusionLFRNS::m_divFC
Array< OneD, Array< OneD, NekDouble > > m_divFC
Definition: DiffusionLFRNS.h:97
Nektar::SolverUtils::DiffusionLFRNS::m_DFC1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Definition: DiffusionLFRNS.h:89
Nektar::SolverUtils::DiffusionLFRNS::m_dGR_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition: DiffusionLFRNS.h:67
Nektar::SolverUtils::DiffusionLFRNS::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: DiffusionLFRNS.h:55
Nektar::SolverUtils::DiffusionLFRNS::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::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
Nektar::SolverUtils::DiffusionLFRNS::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFRNS.h:60
Nektar::SolverUtils::DiffusionLFRNS::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFRNS.h:75
Nektar::SolverUtils::DiffusionLFRNS::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFRNS.h:59
Nektar::SolverUtils::DiffusionLFRNS::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFRNS.h:64
Nektar::SolverUtils::DiffusionLFRNS::v_DivCFlux_2D
virtual void v_DivCFlux_2D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems.
Definition: DiffusionLFRNS.cpp:1999
Nektar::SolverUtils::DiffusionLFRNS::m_viscTensor
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_viscTensor
Definition: DiffusionLFRNS.h:93
Nektar::SolverUtils::DiffusionLFRNS::m_tmp1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Definition: DiffusionLFRNS.h:99
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::v_Diffuse
virtual void v_Diffuse(const std::size_t nConvective, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd)
Calculate FR Diffusion for the Navier-Stokes (NS) equations using an LDG interface flux.
Definition: DiffusionLFRNS.cpp:906
Nektar::SolverUtils::DiffusionLFRNS::m_dGL_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition: DiffusionLFRNS.h:66
Nektar::SolverUtils::DiffusionLFRNS::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFRNS.h:57
Nektar::SolverUtils::DiffusionLFRNS::m_IF1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
Definition: DiffusionLFRNS.h:87
Nektar::SolverUtils::DiffusionLFRNS::m_D1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
Definition: DiffusionLFRNS.h:91
Nektar::SolverUtils::DiffusionLFRNS::DiffusionLFRNS
DiffusionLFRNS(std::string diffType)
DiffusionLFRNS uses the Flux Reconstruction (FR) approach to compute the diffusion term....
Definition: DiffusionLFRNS.cpp:69
Nektar::SolverUtils::DiffusionLFRNS::m_session
LibUtilities::SessionReaderSharedPtr m_session
Definition: DiffusionLFRNS.h:76
Nektar::SolverUtils::DiffusionLFRNS::v_SetupMetrics
virtual void v_SetupMetrics(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform...
Definition: DiffusionLFRNS.cpp:224
Nektar::SolverUtils::DiffusionLFRNS::m_homoDerivs
Array< OneD, Array< OneD, NekDouble > > m_homoDerivs
Definition: DiffusionLFRNS.h:102
Nektar::SolverUtils::DiffusionLFRNS::v_SetupCFunctions
virtual void v_SetupCFunctions(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–72...
Definition: DiffusionLFRNS.cpp:376
Nektar::SolverUtils::DiffusionLFRNS::m_BD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
Definition: DiffusionLFRNS.h:90
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_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFRNS.h:63
Nektar::SolverUtils::DiffusionLFRNS::v_WeakPenaltyO2
virtual void v_WeakPenaltyO2(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux)
Imposes appropriate bcs for the 2nd order derivatives.
Definition: DiffusionLFRNS.cpp:1583
Nektar::SolverUtils::DiffusionLFRNS::create
static DiffusionSharedPtr create(std::string diffType)
Definition: DiffusionLFRNS.h:47
Nektar::SolverUtils::DiffusionLFRNS::m_DU1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Definition: DiffusionLFRNS.h:88
Nektar::SolverUtils::DiffusionLFRNS::m_DD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
Definition: DiffusionLFRNS.h:92
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::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:120
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
Nektar::LibUtilities::eGaussGaussLegendre
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:48
Nektar::SolverUtils::GetDiffusionFactory
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:41
Nektar::SpatialDomains::eDeformed
@ eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:56
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Polylib::jacobd
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition: Polylib.cpp:1134
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:192
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:461
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:513
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:322
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
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:257
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:436
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1226
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
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:372
Generated on Wed Nov 24 2021 10:33:26 for Nektar++ by doxygen 1.9.1