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