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