Nektar++
DiffusionLFR.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: DiffusionLFR.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: LFR diffusion class.
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <LocalRegions/Expansion1D.h>
36 #include <LocalRegions/Expansion2D.h>
37 #include <SolverUtils/Diffusion/DiffusionLFR.h>
38 #include <LibUtilities/Polylib/Polylib.h>
39 #include <MultiRegions/DisContField.h>
40 #include <boost/core/ignore_unused.hpp>
41 #include <boost/math/special_functions/gamma.hpp>
42 #include <iostream>
43 #include <iomanip>
44 
45 using namespace std;
46 
47 namespace Nektar
48 {
49  namespace SolverUtils
50  {
51  std::string DiffusionLFR::type[] = {
52  GetDiffusionFactory().RegisterCreatorFunction(
53  "LFRDG", DiffusionLFR::create),
54  GetDiffusionFactory().RegisterCreatorFunction(
55  "LFRSD", DiffusionLFR::create),
56  GetDiffusionFactory().RegisterCreatorFunction(
57  "LFRHU", DiffusionLFR::create),
58  GetDiffusionFactory().RegisterCreatorFunction(
59  "LFRcmin", DiffusionLFR::create),
60  GetDiffusionFactory().RegisterCreatorFunction(
61  "LFRcinf", DiffusionLFR::create)};
62 
63  /**
64  * @brief DiffusionLFR uses the Flux Reconstruction (FR) approach to
65  * compute the diffusion term. The implementation is only for segments,
66  * quadrilaterals and hexahedra at the moment.
67  *
68  * \todo Extension to triangles, tetrahedra and other shapes.
69  * (Long term objective)
70  */
71  DiffusionLFR::DiffusionLFR(std::string diffType):m_diffType(diffType)
72  {
73  }
74 
75  /**
76  * @brief Initiliase DiffusionLFR objects and store them before starting
77  * the time-stepping.
78  *
79  * This routine calls the virtual functions #v_SetupMetrics and
80  * #v_SetupCFunctions to initialise the objects needed by DiffusionLFR.
81  *
82  * @param pSession Pointer to session reader.
83  * @param pFields Pointer to fields.
84  */
85  void DiffusionLFR::v_InitObject(
86  LibUtilities::SessionReaderSharedPtr pSession,
87  Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
88  {
89  m_session = pSession;
90  v_SetupMetrics(pSession, pFields);
91  v_SetupCFunctions(pSession, pFields);
92 
93  // Initialising arrays
94  int i, j;
95  int nConvectiveFields = pFields.size();
96  int nDim = pFields[0]->GetCoordim(0);
97  int nSolutionPts = pFields[0]->GetTotPoints();
98  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
99 
100  m_IF1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
101  nConvectiveFields);
102  m_DU1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
103  nConvectiveFields);
104  m_DFC1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
105  nConvectiveFields);
106  m_BD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
107  nConvectiveFields);
108  m_D1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
109  nConvectiveFields);
110  m_DD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
111  nConvectiveFields);
112  m_IF2 = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
113  m_DFC2 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
114  nConvectiveFields);
115  m_divFD = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
116  m_divFC = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
117 
118  m_tmp1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
119  nConvectiveFields);
120  m_tmp2 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
121  nConvectiveFields);
122 
123 
124 
125  for (i = 0; i < nConvectiveFields; ++i)
126  {
127  m_IF1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
128  m_DU1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
129  m_DFC1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
130  m_BD1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
131  m_D1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
132  m_DD1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
133  m_IF2[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
134  m_DFC2[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
135  m_divFD[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
136  m_divFC[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
137 
138  m_tmp1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
139  m_tmp2[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
140 
141  for (j = 0; j < nDim; ++j)
142  {
143  m_IF1[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
144  m_DU1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
145  m_DFC1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
146  m_BD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
147  m_D1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
148  m_DD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
149  m_DFC2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
150 
151  m_tmp1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
152  m_tmp2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
153  }
154  }
155 
156  }
157 
158  /**
159  * @brief Setup the metric terms to compute the contravariant
160  * fluxes. (i.e. this special metric terms transform the fluxes
161  * at the interfaces of each element from the physical space to
162  * the standard space).
163  *
164  * This routine calls the function #GetTraceQFactors to compute and
165  * store the metric factors following an anticlockwise conventions
166  * along the edges/faces of each element. Note: for 1D problem
167  * the transformation is not needed.
168  *
169  * @param pSession Pointer to session reader.
170  * @param pFields Pointer to fields.
171  *
172  * \todo Add the metric terms for 3D Hexahedra.
173  */
174  void DiffusionLFR::v_SetupMetrics(
175  LibUtilities::SessionReaderSharedPtr pSession,
176  Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
177  {
178  boost::ignore_unused(pSession);
179 
180  int i, n;
181  int nquad0, nquad1;
182  int phys_offset;
183  int nLocalSolutionPts;
184  int nElements = pFields[0]->GetExpSize();
185  int nDimensions = pFields[0]->GetCoordim(0);
186  int nSolutionPts = pFields[0]->GetTotPoints();
187  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
188 
189  m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(nDimensions);
190  for (i = 0; i < nDimensions; ++i)
191  {
192  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
193  }
194  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
195 
196  Array<TwoD, const NekDouble> gmat;
197  Array<OneD, const NekDouble> jac;
198 
199  m_jac = Array<OneD, NekDouble>(nSolutionPts);
200 
201  Array<OneD, NekDouble> auxArray1;
202  Array<OneD, LibUtilities::BasisSharedPtr> base;
203  LibUtilities::PointsKeyVector ptsKeys;
204 
205  switch (nDimensions)
206  {
207  case 1:
208  {
209  for (n = 0; n < nElements; ++n)
210  {
211  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
212  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
213  phys_offset = pFields[0]->GetPhys_Offset(n);
214  jac = pFields[0]->GetExp(n)
215  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
216  ->GetMetricInfo()->GetJac(ptsKeys);
217  for (i = 0; i < nLocalSolutionPts; ++i)
218  {
219  m_jac[i+phys_offset] = jac[0];
220  }
221  }
222  break;
223  }
224  case 2:
225  {
226  m_gmat = Array<OneD, Array<OneD, NekDouble> >(4);
227  m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
228  m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
229  m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
230  m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
231 
232  m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble> >(nElements);
233  m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
234  m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
235  m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
236 
237  for (n = 0; n < nElements; ++n)
238  {
239  base = pFields[0]->GetExp(n)->GetBase();
240  nquad0 = base[0]->GetNumPoints();
241  nquad1 = base[1]->GetNumPoints();
242 
243  m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
244  m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
245  m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
246  m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
247 
248  // Extract the Q factors at each edge point
249  pFields[0]->GetExp(n)->GetTraceQFactors(
250  0, auxArray1 = m_Q2D_e0[n]);
251  pFields[0]->GetExp(n)->GetTraceQFactors(
252  1, auxArray1 = m_Q2D_e1[n]);
253  pFields[0]->GetExp(n)->GetTraceQFactors(
254  2, auxArray1 = m_Q2D_e2[n]);
255  pFields[0]->GetExp(n)->GetTraceQFactors(
256  3, auxArray1 = m_Q2D_e3[n]);
257 
258  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
259  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
260  phys_offset = pFields[0]->GetPhys_Offset(n);
261 
262  jac = pFields[0]->GetExp(n)
263  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
264  ->GetMetricInfo()->GetJac(ptsKeys);
265  gmat = pFields[0]->GetExp(n)
266  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
267  ->GetMetricInfo()->GetDerivFactors(ptsKeys);
268 
269  if (pFields[0]->GetExp(n)
270  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
271  ->GetMetricInfo()->GetGtype()
272  == SpatialDomains::eDeformed)
273  {
274  for (i = 0; i < nLocalSolutionPts; ++i)
275  {
276  m_jac[i+phys_offset] = jac[i];
277  m_gmat[0][i+phys_offset] = gmat[0][i];
278  m_gmat[1][i+phys_offset] = gmat[1][i];
279  m_gmat[2][i+phys_offset] = gmat[2][i];
280  m_gmat[3][i+phys_offset] = gmat[3][i];
281  }
282  }
283  else
284  {
285  for (i = 0; i < nLocalSolutionPts; ++i)
286  {
287  m_jac[i+phys_offset] = jac[0];
288  m_gmat[0][i+phys_offset] = gmat[0][0];
289  m_gmat[1][i+phys_offset] = gmat[1][0];
290  m_gmat[2][i+phys_offset] = gmat[2][0];
291  m_gmat[3][i+phys_offset] = gmat[3][0];
292  }
293  }
294  }
295  break;
296  }
297  case 3:
298  {
299  ASSERTL0(false,"3DFR Metric terms not implemented yet");
300  break;
301  }
302  default:
303  {
304  ASSERTL0(false, "Expansion dimension not recognised");
305  break;
306  }
307  }
308  }
309 
310  /**
311  * @brief Setup the derivatives of the correction functions. For more
312  * details see J Sci Comput (2011) 47: 50–72
313  *
314  * This routine calls 3 different bases:
315  * #eDG_DG_Left - #eDG_DG_Left which recovers a nodal DG scheme,
316  * #eDG_SD_Left - #eDG_SD_Left which recovers the SD scheme,
317  * #eDG_HU_Left - #eDG_HU_Left which recovers the Huynh scheme.
318  * The values of the derivatives of the correction function are then
319  * stored into global variables and reused into the virtual functions
320  * #v_DivCFlux_1D, #v_DivCFlux_2D, #v_DivCFlux_3D to compute the
321  * the divergence of the correction flux for 1D, 2D or 3D problems.
322  *
323  * @param pSession Pointer to session reader.
324  * @param pFields Pointer to fields.
325  */
326  void DiffusionLFR::v_SetupCFunctions(
327  LibUtilities::SessionReaderSharedPtr pSession,
328  Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
329  {
330  boost::ignore_unused(pSession);
331 
332  int i, n;
333  NekDouble c0 = 0.0;
334  NekDouble c1 = 0.0;
335  NekDouble c2 = 0.0;
336  int nquad0, nquad1, nquad2;
337  int nmodes0, nmodes1, nmodes2;
338  Array<OneD, LibUtilities::BasisSharedPtr> base;
339 
340  int nElements = pFields[0]->GetExpSize();
341  int nDim = pFields[0]->GetCoordim(0);
342 
343  switch (nDim)
344  {
345  case 1:
346  {
347  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
348  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
349 
350  for (n = 0; n < nElements; ++n)
351  {
352  base = pFields[0]->GetExp(n)->GetBase();
353  nquad0 = base[0]->GetNumPoints();
354  nmodes0 = base[0]->GetNumModes();
355  Array<OneD, const NekDouble> z0;
356  Array<OneD, const NekDouble> w0;
357 
358  base[0]->GetZW(z0, w0);
359 
360  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
361  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
362 
363  // Auxiliary vectors to build up the auxiliary
364  // derivatives of the Legendre polynomials
365  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
366  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
367  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
368 
369  // Degree of the correction functions
370  int p0 = nmodes0 - 1;
371 
372  // Function sign to compute left correction function
373  NekDouble sign0 = pow(-1.0, p0);
374 
375  // Factorial factor to build the scheme
376  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
377  / (pow(2.0, p0)
378  * boost::math::tgamma(p0 + 1)
379  * boost::math::tgamma(p0 + 1));
380 
381  // Scalar parameter which recovers the FR schemes
382  if (m_diffType == "LFRDG")
383  {
384  c0 = 0.0;
385  }
386  else if (m_diffType == "LFRSD")
387  {
388  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
389  * (ap0 * boost::math::tgamma(p0 + 1))
390  * (ap0 * boost::math::tgamma(p0 + 1)));
391  }
392  else if (m_diffType == "LFRHU")
393  {
394  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
395  * (ap0 * boost::math::tgamma(p0 + 1))
396  * (ap0 * boost::math::tgamma(p0 + 1)));
397  }
398  else if (m_diffType == "LFRcmin")
399  {
400  c0 = -2.0 / ((2.0 * p0 + 1.0)
401  * (ap0 * boost::math::tgamma(p0 + 1))
402  * (ap0 * boost::math::tgamma(p0 + 1)));
403  }
404  else if (m_diffType == "LFRinf")
405  {
406  c0 = 10000000000000000.0;
407  }
408 
409  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
410  * (ap0 * boost::math::tgamma(p0 + 1))
411  * (ap0 * boost::math::tgamma(p0 + 1));
412 
413  NekDouble overeta0 = 1.0 / (1.0 + etap0);
414 
415  // Derivative of the Legendre polynomials
416  // dLp = derivative of the Legendre polynomial of order p
417  // dLpp = derivative of the Legendre polynomial of order p+1
418  // dLpm = derivative of the Legendre polynomial of order p-1
419  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
420  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
421  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
422 
423  // Building the DG_c_Left
424  for(i = 0; i < nquad0; ++i)
425  {
426  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
427  m_dGL_xi1[n][i] += dLpp0[i];
428  m_dGL_xi1[n][i] *= overeta0;
429  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
430  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
431  }
432 
433  // Building the DG_c_Right
434  for(i = 0; i < nquad0; ++i)
435  {
436  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
437  m_dGR_xi1[n][i] += dLpp0[i];
438  m_dGR_xi1[n][i] *= overeta0;
439  m_dGR_xi1[n][i] += dLp0[i];
440  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
441  }
442  }
443  break;
444  }
445  case 2:
446  {
447  LibUtilities::BasisSharedPtr BasisFR_Left0;
448  LibUtilities::BasisSharedPtr BasisFR_Right0;
449  LibUtilities::BasisSharedPtr BasisFR_Left1;
450  LibUtilities::BasisSharedPtr BasisFR_Right1;
451 
452  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
453  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
454  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
455  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
456 
457  for (n = 0; n < nElements; ++n)
458  {
459  base = pFields[0]->GetExp(n)->GetBase();
460  nquad0 = base[0]->GetNumPoints();
461  nquad1 = base[1]->GetNumPoints();
462  nmodes0 = base[0]->GetNumModes();
463  nmodes1 = base[1]->GetNumModes();
464 
465  Array<OneD, const NekDouble> z0;
466  Array<OneD, const NekDouble> w0;
467  Array<OneD, const NekDouble> z1;
468  Array<OneD, const NekDouble> w1;
469 
470  base[0]->GetZW(z0, w0);
471  base[1]->GetZW(z1, w1);
472 
473  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
474  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
475  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
476  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
477 
478  // Auxiliary vectors to build up the auxiliary
479  // derivatives of the Legendre polynomials
480  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
481  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
482  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
483  Array<OneD,NekDouble> dLp1 (nquad1, 0.0);
484  Array<OneD,NekDouble> dLpp1 (nquad1, 0.0);
485  Array<OneD,NekDouble> dLpm1 (nquad1, 0.0);
486 
487  // Degree of the correction functions
488  int p0 = nmodes0 - 1;
489  int p1 = nmodes1 - 1;
490 
491  // Function sign to compute left correction function
492  NekDouble sign0 = pow(-1.0, p0);
493  NekDouble sign1 = pow(-1.0, p1);
494 
495  // Factorial factor to build the scheme
496  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
497  / (pow(2.0, p0)
498  * boost::math::tgamma(p0 + 1)
499  * boost::math::tgamma(p0 + 1));
500 
501  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
502  / (pow(2.0, p1)
503  * boost::math::tgamma(p1 + 1)
504  * boost::math::tgamma(p1 + 1));
505 
506  // Scalar parameter which recovers the FR schemes
507  if (m_diffType == "LFRDG")
508  {
509  c0 = 0.0;
510  c1 = 0.0;
511  }
512  else if (m_diffType == "LFRSD")
513  {
514  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
515  * (ap0 * boost::math::tgamma(p0 + 1))
516  * (ap0 * boost::math::tgamma(p0 + 1)));
517 
518  c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
519  * (ap1 * boost::math::tgamma(p1 + 1))
520  * (ap1 * boost::math::tgamma(p1 + 1)));
521  }
522  else if (m_diffType == "LFRHU")
523  {
524  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
525  * (ap0 * boost::math::tgamma(p0 + 1))
526  * (ap0 * boost::math::tgamma(p0 + 1)));
527 
528  c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
529  * (ap1 * boost::math::tgamma(p1 + 1))
530  * (ap1 * boost::math::tgamma(p1 + 1)));
531  }
532  else if (m_diffType == "LFRcmin")
533  {
534  c0 = -2.0 / ((2.0 * p0 + 1.0)
535  * (ap0 * boost::math::tgamma(p0 + 1))
536  * (ap0 * boost::math::tgamma(p0 + 1)));
537 
538  c1 = -2.0 / ((2.0 * p1 + 1.0)
539  * (ap1 * boost::math::tgamma(p1 + 1))
540  * (ap1 * boost::math::tgamma(p1 + 1)));
541  }
542  else if (m_diffType == "LFRinf")
543  {
544  c0 = 10000000000000000.0;
545  c1 = 10000000000000000.0;
546  }
547 
548  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
549  * (ap0 * boost::math::tgamma(p0 + 1))
550  * (ap0 * boost::math::tgamma(p0 + 1));
551 
552  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
553  * (ap1 * boost::math::tgamma(p1 + 1))
554  * (ap1 * boost::math::tgamma(p1 + 1));
555 
556  NekDouble overeta0 = 1.0 / (1.0 + etap0);
557  NekDouble overeta1 = 1.0 / (1.0 + etap1);
558 
559  // Derivative of the Legendre polynomials
560  // dLp = derivative of the Legendre polynomial of order p
561  // dLpp = derivative of the Legendre polynomial of order p+1
562  // dLpm = derivative of the Legendre polynomial of order p-1
563  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
564  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
565  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
566 
567  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
568  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1+1, 0.0, 0.0);
569  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1-1, 0.0, 0.0);
570 
571  // Building the DG_c_Left
572  for(i = 0; i < nquad0; ++i)
573  {
574  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
575  m_dGL_xi1[n][i] += dLpp0[i];
576  m_dGL_xi1[n][i] *= overeta0;
577  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
578  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
579  }
580 
581  // Building the DG_c_Left
582  for(i = 0; i < nquad1; ++i)
583  {
584  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
585  m_dGL_xi2[n][i] += dLpp1[i];
586  m_dGL_xi2[n][i] *= overeta1;
587  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
588  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
589  }
590 
591  // Building the DG_c_Right
592  for(i = 0; i < nquad0; ++i)
593  {
594  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
595  m_dGR_xi1[n][i] += dLpp0[i];
596  m_dGR_xi1[n][i] *= overeta0;
597  m_dGR_xi1[n][i] += dLp0[i];
598  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
599  }
600 
601  // Building the DG_c_Right
602  for(i = 0; i < nquad1; ++i)
603  {
604  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
605  m_dGR_xi2[n][i] += dLpp1[i];
606  m_dGR_xi2[n][i] *= overeta1;
607  m_dGR_xi2[n][i] += dLp1[i];
608  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
609  }
610  }
611  break;
612  }
613  case 3:
614  {
615  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
616  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
617  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
618  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
619  m_dGL_xi3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
620  m_dGR_xi3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
621 
622  for (n = 0; n < nElements; ++n)
623  {
624  base = pFields[0]->GetExp(n)->GetBase();
625  nquad0 = base[0]->GetNumPoints();
626  nquad1 = base[1]->GetNumPoints();
627  nquad2 = base[2]->GetNumPoints();
628  nmodes0 = base[0]->GetNumModes();
629  nmodes1 = base[1]->GetNumModes();
630  nmodes2 = base[2]->GetNumModes();
631 
632  Array<OneD, const NekDouble> z0;
633  Array<OneD, const NekDouble> w0;
634  Array<OneD, const NekDouble> z1;
635  Array<OneD, const NekDouble> w1;
636  Array<OneD, const NekDouble> z2;
637  Array<OneD, const NekDouble> w2;
638 
639  base[0]->GetZW(z0, w0);
640  base[1]->GetZW(z1, w1);
641  base[1]->GetZW(z2, w2);
642 
643  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
644  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
645  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
646  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
647  m_dGL_xi3[n] = Array<OneD, NekDouble>(nquad2);
648  m_dGR_xi3[n] = Array<OneD, NekDouble>(nquad2);
649 
650  // Auxiliary vectors to build up the auxiliary
651  // derivatives of the Legendre polynomials
652  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
653  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
654  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
655  Array<OneD,NekDouble> dLp1 (nquad1, 0.0);
656  Array<OneD,NekDouble> dLpp1 (nquad1, 0.0);
657  Array<OneD,NekDouble> dLpm1 (nquad1, 0.0);
658  Array<OneD,NekDouble> dLp2 (nquad2, 0.0);
659  Array<OneD,NekDouble> dLpp2 (nquad2, 0.0);
660  Array<OneD,NekDouble> dLpm2 (nquad2, 0.0);
661 
662  // Degree of the correction functions
663  int p0 = nmodes0 - 1;
664  int p1 = nmodes1 - 1;
665  int p2 = nmodes2 - 1;
666 
667  // Function sign to compute left correction function
668  NekDouble sign0 = pow(-1.0, p0);
669  NekDouble sign1 = pow(-1.0, p1);
670  NekDouble sign2 = pow(-1.0, p2);
671 
672  // Factorial factor to build the scheme
673  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
674  / (pow(2.0, p0)
675  * boost::math::tgamma(p0 + 1)
676  * boost::math::tgamma(p0 + 1));
677 
678  // Factorial factor to build the scheme
679  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
680  / (pow(2.0, p1)
681  * boost::math::tgamma(p1 + 1)
682  * boost::math::tgamma(p1 + 1));
683 
684  // Factorial factor to build the scheme
685  NekDouble ap2 = boost::math::tgamma(2 * p2 + 1)
686  / (pow(2.0, p2)
687  * boost::math::tgamma(p2 + 1)
688  * boost::math::tgamma(p2 + 1));
689 
690  // Scalar parameter which recovers the FR schemes
691  if (m_diffType == "LFRDG")
692  {
693  c0 = 0.0;
694  c1 = 0.0;
695  c2 = 0.0;
696  }
697  else if (m_diffType == "LFRSD")
698  {
699  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
700  * (ap0 * boost::math::tgamma(p0 + 1))
701  * (ap0 * boost::math::tgamma(p0 + 1)));
702 
703  c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
704  * (ap1 * boost::math::tgamma(p1 + 1))
705  * (ap1 * boost::math::tgamma(p1 + 1)));
706 
707  c2 = 2.0 * p2 / ((2.0 * p2 + 1.0) * (p2 + 1.0)
708  * (ap2 * boost::math::tgamma(p2 + 1))
709  * (ap2 * boost::math::tgamma(p2 + 1)));
710  }
711  else if (m_diffType == "LFRHU")
712  {
713  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
714  * (ap0 * boost::math::tgamma(p0 + 1))
715  * (ap0 * boost::math::tgamma(p0 + 1)));
716 
717  c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
718  * (ap1 * boost::math::tgamma(p1 + 1))
719  * (ap1 * boost::math::tgamma(p1 + 1)));
720 
721  c2 = 2.0 * (p2 + 1.0) / ((2.0 * p2 + 1.0) * p2
722  * (ap2 * boost::math::tgamma(p2 + 1))
723  * (ap2 * boost::math::tgamma(p2 + 1)));
724  }
725  else if (m_diffType == "LFRcmin")
726  {
727  c0 = -2.0 / ((2.0 * p0 + 1.0)
728  * (ap0 * boost::math::tgamma(p0 + 1))
729  * (ap0 * boost::math::tgamma(p0 + 1)));
730 
731  c1 = -2.0 / ((2.0 * p1 + 1.0)
732  * (ap1 * boost::math::tgamma(p1 + 1))
733  * (ap1 * boost::math::tgamma(p1 + 1)));
734 
735  c2 = -2.0 / ((2.0 * p2 + 1.0)
736  * (ap2 * boost::math::tgamma(p2 + 1))
737  * (ap2 * boost::math::tgamma(p2 + 1)));
738  }
739  else if (m_diffType == "LFRinf")
740  {
741  c0 = 10000000000000000.0;
742  c1 = 10000000000000000.0;
743  c2 = 10000000000000000.0;
744  }
745 
746  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
747  * (ap0 * boost::math::tgamma(p0 + 1))
748  * (ap0 * boost::math::tgamma(p0 + 1));
749 
750  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
751  * (ap1 * boost::math::tgamma(p1 + 1))
752  * (ap1 * boost::math::tgamma(p1 + 1));
753 
754  NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0)
755  * (ap2 * boost::math::tgamma(p2 + 1))
756  * (ap2 * boost::math::tgamma(p2 + 1));
757 
758  NekDouble overeta0 = 1.0 / (1.0 + etap0);
759  NekDouble overeta1 = 1.0 / (1.0 + etap1);
760  NekDouble overeta2 = 1.0 / (1.0 + etap2);
761 
762  // Derivative of the Legendre polynomials
763  // dLp = derivative of the Legendre polynomial of order p
764  // dLpp = derivative of the Legendre polynomial of order p+1
765  // dLpm = derivative of the Legendre polynomial of order p-1
766  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
767  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
768  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
769 
770  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
771  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1+1, 0.0, 0.0);
772  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1-1, 0.0, 0.0);
773 
774  Polylib::jacobd(nquad2, z2.data(), &(dLp2[0]), p2, 0.0, 0.0);
775  Polylib::jacobd(nquad2, z2.data(), &(dLpp2[0]), p2+1, 0.0, 0.0);
776  Polylib::jacobd(nquad2, z2.data(), &(dLpm2[0]), p2-1, 0.0, 0.0);
777 
778  // Building the DG_c_Left
779  for(i = 0; i < nquad0; ++i)
780  {
781  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
782  m_dGL_xi1[n][i] += dLpp0[i];
783  m_dGL_xi1[n][i] *= overeta0;
784  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
785  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
786  }
787 
788  // Building the DG_c_Left
789  for(i = 0; i < nquad1; ++i)
790  {
791  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
792  m_dGL_xi2[n][i] += dLpp1[i];
793  m_dGL_xi2[n][i] *= overeta1;
794  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
795  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
796  }
797 
798  // Building the DG_c_Left
799  for(i = 0; i < nquad2; ++i)
800  {
801  m_dGL_xi3[n][i] = etap2 * dLpm2[i];
802  m_dGL_xi3[n][i] += dLpp2[i];
803  m_dGL_xi3[n][i] *= overeta2;
804  m_dGL_xi3[n][i] = dLp2[i] - m_dGL_xi3[n][i];
805  m_dGL_xi3[n][i] = 0.5 * sign2 * m_dGL_xi3[n][i];
806  }
807 
808  // Building the DG_c_Right
809  for(i = 0; i < nquad0; ++i)
810  {
811  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
812  m_dGR_xi1[n][i] += dLpp0[i];
813  m_dGR_xi1[n][i] *= overeta0;
814  m_dGR_xi1[n][i] += dLp0[i];
815  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
816  }
817 
818  // Building the DG_c_Right
819  for(i = 0; i < nquad1; ++i)
820  {
821  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
822  m_dGR_xi2[n][i] += dLpp1[i];
823  m_dGR_xi2[n][i] *= overeta1;
824  m_dGR_xi2[n][i] += dLp1[i];
825  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
826  }
827 
828  // Building the DG_c_Right
829  for(i = 0; i < nquad2; ++i)
830  {
831  m_dGR_xi3[n][i] = etap2 * dLpm2[i];
832  m_dGR_xi3[n][i] += dLpp2[i];
833  m_dGR_xi3[n][i] *= overeta2;
834  m_dGR_xi3[n][i] += dLp2[i];
835  m_dGR_xi3[n][i] = 0.5 * m_dGR_xi3[n][i];
836  }
837  }
838  break;
839  }
840  default:
841  {
842  ASSERTL0(false,"Expansion dimension not recognised");
843  break;
844  }
845  }
846  }
847 
848  /**
849  * @brief Calculate FR Diffusion for the linear problems
850  * using an LDG interface flux.
851  *
852  */
853  void DiffusionLFR::v_Diffuse(
854  const std::size_t nConvectiveFields,
855  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
856  const Array<OneD, Array<OneD, NekDouble> > &inarray,
857  Array<OneD, Array<OneD, NekDouble> > &outarray,
858  const Array<OneD, Array<OneD, NekDouble> > &pFwd,
859  const Array<OneD, Array<OneD, NekDouble> > &pBwd)
860  {
861  boost::ignore_unused(pFwd, pBwd);
862 
863  int i, j, n;
864  int phys_offset;
865  //Array<TwoD, const NekDouble> gmat;
866  //Array<OneD, const NekDouble> jac;
867  Array<OneD, NekDouble> auxArray1, auxArray2;
868 
869  Array<OneD, LibUtilities::BasisSharedPtr> Basis;
870  Basis = fields[0]->GetExp(0)->GetBase();
871 
872  int nElements = fields[0]->GetExpSize();
873  int nDim = fields[0]->GetCoordim(0);
874  int nSolutionPts = fields[0]->GetTotPoints();
875  int nCoeffs = fields[0]->GetNcoeffs();
876 
877  Array<OneD, Array<OneD, NekDouble> > outarrayCoeff(nConvectiveFields);
878  for (i = 0; i < nConvectiveFields; ++i)
879  {
880  outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
881  }
882 
883  // Compute interface numerical fluxes for inarray in physical space
884  v_NumFluxforScalar(fields, inarray, m_IF1);
885 
886  switch(nDim)
887  {
888  // 1D problems
889  case 1:
890  {
891  for (i = 0; i < nConvectiveFields; ++i)
892  {
893  // Computing the physical first-order discountinuous
894  // derivative
895  for (n = 0; n < nElements; n++)
896  {
897  phys_offset = fields[0]->GetPhys_Offset(n);
898 
899  fields[i]->GetExp(n)->PhysDeriv(0,
900  auxArray1 = inarray[i] + phys_offset,
901  auxArray2 = m_DU1[i][0] + phys_offset);
902  }
903 
904  // Computing the standard first-order correction
905  // derivative
906  v_DerCFlux_1D(nConvectiveFields, fields, inarray[i],
907  m_IF1[i][0], m_DFC1[i][0]);
908 
909  // Back to the physical space using global operations
910  Vmath::Vdiv(nSolutionPts, &m_DFC1[i][0][0], 1,
911  &m_jac[0], 1, &m_DFC1[i][0][0], 1);
912  Vmath::Vdiv(nSolutionPts, &m_DFC1[i][0][0], 1,
913  &m_jac[0], 1, &m_DFC1[i][0][0], 1);
914 
915  // Computing total first order derivatives
916  Vmath::Vadd(nSolutionPts, &m_DFC1[i][0][0], 1,
917  &m_DU1[i][0][0], 1, &m_D1[i][0][0], 1);
918 
919  Vmath::Vcopy(nSolutionPts, &m_D1[i][0][0], 1,
920  &m_tmp1[i][0][0], 1);
921  }
922 
923  // Computing interface numerical fluxes for m_D1
924  // in physical space
925  v_NumFluxforVector(fields, inarray, m_tmp1, m_IF2);
926 
927  for (i = 0; i < nConvectiveFields; ++i)
928  {
929  // Computing the physical second-order discountinuous
930  // derivative
931  for (n = 0; n < nElements; n++)
932  {
933  phys_offset = fields[0]->GetPhys_Offset(n);
934 
935  fields[i]->GetExp(n)->PhysDeriv(0,
936  auxArray1 = m_D1[i][0] + phys_offset,
937  auxArray2 = m_DD1[i][0] + phys_offset);
938  }
939 
940  // Computing the standard second-order correction
941  // derivative
942  v_DerCFlux_1D(nConvectiveFields, fields, m_D1[i][0],
943  m_IF2[i], m_DFC2[i][0]);
944 
945  // Back to the physical space using global operations
946  Vmath::Vdiv(nSolutionPts, &m_DFC2[i][0][0], 1,
947  &m_jac[0], 1, &m_DFC2[i][0][0], 1);
948  Vmath::Vdiv(nSolutionPts, &m_DFC2[i][0][0], 1,
949  &m_jac[0], 1, &m_DFC2[i][0][0], 1);
950 
951  // Adding the total divergence to outarray (RHS)
952  Vmath::Vadd(nSolutionPts, &m_DFC2[i][0][0], 1,
953  &m_DD1[i][0][0], 1, &outarray[i][0], 1);
954 
955  // Primitive Dealiasing 1D
956  if (!(Basis[0]->Collocation()))
957  {
958  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
959  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
960  }
961  }
962  break;
963  }
964  // 2D problems
965  case 2:
966  {
967  for(i = 0; i < nConvectiveFields; ++i)
968  {
969  for (j = 0; j < nDim; ++j)
970  {
971  // Temporary vectors
972  Array<OneD, NekDouble> u1_hat(nSolutionPts, 0.0);
973  Array<OneD, NekDouble> u2_hat(nSolutionPts, 0.0);
974 
975  if (j == 0)
976  {
977  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
978  &m_gmat[0][0], 1, &u1_hat[0], 1);
979 
980  Vmath::Vmul(nSolutionPts, &u1_hat[0], 1,
981  &m_jac[0], 1, &u1_hat[0], 1);
982 
983  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
984  &m_gmat[1][0], 1, &u2_hat[0], 1);
985 
986  Vmath::Vmul(nSolutionPts, &u2_hat[0], 1,
987  &m_jac[0], 1, &u2_hat[0], 1);
988  }
989  else if (j == 1)
990  {
991  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
992  &m_gmat[2][0], 1, &u1_hat[0], 1);
993 
994  Vmath::Vmul(nSolutionPts, &u1_hat[0], 1,
995  &m_jac[0], 1, &u1_hat[0], 1);
996 
997  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
998  &m_gmat[3][0], 1, &u2_hat[0], 1);
999 
1000  Vmath::Vmul(nSolutionPts, &u2_hat[0], 1,
1001  &m_jac[0], 1, &u2_hat[0], 1);
1002  }
1003 
1004  for (n = 0; n < nElements; n++)
1005  {
1006  phys_offset = fields[0]->GetPhys_Offset(n);
1007 
1008  fields[i]->GetExp(n)->StdPhysDeriv(0,
1009  auxArray1 = u1_hat + phys_offset,
1010  auxArray2 = m_tmp1[i][j] + phys_offset);
1011 
1012  fields[i]->GetExp(n)->StdPhysDeriv(1,
1013  auxArray1 = u2_hat + phys_offset,
1014  auxArray2 = m_tmp2[i][j] + phys_offset);
1015  }
1016 
1017  Vmath::Vadd(nSolutionPts, &m_tmp1[i][j][0], 1,
1018  &m_tmp2[i][j][0], 1,
1019  &m_DU1[i][j][0], 1);
1020 
1021  // Divide by the metric jacobian
1022  Vmath::Vdiv(nSolutionPts, &m_DU1[i][j][0], 1,
1023  &m_jac[0], 1, &m_DU1[i][j][0], 1);
1024 
1025  // Computing the standard first-order correction
1026  // derivatives
1027  v_DerCFlux_2D(nConvectiveFields, j, fields,
1028  inarray[i], m_IF1[i][j],
1029  m_DFC1[i][j]);
1030  }
1031 
1032  // Multiplying derivatives by B matrix to get auxiliary
1033  // variables
1034  for (j = 0; j < nSolutionPts; ++j)
1035  {
1036  // std(q1)
1037  m_BD1[i][0][j] =
1038  (m_gmat[0][j]*m_gmat[0][j] +
1039  m_gmat[2][j]*m_gmat[2][j]) *
1040  m_DFC1[i][0][j] +
1041  (m_gmat[1][j]*m_gmat[0][j] +
1042  m_gmat[3][j]*m_gmat[2][j]) *
1043  m_DFC1[i][1][j];
1044 
1045  // std(q2)
1046  m_BD1[i][1][j] =
1047  (m_gmat[1][j]*m_gmat[0][j] +
1048  m_gmat[3][j]*m_gmat[2][j]) *
1049  m_DFC1[i][0][j] +
1050  (m_gmat[1][j]*m_gmat[1][j] +
1051  m_gmat[3][j]*m_gmat[3][j]) *
1052  m_DFC1[i][1][j];
1053  }
1054 
1055  // Multiplying derivatives by A^(-1) to get back
1056  // into the physical space
1057  for (j = 0; j < nSolutionPts; j++)
1058  {
1059  // q1 = A11^(-1)*std(q1) + A12^(-1)*std(q2)
1060  m_DFC1[i][0][j] =
1061  m_gmat[3][j] * m_BD1[i][0][j] -
1062  m_gmat[2][j] * m_BD1[i][1][j];
1063 
1064  // q2 = A21^(-1)*std(q1) + A22^(-1)*std(q2)
1065  m_DFC1[i][1][j] =
1066  m_gmat[0][j] * m_BD1[i][1][j] -
1067  m_gmat[1][j] * m_BD1[i][0][j];
1068  }
1069 
1070  // Computing the physical first-order derivatives
1071  for (j = 0; j < nDim; ++j)
1072  {
1073  Vmath::Vadd(nSolutionPts, &m_DU1[i][j][0], 1,
1074  &m_DFC1[i][j][0], 1,
1075  &m_D1[i][j][0], 1);
1076  }
1077  }
1078 
1079  // Computing interface numerical fluxes for m_D1
1080  // in physical space
1081  v_NumFluxforVector(fields, inarray, m_D1, m_IF2);
1082 
1083  // Computing the standard second-order discontinuous
1084  // derivatives
1085  for (i = 0; i < nConvectiveFields; ++i)
1086  {
1087  Array<OneD, NekDouble> f_hat(nSolutionPts, 0.0);
1088  Array<OneD, NekDouble> g_hat(nSolutionPts, 0.0);
1089 
1090  for (j = 0; j < nSolutionPts; j++)
1091  {
1092  f_hat[j] = (m_D1[i][0][j] * m_gmat[0][j] +
1093  m_D1[i][1][j] * m_gmat[2][j])*m_jac[j];
1094 
1095  g_hat[j] = (m_D1[i][0][j] * m_gmat[1][j] +
1096  m_D1[i][1][j] * m_gmat[3][j])*m_jac[j];
1097  }
1098 
1099  for (n = 0; n < nElements; n++)
1100  {
1101  phys_offset = fields[0]->GetPhys_Offset(n);
1102 
1103  fields[0]->GetExp(n)->StdPhysDeriv(0,
1104  auxArray1 = f_hat + phys_offset,
1105  auxArray2 = m_DD1[i][0] + phys_offset);
1106 
1107  fields[0]->GetExp(n)->StdPhysDeriv(1,
1108  auxArray1 = g_hat + phys_offset,
1109  auxArray2 = m_DD1[i][1] + phys_offset);
1110  }
1111 
1112  // Divergence of the standard discontinuous flux
1113  Vmath::Vadd(nSolutionPts,
1114  &m_DD1[i][0][0], 1,
1115  &m_DD1[i][1][0], 1,
1116  &m_divFD[i][0], 1);
1117 
1118  // Divergence of the standard correction flux
1119  if (Basis[0]->GetPointsType() ==
1120  LibUtilities::eGaussGaussLegendre &&
1121  Basis[1]->GetPointsType() ==
1122  LibUtilities::eGaussGaussLegendre)
1123  {
1124  v_DivCFlux_2D_Gauss(nConvectiveFields, fields,
1125  f_hat, g_hat,
1126  m_IF2[i], m_divFC[i]);
1127  }
1128  else
1129  {
1130  v_DivCFlux_2D(nConvectiveFields, fields,
1131  m_D1[i][0], m_D1[i][1],
1132  m_IF2[i], m_divFC[i]);
1133  }
1134 
1135  // Divergence of the standard final flux
1136  Vmath::Vadd(nSolutionPts, &m_divFD[i][0], 1,
1137  &m_divFC[i][0], 1, &outarray[i][0], 1);
1138 
1139  // Dividing by the jacobian to get back into
1140  // physical space
1141  Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1,
1142  &m_jac[0], 1, &outarray[i][0], 1);
1143 
1144  // Primitive Dealiasing 2D
1145  if (!(Basis[0]->Collocation()))
1146  {
1147  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1148  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1149  }
1150  }//Close loop on nConvectiveFields
1151  break;
1152  }
1153  // 3D problems
1154  case 3:
1155  {
1156  ASSERTL0(false, "3D FRDG case not implemented yet");
1157  break;
1158  }
1159  }
1160  }
1161 
1162  /**
1163  * @brief Builds the numerical flux for the 1st order derivatives
1164  *
1165  */
1166  void DiffusionLFR::v_NumFluxforScalar(
1167  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1168  const Array<OneD, Array<OneD, NekDouble> > &ufield,
1169  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &uflux)
1170  {
1171  int i, j;
1172  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1173  int nvariables = fields.size();
1174  int nDim = fields[0]->GetCoordim(0);
1175 
1176  Array<OneD, NekDouble > Fwd (nTracePts);
1177  Array<OneD, NekDouble > Bwd (nTracePts);
1178  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1179  Array<OneD, NekDouble > fluxtemp(nTracePts, 0.0);
1180 
1181  // Get the normal velocity Vn
1182  for (i = 0; i < nDim; ++i)
1183  {
1184  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
1185  Vn, 1, Vn, 1);
1186  }
1187 
1188  // Get the sign of (v \cdot n), v = an arbitrary vector
1189  // Evaluate upwind flux:
1190  // uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
1191  for (j = 0; j < nDim; ++j)
1192  {
1193  for (i = 0; i < nvariables ; ++i)
1194  {
1195  // Compute Fwd and Bwd value of ufield of i direction
1196  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1197 
1198  // if Vn >= 0, flux = uFwd, i.e.,
1199  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
1200  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
1201 
1202  // else if Vn < 0, flux = uBwd, i.e.,
1203  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
1204  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
1205 
1206  fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, fluxtemp);
1207 
1208  // Imposing weak boundary condition with flux
1209  // if Vn >= 0, uflux = uBwd at Neumann, i.e.,
1210  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
1211  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
1212 
1213  // if Vn >= 0, uflux = uFwd at Neumann, i.e.,
1214  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
1215  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
1216 
1217  if(fields[0]->GetBndCondExpansions().size())
1218  {
1219  v_WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
1220  }
1221 
1222  // if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
1223  // i.e,
1224  // edge::eForward, uFwd \‍(\tan_{\xi}^Fwd \cdot \vec{n})
1225  // edge::eBackward, uFwd \‍(\tan_{\xi}^Bwd \cdot \vec{n})
1226 
1227  // else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}),
1228  // i.e,
1229  // edge::eForward, uBwd \‍(\tan_{\xi}^Fwd \cdot \vec{n})
1230  // edge::eBackward, uBwd \‍(\tan_{\xi}^Bwd \cdot \vec{n})
1231 
1232  Vmath::Vcopy(nTracePts, fluxtemp, 1, uflux[i][j], 1);
1233  }
1234  }
1235  }
1236 
1237  /**
1238  * @brief Imposes appropriate bcs for the 1st order derivatives
1239  *
1240  */
1241  void DiffusionLFR::v_WeakPenaltyforScalar(
1242  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1243  const int var,
1244  const Array<OneD, const NekDouble> &ufield,
1245  Array<OneD, NekDouble> &penaltyflux)
1246  {
1247  // Variable initialisation
1248  int i, e, id1, id2;
1249  int nBndEdgePts, nBndEdges;
1250  int cnt = 0;
1251  int nBndRegions = fields[var]->GetBndCondExpansions().size();
1252  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1253 
1254  // Extract physical values of the fields at the trace space
1255  Array<OneD, NekDouble > uplus(nTracePts);
1256  fields[var]->ExtractTracePhys(ufield, uplus);
1257 
1258  // Impose boundary conditions
1259  for (i = 0; i < nBndRegions; ++i)
1260  {
1261  // Number of boundary edges related to bcRegion 'i'
1262  nBndEdges = fields[var]->
1263  GetBndCondExpansions()[i]->GetExpSize();
1264 
1265  // Weakly impose boundary conditions by modifying flux values
1266  for (e = 0; e < nBndEdges ; ++e)
1267  {
1268  // Number of points on the expansion
1269  nBndEdgePts = fields[var]->
1270  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
1271 
1272  // Offset of the boundary expansion
1273  id1 = fields[var]->
1274  GetBndCondExpansions()[i]->GetPhys_Offset(e);
1275 
1276  // Offset of the trace space related to boundary expansion
1277  id2 = fields[0]->GetTrace()->
1278  GetPhys_Offset(fields[0]->GetTraceMap()->
1279  GetBndCondIDToGlobalTraceID(cnt++));
1280 
1281  // Dirichlet bcs ==> uflux = gD
1282  if (fields[var]->GetBndConditions()[i]->
1283  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1284  {
1285  Vmath::Vcopy(nBndEdgePts,
1286  &(fields[var]->
1287  GetBndCondExpansions()[i]->
1288  GetPhys())[id1], 1,
1289  &penaltyflux[id2], 1);
1290  }
1291  // Neumann bcs ==> uflux = u+
1292  else if ((fields[var]->GetBndConditions()[i])->
1293  GetBoundaryConditionType() == SpatialDomains::eNeumann)
1294  {
1295  Vmath::Vcopy(nBndEdgePts,
1296  &uplus[id2], 1,
1297  &penaltyflux[id2], 1);
1298  }
1299  }
1300  }
1301  }
1302 
1303  /**
1304  * @brief Build the numerical flux for the 2nd order derivatives
1305  *
1306  */
1307  void DiffusionLFR::v_NumFluxforVector(
1308  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1309  const Array<OneD, Array<OneD, NekDouble> > &ufield,
1310  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &qfield,
1311  Array<OneD, Array<OneD, NekDouble> > &qflux)
1312  {
1313  boost::ignore_unused(ufield);
1314 
1315  int i, j;
1316  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1317  int nvariables = fields.size();
1318  int nDim = fields[0]->GetCoordim(0);
1319 
1320  NekDouble C11 = 0.0;
1321  Array<OneD, NekDouble > Fwd(nTracePts);
1322  Array<OneD, NekDouble > Bwd(nTracePts);
1323  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1324 
1325  Array<OneD, NekDouble > qFwd (nTracePts);
1326  Array<OneD, NekDouble > qBwd (nTracePts);
1327  Array<OneD, NekDouble > qfluxtemp(nTracePts, 0.0);
1328  Array<OneD, NekDouble > uterm(nTracePts);
1329 
1330  // Get the normal velocity Vn
1331  for(i = 0; i < nDim; ++i)
1332  {
1333  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
1334  Vn, 1, Vn, 1);
1335  }
1336 
1337  // Evaulate upwind flux:
1338  // qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
1339  for (i = 0; i < nvariables; ++i)
1340  {
1341  qflux[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
1342  for (j = 0; j < nDim; ++j)
1343  {
1344  // Compute Fwd and Bwd value of ufield of jth direction
1345  fields[i]->GetFwdBwdTracePhys(qfield[i][j], qFwd, qBwd);
1346 
1347  // if Vn >= 0, flux = uFwd, i.e.,
1348  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick
1349  // qflux = qBwd = q+
1350  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick
1351  // qflux = qBwd = q-
1352 
1353  // else if Vn < 0, flux = uBwd, i.e.,
1354  // edge::eForward, if V*n<0 <=> V*n_F<0, pick
1355  // qflux = qFwd = q-
1356  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
1357  // qflux = qFwd = q+
1358 
1359  fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1360 
1361  Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
1362  qfluxtemp, 1);
1363 
1364  /*
1365  // Generate Stability term = - C11 ( u- - u+ )
1366  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1367 
1368  Vmath::Vsub(nTracePts,
1369  Fwd, 1, Bwd, 1,
1370  uterm, 1);
1371 
1372  Vmath::Smul(nTracePts,
1373  -1.0 * C11, uterm, 1,
1374  uterm, 1);
1375 
1376  // Flux = {Fwd, Bwd} * (nx, ny, nz) + uterm * (nx, ny)
1377  Vmath::Vadd(nTracePts,
1378  uterm, 1,
1379  qfluxtemp, 1,
1380  qfluxtemp, 1);
1381  */
1382 
1383  // Imposing weak boundary condition with flux
1384  if (fields[0]->GetBndCondExpansions().size())
1385  {
1386  v_WeakPenaltyforVector(fields, i, j, qfield[i][j],
1387  qfluxtemp, C11);
1388  }
1389 
1390  // q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
1391  // n_xi = n_x * tan_xi_x + n_y * tan_xi_y + n_z * tan_xi_z
1392  // n_xi = n_x * tan_eta_x + n_y * tan_eta_y + n_z*tan_eta_z
1393  Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1,
1394  qflux[i], 1);
1395  }
1396  }
1397  }
1398 
1399 
1400 
1401  /**
1402  * @brief Imposes appropriate bcs for the 2nd order derivatives
1403  *
1404  */
1405  void DiffusionLFR::v_WeakPenaltyforVector(
1406  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1407  const int var,
1408  const int dir,
1409  const Array<OneD, const NekDouble> &qfield,
1410  Array<OneD, NekDouble> &penaltyflux,
1411  NekDouble C11)
1412  {
1413  boost::ignore_unused(C11);
1414 
1415  int i, e, id1, id2;
1416  int nBndEdges, nBndEdgePts;
1417  int nBndRegions = fields[var]->GetBndCondExpansions().size();
1418  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1419 
1420  Array<OneD, NekDouble > uterm(nTracePts);
1421  Array<OneD, NekDouble > qtemp(nTracePts);
1422  int cnt = 0;
1423 
1424  fields[var]->ExtractTracePhys(qfield, qtemp);
1425 
1426  for (i = 0; i < nBndRegions; ++i)
1427  {
1428  nBndEdges = fields[var]->
1429  GetBndCondExpansions()[i]->GetExpSize();
1430 
1431  // Weakly impose boundary conditions by modifying flux values
1432  for (e = 0; e < nBndEdges ; ++e)
1433  {
1434  nBndEdgePts = fields[var]->
1435  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
1436 
1437  id1 = fields[var]->
1438  GetBndCondExpansions()[i]->GetPhys_Offset(e);
1439 
1440  id2 = fields[0]->GetTrace()->
1441  GetPhys_Offset(fields[0]->GetTraceMap()->
1442  GetBndCondIDToGlobalTraceID(cnt++));
1443 
1444  // For Dirichlet boundary condition:
1445  //qflux = q+ - C_11 (u+ - g_D) (nx, ny)
1446  if (fields[var]->GetBndConditions()[i]->
1447  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1448  {
1449  Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1450  &qtemp[id2], 1, &penaltyflux[id2], 1);
1451  }
1452  // For Neumann boundary condition: qflux = g_N
1453  else if ((fields[var]->GetBndConditions()[i])->
1454  GetBoundaryConditionType() == SpatialDomains::eNeumann)
1455  {
1456  Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1457  &(fields[var]->GetBndCondExpansions()[i]->
1458  GetPhys())[id1], 1, &penaltyflux[id2], 1);
1459  }
1460  }
1461  }
1462  }
1463 
1464  /**
1465  * @brief Compute the derivative of the corrective flux for 1D problems.
1466  *
1467  * @param nConvectiveFields Number of fields.
1468  * @param fields Pointer to fields.
1469  * @param flux Volumetric flux in the physical space.
1470  * @param iFlux Numerical interface flux in physical space.
1471  * @param derCFlux Derivative of the corrective flux.
1472  *
1473  */
1474  void DiffusionLFR::v_DerCFlux_1D(
1475  const int nConvectiveFields,
1476  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1477  const Array<OneD, const NekDouble> &flux,
1478  const Array<OneD, const NekDouble> &iFlux,
1479  Array<OneD, NekDouble> &derCFlux)
1480  {
1481  boost::ignore_unused(nConvectiveFields);
1482 
1483  int n;
1484  int nLocalSolutionPts, phys_offset, t_offset;
1485 
1486  Array<OneD, NekDouble> auxArray1, auxArray2;
1487  Array<TwoD, const NekDouble> gmat;
1488  Array<OneD, const NekDouble> jac;
1489 
1490  LibUtilities::PointsKeyVector ptsKeys;
1491  LibUtilities::BasisSharedPtr Basis;
1492  Basis = fields[0]->GetExp(0)->GetBasis(0);
1493 
1494  int nElements = fields[0]->GetExpSize();
1495  int nSolutionPts = fields[0]->GetTotPoints();
1496 
1497  vector<bool> negatedFluxNormal =
1498  std::static_pointer_cast<MultiRegions::DisContField>(
1499  fields[0])->GetNegatedFluxNormal();
1500 
1501  // Arrays to store the derivatives of the correction flux
1502  Array<OneD, NekDouble> DCL(nSolutionPts/nElements, 0.0);
1503  Array<OneD, NekDouble> DCR(nSolutionPts/nElements, 0.0);
1504 
1505  // Arrays to store the intercell numerical flux jumps
1506  Array<OneD, NekDouble> JumpL(nElements);
1507  Array<OneD, NekDouble> JumpR(nElements);
1508 
1509  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1510  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1511 
1512  for (n = 0; n < nElements; ++n)
1513  {
1514  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1515  phys_offset = fields[0]->GetPhys_Offset(n);
1516 
1517  Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1518  Array<OneD, NekDouble> tmpFluxVertex(1,0.0);
1519  Vmath::Vcopy(nLocalSolutionPts,
1520  &flux[phys_offset], 1,
1521  &tmparrayX1[0], 1);
1522 
1523  fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0],
1524  tmparrayX1,
1525  tmpFluxVertex);
1526 
1527  t_offset = fields[0]->GetTrace()
1528  ->GetPhys_Offset(elmtToTrace[n][0]->GetElmtId());
1529 
1530  if(negatedFluxNormal[2*n])
1531  {
1532  JumpL[n] = iFlux[t_offset] - tmpFluxVertex[0];
1533  }
1534  else
1535  {
1536  JumpL[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1537  }
1538 
1539  t_offset = fields[0]->GetTrace()
1540  ->GetPhys_Offset(elmtToTrace[n][1]->GetElmtId());
1541 
1542  fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1],
1543  tmparrayX1,
1544  tmpFluxVertex);
1545  if(negatedFluxNormal[2*n+1])
1546  {
1547  JumpR[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1548  }
1549  else
1550  {
1551  JumpR[n] = iFlux[t_offset] - tmpFluxVertex[0];
1552  }
1553  }
1554 
1555  for (n = 0; n < nElements; ++n)
1556  {
1557  ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1558  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1559  phys_offset = fields[0]->GetPhys_Offset(n);
1560  jac = fields[0]->GetExp(n)
1561  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
1562  ->GetMetricInfo()->GetJac(ptsKeys);
1563 
1564  JumpL[n] = JumpL[n] * jac[0];
1565  JumpR[n] = JumpR[n] * jac[0];
1566 
1567  // Left jump multiplied by left derivative of C function
1568  Vmath::Smul(nLocalSolutionPts, JumpL[n],
1569  &m_dGL_xi1[n][0], 1, &DCL[0], 1);
1570 
1571  // Right jump multiplied by right derivative of C function
1572  Vmath::Smul(nLocalSolutionPts, JumpR[n],
1573  &m_dGR_xi1[n][0], 1, &DCR[0], 1);
1574 
1575  // Assembling derivative of the correction flux
1576  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1577  &derCFlux[phys_offset], 1);
1578  }
1579  }
1580 
1581  /**
1582  * @brief Compute the derivative of the corrective flux wrt a given
1583  * coordinate for 2D problems.
1584  *
1585  * There could be a bug for deformed elements since first derivatives
1586  * are not exactly the same results of DiffusionLDG as expected
1587  *
1588  * @param nConvectiveFields Number of fields.
1589  * @param fields Pointer to fields.
1590  * @param flux Volumetric flux in the physical space.
1591  * @param numericalFlux Numerical interface flux in physical space.
1592  * @param derCFlux Derivative of the corrective flux.
1593  *
1594  * \todo: Switch on shapes eventually here.
1595  */
1596  void DiffusionLFR::v_DerCFlux_2D(
1597  const int nConvectiveFields,
1598  const int direction,
1599  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1600  const Array<OneD, const NekDouble> &flux,
1601  const Array<OneD, NekDouble> &iFlux,
1602  Array<OneD, NekDouble> &derCFlux)
1603  {
1604  boost::ignore_unused(nConvectiveFields);
1605 
1606  int n, e, i, j, cnt;
1607 
1608  Array<TwoD, const NekDouble> gmat;
1609  Array<OneD, const NekDouble> jac;
1610 
1611  int nElements = fields[0]->GetExpSize();
1612 
1613  int trace_offset, phys_offset;
1614  int nLocalSolutionPts;
1615  int nquad0, nquad1;
1616  int nEdgePts;
1617 
1618  LibUtilities::PointsKeyVector ptsKeys;
1619  Array<OneD, NekDouble> auxArray1, auxArray2;
1620  Array<OneD, LibUtilities::BasisSharedPtr> base;
1621  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1622  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1623 
1624  // Loop on the elements
1625  for (n = 0; n < nElements; ++n)
1626  {
1627  // Offset of the element on the global vector
1628  phys_offset = fields[0]->GetPhys_Offset(n);
1629  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1630  ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1631 
1632  jac = fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1633  ->GetGeom2D()->GetMetricInfo()->GetJac(ptsKeys);
1634 
1635  base = fields[0]->GetExp(n)->GetBase();
1636  nquad0 = base[0]->GetNumPoints();
1637  nquad1 = base[1]->GetNumPoints();
1638 
1639  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1640  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1641  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1642  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1643 
1644  // Loop on the edges
1645  for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1646  {
1647  // Number of edge points of edge 'e'
1648  nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1649 
1650  // Array for storing volumetric fluxes on each edge
1651  Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1652 
1653  // Offset of the trace space correspondent to edge 'e'
1654  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1655  elmtToTrace[n][e]->GetElmtId());
1656 
1657  // Get the normals of edge 'e'
1658  //const Array<OneD, const Array<OneD, NekDouble> > &normals =
1659  //fields[0]->GetExp(n)->GetTraceNormal(e);
1660 
1661  // Extract the edge values of the volumetric fluxes
1662  // on edge 'e' and order them accordingly to the order
1663  // of the trace space
1664  fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1665  flux + phys_offset,
1666  auxArray1 = tmparray);
1667 
1668  // Compute the fluxJumps per each edge 'e' and each
1669  // flux point
1670  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1671  Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1,
1672  &tmparray[0], 1, &fluxJumps[0], 1);
1673 
1674  // Check the ordering of the fluxJumps and reverse
1675  // it in case of backward definition of edge 'e'
1676  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1677  StdRegions::eBackwards)
1678  {
1679  Vmath::Reverse(nEdgePts,
1680  &fluxJumps[0], 1,
1681  &fluxJumps[0], 1);
1682  }
1683 
1684  // Deformed elements
1685  if (fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1686  ->GetGeom2D()->GetMetricInfo()->GetGtype()
1687  == SpatialDomains::eDeformed)
1688  {
1689  // Extract the Jacobians along edge 'e'
1690  Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1691 
1692  fields[0]->GetExp(n)->GetTracePhysVals(
1693  e, elmtToTrace[n][e],
1694  jac, auxArray1 = jacEdge);
1695 
1696  // Check the ordering of the fluxJumps and reverse
1697  // it in case of backward definition of edge 'e'
1698  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1699  StdRegions::eBackwards)
1700  {
1701  Vmath::Reverse(nEdgePts,
1702  &jacEdge[0], 1,
1703  &jacEdge[0], 1);
1704  }
1705 
1706  // Multiply the fluxJumps by the edge 'e' Jacobians
1707  // to bring the fluxJumps into the standard space
1708  for (j = 0; j < nEdgePts; j++)
1709  {
1710  fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1711  }
1712  }
1713  // Non-deformed elements
1714  else
1715  {
1716  // Multiply the fluxJumps by the edge 'e' Jacobians
1717  // to bring the fluxJumps into the standard space
1718  Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1,
1719  fluxJumps, 1);
1720  }
1721 
1722  // Multiply jumps by derivatives of the correction functions
1723  // All the quntities at this level should be defined into
1724  // the standard space
1725  switch (e)
1726  {
1727  case 0:
1728  for (i = 0; i < nquad0; ++i)
1729  {
1730  // Multiply fluxJumps by Q factors
1731  //fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1732 
1733  for (j = 0; j < nquad1; ++j)
1734  {
1735  cnt = i + j*nquad0;
1736  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1737  }
1738  }
1739  break;
1740  case 1:
1741  for (i = 0; i < nquad1; ++i)
1742  {
1743  // Multiply fluxJumps by Q factors
1744  //fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1745 
1746  for (j = 0; j < nquad0; ++j)
1747  {
1748  cnt = (nquad0)*i + j;
1749  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1750  }
1751  }
1752  break;
1753  case 2:
1754  for (i = 0; i < nquad0; ++i)
1755  {
1756  // Multiply fluxJumps by Q factors
1757  //fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1758 
1759  for (j = 0; j < nquad1; ++j)
1760  {
1761  cnt = j*nquad0 + i;
1762  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1763  }
1764  }
1765  break;
1766  case 3:
1767  for (i = 0; i < nquad1; ++i)
1768  {
1769  // Multiply fluxJumps by Q factors
1770  //fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1771  for (j = 0; j < nquad0; ++j)
1772  {
1773  cnt = j + i*nquad0;
1774  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1775  }
1776  }
1777  break;
1778 
1779  default:
1780  ASSERTL0(false, "edge value (< 3) is out of range");
1781  break;
1782  }
1783  }
1784 
1785  // Summing the component of the flux for calculating the
1786  // derivatives per each direction
1787  if (direction == 0)
1788  {
1789  Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1,
1790  &divCFluxE3[0], 1, &derCFlux[phys_offset], 1);
1791  }
1792  else if (direction == 1)
1793  {
1794  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1795  &divCFluxE2[0], 1, &derCFlux[phys_offset], 1);
1796  }
1797  }
1798  }
1799 
1800  /**
1801  * @brief Compute the divergence of the corrective flux for 2D problems.
1802  *
1803  * @param nConvectiveFields Number of fields.
1804  * @param fields Pointer to fields.
1805  * @param fluxX1 X1 - volumetric flux in physical space.
1806  * @param fluxX2 X2 - volumetric flux in physical space.
1807  * @param numericalFlux Interface flux in physical space.
1808  * @param divCFlux Divergence of the corrective flux for
1809  * 2D problems.
1810  *
1811  * \todo: Switch on shapes eventually here.
1812  */
1813  void DiffusionLFR::v_DivCFlux_2D(
1814  const int nConvectiveFields,
1815  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1816  const Array<OneD, const NekDouble> &fluxX1,
1817  const Array<OneD, const NekDouble> &fluxX2,
1818  const Array<OneD, const NekDouble> &numericalFlux,
1819  Array<OneD, NekDouble> &divCFlux)
1820  {
1821  boost::ignore_unused(nConvectiveFields);
1822 
1823  int n, e, i, j, cnt;
1824 
1825  int nElements = fields[0]->GetExpSize();
1826 
1827  int nLocalSolutionPts;
1828  int nEdgePts;
1829  int trace_offset;
1830  int phys_offset;
1831  int nquad0;
1832  int nquad1;
1833 
1834  Array<OneD, NekDouble> auxArray1, auxArray2;
1835  Array<OneD, LibUtilities::BasisSharedPtr> base;
1836 
1837  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1838  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1839 
1840  // Loop on the elements
1841  for(n = 0; n < nElements; ++n)
1842  {
1843  // Offset of the element on the global vector
1844  phys_offset = fields[0]->GetPhys_Offset(n);
1845  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1846 
1847  base = fields[0]->GetExp(n)->GetBase();
1848  nquad0 = base[0]->GetNumPoints();
1849  nquad1 = base[1]->GetNumPoints();
1850 
1851  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1852  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1853  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1854  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1855 
1856  // Loop on the edges
1857  for(e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1858  {
1859  // Number of edge points of edge e
1860  nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1861 
1862  Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1863  Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1864  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
1865  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
1866  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1867 
1868  // Offset of the trace space correspondent to edge e
1869  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1870  elmtToTrace[n][e]->GetElmtId());
1871 
1872  // Get the normals of edge e
1873  const Array<OneD, const Array<OneD, NekDouble> > &normals =
1874  fields[0]->GetExp(n)->GetTraceNormal(e);
1875 
1876  // Extract the edge values of flux-x on edge e and order
1877  // them accordingly to the order of the trace space
1878  fields[0]->GetExp(n)->GetTracePhysVals(
1879  e, elmtToTrace[n][e],
1880  fluxX1 + phys_offset,
1881  auxArray1 = tmparrayX1);
1882 
1883  // Extract the edge values of flux-y on edge e and order
1884  // them accordingly to the order of the trace space
1885  fields[0]->GetExp(n)->GetTracePhysVals(
1886  e, elmtToTrace[n][e],
1887  fluxX2 + phys_offset,
1888  auxArray1 = tmparrayX2);
1889 
1890  // Multiply the edge components of the flux by the normal
1891  for (i = 0; i < nEdgePts; ++i)
1892  {
1893  fluxN[i] =
1894  tmparrayX1[i]*m_traceNormals[0][trace_offset+i] +
1895  tmparrayX2[i]*m_traceNormals[1][trace_offset+i];
1896  }
1897 
1898  // Subtract to the Riemann flux the discontinuous flux
1899  Vmath::Vsub(nEdgePts,
1900  &numericalFlux[trace_offset], 1,
1901  &fluxN[0], 1, &fluxJumps[0], 1);
1902 
1903  // Check the ordering of the jump vectors
1904  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1905  StdRegions::eBackwards)
1906  {
1907  Vmath::Reverse(nEdgePts,
1908  auxArray1 = fluxJumps, 1,
1909  auxArray2 = fluxJumps, 1);
1910  }
1911 
1912  for (i = 0; i < nEdgePts; ++i)
1913  {
1914  if (m_traceNormals[0][trace_offset+i] != normals[0][i]
1915  || m_traceNormals[1][trace_offset+i] != normals[1][i])
1916  {
1917  fluxJumps[i] = -fluxJumps[i];
1918  }
1919  }
1920 
1921  // Multiply jumps by derivatives of the correction functions
1922  switch (e)
1923  {
1924  case 0:
1925  for (i = 0; i < nquad0; ++i)
1926  {
1927  // Multiply fluxJumps by Q factors
1928  fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1929 
1930  for (j = 0; j < nquad1; ++j)
1931  {
1932  cnt = i + j*nquad0;
1933  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1934  }
1935  }
1936  break;
1937  case 1:
1938  for (i = 0; i < nquad1; ++i)
1939  {
1940  // Multiply fluxJumps by Q factors
1941  fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1942 
1943  for (j = 0; j < nquad0; ++j)
1944  {
1945  cnt = (nquad0)*i + j;
1946  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1947  }
1948  }
1949  break;
1950  case 2:
1951  for (i = 0; i < nquad0; ++i)
1952  {
1953  // Multiply fluxJumps by Q factors
1954  fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1955 
1956  for (j = 0; j < nquad1; ++j)
1957  {
1958  cnt = j*nquad0 + i;
1959  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1960  }
1961  }
1962  break;
1963  case 3:
1964  for (i = 0; i < nquad1; ++i)
1965  {
1966  // Multiply fluxJumps by Q factors
1967  fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1968  for (j = 0; j < nquad0; ++j)
1969  {
1970  cnt = j + i*nquad0;
1971  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1972  }
1973  }
1974  break;
1975 
1976  default:
1977  ASSERTL0(false,"edge value (< 3) is out of range");
1978  break;
1979  }
1980  }
1981 
1982  // Sum each edge contribution
1983  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1984  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1985 
1986  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1987  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1988 
1989  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1990  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1991  }
1992  }
1993 
1994  /**
1995  * @brief Compute the divergence of the corrective flux for 2D problems
1996  * where POINTSTYPE="GaussGaussLegendre"
1997  *
1998  * @param nConvectiveFields Number of fields.
1999  * @param fields Pointer to fields.
2000  * @param fluxX1 X1-volumetric flux in physical space.
2001  * @param fluxX2 X2-volumetric flux in physical space.
2002  * @param numericalFlux Interface flux in physical space.
2003  * @param divCFlux Divergence of the corrective flux.
2004  *
2005  * \todo: Switch on shapes eventually here.
2006  */
2007 
2008  void DiffusionLFR::v_DivCFlux_2D_Gauss(
2009  const int nConvectiveFields,
2010  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
2011  const Array<OneD, const NekDouble> &fluxX1,
2012  const Array<OneD, const NekDouble> &fluxX2,
2013  const Array<OneD, const NekDouble> &numericalFlux,
2014  Array<OneD, NekDouble> &divCFlux)
2015  {
2016  boost::ignore_unused(nConvectiveFields);
2017 
2018  int n, e, i, j, cnt;
2019 
2020  int nElements = fields[0]->GetExpSize();
2021  int nLocalSolutionPts;
2022  int nEdgePts;
2023  int trace_offset;
2024  int phys_offset;
2025  int nquad0;
2026  int nquad1;
2027 
2028  Array<OneD, NekDouble> auxArray1, auxArray2;
2029  Array<OneD, LibUtilities::BasisSharedPtr> base;
2030 
2031  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
2032  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
2033 
2034  // Loop on the elements
2035  for(n = 0; n < nElements; ++n)
2036  {
2037  // Offset of the element on the global vector
2038  phys_offset = fields[0]->GetPhys_Offset(n);
2039  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
2040 
2041  base = fields[0]->GetExp(n)->GetBase();
2042  nquad0 = base[0]->GetNumPoints();
2043  nquad1 = base[1]->GetNumPoints();
2044 
2045  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
2046  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
2047  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
2048  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
2049 
2050  // Loop on the edges
2051  for(e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
2052  {
2053  // Number of edge points of edge e
2054  nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
2055 
2056  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
2057  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
2058  Array<OneD, NekDouble> fluxN_R (nEdgePts, 0.0);
2059  Array<OneD, NekDouble> fluxN_D (nEdgePts, 0.0);
2060  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
2061 
2062  // Offset of the trace space correspondent to edge e
2063  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
2064  elmtToTrace[n][e]->GetElmtId());
2065 
2066  // Get the normals of edge e
2067  const Array<OneD, const Array<OneD, NekDouble> > &normals =
2068  fields[0]->GetExp(n)->GetTraceNormal(e);
2069 
2070  // Extract the trasformed normal flux at each edge
2071  switch (e)
2072  {
2073  case 0:
2074  // Extract the edge values of transformed flux-y on
2075  // edge e and order them accordingly to the order of
2076  // the trace space
2077  fields[0]->GetExp(n)->GetTracePhysVals(
2078  e, elmtToTrace[n][e],
2079  fluxX2 + phys_offset,
2080  auxArray1 = fluxN_D);
2081 
2082  Vmath::Neg (nEdgePts, fluxN_D, 1);
2083 
2084  // Extract the physical Rieamann flux at each edge
2085  Vmath::Vcopy(nEdgePts,
2086  &numericalFlux[trace_offset], 1,
2087  &fluxN[0], 1);
2088 
2089  // Check the ordering of vectors
2090  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2091  StdRegions::eBackwards)
2092  {
2093  Vmath::Reverse(nEdgePts,
2094  auxArray1 = fluxN, 1,
2095  auxArray2 = fluxN, 1);
2096 
2097  Vmath::Reverse(nEdgePts,
2098  auxArray1 = fluxN_D, 1,
2099  auxArray2 = fluxN_D, 1);
2100  }
2101 
2102  // Transform Riemann Fluxes in the standard element
2103  for (i = 0; i < nquad0; ++i)
2104  {
2105  // Multiply Riemann Flux by Q factors
2106  fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
2107  }
2108 
2109  for (i = 0; i < nEdgePts; ++i)
2110  {
2111  if (m_traceNormals[0][trace_offset+i] !=
2112  normals[0][i] ||
2113  m_traceNormals[1][trace_offset+i] !=
2114  normals[1][i])
2115  {
2116  fluxN_R[i] = -fluxN_R[i];
2117  }
2118  }
2119 
2120  // Subtract to the Riemann flux the discontinuous
2121  // flux
2122  Vmath::Vsub(nEdgePts,
2123  &fluxN_R[0], 1,
2124  &fluxN_D[0], 1, &fluxJumps[0], 1);
2125 
2126  // Multiplicate the flux jumps for the correction
2127  // function
2128  for (i = 0; i < nquad0; ++i)
2129  {
2130  for (j = 0; j < nquad1; ++j)
2131  {
2132  cnt = i + j*nquad0;
2133  divCFluxE0[cnt] =
2134  -fluxJumps[i] * m_dGL_xi2[n][j];
2135  }
2136  }
2137  break;
2138  case 1:
2139  // Extract the edge values of transformed flux-x on
2140  // edge e and order them accordingly to the order of
2141  // the trace space
2142  fields[0]->GetExp(n)->GetTracePhysVals(
2143  e, elmtToTrace[n][e],
2144  fluxX1 + phys_offset,
2145  auxArray1 = fluxN_D);
2146 
2147  // Extract the physical Rieamann flux at each edge
2148  Vmath::Vcopy(nEdgePts,
2149  &numericalFlux[trace_offset], 1,
2150  &fluxN[0], 1);
2151 
2152  // Check the ordering of vectors
2153  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2154  StdRegions::eBackwards)
2155  {
2156  Vmath::Reverse(nEdgePts,
2157  auxArray1 = fluxN, 1,
2158  auxArray2 = fluxN, 1);
2159 
2160  Vmath::Reverse(nEdgePts,
2161  auxArray1 = fluxN_D, 1,
2162  auxArray2 = fluxN_D, 1);
2163  }
2164 
2165  // Transform Riemann Fluxes in the standard element
2166  for (i = 0; i < nquad1; ++i)
2167  {
2168  // Multiply Riemann Flux by Q factors
2169  fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
2170  }
2171 
2172  for (i = 0; i < nEdgePts; ++i)
2173  {
2174  if (m_traceNormals[0][trace_offset+i] !=
2175  normals[0][i] ||
2176  m_traceNormals[1][trace_offset+i] !=
2177  normals[1][i])
2178  {
2179  fluxN_R[i] = -fluxN_R[i];
2180  }
2181  }
2182 
2183  // Subtract to the Riemann flux the discontinuous
2184  // flux
2185  Vmath::Vsub(nEdgePts,
2186  &fluxN_R[0], 1,
2187  &fluxN_D[0], 1, &fluxJumps[0], 1);
2188 
2189  // Multiplicate the flux jumps for the correction
2190  // function
2191  for (i = 0; i < nquad1; ++i)
2192  {
2193  for (j = 0; j < nquad0; ++j)
2194  {
2195  cnt = (nquad0)*i + j;
2196  divCFluxE1[cnt] =
2197  fluxJumps[i] * m_dGR_xi1[n][j];
2198  }
2199  }
2200  break;
2201  case 2:
2202 
2203  // Extract the edge values of transformed flux-y on
2204  // edge e and order them accordingly to the order of
2205  // the trace space
2206 
2207  fields[0]->GetExp(n)->GetTracePhysVals(
2208  e, elmtToTrace[n][e],
2209  fluxX2 + phys_offset,
2210  auxArray1 = fluxN_D);
2211 
2212  // Extract the physical Rieamann flux at each edge
2213  Vmath::Vcopy(nEdgePts,
2214  &numericalFlux[trace_offset], 1,
2215  &fluxN[0], 1);
2216 
2217  // Check the ordering of vectors
2218  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2219  StdRegions::eBackwards)
2220  {
2221  Vmath::Reverse(nEdgePts,
2222  auxArray1 = fluxN, 1,
2223  auxArray2 = fluxN, 1);
2224 
2225  Vmath::Reverse(nEdgePts,
2226  auxArray1 = fluxN_D, 1,
2227  auxArray2 = fluxN_D, 1);
2228  }
2229 
2230  // Transform Riemann Fluxes in the standard element
2231  for (i = 0; i < nquad0; ++i)
2232  {
2233  // Multiply Riemann Flux by Q factors
2234  fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
2235  }
2236 
2237  for (i = 0; i < nEdgePts; ++i)
2238  {
2239  if (m_traceNormals[0][trace_offset+i] !=
2240  normals[0][i] ||
2241  m_traceNormals[1][trace_offset+i] !=
2242  normals[1][i])
2243  {
2244  fluxN_R[i] = -fluxN_R[i];
2245  }
2246  }
2247 
2248  // Subtract to the Riemann flux the discontinuous
2249  // flux
2250 
2251  Vmath::Vsub(nEdgePts,
2252  &fluxN_R[0], 1,
2253  &fluxN_D[0], 1, &fluxJumps[0], 1);
2254 
2255  // Multiplicate the flux jumps for the correction
2256  // function
2257  for (i = 0; i < nquad0; ++i)
2258  {
2259  for (j = 0; j < nquad1; ++j)
2260  {
2261  cnt = j*nquad0 + i;
2262  divCFluxE2[cnt] =
2263  fluxJumps[i] * m_dGR_xi2[n][j];
2264  }
2265  }
2266  break;
2267  case 3:
2268  // Extract the edge values of transformed flux-x on
2269  // edge e and order them accordingly to the order of
2270  // the trace space
2271 
2272  fields[0]->GetExp(n)->GetTracePhysVals(
2273  e, elmtToTrace[n][e],
2274  fluxX1 + phys_offset,
2275  auxArray1 = fluxN_D);
2276  Vmath::Neg (nEdgePts, fluxN_D, 1);
2277 
2278  // Extract the physical Rieamann flux at each edge
2279  Vmath::Vcopy(nEdgePts,
2280  &numericalFlux[trace_offset], 1,
2281  &fluxN[0], 1);
2282 
2283  // Check the ordering of vectors
2284  if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2285  StdRegions::eBackwards)
2286  {
2287  Vmath::Reverse(nEdgePts,
2288  auxArray1 = fluxN, 1,
2289  auxArray2 = fluxN, 1);
2290 
2291  Vmath::Reverse(nEdgePts,
2292  auxArray1 = fluxN_D, 1,
2293  auxArray2 = fluxN_D, 1);
2294  }
2295 
2296  // Transform Riemann Fluxes in the standard element
2297  for (i = 0; i < nquad1; ++i)
2298  {
2299  // Multiply Riemann Flux by Q factors
2300  fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
2301  }
2302 
2303  for (i = 0; i < nEdgePts; ++i)
2304  {
2305  if (m_traceNormals[0][trace_offset+i] !=
2306  normals[0][i] ||
2307  m_traceNormals[1][trace_offset+i] !=
2308  normals[1][i])
2309  {
2310  fluxN_R[i] = -fluxN_R[i];
2311  }
2312  }
2313 
2314  // Subtract to the Riemann flux the discontinuous
2315  // flux
2316 
2317  Vmath::Vsub(nEdgePts,
2318  &fluxN_R[0], 1,
2319  &fluxN_D[0], 1, &fluxJumps[0], 1);
2320 
2321  // Multiplicate the flux jumps for the correction
2322  // function
2323  for (i = 0; i < nquad1; ++i)
2324  {
2325  for (j = 0; j < nquad0; ++j)
2326  {
2327  cnt = j + i*nquad0;
2328  divCFluxE3[cnt] =
2329  -fluxJumps[i] * m_dGL_xi1[n][j];
2330  }
2331  }
2332  break;
2333  default:
2334  ASSERTL0(false,"edge value (< 3) is out of range");
2335  break;
2336  }
2337  }
2338 
2339 
2340  // Sum each edge contribution
2341  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
2342  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
2343 
2344  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2345  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2346 
2347  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2348  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2349  }
2350  }
2351 
2352  }// close namespace SolverUtils
2353 }// close namespace nektar++
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::LibUtilities::Basis
Represents a basis of a given type.
Definition: Basis.h:212
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:200
Nektar::LocalRegions::Expansion::GetMetricInfo
const SpatialDomains::GeomFactorsSharedPtr & GetMetricInfo() const
Definition: Expansion.cpp:251
Nektar::SolverUtils::DiffusionLFR::v_NumFluxforVector
virtual void v_NumFluxforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
Build the numerical flux for the 2nd order derivatives.
Definition: DiffusionLFR.cpp:1307
Nektar::SolverUtils::DiffusionLFR::m_tmp2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Definition: DiffusionLFR.h:89
Nektar::SolverUtils::DiffusionLFR::v_DerCFlux_1D
virtual void v_DerCFlux_1D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, const NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux for 1D problems.
Definition: DiffusionLFR.cpp:1474
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFR.h:64
Nektar::SolverUtils::DiffusionLFR::v_Diffuse
virtual void v_Diffuse(const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd)
Calculate FR Diffusion for the linear problems using an LDG interface flux.
Definition: DiffusionLFR.cpp:853
Nektar::SolverUtils::DiffusionLFR::m_IF1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
Definition: DiffusionLFR.h:77
Nektar::SolverUtils::DiffusionLFR::m_jac
Array< OneD, NekDouble > m_jac
Definition: DiffusionLFR.h:54
Nektar::SolverUtils::DiffusionLFR::v_NumFluxforScalar
virtual void v_NumFluxforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
Builds the numerical flux for the 1st order derivatives.
Definition: DiffusionLFR.cpp:1166
Nektar::SolverUtils::DiffusionLFR::v_WeakPenaltyforVector
virtual void v_WeakPenaltyforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux, NekDouble C11)
Imposes appropriate bcs for the 2nd order derivatives.
Definition: DiffusionLFR.cpp:1405
Nektar::SolverUtils::DiffusionLFR::m_DFC1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Definition: DiffusionLFR.h:79
Nektar::SolverUtils::DiffusionLFR::m_DD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
Definition: DiffusionLFR.h:82
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:63
Nektar::SolverUtils::DiffusionLFR::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: DiffusionLFR.cpp:174
Nektar::SolverUtils::DiffusionLFR::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: DiffusionLFR.cpp:2008
Nektar::SolverUtils::DiffusionLFR::m_divFC
Array< OneD, Array< OneD, NekDouble > > m_divFC
Definition: DiffusionLFR.h:86
Nektar::SolverUtils::DiffusionLFR::v_InitObject
virtual void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initiliase DiffusionLFR objects and store them before starting the time-stepping.
Definition: DiffusionLFR.cpp:85
Nektar::SolverUtils::DiffusionLFR::m_session
LibUtilities::SessionReaderSharedPtr m_session
Definition: DiffusionLFR.h:75
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFR.h:65
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:62
Nektar::SolverUtils::DiffusionLFR::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: DiffusionLFR.cpp:1813
Nektar::SolverUtils::DiffusionLFR::v_DerCFlux_2D
virtual void v_DerCFlux_2D(const int nConvectiveFields, const int direction, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux wrt a given coordinate for 2D problems.
Definition: DiffusionLFR.cpp:1596
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition: DiffusionLFR.h:67
Nektar::SolverUtils::DiffusionLFR::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: DiffusionLFR.h:55
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:74
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFR.h:57
Nektar::SolverUtils::DiffusionLFR::m_DU1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Definition: DiffusionLFR.h:78
Nektar::SolverUtils::DiffusionLFR::m_DFC2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Definition: DiffusionLFR.h:84
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFR.h:60
Nektar::SolverUtils::DiffusionLFR::m_tmp1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Definition: DiffusionLFR.h:88
Nektar::SolverUtils::DiffusionLFR::v_SetupCFunctions
virtual void v_SetupCFunctions(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–72...
Definition: DiffusionLFR.cpp:326
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFR.h:59
Nektar::SolverUtils::DiffusionLFR::m_D1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
Definition: DiffusionLFR.h:81
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition: DiffusionLFR.h:66
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFR.h:58
Nektar::SolverUtils::DiffusionLFR::m_IF2
Array< OneD, Array< OneD, NekDouble > > m_IF2
Definition: DiffusionLFR.h:83
Nektar::SolverUtils::DiffusionLFR::m_divFD
Array< OneD, Array< OneD, NekDouble > > m_divFD
Definition: DiffusionLFR.h:85
Nektar::SolverUtils::DiffusionLFR::v_WeakPenaltyforScalar
virtual void v_WeakPenaltyforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const Array< OneD, const NekDouble > &ufield, Array< OneD, NekDouble > &penaltyflux)
Imposes appropriate bcs for the 1st order derivatives.
Definition: DiffusionLFR.cpp:1241
Nektar::SolverUtils::DiffusionLFR::m_BD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
Definition: DiffusionLFR.h:80
Nektar::LibUtilities::BasisSharedPtr
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:74
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:120
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
Nektar::LibUtilities::eGaussGaussLegendre
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:48
Nektar::SolverUtils::GetDiffusionFactory
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:41
Nektar::SpatialDomains::eDeformed
@ eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:56
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Polylib::jacobd
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition: Polylib.cpp:1134
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:192
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:565
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:461
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:322
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:257
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1226
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:372
Generated on Wed Nov 24 2021 10:33:26 for Nektar++ by doxygen 1.9.1