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