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