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/DisContField1D.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)->GetEdgeQFactors(
187  0, auxArray1 = m_Q2D_e0[n]);
188  pFields[0]->GetExp(n)->GetEdgeQFactors(
189  1, auxArray1 = m_Q2D_e1[n]);
190  pFields[0]->GetExp(n)->GetEdgeQFactors(
191  2, auxArray1 = m_Q2D_e2[n]);
192  pFields[0]->GetExp(n)->GetEdgeQFactors(
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::DisContField1D>(
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  NekDouble tmpFluxVertex = 0;
1060  Vmath::Vcopy(nLocalSolutionPts,
1061  &fluxX1[phys_offset], 1,
1062  &tmparrayX1[0], 1);
1063 
1064  fields[0]->GetExp(n)->GetVertexPhysVals(0, tmparrayX1,
1065  tmpFluxVertex);
1066 
1067  t_offset = fields[0]->GetTrace()
1068  ->GetPhys_Offset(elmtToTrace[n][0]->GetElmtId());
1069 
1070  if(negatedFluxNormal[2*n])
1071  {
1072  JumpL[n] = numericalFlux[t_offset] - tmpFluxVertex;
1073  }
1074  else
1075  {
1076  JumpL[n] = -numericalFlux[t_offset] - tmpFluxVertex;
1077  }
1078 
1079  fields[0]->GetExp(n)->GetVertexPhysVals(1, tmparrayX1,
1080  tmpFluxVertex);
1081 
1082  t_offset = fields[0]->GetTrace()
1083  ->GetPhys_Offset(elmtToTrace[n][1]->GetElmtId());
1084 
1085  if(negatedFluxNormal[2*n+1])
1086  {
1087  JumpR[n] = -numericalFlux[t_offset] - tmpFluxVertex;
1088  }
1089  else
1090  {
1091  JumpR[n] = numericalFlux[t_offset] - tmpFluxVertex;
1092  }
1093  }
1094 
1095  for (n = 0; n < nElements; ++n)
1096  {
1097  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1098  phys_offset = fields[0]->GetPhys_Offset(n);
1099 
1100  // Left jump multiplied by left derivative of C function
1101  Vmath::Smul(nLocalSolutionPts, JumpL[n], &m_dGL_xi1[n][0], 1,
1102  &DCL[0], 1);
1103 
1104  // Right jump multiplied by right derivative of C function
1105  Vmath::Smul(nLocalSolutionPts, JumpR[n], &m_dGR_xi1[n][0], 1,
1106  &DCR[0], 1);
1107 
1108  // Assembling divergence of the correction flux
1109  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1110  &divCFlux[phys_offset], 1);
1111  }
1112  }
1113 
1114  /**
1115  * @brief Compute the divergence of the corrective flux for 2D problems.
1116  *
1117  * @param nConvectiveFields Number of fields.
1118  * @param fields Pointer to fields.
1119  * @param fluxX1 X1-volumetric flux in physical space.
1120  * @param fluxX2 X2-volumetric flux in physical space.
1121  * @param numericalFlux Interface flux in physical space.
1122  * @param divCFlux Divergence of the corrective flux.
1123  *
1124  * \todo: Switch on shapes eventually here.
1125  */
1126  void AdvectionFR::v_DivCFlux_2D(
1127  const int nConvectiveFields,
1128  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1129  const Array<OneD, const NekDouble> &fluxX1,
1130  const Array<OneD, const NekDouble> &fluxX2,
1131  const Array<OneD, const NekDouble> &numericalFlux,
1132  Array<OneD, NekDouble> &divCFlux)
1133  {
1134  boost::ignore_unused(nConvectiveFields);
1135 
1136  int n, e, i, j, cnt;
1137 
1138  int nElements = fields[0]->GetExpSize();
1139 
1140  int nLocalSolutionPts, nEdgePts;
1141  int trace_offset, phys_offset;
1142  int nquad0, nquad1;
1143 
1144  Array<OneD, NekDouble> auxArray1;
1145  Array<OneD, LibUtilities::BasisSharedPtr> base;
1146 
1147  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1148  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1149 
1150  // Loop on the elements
1151  for(n = 0; n < nElements; ++n)
1152  {
1153  // Offset of the element on the global vector
1154  phys_offset = fields[0]->GetPhys_Offset(n);
1155  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1156 
1157  base = fields[0]->GetExp(n)->GetBase();
1158  nquad0 = base[0]->GetNumPoints();
1159  nquad1 = base[1]->GetNumPoints();
1160 
1161  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1162  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1163  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1164  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1165 
1166  // Loop on the edges
1167  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1168  {
1169  // Number of edge points of edge e
1170  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1171 
1172  Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1173  Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1174  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
1175  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
1176  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
1177 
1178  // Offset of the trace space correspondent to edge e
1179  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1180  elmtToTrace[n][e]->GetElmtId());
1181 
1182  // Get the normals of edge e
1183  const Array<OneD, const Array<OneD, NekDouble> > &normals =
1184  fields[0]->GetExp(n)->GetEdgeNormal(e);
1185 
1186  // Extract the edge values of flux-x on edge e and order
1187  // them accordingly to the order of the trace space
1188  fields[0]->GetExp(n)->GetEdgePhysVals(
1189  e, elmtToTrace[n][e],
1190  fluxX1 + phys_offset,
1191  auxArray1 = tmparrayX1);
1192 
1193  // Extract the edge values of flux-y on edge e and order
1194  // them accordingly to the order of the trace space
1195  fields[0]->GetExp(n)->GetEdgePhysVals(
1196  e, elmtToTrace[n][e],
1197  fluxX2 + phys_offset,
1198  auxArray1 = tmparrayX2);
1199 
1200  // Multiply the edge components of the flux by the normal
1201  Vmath::Vvtvvtp(nEdgePts, &tmparrayX1[0], 1,
1202  &m_traceNormals[0][trace_offset], 1,
1203  &tmparrayX2[0], 1,
1204  &m_traceNormals[1][trace_offset], 1,
1205  &fluxN[0], 1);
1206 
1207  // Subtract to the Riemann flux the discontinuous flux
1208  Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1,
1209  &fluxN[0], 1, &fluxJumps[0], 1);
1210 
1211  // Check the ordering of the jump vectors
1212  if (fields[0]->GetExp(n)->GetEorient(e) ==
1213  StdRegions::eBackwards)
1214  {
1215  Vmath::Reverse(nEdgePts, &fluxJumps[0], 1,
1216  &fluxJumps[0], 1);
1217  }
1218 
1219  NekDouble fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1220  -1.0 : 1.0;
1221 
1222  for (i = 0; i < nEdgePts; ++i)
1223  {
1224  if (m_traceNormals[0][trace_offset+i] != fac*normals[0][i]
1225  || m_traceNormals[1][trace_offset+i] != fac*normals[1][i])
1226  {
1227  fluxJumps[i] = -fluxJumps[i];
1228  }
1229  }
1230 
1231  // Multiply jumps by derivatives of the correction functions
1232  switch (e)
1233  {
1234  case 0:
1235  for (i = 0; i < nquad0; ++i)
1236  {
1237  // Multiply fluxJumps by Q factors
1238  fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1239 
1240  for (j = 0; j < nquad1; ++j)
1241  {
1242  cnt = i + j*nquad0;
1243  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1244  }
1245  }
1246  break;
1247  case 1:
1248  for (i = 0; i < nquad1; ++i)
1249  {
1250  // Multiply fluxJumps by Q factors
1251  fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1252 
1253  for (j = 0; j < nquad0; ++j)
1254  {
1255  cnt = (nquad0)*i + j;
1256  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1257  }
1258  }
1259  break;
1260  case 2:
1261  for (i = 0; i < nquad0; ++i)
1262  {
1263  // Multiply fluxJumps by Q factors
1264  fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1265 
1266  for (j = 0; j < nquad1; ++j)
1267  {
1268  cnt = j*nquad0 + i;
1269  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1270  }
1271  }
1272  break;
1273  case 3:
1274  for (i = 0; i < nquad1; ++i)
1275  {
1276  // Multiply fluxJumps by Q factors
1277  fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1278  for (j = 0; j < nquad0; ++j)
1279  {
1280  cnt = j + i*nquad0;
1281  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1282  }
1283  }
1284  break;
1285  default:
1286  ASSERTL0(false,"edge value (< 3) is out of range");
1287  break;
1288  }
1289  }
1290 
1291  // Sum each edge contribution
1292  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1293  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1294 
1295  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1296  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1297 
1298  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1299  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1300  }
1301  }
1302 
1303  /**
1304  * @brief Compute the divergence of the corrective flux for 2D problems
1305  * where POINTSTYPE="GaussGaussLegendre"
1306  *
1307  * @param nConvectiveFields Number of fields.
1308  * @param fields Pointer to fields.
1309  * @param fluxX1 X1-volumetric flux in physical space.
1310  * @param fluxX2 X2-volumetric flux in physical space.
1311  * @param numericalFlux Interface flux in physical space.
1312  * @param divCFlux Divergence of the corrective flux.
1313  *
1314  * \todo: Switch on shapes eventually here.
1315  */
1316 
1317  void AdvectionFR::v_DivCFlux_2D_Gauss(
1318  const int nConvectiveFields,
1319  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1320  const Array<OneD, const NekDouble> &fluxX1,
1321  const Array<OneD, const NekDouble> &fluxX2,
1322  const Array<OneD, const NekDouble> &numericalFlux,
1323  Array<OneD, NekDouble> &divCFlux)
1324  {
1325  boost::ignore_unused(nConvectiveFields);
1326 
1327  int n, e, i, j, cnt;
1328 
1329  int nElements = fields[0]->GetExpSize();
1330  int nLocalSolutionPts;
1331  int nEdgePts;
1332  int trace_offset;
1333  int phys_offset;
1334  int nquad0;
1335  int nquad1;
1336 
1337  NekDouble fac;
1338 
1339  Array<OneD, NekDouble> auxArray1, auxArray2;
1340  Array<OneD, LibUtilities::BasisSharedPtr> base;
1341 
1342  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1343  &elmtToTrace = 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)->GetNedges(); ++e)
1363  {
1364  // Number of edge points of edge e
1365  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(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)->GetEdgeNormal(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)->GetEdgePhysVals(
1389  e, elmtToTrace[n][e],
1390  fluxX2 + phys_offset,
1391  auxArray1 = fluxN_D);
1392 
1393  Vmath::Neg (nEdgePts, fluxN_D, 1);
1394 
1395  // Extract the physical Rieamann flux at each edge
1396  Vmath::Vcopy(nEdgePts,
1397  &numericalFlux[trace_offset], 1,
1398  &fluxN[0], 1);
1399 
1400  // Check the ordering of vectors
1401  if (fields[0]->GetExp(n)->GetEorient(e) ==
1402  StdRegions::eBackwards)
1403  {
1404  Vmath::Reverse(nEdgePts,
1405  auxArray1 = fluxN, 1,
1406  auxArray2 = fluxN, 1);
1407 
1408  Vmath::Reverse(nEdgePts,
1409  auxArray1 = fluxN_D, 1,
1410  auxArray2 = fluxN_D, 1);
1411  }
1412 
1413  // Transform Riemann Fluxes in the standard element
1414  for (i = 0; i < nquad0; ++i)
1415  {
1416  // Multiply Riemann Flux by Q factors
1417  fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
1418  }
1419 
1420  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1421  -1.0 : 1.0;
1422 
1423  for (i = 0; i < nEdgePts; ++i)
1424  {
1425  if (m_traceNormals[0][trace_offset+i]
1426  != fac*normals[0][i] ||
1427  m_traceNormals[1][trace_offset+i]
1428  != fac*normals[1][i])
1429  {
1430  fluxN_R[i] = -fluxN_R[i];
1431  }
1432  }
1433 
1434  // Subtract to the Riemann flux the discontinuous
1435  // flux
1436  Vmath::Vsub(nEdgePts,
1437  &fluxN_R[0], 1,
1438  &fluxN_D[0], 1, &fluxJumps[0], 1);
1439 
1440  // Multiplicate the flux jumps for the correction
1441  // function
1442  for (i = 0; i < nquad0; ++i)
1443  {
1444  for (j = 0; j < nquad1; ++j)
1445  {
1446  cnt = i + j*nquad0;
1447  divCFluxE0[cnt] =
1448  -fluxJumps[i] * m_dGL_xi2[n][j];
1449  }
1450  }
1451  break;
1452  case 1:
1453  // Extract the edge values of transformed flux-x on
1454  // edge e and order them accordingly to the order of
1455  // the trace space
1456  fields[0]->GetExp(n)->GetEdgePhysVals(
1457  e, elmtToTrace[n][e],
1458  fluxX1 + phys_offset,
1459  auxArray1 = fluxN_D);
1460 
1461  // Extract the physical Rieamann flux at each edge
1462  Vmath::Vcopy(nEdgePts,
1463  &numericalFlux[trace_offset], 1,
1464  &fluxN[0], 1);
1465 
1466  // Check the ordering of vectors
1467  if (fields[0]->GetExp(n)->GetEorient(e) ==
1468  StdRegions::eBackwards)
1469  {
1470  Vmath::Reverse(nEdgePts,
1471  auxArray1 = fluxN, 1,
1472  auxArray2 = fluxN, 1);
1473 
1474  Vmath::Reverse(nEdgePts,
1475  auxArray1 = fluxN_D, 1,
1476  auxArray2 = fluxN_D, 1);
1477  }
1478 
1479  // Transform Riemann Fluxes in the standard element
1480  for (i = 0; i < nquad1; ++i)
1481  {
1482  // Multiply Riemann Flux by Q factors
1483  fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
1484  }
1485 
1486  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1487  -1.0 : 1.0;
1488 
1489  for (i = 0; i < nEdgePts; ++i)
1490  {
1491  if (m_traceNormals[0][trace_offset+i]
1492  != fac*normals[0][i] ||
1493  m_traceNormals[1][trace_offset+i]
1494  != fac*normals[1][i])
1495  {
1496  fluxN_R[i] = -fluxN_R[i];
1497  }
1498  }
1499 
1500  // Subtract to the Riemann flux the discontinuous
1501  // flux
1502  Vmath::Vsub(nEdgePts,
1503  &fluxN_R[0], 1,
1504  &fluxN_D[0], 1, &fluxJumps[0], 1);
1505 
1506  // Multiplicate the flux jumps for the correction
1507  // function
1508  for (i = 0; i < nquad1; ++i)
1509  {
1510  for (j = 0; j < nquad0; ++j)
1511  {
1512  cnt = (nquad0)*i + j;
1513  divCFluxE1[cnt] =
1514  fluxJumps[i] * m_dGR_xi1[n][j];
1515  }
1516  }
1517  break;
1518  case 2:
1519 
1520  // Extract the edge values of transformed flux-y on
1521  // edge e and order them accordingly to the order of
1522  // the trace space
1523 
1524  fields[0]->GetExp(n)->GetEdgePhysVals(
1525  e, elmtToTrace[n][e],
1526  fluxX2 + phys_offset,
1527  auxArray1 = fluxN_D);
1528 
1529  // Extract the physical Rieamann flux at each edge
1530  Vmath::Vcopy(nEdgePts,
1531  &numericalFlux[trace_offset], 1,
1532  &fluxN[0], 1);
1533 
1534  // Check the ordering of vectors
1535  if (fields[0]->GetExp(n)->GetEorient(e) ==
1536  StdRegions::eBackwards)
1537  {
1538  Vmath::Reverse(nEdgePts,
1539  auxArray1 = fluxN, 1,
1540  auxArray2 = fluxN, 1);
1541 
1542  Vmath::Reverse(nEdgePts,
1543  auxArray1 = fluxN_D, 1,
1544  auxArray2 = fluxN_D, 1);
1545  }
1546 
1547  // Transform Riemann Fluxes in the standard element
1548  for (i = 0; i < nquad0; ++i)
1549  {
1550  // Multiply Riemann Flux by Q factors
1551  fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
1552  }
1553 
1554  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1555  -1.0 : 1.0;
1556 
1557  for (i = 0; i < nEdgePts; ++i)
1558  {
1559  if (m_traceNormals[0][trace_offset+i]
1560  != fac*normals[0][i] ||
1561  m_traceNormals[1][trace_offset+i]
1562  != fac*normals[1][i])
1563  {
1564  fluxN_R[i] = -fluxN_R[i];
1565  }
1566  }
1567 
1568  // Subtract to the Riemann flux the discontinuous
1569  // flux
1570 
1571  Vmath::Vsub(nEdgePts,
1572  &fluxN_R[0], 1,
1573  &fluxN_D[0], 1, &fluxJumps[0], 1);
1574 
1575  // Multiplicate the flux jumps for the correction
1576  // function
1577  for (i = 0; i < nquad0; ++i)
1578  {
1579  for (j = 0; j < nquad1; ++j)
1580  {
1581  cnt = j*nquad0 + i;
1582  divCFluxE2[cnt] =
1583  fluxJumps[i] * m_dGR_xi2[n][j];
1584  }
1585  }
1586  break;
1587  case 3:
1588  // Extract the edge values of transformed flux-x on
1589  // edge e and order them accordingly to the order of
1590  // the trace space
1591 
1592  fields[0]->GetExp(n)->GetEdgePhysVals(
1593  e, elmtToTrace[n][e],
1594  fluxX1 + phys_offset,
1595  auxArray1 = fluxN_D);
1596  Vmath::Neg (nEdgePts, fluxN_D, 1);
1597 
1598  // Extract the physical Rieamann flux at each edge
1599  Vmath::Vcopy(nEdgePts,
1600  &numericalFlux[trace_offset], 1,
1601  &fluxN[0], 1);
1602 
1603  // Check the ordering of vectors
1604  if (fields[0]->GetExp(n)->GetEorient(e) ==
1605  StdRegions::eBackwards)
1606  {
1607  Vmath::Reverse(nEdgePts,
1608  auxArray1 = fluxN, 1,
1609  auxArray2 = fluxN, 1);
1610 
1611  Vmath::Reverse(nEdgePts,
1612  auxArray1 = fluxN_D, 1,
1613  auxArray2 = fluxN_D, 1);
1614  }
1615 
1616  // Transform Riemann Fluxes in the standard element
1617  for (i = 0; i < nquad1; ++i)
1618  {
1619  // Multiply Riemann Flux by Q factors
1620  fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
1621  }
1622 
1623  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1624  -1.0 : 1.0;
1625 
1626  for (i = 0; i < nEdgePts; ++i)
1627  {
1628  if (m_traceNormals[0][trace_offset+i]
1629  != fac*normals[0][i] ||
1630  m_traceNormals[1][trace_offset+i]
1631  != fac*normals[1][i])
1632  {
1633  fluxN_R[i] = -fluxN_R[i];
1634  }
1635  }
1636 
1637  // Subtract to the Riemann flux the discontinuous
1638  // flux
1639 
1640  Vmath::Vsub(nEdgePts,
1641  &fluxN_R[0], 1,
1642  &fluxN_D[0], 1, &fluxJumps[0], 1);
1643 
1644  // Multiplicate the flux jumps for the correction
1645  // function
1646  for (i = 0; i < nquad1; ++i)
1647  {
1648  for (j = 0; j < nquad0; ++j)
1649  {
1650  cnt = j + i*nquad0;
1651  divCFluxE3[cnt] =
1652  -fluxJumps[i] * m_dGL_xi1[n][j];
1653  }
1654  }
1655  break;
1656  default:
1657  ASSERTL0(false,"edge value (< 3) is out of range");
1658  break;
1659  }
1660  }
1661 
1662 
1663  // Sum each edge contribution
1664  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1665  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1666 
1667  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1668  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1669 
1670  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1671  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1672  }
1673  }
1674 
1675 
1676  /**
1677  * @brief Compute the divergence of the corrective flux for 3D problems.
1678  *
1679  * @param nConvectiveFields Number of fields.
1680  * @param fields Pointer to fields.
1681  * @param fluxX1 X1-volumetric flux in physical space.
1682  * @param fluxX2 X2-volumetric flux in physical space.
1683  * @param fluxX3 X3-volumetric flux in physical space.
1684  * @param numericalFlux Interface flux in physical space.
1685  * @param divCFlux Divergence of the corrective flux.
1686  *
1687  * \todo: To be implemented. Switch on shapes eventually here.
1688  */
1689  void AdvectionFR::v_DivCFlux_3D(
1690  const int nConvectiveFields,
1691  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1692  const Array<OneD, const NekDouble> &fluxX1,
1693  const Array<OneD, const NekDouble> &fluxX2,
1694  const Array<OneD, const NekDouble> &fluxX3,
1695  const Array<OneD, const NekDouble> &numericalFlux,
1696  Array<OneD, NekDouble> &divCFlux)
1697  {
1698  boost::ignore_unused(nConvectiveFields, fields, fluxX1, fluxX2,
1699  fluxX3, numericalFlux, divCFlux);
1700  }
1701 
1702  }
1703 }
Nektar::SolverUtils::AdvectionFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: AdvectionFR.h:60
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::SolverUtils::AdvectionFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: AdvectionFR.h:64
Nektar::SolverUtils::AdvectionFR::m_dGL_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition: AdvectionFR.h:68
Nektar::LibUtilities::Basis
Represents a basis of a given type.
Definition: Basis.h:211
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
Nektar::SolverUtils::AdvectionFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: AdvectionFR.h:65
Nektar::SolverUtils::AdvectionFR::m_dGR_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition: AdvectionFR.h:69
Nektar::SolverUtils::AdvectionFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: AdvectionFR.h:59
Nektar::SolverUtils::AdvectionFR::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: AdvectionFR.h:57
Nektar::SolverUtils::AdvectionFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: AdvectionFR.h:67
Expansion2D.h
std
STL namespace.
Nektar::SolverUtils::Advection::m_riemann
RiemannSolverSharedPtr m_riemann
Riemann solver for DG-type schemes.
Definition: Advection.h:143
Nektar::SolverUtils::AdvectionFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: AdvectionFR.h:62
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:244
Nektar::LibUtilities::BasisSharedPtr
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:74
Points.h
Nektar::LibUtilities::eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:48
Nektar::SolverUtils::AdvectionFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: AdvectionFR.h:61
StdSegExp.h
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1088
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
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:1131
Basis.h
ManagerAccess.h
AdvectionFR.h
Nektar::MultiRegions::DisContField1D
This class is the abstraction of a global discontinuous two- dimensional spectral/hp element expansio...
Definition: DisContField1D.h:54
Nektar::SolverUtils::GetAdvectionFactory
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:47
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:399
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
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:1689
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:1317
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:1126
Nektar::LocalRegions::Expansion::GetMetricInfo
const SpatialDomains::GeomFactorsSharedPtr & GetMetricInfo() const
Definition: Expansion.cpp:181
Nektar::SolverUtils::AdvectionFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: AdvectionFR.h:66
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:346
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
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:540
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
Polylib.h
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::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_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: AdvectionFR.h:76
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:199
Nektar::SolverUtils::Advection::m_fluxVector
AdvectionFluxVecCB m_fluxVector
Callback function to the flux vector (set when advection is in conservative form).
Definition: Advection.h:141
Expansion1D.h
AssemblyMapDG.h
Nektar::SolverUtils::Advection::m_spaceDim
int m_spaceDim
Storage for space dimension. Used for homogeneous extension.
Definition: Advection.h:145
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Nektar::SpatialDomains::eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:56
DisContField1D.h
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:112
Nektar::SolverUtils::AdvectionFR::m_jac
Array< OneD, NekDouble > m_jac
Definition: AdvectionFR.h:56
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:98
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:302
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:186
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–7...
Definition: AdvectionFR.cpp:272
Generated on Wed Jan 20 2021 03:33:31 for Nektar++ by   doxygen 1.8.13