Nektar++
AdvectionFR.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: AdvectionFR.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: FR advection class.
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
37 #include <LibUtilities/Foundations/Basis.h>
38 #include <LibUtilities/Foundations/ManagerAccess.h>
39 #include <LibUtilities/Foundations/Points.h>
40 #include <LibUtilities/Polylib/Polylib.h>
41 #include <LocalRegions/Expansion1D.h>
42 #include <LocalRegions/Expansion2D.h>
43 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
44 #include <MultiRegions/DisContField.h>
45 #include <SolverUtils/Advection/AdvectionFR.h>
46 #include <StdRegions/StdSegExp.h>
47 #include <boost/math/special_functions/gamma.hpp>
48 
49 #include <iomanip>
50 #include <iostream>
51 
52 using namespace std;
53 
54 namespace Nektar
55 {
56 namespace SolverUtils
57 {
58 std::string AdvectionFR::type[] = {
59  GetAdvectionFactory().RegisterCreatorFunction("FRDG", AdvectionFR::create),
60  GetAdvectionFactory().RegisterCreatorFunction("FRSD", AdvectionFR::create),
61  GetAdvectionFactory().RegisterCreatorFunction("FRHU", AdvectionFR::create),
62  GetAdvectionFactory().RegisterCreatorFunction("FRcmin",
63  AdvectionFR::create),
64  GetAdvectionFactory().RegisterCreatorFunction("FRcinf",
65  AdvectionFR::create)};
66 
67 /**
68  * @brief AdvectionFR uses the Flux Reconstruction (FR) approach to
69  * compute the advection term. The implementation is only for segments,
70  * quadrilaterals and hexahedra at the moment.
71  *
72  * \todo Extension to triangles, tetrahedra and other shapes.
73  * (Long term objective)
74  */
75 AdvectionFR::AdvectionFR(std::string advType) : m_advType(advType)
76 {
77 }
78 
79 /**
80  * @brief Initiliase AdvectionFR objects and store them before starting
81  * the time-stepping.
82  *
83  * This routine calls the functions #SetupMetrics and
84  * #SetupCFunctions to initialise the objects needed by AdvectionFR.
85  *
86  * @param pSession Pointer to session reader.
87  * @param pFields Pointer to fields.
88  */
89 void AdvectionFR::v_InitObject(
90  LibUtilities::SessionReaderSharedPtr pSession,
91  Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
92 {
93  Advection::v_InitObject(pSession, pFields);
94  SetupMetrics(pSession, pFields);
95  SetupCFunctions(pSession, pFields);
96 }
97 
98 /**
99  * @brief Setup the metric terms to compute the contravariant
100  * fluxes. (i.e. this special metric terms transform the fluxes
101  * at the interfaces of each element from the physical space to
102  * the standard space).
103  *
104  * This routine calls the function #GetEdgeQFactors to compute and
105  * store the metric factors following an anticlockwise conventions
106  * along the edges/faces of each element. Note: for 1D problem
107  * the transformation is not needed.
108  *
109  * @param pSession Pointer to session reader.
110  * @param pFields Pointer to fields.
111  *
112  * \todo Add the metric terms for 3D Hexahedra.
113  */
114 void AdvectionFR::SetupMetrics(
115  LibUtilities::SessionReaderSharedPtr pSession,
116  Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
117 {
118  boost::ignore_unused(pSession);
119  int i, n;
120  int nquad0, nquad1;
121  int phys_offset;
122  int nLocalSolutionPts;
123  int nElements = pFields[0]->GetExpSize();
124  int nDimensions = pFields[0]->GetCoordim(0);
125  int nSolutionPts = pFields[0]->GetTotPoints();
126  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
127 
128  Array<TwoD, const NekDouble> gmat;
129  Array<OneD, const NekDouble> jac;
130 
131  m_jac = Array<OneD, NekDouble>(nSolutionPts);
132 
133  Array<OneD, NekDouble> auxArray1;
134  Array<OneD, LibUtilities::BasisSharedPtr> base;
135  LibUtilities::PointsKeyVector ptsKeys;
136 
137  switch (nDimensions)
138  {
139  case 1:
140  {
141  for (n = 0; n < nElements; ++n)
142  {
143  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
144  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
145  phys_offset = pFields[0]->GetPhys_Offset(n);
146  jac = pFields[0]
147  ->GetExp(n)
148  ->as<LocalRegions::Expansion1D>()
149  ->GetGeom1D()
150  ->GetMetricInfo()
151  ->GetJac(ptsKeys);
152  for (i = 0; i < nLocalSolutionPts; ++i)
153  {
154  m_jac[i + phys_offset] = jac[0];
155  }
156  }
157  break;
158  }
159  case 2:
160  {
161  m_gmat = Array<OneD, Array<OneD, NekDouble>>(4);
162  m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
163  m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
164  m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
165  m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
166 
167  m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble>>(nElements);
168  m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
169  m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
170  m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
171 
172  LibUtilities::PointsKeyVector ptsKeys;
173 
174  for (n = 0; n < nElements; ++n)
175  {
176  base = pFields[0]->GetExp(n)->GetBase();
177  nquad0 = base[0]->GetNumPoints();
178  nquad1 = base[1]->GetNumPoints();
179 
180  m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
181  m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
182  m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
183  m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
184 
185  // Extract the Q factors at each edge point
186  pFields[0]->GetExp(n)->GetTraceQFactors(0, auxArray1 =
187  m_Q2D_e0[n]);
188  pFields[0]->GetExp(n)->GetTraceQFactors(1, auxArray1 =
189  m_Q2D_e1[n]);
190  pFields[0]->GetExp(n)->GetTraceQFactors(2, auxArray1 =
191  m_Q2D_e2[n]);
192  pFields[0]->GetExp(n)->GetTraceQFactors(3, auxArray1 =
193  m_Q2D_e3[n]);
194 
195  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
196  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
197  phys_offset = pFields[0]->GetPhys_Offset(n);
198 
199  jac = pFields[0]
200  ->GetExp(n)
201  ->as<LocalRegions::Expansion2D>()
202  ->GetGeom2D()
203  ->GetMetricInfo()
204  ->GetJac(ptsKeys);
205  gmat = pFields[0]
206  ->GetExp(n)
207  ->as<LocalRegions::Expansion2D>()
208  ->GetGeom2D()
209  ->GetMetricInfo()
210  ->GetDerivFactors(ptsKeys);
211 
212  if (pFields[0]
213  ->GetExp(n)
214  ->as<LocalRegions::Expansion2D>()
215  ->GetGeom2D()
216  ->GetMetricInfo()
217  ->GetGtype() == SpatialDomains::eDeformed)
218  {
219  for (i = 0; i < nLocalSolutionPts; ++i)
220  {
221  m_jac[i + phys_offset] = jac[i];
222  m_gmat[0][i + phys_offset] = gmat[0][i];
223  m_gmat[1][i + phys_offset] = gmat[1][i];
224  m_gmat[2][i + phys_offset] = gmat[2][i];
225  m_gmat[3][i + phys_offset] = gmat[3][i];
226  }
227  }
228  else
229  {
230  for (i = 0; i < nLocalSolutionPts; ++i)
231  {
232  m_jac[i + phys_offset] = jac[0];
233  m_gmat[0][i + phys_offset] = gmat[0][0];
234  m_gmat[1][i + phys_offset] = gmat[1][0];
235  m_gmat[2][i + phys_offset] = gmat[2][0];
236  m_gmat[3][i + phys_offset] = gmat[3][0];
237  }
238  }
239  }
240 
241  m_traceNormals = Array<OneD, Array<OneD, NekDouble>>(nDimensions);
242  for (i = 0; i < nDimensions; ++i)
243  {
244  m_traceNormals[i] = Array<OneD, NekDouble>(nTracePts);
245  }
246  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
247 
248  break;
249  }
250  case 3:
251  {
252  ASSERTL0(false, "3DFR Metric terms not implemented yet");
253  break;
254  }
255  default:
256  {
257  ASSERTL0(false, "Expansion dimension not recognised");
258  break;
259  }
260  }
261 }
262 
263 /**
264  * @brief Setup the derivatives of the correction functions. For more
265  * details see J Sci Comput (2011) 47: 50–72
266  *
267  * This routine calls 3 different bases:
268  * #eDG_DG_Left - #eDG_DG_Left which recovers a nodal DG scheme,
269  * #eDG_SD_Left - #eDG_SD_Left which recovers the SD scheme,
270  * #eDG_HU_Left - #eDG_HU_Left which recovers the Huynh scheme.
271  * The values of the derivatives of the correction function are then
272  * stored into global variables and reused into the functions
273  * #DivCFlux_1D, #DivCFlux_2D, #DivCFlux_3D to compute the
274  * the divergence of the correction flux for 1D, 2D or 3D problems.
275  *
276  * @param pSession Pointer to session reader.
277  * @param pFields Pointer to fields.
278  */
279 void AdvectionFR::SetupCFunctions(
280  LibUtilities::SessionReaderSharedPtr pSession,
281  Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
282 {
283  boost::ignore_unused(pSession);
284 
285  int i, n;
286  NekDouble c0 = 0.0;
287  NekDouble c1 = 0.0;
288  NekDouble c2 = 0.0;
289  int nquad0, nquad1, nquad2;
290  int nmodes0, nmodes1, nmodes2;
291  Array<OneD, LibUtilities::BasisSharedPtr> base;
292 
293  int nElements = pFields[0]->GetExpSize();
294  int nDimensions = pFields[0]->GetCoordim(0);
295 
296  switch (nDimensions)
297  {
298  case 1:
299  {
300  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
301  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
302 
303  for (n = 0; n < nElements; ++n)
304  {
305  base = pFields[0]->GetExp(n)->GetBase();
306  nquad0 = base[0]->GetNumPoints();
307  nmodes0 = base[0]->GetNumModes();
308  Array<OneD, const NekDouble> z0;
309  Array<OneD, const NekDouble> w0;
310 
311  base[0]->GetZW(z0, w0);
312 
313  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
314  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
315 
316  // Auxiliary vectors to build up the auxiliary
317  // derivatives of the Legendre polynomials
318  Array<OneD, NekDouble> dLp0(nquad0, 0.0);
319  Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
320  Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
321 
322  // Degree of the correction functions
323  int p0 = nmodes0 - 1;
324 
325  // Function sign to compute left correction function
326  NekDouble sign0 = pow(-1.0, p0);
327 
328  // Factorial factor to build the scheme
329  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1) /
330  (pow(2.0, p0) * boost::math::tgamma(p0 + 1) *
331  boost::math::tgamma(p0 + 1));
332 
333  // Scalar parameter which recovers the FR schemes
334  if (m_advType == "FRDG")
335  {
336  c0 = 0.0;
337  }
338  else if (m_advType == "FRSD")
339  {
340  c0 = 2.0 * p0 /
341  ((2.0 * p0 + 1.0) * (p0 + 1.0) *
342  (ap0 * boost::math::tgamma(p0 + 1)) *
343  (ap0 * boost::math::tgamma(p0 + 1)));
344  }
345  else if (m_advType == "FRHU")
346  {
347  c0 = 2.0 * (p0 + 1.0) /
348  ((2.0 * p0 + 1.0) * p0 *
349  (ap0 * boost::math::tgamma(p0 + 1)) *
350  (ap0 * boost::math::tgamma(p0 + 1)));
351  }
352  else if (m_advType == "FRcmin")
353  {
354  c0 = -2.0 / ((2.0 * p0 + 1.0) *
355  (ap0 * boost::math::tgamma(p0 + 1)) *
356  (ap0 * boost::math::tgamma(p0 + 1)));
357  }
358  else if (m_advType == "FRcinf")
359  {
360  c0 = 10000000000000000.0;
361  }
362 
363  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
364  (ap0 * boost::math::tgamma(p0 + 1)) *
365  (ap0 * boost::math::tgamma(p0 + 1));
366 
367  NekDouble overeta0 = 1.0 / (1.0 + etap0);
368 
369  // Derivative of the Legendre polynomials
370  // dLp = derivative of the Legendre polynomial of order p
371  // dLpp = derivative of the Legendre polynomial of order p+1
372  // dLpm = derivative of the Legendre polynomial of order p-1
373  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
374  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
375  0.0);
376  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
377  0.0);
378 
379  // Building the DG_c_Left
380  for (i = 0; i < nquad0; ++i)
381  {
382  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
383  m_dGL_xi1[n][i] += dLpp0[i];
384  m_dGL_xi1[n][i] *= overeta0;
385  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
386  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
387  }
388 
389  // Building the DG_c_Right
390  for (i = 0; i < nquad0; ++i)
391  {
392  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
393  m_dGR_xi1[n][i] += dLpp0[i];
394  m_dGR_xi1[n][i] *= overeta0;
395  m_dGR_xi1[n][i] += dLp0[i];
396  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
397  }
398  }
399  break;
400  }
401  case 2:
402  {
403  LibUtilities::BasisSharedPtr BasisFR_Left0;
404  LibUtilities::BasisSharedPtr BasisFR_Right0;
405  LibUtilities::BasisSharedPtr BasisFR_Left1;
406  LibUtilities::BasisSharedPtr BasisFR_Right1;
407 
408  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
409  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
410  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
411  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
412 
413  for (n = 0; n < nElements; ++n)
414  {
415  base = pFields[0]->GetExp(n)->GetBase();
416  nquad0 = base[0]->GetNumPoints();
417  nquad1 = base[1]->GetNumPoints();
418  nmodes0 = base[0]->GetNumModes();
419  nmodes1 = base[1]->GetNumModes();
420 
421  Array<OneD, const NekDouble> z0;
422  Array<OneD, const NekDouble> w0;
423  Array<OneD, const NekDouble> z1;
424  Array<OneD, const NekDouble> w1;
425 
426  base[0]->GetZW(z0, w0);
427  base[1]->GetZW(z1, w1);
428 
429  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
430  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
431  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
432  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
433 
434  // Auxiliary vectors to build up the auxiliary
435  // derivatives of the Legendre polynomials
436  Array<OneD, NekDouble> dLp0(nquad0, 0.0);
437  Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
438  Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
439  Array<OneD, NekDouble> dLp1(nquad1, 0.0);
440  Array<OneD, NekDouble> dLpp1(nquad1, 0.0);
441  Array<OneD, NekDouble> dLpm1(nquad1, 0.0);
442 
443  // Degree of the correction functions
444  int p0 = nmodes0 - 1;
445  int p1 = nmodes1 - 1;
446 
447  // Function sign to compute left correction function
448  NekDouble sign0 = pow(-1.0, p0);
449  NekDouble sign1 = pow(-1.0, p1);
450 
451  // Factorial factor to build the scheme
452  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1) /
453  (pow(2.0, p0) * boost::math::tgamma(p0 + 1) *
454  boost::math::tgamma(p0 + 1));
455 
456  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1) /
457  (pow(2.0, p1) * boost::math::tgamma(p1 + 1) *
458  boost::math::tgamma(p1 + 1));
459 
460  // Scalar parameter which recovers the FR schemes
461  if (m_advType == "FRDG")
462  {
463  c0 = 0.0;
464  c1 = 0.0;
465  }
466  else if (m_advType == "FRSD")
467  {
468  c0 = 2.0 * p0 /
469  ((2.0 * p0 + 1.0) * (p0 + 1.0) *
470  (ap0 * boost::math::tgamma(p0 + 1)) *
471  (ap0 * boost::math::tgamma(p0 + 1)));
472 
473  c1 = 2.0 * p1 /
474  ((2.0 * p1 + 1.0) * (p1 + 1.0) *
475  (ap1 * boost::math::tgamma(p1 + 1)) *
476  (ap1 * boost::math::tgamma(p1 + 1)));
477  }
478  else if (m_advType == "FRHU")
479  {
480  c0 = 2.0 * (p0 + 1.0) /
481  ((2.0 * p0 + 1.0) * p0 *
482  (ap0 * boost::math::tgamma(p0 + 1)) *
483  (ap0 * boost::math::tgamma(p0 + 1)));
484 
485  c1 = 2.0 * (p1 + 1.0) /
486  ((2.0 * p1 + 1.0) * p1 *
487  (ap1 * boost::math::tgamma(p1 + 1)) *
488  (ap1 * boost::math::tgamma(p1 + 1)));
489  }
490  else if (m_advType == "FRcmin")
491  {
492  c0 = -2.0 / ((2.0 * p0 + 1.0) *
493  (ap0 * boost::math::tgamma(p0 + 1)) *
494  (ap0 * boost::math::tgamma(p0 + 1)));
495 
496  c1 = -2.0 / ((2.0 * p1 + 1.0) *
497  (ap1 * boost::math::tgamma(p1 + 1)) *
498  (ap1 * boost::math::tgamma(p1 + 1)));
499  }
500  else if (m_advType == "FRcinf")
501  {
502  c0 = 10000000000000000.0;
503  c1 = 10000000000000000.0;
504  }
505 
506  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
507  (ap0 * boost::math::tgamma(p0 + 1)) *
508  (ap0 * boost::math::tgamma(p0 + 1));
509 
510  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
511  (ap1 * boost::math::tgamma(p1 + 1)) *
512  (ap1 * boost::math::tgamma(p1 + 1));
513 
514  NekDouble overeta0 = 1.0 / (1.0 + etap0);
515  NekDouble overeta1 = 1.0 / (1.0 + etap1);
516 
517  // Derivative of the Legendre polynomials
518  // dLp = derivative of the Legendre polynomial of order p
519  // dLpp = derivative of the Legendre polynomial of order p+1
520  // dLpm = derivative of the Legendre polynomial of order p-1
521  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
522  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
523  0.0);
524  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
525  0.0);
526 
527  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
528  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1 + 1, 0.0,
529  0.0);
530  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1 - 1, 0.0,
531  0.0);
532 
533  // Building the DG_c_Left
534  for (i = 0; i < nquad0; ++i)
535  {
536  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
537  m_dGL_xi1[n][i] += dLpp0[i];
538  m_dGL_xi1[n][i] *= overeta0;
539  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
540  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
541  }
542 
543  // Building the DG_c_Left
544  for (i = 0; i < nquad1; ++i)
545  {
546  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
547  m_dGL_xi2[n][i] += dLpp1[i];
548  m_dGL_xi2[n][i] *= overeta1;
549  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
550  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
551  }
552 
553  // Building the DG_c_Right
554  for (i = 0; i < nquad0; ++i)
555  {
556  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
557  m_dGR_xi1[n][i] += dLpp0[i];
558  m_dGR_xi1[n][i] *= overeta0;
559  m_dGR_xi1[n][i] += dLp0[i];
560  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
561  }
562 
563  // Building the DG_c_Right
564  for (i = 0; i < nquad1; ++i)
565  {
566  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
567  m_dGR_xi2[n][i] += dLpp1[i];
568  m_dGR_xi2[n][i] *= overeta1;
569  m_dGR_xi2[n][i] += dLp1[i];
570  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
571  }
572  }
573  break;
574  }
575  case 3:
576  {
577  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
578  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
579  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
580  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
581  m_dGL_xi3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
582  m_dGR_xi3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
583 
584  for (n = 0; n < nElements; ++n)
585  {
586  base = pFields[0]->GetExp(n)->GetBase();
587  nquad0 = base[0]->GetNumPoints();
588  nquad1 = base[1]->GetNumPoints();
589  nquad2 = base[2]->GetNumPoints();
590  nmodes0 = base[0]->GetNumModes();
591  nmodes1 = base[1]->GetNumModes();
592  nmodes2 = base[2]->GetNumModes();
593 
594  Array<OneD, const NekDouble> z0;
595  Array<OneD, const NekDouble> w0;
596  Array<OneD, const NekDouble> z1;
597  Array<OneD, const NekDouble> w1;
598  Array<OneD, const NekDouble> z2;
599  Array<OneD, const NekDouble> w2;
600 
601  base[0]->GetZW(z0, w0);
602  base[1]->GetZW(z1, w1);
603  base[1]->GetZW(z2, w2);
604 
605  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
606  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
607  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
608  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
609  m_dGL_xi3[n] = Array<OneD, NekDouble>(nquad2);
610  m_dGR_xi3[n] = Array<OneD, NekDouble>(nquad2);
611 
612  // Auxiliary vectors to build up the auxiliary
613  // derivatives of the Legendre polynomials
614  Array<OneD, NekDouble> dLp0(nquad0, 0.0);
615  Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
616  Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
617  Array<OneD, NekDouble> dLp1(nquad1, 0.0);
618  Array<OneD, NekDouble> dLpp1(nquad1, 0.0);
619  Array<OneD, NekDouble> dLpm1(nquad1, 0.0);
620  Array<OneD, NekDouble> dLp2(nquad2, 0.0);
621  Array<OneD, NekDouble> dLpp2(nquad2, 0.0);
622  Array<OneD, NekDouble> dLpm2(nquad2, 0.0);
623 
624  // Degree of the correction functions
625  int p0 = nmodes0 - 1;
626  int p1 = nmodes1 - 1;
627  int p2 = nmodes2 - 1;
628 
629  // Function sign to compute left correction function
630  NekDouble sign0 = pow(-1.0, p0);
631  NekDouble sign1 = pow(-1.0, p1);
632 
633  // Factorial factor to build the scheme
634  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1) /
635  (pow(2.0, p0) * boost::math::tgamma(p0 + 1) *
636  boost::math::tgamma(p0 + 1));
637 
638  // Factorial factor to build the scheme
639  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1) /
640  (pow(2.0, p1) * boost::math::tgamma(p1 + 1) *
641  boost::math::tgamma(p1 + 1));
642 
643  // Factorial factor to build the scheme
644  NekDouble ap2 = boost::math::tgamma(2 * p2 + 1) /
645  (pow(2.0, p2) * boost::math::tgamma(p2 + 1) *
646  boost::math::tgamma(p2 + 1));
647 
648  // Scalar parameter which recovers the FR schemes
649  if (m_advType == "FRDG")
650  {
651  c0 = 0.0;
652  c1 = 0.0;
653  c2 = 0.0;
654  }
655  else if (m_advType == "FRSD")
656  {
657  c0 = 2.0 * p0 /
658  ((2.0 * p0 + 1.0) * (p0 + 1.0) *
659  (ap0 * boost::math::tgamma(p0 + 1)) *
660  (ap0 * boost::math::tgamma(p0 + 1)));
661 
662  c1 = 2.0 * p1 /
663  ((2.0 * p1 + 1.0) * (p1 + 1.0) *
664  (ap1 * boost::math::tgamma(p1 + 1)) *
665  (ap1 * boost::math::tgamma(p1 + 1)));
666 
667  c2 = 2.0 * p2 /
668  ((2.0 * p2 + 1.0) * (p2 + 1.0) *
669  (ap2 * boost::math::tgamma(p2 + 1)) *
670  (ap2 * boost::math::tgamma(p2 + 1)));
671  }
672  else if (m_advType == "FRHU")
673  {
674  c0 = 2.0 * (p0 + 1.0) /
675  ((2.0 * p0 + 1.0) * p0 *
676  (ap0 * boost::math::tgamma(p0 + 1)) *
677  (ap0 * boost::math::tgamma(p0 + 1)));
678 
679  c1 = 2.0 * (p1 + 1.0) /
680  ((2.0 * p1 + 1.0) * p1 *
681  (ap1 * boost::math::tgamma(p1 + 1)) *
682  (ap1 * boost::math::tgamma(p1 + 1)));
683 
684  c2 = 2.0 * (p2 + 1.0) /
685  ((2.0 * p2 + 1.0) * p2 *
686  (ap2 * boost::math::tgamma(p2 + 1)) *
687  (ap2 * boost::math::tgamma(p2 + 1)));
688  }
689  else if (m_advType == "FRcmin")
690  {
691  c0 = -2.0 / ((2.0 * p0 + 1.0) *
692  (ap0 * boost::math::tgamma(p0 + 1)) *
693  (ap0 * boost::math::tgamma(p0 + 1)));
694 
695  c1 = -2.0 / ((2.0 * p1 + 1.0) *
696  (ap1 * boost::math::tgamma(p1 + 1)) *
697  (ap1 * boost::math::tgamma(p1 + 1)));
698 
699  c2 = -2.0 / ((2.0 * p2 + 1.0) *
700  (ap2 * boost::math::tgamma(p2 + 1)) *
701  (ap2 * boost::math::tgamma(p2 + 1)));
702  }
703  else if (m_advType == "FRcinf")
704  {
705  c0 = 10000000000000000.0;
706  c1 = 10000000000000000.0;
707  c2 = 10000000000000000.0;
708  }
709 
710  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
711  (ap0 * boost::math::tgamma(p0 + 1)) *
712  (ap0 * boost::math::tgamma(p0 + 1));
713 
714  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
715  (ap1 * boost::math::tgamma(p1 + 1)) *
716  (ap1 * boost::math::tgamma(p1 + 1));
717 
718  NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0) *
719  (ap2 * boost::math::tgamma(p2 + 1)) *
720  (ap2 * boost::math::tgamma(p2 + 1));
721 
722  NekDouble overeta0 = 1.0 / (1.0 + etap0);
723  NekDouble overeta1 = 1.0 / (1.0 + etap1);
724  NekDouble overeta2 = 1.0 / (1.0 + etap2);
725 
726  // Derivative of the Legendre polynomials
727  // dLp = derivative of the Legendre polynomial of order p
728  // dLpp = derivative of the Legendre polynomial of order p+1
729  // dLpm = derivative of the Legendre polynomial of order p-1
730  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
731  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
732  0.0);
733  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
734  0.0);
735 
736  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
737  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1 + 1, 0.0,
738  0.0);
739  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1 - 1, 0.0,
740  0.0);
741 
742  Polylib::jacobd(nquad2, z2.data(), &(dLp2[0]), p2, 0.0, 0.0);
743  Polylib::jacobd(nquad2, z2.data(), &(dLpp2[0]), p2 + 1, 0.0,
744  0.0);
745  Polylib::jacobd(nquad2, z2.data(), &(dLpm2[0]), p2 - 1, 0.0,
746  0.0);
747 
748  // Building the DG_c_Left
749  for (i = 0; i < nquad0; ++i)
750  {
751  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
752  m_dGL_xi1[n][i] += dLpp0[i];
753  m_dGL_xi1[n][i] *= overeta0;
754  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
755  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
756  }
757 
758  // Building the DG_c_Left
759  for (i = 0; i < nquad1; ++i)
760  {
761  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
762  m_dGL_xi2[n][i] += dLpp1[i];
763  m_dGL_xi2[n][i] *= overeta1;
764  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
765  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
766  }
767 
768  // Building the DG_c_Left
769  for (i = 0; i < nquad2; ++i)
770  {
771  m_dGL_xi3[n][i] = etap2 * dLpm2[i];
772  m_dGL_xi3[n][i] += dLpp2[i];
773  m_dGL_xi3[n][i] *= overeta2;
774  m_dGL_xi3[n][i] = dLp2[i] - m_dGL_xi3[n][i];
775  m_dGL_xi3[n][i] = 0.5 * sign1 * m_dGL_xi3[n][i];
776  }
777 
778  // Building the DG_c_Right
779  for (i = 0; i < nquad0; ++i)
780  {
781  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
782  m_dGR_xi1[n][i] += dLpp0[i];
783  m_dGR_xi1[n][i] *= overeta0;
784  m_dGR_xi1[n][i] += dLp0[i];
785  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
786  }
787 
788  // Building the DG_c_Right
789  for (i = 0; i < nquad1; ++i)
790  {
791  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
792  m_dGR_xi2[n][i] += dLpp1[i];
793  m_dGR_xi2[n][i] *= overeta1;
794  m_dGR_xi2[n][i] += dLp1[i];
795  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
796  }
797 
798  // Building the DG_c_Right
799  for (i = 0; i < nquad2; ++i)
800  {
801  m_dGR_xi3[n][i] = etap2 * dLpm2[i];
802  m_dGR_xi3[n][i] += dLpp2[i];
803  m_dGR_xi3[n][i] *= overeta2;
804  m_dGR_xi3[n][i] += dLp2[i];
805  m_dGR_xi3[n][i] = 0.5 * m_dGR_xi3[n][i];
806  }
807  }
808  break;
809  }
810  default:
811  {
812  ASSERTL0(false, "Expansion dimension not recognised");
813  break;
814  }
815  }
816 }
817 
818 /**
819  * @brief Compute the advection term at each time-step using the Flux
820  * Reconstruction approach (FR).
821  *
822  * @param nConvectiveFields Number of fields.
823  * @param fields Pointer to fields.
824  * @param advVel Advection velocities.
825  * @param inarray Solution at the previous time-step.
826  * @param outarray Advection term to be passed at the
827  * time integration class.
828  *
829  */
830 void AdvectionFR::v_Advect(
831  const int nConvectiveFields,
832  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
833  const Array<OneD, Array<OneD, NekDouble>> &advVel,
834  const Array<OneD, Array<OneD, NekDouble>> &inarray,
835  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble &time,
836  const Array<OneD, Array<OneD, NekDouble>> &pFwd,
837  const Array<OneD, Array<OneD, NekDouble>> &pBwd)
838 {
839  boost::ignore_unused(advVel, time, pFwd, pBwd);
840 
841  int i, j, n;
842  int phys_offset;
843 
844  Array<OneD, NekDouble> auxArray1, auxArray2;
845 
846  Array<OneD, LibUtilities::BasisSharedPtr> Basis;
847  Basis = fields[0]->GetExp(0)->GetBase();
848 
849  int nElements = fields[0]->GetExpSize();
850  int nDimensions = fields[0]->GetCoordim(0);
851  int nSolutionPts = fields[0]->GetTotPoints();
852  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
853  int nCoeffs = fields[0]->GetNcoeffs();
854 
855  // Storage for flux vector.
856  Array<OneD, Array<OneD, Array<OneD, NekDouble>>> fluxvector(
857  nConvectiveFields);
858  // Outarray for Galerkin projection in case of primitive dealising
859  Array<OneD, Array<OneD, NekDouble>> outarrayCoeff(nConvectiveFields);
860  // Store forwards/backwards space along trace space
861  Array<OneD, Array<OneD, NekDouble>> Fwd(nConvectiveFields);
862  Array<OneD, Array<OneD, NekDouble>> Bwd(nConvectiveFields);
863  Array<OneD, Array<OneD, NekDouble>> numflux(nConvectiveFields);
864 
865  // Set up storage for flux vector.
866  for (i = 0; i < nConvectiveFields; ++i)
867  {
868  fluxvector[i] = Array<OneD, Array<OneD, NekDouble>>(m_spaceDim);
869 
870  for (j = 0; j < m_spaceDim; ++j)
871  {
872  fluxvector[i][j] = Array<OneD, NekDouble>(nSolutionPts);
873  }
874  }
875 
876  for (i = 0; i < nConvectiveFields; ++i)
877  {
878  outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
879  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
880  Bwd[i] = Array<OneD, NekDouble>(nTracePts);
881  numflux[i] = Array<OneD, NekDouble>(nTracePts);
882  fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
883  }
884 
885  // Computing the interface flux at each trace point
886  m_riemann->Solve(m_spaceDim, Fwd, Bwd, numflux);
887 
888  // Divergence of the flux (computing the RHS)
889  switch (nDimensions)
890  {
891  // 1D problems
892  case 1:
893  {
894  Array<OneD, NekDouble> DfluxvectorX1(nSolutionPts);
895  Array<OneD, NekDouble> divFC(nSolutionPts);
896 
897  // Get the advection flux vector (solver specific)
898  m_fluxVector(inarray, fluxvector);
899 
900  // Get the discontinuous flux divergence
901  for (i = 0; i < nConvectiveFields; ++i)
902  {
903  for (n = 0; n < nElements; n++)
904  {
905  phys_offset = fields[0]->GetPhys_Offset(n);
906 
907  fields[i]->GetExp(n)->PhysDeriv(
908  0, fluxvector[i][0] + phys_offset,
909  auxArray2 = DfluxvectorX1 + phys_offset);
910  }
911 
912  // Get the correction flux divergence
913  DivCFlux_1D(nConvectiveFields, fields, fluxvector[i][0],
914  numflux[i], divFC);
915 
916  // Back to the physical space using global operations
917  Vmath::Vdiv(nSolutionPts, &divFC[0], 1, &m_jac[0], 1,
918  &outarray[i][0], 1);
919 
920  // Adding the total divergence to outarray (RHS)
921  Vmath::Vadd(nSolutionPts, &outarray[i][0], 1, &DfluxvectorX1[0],
922  1, &outarray[i][0], 1);
923 
924  // Primitive Dealiasing 1D
925  if (!(Basis[0]->Collocation()))
926  {
927  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
928  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
929  }
930  }
931  break;
932  }
933  // 2D problems
934  case 2:
935  {
936  Array<OneD, NekDouble> DfluxvectorX1(nSolutionPts);
937  Array<OneD, NekDouble> DfluxvectorX2(nSolutionPts);
938  Array<OneD, NekDouble> divFD(nSolutionPts);
939  Array<OneD, NekDouble> divFC(nSolutionPts);
940 
941  // Get the advection flux vector (solver specific)
942  m_fluxVector(inarray, fluxvector);
943 
944  for (i = 0; i < nConvectiveFields; ++i)
945  {
946  // Temporary vectors
947  Array<OneD, NekDouble> f_hat(nSolutionPts);
948  Array<OneD, NekDouble> g_hat(nSolutionPts);
949 
950  Vmath::Vvtvvtp(nSolutionPts, &m_gmat[0][0], 1,
951  &fluxvector[i][0][0], 1, &m_gmat[2][0], 1,
952  &fluxvector[i][1][0], 1, &f_hat[0], 1);
953 
954  Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &f_hat[0], 1, &f_hat[0],
955  1);
956 
957  Vmath::Vvtvvtp(nSolutionPts, &m_gmat[1][0], 1,
958  &fluxvector[i][0][0], 1, &m_gmat[3][0], 1,
959  &fluxvector[i][1][0], 1, &g_hat[0], 1);
960 
961  Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &g_hat[0], 1, &g_hat[0],
962  1);
963 
964  // Get the discontinuous flux derivatives
965  for (n = 0; n < nElements; n++)
966  {
967  phys_offset = fields[0]->GetPhys_Offset(n);
968  fields[0]->GetExp(n)->StdPhysDeriv(
969  0, auxArray1 = f_hat + phys_offset,
970  auxArray2 = DfluxvectorX1 + phys_offset);
971  fields[0]->GetExp(n)->StdPhysDeriv(
972  1, auxArray1 = g_hat + phys_offset,
973  auxArray2 = DfluxvectorX2 + phys_offset);
974  }
975 
976  // Divergence of the discontinuous flux
977  Vmath::Vadd(nSolutionPts, DfluxvectorX1, 1, DfluxvectorX2, 1,
978  divFD, 1);
979 
980  // Divergence of the correction flux
981  if (Basis[0]->GetPointsType() ==
982  LibUtilities::eGaussGaussLegendre &&
983  Basis[1]->GetPointsType() ==
984  LibUtilities::eGaussGaussLegendre)
985  {
986  DivCFlux_2D_Gauss(nConvectiveFields, fields, f_hat, g_hat,
987  numflux[i], divFC);
988  }
989  else
990  {
991  DivCFlux_2D(nConvectiveFields, fields, fluxvector[i][0],
992  fluxvector[i][1], numflux[i], divFC);
993  }
994 
995  // Divergence of the final flux
996  Vmath::Vadd(nSolutionPts, divFD, 1, divFC, 1, outarray[i], 1);
997 
998  // Back to the physical space using a global operation
999  Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1, &m_jac[0], 1,
1000  &outarray[i][0], 1);
1001 
1002  // Primitive Dealiasing 2D
1003  if (!(Basis[0]->Collocation()))
1004  {
1005  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1006  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1007  }
1008  } // close nConvectiveFields loop
1009  break;
1010  }
1011  // 3D problems
1012  case 3:
1013  {
1014  ASSERTL0(false, "3D FRDG case not implemented yet");
1015  break;
1016  }
1017  }
1018 }
1019 
1020 /**
1021  * @brief Compute the divergence of the corrective flux for 1D problems.
1022  *
1023  * @param nConvectiveFields Number of fields.
1024  * @param fields Pointer to fields.
1025  * @param fluxX1 X1-volumetric flux in physical space.
1026  * @param numericalFlux Interface flux in physical space.
1027  * @param divCFlux Divergence of the corrective flux.
1028  *
1029  */
1030 void AdvectionFR::DivCFlux_1D(
1031  const int nConvectiveFields,
1032  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1033  const Array<OneD, const NekDouble> &fluxX1,
1034  const Array<OneD, const NekDouble> &numericalFlux,
1035  Array<OneD, NekDouble> &divCFlux)
1036 {
1037  boost::ignore_unused(nConvectiveFields);
1038 
1039  int n;
1040  int nLocalSolutionPts, phys_offset, t_offset;
1041 
1042  LibUtilities::BasisSharedPtr Basis;
1043  Basis = fields[0]->GetExp(0)->GetBasis(0);
1044 
1045  int nElements = fields[0]->GetExpSize();
1046  int nSolutionPts = fields[0]->GetTotPoints();
1047 
1048  vector<bool> leftAdjacentTraces =
1049  std::static_pointer_cast<MultiRegions::DisContField>(fields[0])
1050  ->GetLeftAdjacentTraces();
1051 
1052  // Arrays to store the derivatives of the correction flux
1053  Array<OneD, NekDouble> DCL(nSolutionPts / nElements, 0.0);
1054  Array<OneD, NekDouble> DCR(nSolutionPts / nElements, 0.0);
1055 
1056  // Arrays to store the intercell numerical flux jumps
1057  Array<OneD, NekDouble> JumpL(nElements);
1058  Array<OneD, NekDouble> JumpR(nElements);
1059 
1060  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1061  fields[0]->GetTraceMap()->GetElmtToTrace();
1062 
1063  for (n = 0; n < nElements; ++n)
1064  {
1065  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1066  phys_offset = fields[0]->GetPhys_Offset(n);
1067 
1068  Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1069  Array<OneD, NekDouble> tmpFluxVertex(1, 0.0);
1070  Vmath::Vcopy(nLocalSolutionPts, &fluxX1[phys_offset], 1, &tmparrayX1[0],
1071  1);
1072 
1073  fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0], tmparrayX1,
1074  tmpFluxVertex);
1075 
1076  t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1077  elmtToTrace[n][0]->GetElmtId());
1078 
1079  if (leftAdjacentTraces[2 * n])
1080  {
1081  JumpL[n] = -numericalFlux[t_offset] - tmpFluxVertex[0];
1082  }
1083  else
1084  {
1085  JumpL[n] = numericalFlux[t_offset] - tmpFluxVertex[0];
1086  }
1087 
1088  fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1], tmparrayX1,
1089  tmpFluxVertex);
1090 
1091  t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1092  elmtToTrace[n][1]->GetElmtId());
1093 
1094  if (leftAdjacentTraces[2 * n + 1])
1095  {
1096  JumpR[n] = numericalFlux[t_offset] - tmpFluxVertex[0];
1097  }
1098  else
1099  {
1100  JumpR[n] = -numericalFlux[t_offset] - tmpFluxVertex[0];
1101  }
1102  }
1103 
1104  for (n = 0; n < nElements; ++n)
1105  {
1106  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1107  phys_offset = fields[0]->GetPhys_Offset(n);
1108 
1109  // Left jump multiplied by left derivative of C function
1110  Vmath::Smul(nLocalSolutionPts, JumpL[n], &m_dGL_xi1[n][0], 1, &DCL[0],
1111  1);
1112 
1113  // Right jump multiplied by right derivative of C function
1114  Vmath::Smul(nLocalSolutionPts, JumpR[n], &m_dGR_xi1[n][0], 1, &DCR[0],
1115  1);
1116 
1117  // Assembling divergence of the correction flux
1118  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1119  &divCFlux[phys_offset], 1);
1120  }
1121 }
1122 
1123 /**
1124  * @brief Compute the divergence of the corrective flux for 2D problems.
1125  *
1126  * @param nConvectiveFields Number of fields.
1127  * @param fields Pointer to fields.
1128  * @param fluxX1 X1-volumetric flux in physical space.
1129  * @param fluxX2 X2-volumetric flux in physical space.
1130  * @param numericalFlux Interface flux in physical space.
1131  * @param divCFlux Divergence of the corrective flux.
1132  *
1133  * \todo: Switch on shapes eventually here.
1134  */
1135 void AdvectionFR::DivCFlux_2D(
1136  const int nConvectiveFields,
1137  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1138  const Array<OneD, const NekDouble> &fluxX1,
1139  const Array<OneD, const NekDouble> &fluxX2,
1140  const Array<OneD, const NekDouble> &numericalFlux,
1141  Array<OneD, NekDouble> &divCFlux)
1142 {
1143  boost::ignore_unused(nConvectiveFields);
1144 
1145  int n, e, i, j, cnt;
1146 
1147  int nElements = fields[0]->GetExpSize();
1148 
1149  int nLocalSolutionPts, nEdgePts;
1150  int trace_offset, phys_offset;
1151  int nquad0, nquad1;
1152 
1153  Array<OneD, NekDouble> auxArray1;
1154  Array<OneD, LibUtilities::BasisSharedPtr> base;
1155 
1156  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1157  fields[0]->GetTraceMap()->GetElmtToTrace();
1158 
1159  // Loop on the elements
1160  for (n = 0; n < nElements; ++n)
1161  {
1162  // Offset of the element on the global vector
1163  phys_offset = fields[0]->GetPhys_Offset(n);
1164  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1165 
1166  base = fields[0]->GetExp(n)->GetBase();
1167  nquad0 = base[0]->GetNumPoints();
1168  nquad1 = base[1]->GetNumPoints();
1169 
1170  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1171  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1172  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1173  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1174 
1175  // Loop on the edges
1176  for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1177  {
1178  // Number of edge points of edge e
1179  nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1180 
1181  Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1182  Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1183  Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1184  Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1185  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1186 
1187  // Offset of the trace space correspondent to edge e
1188  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1189  elmtToTrace[n][e]->GetElmtId());
1190 
1191  // Get the normals of edge e
1192  const Array<OneD, const Array<OneD, NekDouble>> &normals =
1193  fields[0]->GetExp(n)->GetTraceNormal(e);
1194 
1195  // Extract the edge values of flux-x on edge e and order
1196  // them accordingly to the order of the trace space
1197  fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1198  fluxX1 + phys_offset,
1199  auxArray1 = tmparrayX1);
1200 
1201  // Extract the edge values of flux-y on edge e and order
1202  // them accordingly to the order of the trace space
1203  fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1204  fluxX2 + phys_offset,
1205  auxArray1 = tmparrayX2);
1206 
1207  // Multiply the edge components of the flux by the normal
1208  Vmath::Vvtvvtp(nEdgePts, &tmparrayX1[0], 1,
1209  &m_traceNormals[0][trace_offset], 1, &tmparrayX2[0],
1210  1, &m_traceNormals[1][trace_offset], 1, &fluxN[0],
1211  1);
1212 
1213  // Subtract to the Riemann flux the discontinuous flux
1214  Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1, &fluxN[0], 1,
1215  &fluxJumps[0], 1);
1216 
1217  // Check the ordering of the jump vectors
1218  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1219  StdRegions::eBackwards)
1220  {
1221  Vmath::Reverse(nEdgePts, &fluxJumps[0], 1, &fluxJumps[0], 1);
1222  }
1223 
1224  for (i = 0; i < nEdgePts; ++i)
1225  {
1226  if (m_traceNormals[0][trace_offset + i] != normals[0][i] ||
1227  m_traceNormals[1][trace_offset + i] != normals[1][i])
1228  {
1229  fluxJumps[i] = -fluxJumps[i];
1230  }
1231  }
1232 
1233  // Multiply jumps by derivatives of the correction functions
1234  switch (e)
1235  {
1236  case 0:
1237  for (i = 0; i < nquad0; ++i)
1238  {
1239  // Multiply fluxJumps by Q factors
1240  fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1241 
1242  for (j = 0; j < nquad1; ++j)
1243  {
1244  cnt = i + j * nquad0;
1245  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1246  }
1247  }
1248  break;
1249  case 1:
1250  for (i = 0; i < nquad1; ++i)
1251  {
1252  // Multiply fluxJumps by Q factors
1253  fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1254 
1255  for (j = 0; j < nquad0; ++j)
1256  {
1257  cnt = (nquad0)*i + j;
1258  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1259  }
1260  }
1261  break;
1262  case 2:
1263  for (i = 0; i < nquad0; ++i)
1264  {
1265  // Multiply fluxJumps by Q factors
1266  fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1267 
1268  for (j = 0; j < nquad1; ++j)
1269  {
1270  cnt = j * nquad0 + i;
1271  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1272  }
1273  }
1274  break;
1275  case 3:
1276  for (i = 0; i < nquad1; ++i)
1277  {
1278  // Multiply fluxJumps by Q factors
1279  fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1280  for (j = 0; j < nquad0; ++j)
1281  {
1282  cnt = j + i * nquad0;
1283  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1284  }
1285  }
1286  break;
1287  default:
1288  ASSERTL0(false, "edge value (< 3) is out of range");
1289  break;
1290  }
1291  }
1292 
1293  // Sum each edge contribution
1294  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
1295  &divCFlux[phys_offset], 1);
1296 
1297  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1298  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1299 
1300  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1301  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1302  }
1303 }
1304 
1305 /**
1306  * @brief Compute the divergence of the corrective flux for 2D problems
1307  * where POINTSTYPE="GaussGaussLegendre"
1308  *
1309  * @param nConvectiveFields Number of fields.
1310  * @param fields Pointer to fields.
1311  * @param fluxX1 X1-volumetric flux in physical space.
1312  * @param fluxX2 X2-volumetric flux in physical space.
1313  * @param numericalFlux Interface flux in physical space.
1314  * @param divCFlux Divergence of the corrective flux.
1315  *
1316  * \todo: Switch on shapes eventually here.
1317  */
1318 
1319 void AdvectionFR::DivCFlux_2D_Gauss(
1320  const int nConvectiveFields,
1321  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1322  const Array<OneD, const NekDouble> &fluxX1,
1323  const Array<OneD, const NekDouble> &fluxX2,
1324  const Array<OneD, const NekDouble> &numericalFlux,
1325  Array<OneD, NekDouble> &divCFlux)
1326 {
1327  boost::ignore_unused(nConvectiveFields);
1328 
1329  int n, e, i, j, cnt;
1330 
1331  int nElements = fields[0]->GetExpSize();
1332  int nLocalSolutionPts;
1333  int nEdgePts;
1334  int trace_offset;
1335  int phys_offset;
1336  int nquad0;
1337  int nquad1;
1338 
1339  Array<OneD, NekDouble> auxArray1, auxArray2;
1340  Array<OneD, LibUtilities::BasisSharedPtr> base;
1341 
1342  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1343  fields[0]->GetTraceMap()->GetElmtToTrace();
1344 
1345  // Loop on the elements
1346  for (n = 0; n < nElements; ++n)
1347  {
1348  // Offset of the element on the global vector
1349  phys_offset = fields[0]->GetPhys_Offset(n);
1350  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1351 
1352  base = fields[0]->GetExp(n)->GetBase();
1353  nquad0 = base[0]->GetNumPoints();
1354  nquad1 = base[1]->GetNumPoints();
1355 
1356  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1357  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1358  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1359  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1360 
1361  // Loop on the edges
1362  for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1363  {
1364  // Number of edge points of edge e
1365  nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1366 
1367  Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1368  Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1369  Array<OneD, NekDouble> fluxN_R(nEdgePts, 0.0);
1370  Array<OneD, NekDouble> fluxN_D(nEdgePts, 0.0);
1371  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1372 
1373  // Offset of the trace space correspondent to edge e
1374  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1375  elmtToTrace[n][e]->GetElmtId());
1376 
1377  // Get the normals of edge e
1378  const Array<OneD, const Array<OneD, NekDouble>> &normals =
1379  fields[0]->GetExp(n)->GetTraceNormal(e);
1380 
1381  // Extract the trasformed normal flux at each edge
1382  switch (e)
1383  {
1384  case 0:
1385  // Extract the edge values of transformed flux-y on
1386  // edge e and order them accordingly to the order of
1387  // the trace space
1388  fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1389  fluxX2 + phys_offset,
1390  auxArray1 = fluxN_D);
1391 
1392  Vmath::Neg(nEdgePts, fluxN_D, 1);
1393 
1394  // Extract the physical Rieamann flux at each edge
1395  Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
1396  &fluxN[0], 1);
1397 
1398  // Check the ordering of vectors
1399  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1400  StdRegions::eBackwards)
1401  {
1402  Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
1403  auxArray2 = fluxN, 1);
1404 
1405  Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
1406  auxArray2 = fluxN_D, 1);
1407  }
1408 
1409  // Transform Riemann Fluxes in the standard element
1410  for (i = 0; i < nquad0; ++i)
1411  {
1412  // Multiply Riemann Flux by Q factors
1413  fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
1414  }
1415 
1416  for (i = 0; i < nEdgePts; ++i)
1417  {
1418  if (m_traceNormals[0][trace_offset + i] !=
1419  normals[0][i] ||
1420  m_traceNormals[1][trace_offset + i] !=
1421  normals[1][i])
1422  {
1423  fluxN_R[i] = -fluxN_R[i];
1424  }
1425  }
1426 
1427  // Subtract to the Riemann flux the discontinuous
1428  // flux
1429  Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
1430  &fluxJumps[0], 1);
1431 
1432  // Multiplicate the flux jumps for the correction
1433  // function
1434  for (i = 0; i < nquad0; ++i)
1435  {
1436  for (j = 0; j < nquad1; ++j)
1437  {
1438  cnt = i + j * nquad0;
1439  divCFluxE0[cnt] = -fluxJumps[i] * m_dGL_xi2[n][j];
1440  }
1441  }
1442  break;
1443  case 1:
1444  // Extract the edge values of transformed flux-x on
1445  // edge e and order them accordingly to the order of
1446  // the trace space
1447  fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1448  fluxX1 + phys_offset,
1449  auxArray1 = fluxN_D);
1450 
1451  // Extract the physical Rieamann flux at each edge
1452  Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
1453  &fluxN[0], 1);
1454 
1455  // Check the ordering of vectors
1456  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1457  StdRegions::eBackwards)
1458  {
1459  Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
1460  auxArray2 = fluxN, 1);
1461 
1462  Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
1463  auxArray2 = fluxN_D, 1);
1464  }
1465 
1466  // Transform Riemann Fluxes in the standard element
1467  for (i = 0; i < nquad1; ++i)
1468  {
1469  // Multiply Riemann Flux by Q factors
1470  fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
1471  }
1472 
1473  for (i = 0; i < nEdgePts; ++i)
1474  {
1475  if (m_traceNormals[0][trace_offset + i] !=
1476  normals[0][i] ||
1477  m_traceNormals[1][trace_offset + i] !=
1478  normals[1][i])
1479  {
1480  fluxN_R[i] = -fluxN_R[i];
1481  }
1482  }
1483 
1484  // Subtract to the Riemann flux the discontinuous
1485  // flux
1486  Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
1487  &fluxJumps[0], 1);
1488 
1489  // Multiplicate the flux jumps for the correction
1490  // function
1491  for (i = 0; i < nquad1; ++i)
1492  {
1493  for (j = 0; j < nquad0; ++j)
1494  {
1495  cnt = (nquad0)*i + j;
1496  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1497  }
1498  }
1499  break;
1500  case 2:
1501 
1502  // Extract the edge values of transformed flux-y on
1503  // edge e and order them accordingly to the order of
1504  // the trace space
1505 
1506  fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1507  fluxX2 + phys_offset,
1508  auxArray1 = fluxN_D);
1509 
1510  // Extract the physical Rieamann flux at each edge
1511  Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
1512  &fluxN[0], 1);
1513 
1514  // Check the ordering of vectors
1515  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1516  StdRegions::eBackwards)
1517  {
1518  Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
1519  auxArray2 = fluxN, 1);
1520 
1521  Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
1522  auxArray2 = fluxN_D, 1);
1523  }
1524 
1525  // Transform Riemann Fluxes in the standard element
1526  for (i = 0; i < nquad0; ++i)
1527  {
1528  // Multiply Riemann Flux by Q factors
1529  fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
1530  }
1531 
1532  for (i = 0; i < nEdgePts; ++i)
1533  {
1534  if (m_traceNormals[0][trace_offset + i] !=
1535  normals[0][i] ||
1536  m_traceNormals[1][trace_offset + i] !=
1537  normals[1][i])
1538  {
1539  fluxN_R[i] = -fluxN_R[i];
1540  }
1541  }
1542 
1543  // Subtract to the Riemann flux the discontinuous
1544  // flux
1545 
1546  Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
1547  &fluxJumps[0], 1);
1548 
1549  // Multiplicate the flux jumps for the correction
1550  // function
1551  for (i = 0; i < nquad0; ++i)
1552  {
1553  for (j = 0; j < nquad1; ++j)
1554  {
1555  cnt = j * nquad0 + i;
1556  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1557  }
1558  }
1559  break;
1560  case 3:
1561  // Extract the edge values of transformed flux-x on
1562  // edge e and order them accordingly to the order of
1563  // the trace space
1564 
1565  fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1566  fluxX1 + phys_offset,
1567  auxArray1 = fluxN_D);
1568  Vmath::Neg(nEdgePts, fluxN_D, 1);
1569 
1570  // Extract the physical Rieamann flux at each edge
1571  Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
1572  &fluxN[0], 1);
1573 
1574  // Check the ordering of vectors
1575  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1576  StdRegions::eBackwards)
1577  {
1578  Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
1579  auxArray2 = fluxN, 1);
1580 
1581  Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
1582  auxArray2 = fluxN_D, 1);
1583  }
1584 
1585  // Transform Riemann Fluxes in the standard element
1586  for (i = 0; i < nquad1; ++i)
1587  {
1588  // Multiply Riemann Flux by Q factors
1589  fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
1590  }
1591 
1592  for (i = 0; i < nEdgePts; ++i)
1593  {
1594  if (m_traceNormals[0][trace_offset + i] !=
1595  normals[0][i] ||
1596  m_traceNormals[1][trace_offset + i] !=
1597  normals[1][i])
1598  {
1599  fluxN_R[i] = -fluxN_R[i];
1600  }
1601  }
1602 
1603  // Subtract to the Riemann flux the discontinuous
1604  // flux
1605 
1606  Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
1607  &fluxJumps[0], 1);
1608 
1609  // Multiplicate the flux jumps for the correction
1610  // function
1611  for (i = 0; i < nquad1; ++i)
1612  {
1613  for (j = 0; j < nquad0; ++j)
1614  {
1615  cnt = j + i * nquad0;
1616  divCFluxE3[cnt] = -fluxJumps[i] * m_dGL_xi1[n][j];
1617  }
1618  }
1619  break;
1620  default:
1621  ASSERTL0(false, "edge value (< 3) is out of range");
1622  break;
1623  }
1624  }
1625 
1626  // Sum each edge contribution
1627  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
1628  &divCFlux[phys_offset], 1);
1629 
1630  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1631  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1632 
1633  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1634  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1635  }
1636 }
1637 
1638 /**
1639  * @brief Compute the divergence of the corrective flux for 3D problems.
1640  *
1641  * @param nConvectiveFields Number of fields.
1642  * @param fields Pointer to fields.
1643  * @param fluxX1 X1-volumetric flux in physical space.
1644  * @param fluxX2 X2-volumetric flux in physical space.
1645  * @param fluxX3 X3-volumetric flux in physical space.
1646  * @param numericalFlux Interface flux in physical space.
1647  * @param divCFlux Divergence of the corrective flux.
1648  *
1649  * \todo: To be implemented. Switch on shapes eventually here.
1650  */
1651 void AdvectionFR::DivCFlux_3D(
1652  const int nConvectiveFields,
1653  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1654  const Array<OneD, const NekDouble> &fluxX1,
1655  const Array<OneD, const NekDouble> &fluxX2,
1656  const Array<OneD, const NekDouble> &fluxX3,
1657  const Array<OneD, const NekDouble> &numericalFlux,
1658  Array<OneD, NekDouble> &divCFlux)
1659 {
1660  boost::ignore_unused(nConvectiveFields, fields, fluxX1, fluxX2, fluxX3,
1661  numericalFlux, divCFlux);
1662 }
1663 
1664 } // namespace SolverUtils
1665 } // 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::AdvectionFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: AdvectionFR.h:64
Nektar::SolverUtils::AdvectionFR::m_dGR_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition: AdvectionFR.h:69
Nektar::SolverUtils::AdvectionFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: AdvectionFR.h:65
Nektar::SolverUtils::AdvectionFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: AdvectionFR.h:67
Nektar::SolverUtils::AdvectionFR::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: AdvectionFR.cpp:114
Nektar::SolverUtils::AdvectionFR::m_dGL_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition: AdvectionFR.h:68
Nektar::SolverUtils::AdvectionFR::m_jac
Array< OneD, NekDouble > m_jac
Definition: AdvectionFR.h:56
Nektar::SolverUtils::AdvectionFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: AdvectionFR.h:59
Nektar::SolverUtils::AdvectionFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: AdvectionFR.h:62
Nektar::SolverUtils::AdvectionFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: AdvectionFR.h:60
Nektar::SolverUtils::AdvectionFR::v_InitObject
virtual void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
Initiliase AdvectionFR objects and store them before starting the time-stepping.
Definition: AdvectionFR.cpp:89
Nektar::SolverUtils::AdvectionFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: AdvectionFR.h:66
Nektar::SolverUtils::AdvectionFR::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: AdvectionFR.cpp:1319
Nektar::SolverUtils::AdvectionFR::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: AdvectionFR.cpp:279
Nektar::SolverUtils::AdvectionFR::v_Advect
virtual void v_Advect(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &advVel, const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble &time, const Array< OneD, Array< OneD, NekDouble >> &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble >> &pBwd=NullNekDoubleArrayOfArray) override
Compute the advection term at each time-step using the Flux Reconstruction approach (FR).
Definition: AdvectionFR.cpp:830
Nektar::SolverUtils::AdvectionFR::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: AdvectionFR.cpp:1135
Nektar::SolverUtils::AdvectionFR::DivCFlux_3D
void DivCFlux_3D(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 > &fluxX3, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 3D problems.
Definition: AdvectionFR.cpp:1651
Nektar::SolverUtils::AdvectionFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: AdvectionFR.h:76
Nektar::SolverUtils::AdvectionFR::DivCFlux_1D
void DivCFlux_1D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 1D problems.
Definition: AdvectionFR.cpp:1030
Nektar::SolverUtils::AdvectionFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: AdvectionFR.h:61
Nektar::SolverUtils::AdvectionFR::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: AdvectionFR.h:57
Nektar::SolverUtils::Advection::m_spaceDim
int m_spaceDim
Storage for space dimension. Used for homogeneous extension.
Definition: Advection.h:215
Nektar::SolverUtils::Advection::m_fluxVector
AdvectionFluxVecCB m_fluxVector
Callback function to the flux vector (set when advection is in conservative form).
Definition: Advection.h:211
Nektar::SolverUtils::Advection::v_InitObject
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Definition: Advection.cpp:353
Nektar::SolverUtils::Advection::m_riemann
RiemannSolverSharedPtr m_riemann
Riemann solver for DG-type schemes.
Definition: Advection.h:213
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::GetAdvectionFactory
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:47
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::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::Vvtvvtp
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:692
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