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