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