Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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  {
809  int i, j, n;
810  int phys_offset;
811 
812  Array<OneD, NekDouble> auxArray1, auxArray2;
813 
814  Array<OneD, LibUtilities::BasisSharedPtr> Basis;
815  Basis = fields[0]->GetExp(0)->GetBase();
816 
817  int nElements = fields[0]->GetExpSize();
818  int nDimensions = fields[0]->GetCoordim(0);
819  int nSolutionPts = fields[0]->GetTotPoints();
820  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
821  int nCoeffs = fields[0]->GetNcoeffs();
822 
823  // Storage for flux vector.
824  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > fluxvector
825  (nConvectiveFields);
826  // Outarray for Galerkin projection in case of primitive dealising
827  Array<OneD, Array<OneD, NekDouble> > outarrayCoeff
828  (nConvectiveFields);
829  // Store forwards/backwards space along trace space
830  Array<OneD, Array<OneD, NekDouble> > Fwd (nConvectiveFields);
831  Array<OneD, Array<OneD, NekDouble> > Bwd (nConvectiveFields);
832  Array<OneD, Array<OneD, NekDouble> > numflux(nConvectiveFields);
833 
834  // Set up storage for flux vector.
835  for (i = 0; i < nConvectiveFields; ++i)
836  {
837  fluxvector[i] =
838  Array<OneD, Array<OneD, NekDouble> >(m_spaceDim);
839 
840  for (j = 0; j < m_spaceDim; ++j)
841  {
842  fluxvector[i][j] = Array<OneD, NekDouble>(nSolutionPts);
843  }
844  }
845 
846  for (i = 0; i < nConvectiveFields; ++i)
847  {
848  outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
849  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
850  Bwd[i] = Array<OneD, NekDouble>(nTracePts);
851  numflux[i] = Array<OneD, NekDouble>(nTracePts);
852  fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
853  }
854 
855  // Computing the interface flux at each trace point
856  m_riemann->Solve(m_spaceDim, Fwd, Bwd, numflux);
857 
858  // Divergence of the flux (computing the RHS)
859  switch(nDimensions)
860  {
861  // 1D problems
862  case 1:
863  {
864  Array<OneD, NekDouble> DfluxvectorX1(nSolutionPts);
865  Array<OneD, NekDouble> divFC (nSolutionPts);
866 
867  // Get the advection flux vector (solver specific)
868  m_fluxVector(inarray, fluxvector);
869 
870  // Get the discontinuous flux divergence
871  for(i = 0; i < nConvectiveFields; ++i)
872  {
873  for (n = 0; n < nElements; n++)
874  {
875  phys_offset = fields[0]->GetPhys_Offset(n);
876 
877  fields[i]->GetExp(n)->PhysDeriv(
878  0, fluxvector[i][0] + phys_offset,
879  auxArray2 = DfluxvectorX1 + phys_offset);
880  }
881 
882  // Get the correction flux divergence
883  v_DivCFlux_1D(nConvectiveFields, fields,
884  fluxvector[i][0], numflux[i], divFC);
885 
886  // Back to the physical space using global operations
887  Vmath::Vdiv(nSolutionPts, &divFC[0], 1, &m_jac[0], 1,
888  &outarray[i][0], 1);
889 
890  // Adding the total divergence to outarray (RHS)
891  Vmath::Vadd(nSolutionPts, &outarray[i][0], 1,
892  &DfluxvectorX1[0], 1, &outarray[i][0], 1);
893 
894  // Primitive Dealiasing 1D
895  if (!(Basis[0]->Collocation()))
896  {
897  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
898  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
899  }
900  }
901  break;
902  }
903  // 2D problems
904  case 2:
905  {
906  Array<OneD, NekDouble> DfluxvectorX1(nSolutionPts);
907  Array<OneD, NekDouble> DfluxvectorX2(nSolutionPts);
908  Array<OneD, NekDouble> divFD(nSolutionPts);
909  Array<OneD, NekDouble> divFC(nSolutionPts);
910 
911  // Get the advection flux vector (solver specific)
912  m_fluxVector(inarray, fluxvector);
913 
914  for (i = 0; i < nConvectiveFields; ++i)
915  {
916  // Temporary vectors
917  Array<OneD, NekDouble> f_hat(nSolutionPts);
918  Array<OneD, NekDouble> g_hat(nSolutionPts);
919 
920  Vmath::Vvtvvtp(nSolutionPts,
921  &m_gmat[0][0], 1,
922  &fluxvector[i][0][0], 1,
923  &m_gmat[2][0], 1,
924  &fluxvector[i][1][0], 1,
925  &f_hat[0], 1);
926 
927  Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &f_hat[0], 1,
928  &f_hat[0], 1);
929 
930  Vmath::Vvtvvtp(nSolutionPts,
931  &m_gmat[1][0], 1,
932  &fluxvector[i][0][0], 1,
933  &m_gmat[3][0], 1,
934  &fluxvector[i][1][0], 1,
935  &g_hat[0], 1);
936 
937  Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &g_hat[0], 1,
938  &g_hat[0], 1);
939 
940  // Get the discontinuous flux derivatives
941  for (n = 0; n < nElements; n++)
942  {
943  phys_offset = fields[0]->GetPhys_Offset(n);
944  fields[0]->GetExp(n)->StdPhysDeriv(0,
945  auxArray1 = f_hat + phys_offset,
946  auxArray2 = DfluxvectorX1 + phys_offset);
947  fields[0]->GetExp(n)->StdPhysDeriv(1,
948  auxArray1 = g_hat + phys_offset,
949  auxArray2 = DfluxvectorX2 + phys_offset);
950  }
951 
952  // Divergence of the discontinuous flux
953  Vmath::Vadd(nSolutionPts, DfluxvectorX1, 1,
954  DfluxvectorX2, 1, divFD, 1);
955 
956  // Divergence of the correction flux
957  if (Basis[0]->GetPointsType() ==
958  LibUtilities::eGaussGaussLegendre &&
959  Basis[1]->GetPointsType() ==
960  LibUtilities::eGaussGaussLegendre)
961  {
962  v_DivCFlux_2D_Gauss(
963  nConvectiveFields, fields, f_hat, g_hat,
964  numflux[i], divFC);
965  }
966  else
967  {
968  v_DivCFlux_2D(
969  nConvectiveFields, fields,
970  fluxvector[i][0], fluxvector[i][1],
971  numflux[i], divFC);
972 
973  }
974 
975  // Divergence of the final flux
976  Vmath::Vadd(nSolutionPts, divFD, 1, divFC, 1,
977  outarray[i], 1);
978 
979  // Back to the physical space using a global operation
980  Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1,
981  &m_jac[0], 1, &outarray[i][0], 1);
982 
983  // Primitive Dealiasing 2D
984  if (!(Basis[0]->Collocation()))
985  {
986  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
987  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
988  }
989  } // close nConvectiveFields loop
990  break;
991  }
992  // 3D problems
993  case 3:
994  {
995  ASSERTL0(false,"3D FRDG case not implemented yet");
996  break;
997  }
998  }
999  }
1000 
1001  /**
1002  * @brief Compute the divergence of the corrective flux for 1D problems.
1003  *
1004  * @param nConvectiveFields Number of fields.
1005  * @param fields Pointer to fields.
1006  * @param fluxX1 X1-volumetric flux in physical space.
1007  * @param numericalFlux Interface flux in physical space.
1008  * @param divCFlux Divergence of the corrective flux.
1009  *
1010  */
1011  void AdvectionFR::v_DivCFlux_1D(
1012  const int nConvectiveFields,
1013  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1014  const Array<OneD, const NekDouble> &fluxX1,
1015  const Array<OneD, const NekDouble> &numericalFlux,
1016  Array<OneD, NekDouble> &divCFlux)
1017  {
1018  int n;
1019  int nLocalSolutionPts, phys_offset, t_offset;
1020 
1021  LibUtilities::BasisSharedPtr Basis;
1022  Basis = fields[0]->GetExp(0)->GetBasis(0);
1023 
1024  int nElements = fields[0]->GetExpSize();
1025  int nSolutionPts = fields[0]->GetTotPoints();
1026 
1027 
1028  vector<bool> negatedFluxNormal = (boost::static_pointer_cast<MultiRegions::DisContField1D>(fields[0]))->GetNegatedFluxNormal();
1029 
1030  // Arrays to store the derivatives of the correction flux
1031  Array<OneD, NekDouble> DCL(nSolutionPts/nElements, 0.0);
1032  Array<OneD, NekDouble> DCR(nSolutionPts/nElements, 0.0);
1033 
1034  // Arrays to store the intercell numerical flux jumps
1035  Array<OneD, NekDouble> JumpL(nElements);
1036  Array<OneD, NekDouble> JumpR(nElements);
1037 
1038  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1039  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1040 
1041  for (n = 0; n < nElements; ++n)
1042  {
1043  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1044  phys_offset = fields[0]->GetPhys_Offset(n);
1045 
1046  Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1047  NekDouble tmpFluxVertex = 0;
1048  Vmath::Vcopy(nLocalSolutionPts,
1049  &fluxX1[phys_offset], 1,
1050  &tmparrayX1[0], 1);
1051 
1052  fields[0]->GetExp(n)->GetVertexPhysVals(0, tmparrayX1,
1053  tmpFluxVertex);
1054 
1055  t_offset = fields[0]->GetTrace()
1056  ->GetPhys_Offset(elmtToTrace[n][0]->GetElmtId());
1057 
1058  if(negatedFluxNormal[2*n])
1059  {
1060  JumpL[n] = numericalFlux[t_offset] - tmpFluxVertex;
1061  }
1062  else
1063  {
1064  JumpL[n] = -numericalFlux[t_offset] - tmpFluxVertex;
1065  }
1066 
1067  fields[0]->GetExp(n)->GetVertexPhysVals(1, tmparrayX1,
1068  tmpFluxVertex);
1069 
1070  t_offset = fields[0]->GetTrace()
1071  ->GetPhys_Offset(elmtToTrace[n][1]->GetElmtId());
1072 
1073  if(negatedFluxNormal[2*n+1])
1074  {
1075  JumpR[n] = -numericalFlux[t_offset] - tmpFluxVertex;
1076  }
1077  else
1078  {
1079  JumpR[n] = numericalFlux[t_offset] - tmpFluxVertex;
1080  }
1081  }
1082 
1083  for (n = 0; n < nElements; ++n)
1084  {
1085  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1086  phys_offset = fields[0]->GetPhys_Offset(n);
1087 
1088  // Left jump multiplied by left derivative of C function
1089  Vmath::Smul(nLocalSolutionPts, JumpL[n], &m_dGL_xi1[n][0], 1,
1090  &DCL[0], 1);
1091 
1092  // Right jump multiplied by right derivative of C function
1093  Vmath::Smul(nLocalSolutionPts, JumpR[n], &m_dGR_xi1[n][0], 1,
1094  &DCR[0], 1);
1095 
1096  // Assembling divergence of the correction flux
1097  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1098  &divCFlux[phys_offset], 1);
1099  }
1100  }
1101 
1102  /**
1103  * @brief Compute the divergence of the corrective flux for 2D problems.
1104  *
1105  * @param nConvectiveFields Number of fields.
1106  * @param fields Pointer to fields.
1107  * @param fluxX1 X1-volumetric flux in physical space.
1108  * @param fluxX2 X2-volumetric flux in physical space.
1109  * @param numericalFlux Interface flux in physical space.
1110  * @param divCFlux Divergence of the corrective flux.
1111  *
1112  * \todo: Switch on shapes eventually here.
1113  */
1114  void AdvectionFR::v_DivCFlux_2D(
1115  const int nConvectiveFields,
1116  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1117  const Array<OneD, const NekDouble> &fluxX1,
1118  const Array<OneD, const NekDouble> &fluxX2,
1119  const Array<OneD, const NekDouble> &numericalFlux,
1120  Array<OneD, NekDouble> &divCFlux)
1121  {
1122  int n, e, i, j, cnt;
1123 
1124  int nElements = fields[0]->GetExpSize();
1125 
1126  int nLocalSolutionPts, nEdgePts;
1127  int trace_offset, phys_offset;
1128  int nquad0, nquad1;
1129 
1130  Array<OneD, NekDouble> auxArray1;
1131  Array<OneD, LibUtilities::BasisSharedPtr> base;
1132 
1133  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1134  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1135 
1136  // Loop on the elements
1137  for(n = 0; n < nElements; ++n)
1138  {
1139  // Offset of the element on the global vector
1140  phys_offset = fields[0]->GetPhys_Offset(n);
1141  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1142 
1143  base = fields[0]->GetExp(n)->GetBase();
1144  nquad0 = base[0]->GetNumPoints();
1145  nquad1 = base[1]->GetNumPoints();
1146 
1147  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1148  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1149  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1150  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1151 
1152  // Loop on the edges
1153  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1154  {
1155  // Number of edge points of edge e
1156  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1157 
1158  Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1159  Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1160  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
1161  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
1162  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
1163 
1164  // Offset of the trace space correspondent to edge e
1165  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1166  elmtToTrace[n][e]->GetElmtId());
1167 
1168  // Get the normals of edge e
1169  const Array<OneD, const Array<OneD, NekDouble> > &normals =
1170  fields[0]->GetExp(n)->GetEdgeNormal(e);
1171 
1172  // Extract the edge values of flux-x on edge e and order
1173  // them accordingly to the order of the trace space
1174  fields[0]->GetExp(n)->GetEdgePhysVals(
1175  e, elmtToTrace[n][e],
1176  fluxX1 + phys_offset,
1177  auxArray1 = tmparrayX1);
1178 
1179  // Extract the edge values of flux-y on edge e and order
1180  // them accordingly to the order of the trace space
1181  fields[0]->GetExp(n)->GetEdgePhysVals(
1182  e, elmtToTrace[n][e],
1183  fluxX2 + phys_offset,
1184  auxArray1 = tmparrayX2);
1185 
1186  // Multiply the edge components of the flux by the normal
1187  Vmath::Vvtvvtp(nEdgePts, &tmparrayX1[0], 1,
1188  &m_traceNormals[0][trace_offset], 1,
1189  &tmparrayX2[0], 1,
1190  &m_traceNormals[1][trace_offset], 1,
1191  &fluxN[0], 1);
1192 
1193  // Subtract to the Riemann flux the discontinuous flux
1194  Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1,
1195  &fluxN[0], 1, &fluxJumps[0], 1);
1196 
1197  // Check the ordering of the jump vectors
1198  if (fields[0]->GetExp(n)->GetEorient(e) ==
1199  StdRegions::eBackwards)
1200  {
1201  Vmath::Reverse(nEdgePts, &fluxJumps[0], 1,
1202  &fluxJumps[0], 1);
1203  }
1204 
1205  NekDouble fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1206  -1.0 : 1.0;
1207 
1208  for (i = 0; i < nEdgePts; ++i)
1209  {
1210  if (m_traceNormals[0][trace_offset+i] != fac*normals[0][i]
1211  || m_traceNormals[1][trace_offset+i] != fac*normals[1][i])
1212  {
1213  fluxJumps[i] = -fluxJumps[i];
1214  }
1215  }
1216 
1217  // Multiply jumps by derivatives of the correction functions
1218  switch (e)
1219  {
1220  case 0:
1221  for (i = 0; i < nquad0; ++i)
1222  {
1223  // Multiply fluxJumps by Q factors
1224  fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1225 
1226  for (j = 0; j < nquad1; ++j)
1227  {
1228  cnt = i + j*nquad0;
1229  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1230  }
1231  }
1232  break;
1233  case 1:
1234  for (i = 0; i < nquad1; ++i)
1235  {
1236  // Multiply fluxJumps by Q factors
1237  fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1238 
1239  for (j = 0; j < nquad0; ++j)
1240  {
1241  cnt = (nquad0)*i + j;
1242  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1243  }
1244  }
1245  break;
1246  case 2:
1247  for (i = 0; i < nquad0; ++i)
1248  {
1249  // Multiply fluxJumps by Q factors
1250  fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1251 
1252  for (j = 0; j < nquad1; ++j)
1253  {
1254  cnt = (nquad0 - 1) + j*nquad0 - i;
1255  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1256  }
1257  }
1258  break;
1259  case 3:
1260  for (i = 0; i < nquad1; ++i)
1261  {
1262  // Multiply fluxJumps by Q factors
1263  fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1264  for (j = 0; j < nquad0; ++j)
1265  {
1266  cnt = (nquad0*nquad1 - nquad0) + j - i*nquad0;
1267  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1268  }
1269  }
1270  break;
1271  default:
1272  ASSERTL0(false,"edge value (< 3) is out of range");
1273  break;
1274  }
1275  }
1276 
1277  // Sum each edge contribution
1278  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1279  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1280 
1281  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1282  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1283 
1284  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1285  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1286  }
1287  }
1288 
1289  /**
1290  * @brief Compute the divergence of the corrective flux for 2D problems
1291  * where POINTSTYPE="GaussGaussLegendre"
1292  *
1293  * @param nConvectiveFields Number of fields.
1294  * @param fields Pointer to fields.
1295  * @param fluxX1 X1-volumetric flux in physical space.
1296  * @param fluxX2 X2-volumetric flux in physical space.
1297  * @param numericalFlux Interface flux in physical space.
1298  * @param divCFlux Divergence of the corrective flux.
1299  *
1300  * \todo: Switch on shapes eventually here.
1301  */
1302 
1303  void AdvectionFR::v_DivCFlux_2D_Gauss(
1304  const int nConvectiveFields,
1305  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1306  const Array<OneD, const NekDouble> &fluxX1,
1307  const Array<OneD, const NekDouble> &fluxX2,
1308  const Array<OneD, const NekDouble> &numericalFlux,
1309  Array<OneD, NekDouble> &divCFlux)
1310  {
1311  int n, e, i, j, cnt;
1312 
1313  int nElements = fields[0]->GetExpSize();
1314  int nLocalSolutionPts;
1315  int nEdgePts;
1316  int trace_offset;
1317  int phys_offset;
1318  int nquad0;
1319  int nquad1;
1320 
1321  NekDouble fac;
1322 
1323  Array<OneD, NekDouble> auxArray1, auxArray2;
1324  Array<OneD, LibUtilities::BasisSharedPtr> base;
1325 
1326  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1327  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1328 
1329  // Loop on the elements
1330  for(n = 0; n < nElements; ++n)
1331  {
1332  // Offset of the element on the global vector
1333  phys_offset = fields[0]->GetPhys_Offset(n);
1334  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1335 
1336  base = fields[0]->GetExp(n)->GetBase();
1337  nquad0 = base[0]->GetNumPoints();
1338  nquad1 = base[1]->GetNumPoints();
1339 
1340  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1341  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1342  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1343  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1344 
1345  // Loop on the edges
1346  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1347  {
1348  // Number of edge points of edge e
1349  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1350 
1351  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
1352  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
1353  Array<OneD, NekDouble> fluxN_R (nEdgePts, 0.0);
1354  Array<OneD, NekDouble> fluxN_D (nEdgePts, 0.0);
1355  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
1356 
1357  // Offset of the trace space correspondent to edge e
1358  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1359  elmtToTrace[n][e]->GetElmtId());
1360 
1361  // Get the normals of edge e
1362  const Array<OneD, const Array<OneD, NekDouble> > &normals =
1363  fields[0]->GetExp(n)->GetEdgeNormal(e);
1364 
1365  // Extract the trasformed normal flux at each edge
1366  switch (e)
1367  {
1368  case 0:
1369  // Extract the edge values of transformed flux-y on
1370  // edge e and order them accordingly to the order of
1371  // the trace space
1372  fields[0]->GetExp(n)->GetEdgePhysVals(
1373  e, elmtToTrace[n][e],
1374  fluxX2 + phys_offset,
1375  auxArray1 = fluxN_D);
1376 
1377  Vmath::Neg (nEdgePts, fluxN_D, 1);
1378 
1379  // Extract the physical Rieamann flux at each edge
1380  Vmath::Vcopy(nEdgePts,
1381  &numericalFlux[trace_offset], 1,
1382  &fluxN[0], 1);
1383 
1384  // Check the ordering of vectors
1385  if (fields[0]->GetExp(n)->GetEorient(e) ==
1386  StdRegions::eBackwards)
1387  {
1388  Vmath::Reverse(nEdgePts,
1389  auxArray1 = fluxN, 1,
1390  auxArray2 = fluxN, 1);
1391 
1392  Vmath::Reverse(nEdgePts,
1393  auxArray1 = fluxN_D, 1,
1394  auxArray2 = fluxN_D, 1);
1395  }
1396 
1397  // Transform Riemann Fluxes in the standard element
1398  for (i = 0; i < nquad0; ++i)
1399  {
1400  // Multiply Riemann Flux by Q factors
1401  fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
1402  }
1403 
1404  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1405  -1.0 : 1.0;
1406 
1407  for (i = 0; i < nEdgePts; ++i)
1408  {
1409  if (m_traceNormals[0][trace_offset+i]
1410  != fac*normals[0][i] ||
1411  m_traceNormals[1][trace_offset+i]
1412  != fac*normals[1][i])
1413  {
1414  fluxN_R[i] = -fluxN_R[i];
1415  }
1416  }
1417 
1418  // Subtract to the Riemann flux the discontinuous
1419  // flux
1420  Vmath::Vsub(nEdgePts,
1421  &fluxN_R[0], 1,
1422  &fluxN_D[0], 1, &fluxJumps[0], 1);
1423 
1424  // Multiplicate the flux jumps for the correction
1425  // function
1426  for (i = 0; i < nquad0; ++i)
1427  {
1428  for (j = 0; j < nquad1; ++j)
1429  {
1430  cnt = i + j*nquad0;
1431  divCFluxE0[cnt] =
1432  -fluxJumps[i] * m_dGL_xi2[n][j];
1433  }
1434  }
1435  break;
1436  case 1:
1437  // Extract the edge values of transformed flux-x on
1438  // edge e and order them accordingly to the order of
1439  // the trace space
1440  fields[0]->GetExp(n)->GetEdgePhysVals(
1441  e, elmtToTrace[n][e],
1442  fluxX1 + phys_offset,
1443  auxArray1 = fluxN_D);
1444 
1445  // Extract the physical Rieamann flux at each edge
1446  Vmath::Vcopy(nEdgePts,
1447  &numericalFlux[trace_offset], 1,
1448  &fluxN[0], 1);
1449 
1450  // Check the ordering of vectors
1451  if (fields[0]->GetExp(n)->GetEorient(e) ==
1452  StdRegions::eBackwards)
1453  {
1454  Vmath::Reverse(nEdgePts,
1455  auxArray1 = fluxN, 1,
1456  auxArray2 = fluxN, 1);
1457 
1458  Vmath::Reverse(nEdgePts,
1459  auxArray1 = fluxN_D, 1,
1460  auxArray2 = fluxN_D, 1);
1461  }
1462 
1463  // Transform Riemann Fluxes in the standard element
1464  for (i = 0; i < nquad1; ++i)
1465  {
1466  // Multiply Riemann Flux by Q factors
1467  fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
1468  }
1469 
1470  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1471  -1.0 : 1.0;
1472 
1473  for (i = 0; i < nEdgePts; ++i)
1474  {
1475  if (m_traceNormals[0][trace_offset+i]
1476  != fac*normals[0][i] ||
1477  m_traceNormals[1][trace_offset+i]
1478  != fac*normals[1][i])
1479  {
1480  fluxN_R[i] = -fluxN_R[i];
1481  }
1482  }
1483 
1484  // Subtract to the Riemann flux the discontinuous
1485  // flux
1486  Vmath::Vsub(nEdgePts,
1487  &fluxN_R[0], 1,
1488  &fluxN_D[0], 1, &fluxJumps[0], 1);
1489 
1490  // Multiplicate the flux jumps for the correction
1491  // function
1492  for (i = 0; i < nquad1; ++i)
1493  {
1494  for (j = 0; j < nquad0; ++j)
1495  {
1496  cnt = (nquad0)*i + j;
1497  divCFluxE1[cnt] =
1498  fluxJumps[i] * m_dGR_xi1[n][j];
1499  }
1500  }
1501  break;
1502  case 2:
1503 
1504  // Extract the edge values of transformed flux-y on
1505  // edge e and order them accordingly to the order of
1506  // the trace space
1507 
1508  fields[0]->GetExp(n)->GetEdgePhysVals(
1509  e, elmtToTrace[n][e],
1510  fluxX2 + phys_offset,
1511  auxArray1 = fluxN_D);
1512 
1513  // Extract the physical Rieamann flux at each edge
1514  Vmath::Vcopy(nEdgePts,
1515  &numericalFlux[trace_offset], 1,
1516  &fluxN[0], 1);
1517 
1518  // Check the ordering of vectors
1519  if (fields[0]->GetExp(n)->GetEorient(e) ==
1520  StdRegions::eBackwards)
1521  {
1522  Vmath::Reverse(nEdgePts,
1523  auxArray1 = fluxN, 1,
1524  auxArray2 = fluxN, 1);
1525 
1526  Vmath::Reverse(nEdgePts,
1527  auxArray1 = fluxN_D, 1,
1528  auxArray2 = fluxN_D, 1);
1529  }
1530 
1531  // Transform Riemann Fluxes in the standard element
1532  for (i = 0; i < nquad0; ++i)
1533  {
1534  // Multiply Riemann Flux by Q factors
1535  fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
1536  }
1537 
1538  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1539  -1.0 : 1.0;
1540 
1541  for (i = 0; i < nEdgePts; ++i)
1542  {
1543  if (m_traceNormals[0][trace_offset+i]
1544  != fac*normals[0][i] ||
1545  m_traceNormals[1][trace_offset+i]
1546  != fac*normals[1][i])
1547  {
1548  fluxN_R[i] = -fluxN_R[i];
1549  }
1550  }
1551 
1552  // Subtract to the Riemann flux the discontinuous
1553  // flux
1554 
1555  Vmath::Vsub(nEdgePts,
1556  &fluxN_R[0], 1,
1557  &fluxN_D[0], 1, &fluxJumps[0], 1);
1558 
1559  // Multiplicate the flux jumps for the correction
1560  // function
1561  for (i = 0; i < nquad0; ++i)
1562  {
1563  for (j = 0; j < nquad1; ++j)
1564  {
1565  cnt = (nquad0 - 1) + j*nquad0 - i;
1566  divCFluxE2[cnt] =
1567  fluxJumps[i] * m_dGR_xi2[n][j];
1568  }
1569  }
1570  break;
1571  case 3:
1572  // Extract the edge values of transformed flux-x on
1573  // edge e and order them accordingly to the order of
1574  // the trace space
1575 
1576  fields[0]->GetExp(n)->GetEdgePhysVals(
1577  e, elmtToTrace[n][e],
1578  fluxX1 + phys_offset,
1579  auxArray1 = fluxN_D);
1580  Vmath::Neg (nEdgePts, fluxN_D, 1);
1581 
1582  // Extract the physical Rieamann flux at each edge
1583  Vmath::Vcopy(nEdgePts,
1584  &numericalFlux[trace_offset], 1,
1585  &fluxN[0], 1);
1586 
1587  // Check the ordering of vectors
1588  if (fields[0]->GetExp(n)->GetEorient(e) ==
1589  StdRegions::eBackwards)
1590  {
1591  Vmath::Reverse(nEdgePts,
1592  auxArray1 = fluxN, 1,
1593  auxArray2 = fluxN, 1);
1594 
1595  Vmath::Reverse(nEdgePts,
1596  auxArray1 = fluxN_D, 1,
1597  auxArray2 = fluxN_D, 1);
1598  }
1599 
1600  // Transform Riemann Fluxes in the standard element
1601  for (i = 0; i < nquad1; ++i)
1602  {
1603  // Multiply Riemann Flux by Q factors
1604  fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
1605  }
1606 
1607  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1608  -1.0 : 1.0;
1609 
1610  for (i = 0; i < nEdgePts; ++i)
1611  {
1612  if (m_traceNormals[0][trace_offset+i]
1613  != fac*normals[0][i] ||
1614  m_traceNormals[1][trace_offset+i]
1615  != fac*normals[1][i])
1616  {
1617  fluxN_R[i] = -fluxN_R[i];
1618  }
1619  }
1620 
1621  // Subtract to the Riemann flux the discontinuous
1622  // flux
1623 
1624  Vmath::Vsub(nEdgePts,
1625  &fluxN_R[0], 1,
1626  &fluxN_D[0], 1, &fluxJumps[0], 1);
1627 
1628  // Multiplicate the flux jumps for the correction
1629  // function
1630  for (i = 0; i < nquad1; ++i)
1631  {
1632  for (j = 0; j < nquad0; ++j)
1633  {
1634  cnt = (nquad0*nquad1 - nquad0) + j
1635  - i*nquad0;
1636  divCFluxE3[cnt] =
1637  -fluxJumps[i] * m_dGL_xi1[n][j];
1638  }
1639  }
1640  break;
1641  default:
1642  ASSERTL0(false,"edge value (< 3) is out of range");
1643  break;
1644  }
1645  }
1646 
1647 
1648  // Sum each edge contribution
1649  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1650  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1651 
1652  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1653  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1654 
1655  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1656  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1657  }
1658  }
1659 
1660 
1661  /**
1662  * @brief Compute the divergence of the corrective flux for 3D problems.
1663  *
1664  * @param nConvectiveFields Number of fields.
1665  * @param fields Pointer to fields.
1666  * @param fluxX1 X1-volumetric flux in physical space.
1667  * @param fluxX2 X2-volumetric flux in physical space.
1668  * @param fluxX3 X3-volumetric flux in physical space.
1669  * @param numericalFlux Interface flux in physical space.
1670  * @param divCFlux Divergence of the corrective flux.
1671  *
1672  * \todo: To be implemented. Switch on shapes eventually here.
1673  */
1674  void AdvectionFR::v_DivCFlux_3D(
1675  const int nConvectiveFields,
1676  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1677  const Array<OneD, const NekDouble> &fluxX1,
1678  const Array<OneD, const NekDouble> &fluxX2,
1679  const Array<OneD, const NekDouble> &fluxX3,
1680  const Array<OneD, const NekDouble> &numericalFlux,
1681  Array<OneD, NekDouble> &divCFlux)
1682  {
1683 
1684  }
1685 
1686  }
1687 }
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:161
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:220
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:134
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:227
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:47
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:1071
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:199
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:2146
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:382
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:1674
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:1303
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:1114
Nektar::SolverUtils::AdvectionFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: AdvectionFR.h:67
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)
Compute the advection term at each time-step using the Flux Reconstruction approach (FR)...
Definition: AdvectionFR.cpp:801
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:329
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:523
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:1179
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:1011
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:132
Expansion1D.h
AssemblyMapDG.h
Nektar::LibUtilities::BasisSharedPtr
boost::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:79
Nektar::SolverUtils::Advection::m_spaceDim
int m_spaceDim
Storage for space dimension. Used for homogeneous extension.
Definition: Advection.h:136
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
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:97
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:285
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:169
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 Wed May 4 2016 17:54:10 for Nektar++ by   doxygen 1.8.8