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