Nektar++
DiffusionLDGNS.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: DiffusionLDGNS.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: LDGNS diffusion class.
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include "DiffusionLDGNS.h"
36 #include <iomanip>
37 #include <iostream>
38 
39 #include <boost/algorithm/string/predicate.hpp>
40 #include <boost/core/ignore_unused.hpp>
41 
43 
44 namespace Nektar
45 {
46 std::string DiffusionLDGNS::type =
48  "LDGNS", DiffusionLDGNS::create, "LDG for Navier-Stokes");
49 
51 {
52 }
53 
57 {
58  m_session = pSession;
59  m_session->LoadParameter("Twall", m_Twall, 300.15);
60 
61  // Setting up the normals
62  std::size_t nDim = pFields[0]->GetCoordim(0);
63  std::size_t nTracePts = pFields[0]->GetTrace()->GetTotPoints();
64 
65  m_spaceDim = nDim;
66  if (pSession->DefinesSolverInfo("HOMOGENEOUS"))
67  {
68  m_spaceDim = 3;
69  }
70 
71  m_diffDim = m_spaceDim - nDim;
72 
75  for (std::size_t i = 0; i < m_spaceDim; ++i)
76  {
77  m_traceVel[i] = Array<OneD, NekDouble>{nTracePts, 0.0};
78  m_traceNormals[i] = Array<OneD, NekDouble>{nTracePts};
79  }
80  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
81 
82  // Create equation of state object
83  std::string eosType;
84  m_session->LoadSolverInfo("EquationOfState", eosType, "IdealGas");
86 
87  // Set up {h} reference on the trace for penalty term
88  //
89  // Note, this shold be replaced with something smarter when merging
90  // LDG with IP
91 
92  // Get min h per element
93  std::size_t nElements = pFields[0]->GetExpSize();
94  Array<OneD, NekDouble> hEle{nElements, 1.0};
95  for (std::size_t e = 0; e < nElements; e++)
96  {
97  NekDouble h{1.0e+10};
98  std::size_t expDim = pFields[0]->GetShapeDimension();
99  switch (expDim)
100  {
101  case 3:
102  {
104  exp3D = pFields[0]->GetExp(e)->as<LocalRegions::Expansion3D>();
105  for (std::size_t i = 0; i < exp3D->GetNtraces(); ++i)
106  {
107  h = std::min(
108  h, exp3D->GetGeom3D()->GetEdge(i)->GetVertex(0)->dist(*(
109  exp3D->GetGeom3D()->GetEdge(i)->GetVertex(1))));
110  }
111  break;
112  }
113 
114  case 2:
115  {
117  exp2D = pFields[0]->GetExp(e)->as<LocalRegions::Expansion2D>();
118  for (std::size_t i = 0; i < exp2D->GetNtraces(); ++i)
119  {
120  h = std::min(
121  h, exp2D->GetGeom2D()->GetEdge(i)->GetVertex(0)->dist(*(
122  exp2D->GetGeom2D()->GetEdge(i)->GetVertex(1))));
123  }
124  break;
125  }
126  case 1:
127  {
129  exp1D = pFields[0]->GetExp(e)->as<LocalRegions::Expansion1D>();
130 
131  h = std::min(h, exp1D->GetGeom1D()->GetVertex(0)->dist(
132  *(exp1D->GetGeom1D()->GetVertex(1))));
133  break;
134  }
135  default:
136  {
137  ASSERTL0(false, "Dimension out of bound.")
138  }
139  }
140 
141  // Store scaling
142  hEle[e] = h;
143  }
144  // Expand h from elements to points
145  std::size_t nPts = pFields[0]->GetTotPoints();
146  Array<OneD, NekDouble> hElePts{nPts, 0.0};
148  for (std::size_t e = 0; e < pFields[0]->GetExpSize(); e++)
149  {
150  std::size_t nElmtPoints = pFields[0]->GetExp(e)->GetTotPoints();
151  std::size_t physOffset = pFields[0]->GetPhys_Offset(e);
152  Vmath::Fill(nElmtPoints, hEle[e], tmp = hElePts + physOffset, 1);
153  }
154  // Get Fwd and Bwd traces
155  Array<OneD, NekDouble> Fwd{nTracePts, 0.0};
156  Array<OneD, NekDouble> Bwd{nTracePts, 0.0};
157  pFields[0]->GetFwdBwdTracePhys(hElePts, Fwd, Bwd);
158  // Fix boundaries
159  std::size_t cnt = 0;
160  std::size_t nBndRegions = pFields[0]->GetBndCondExpansions().size();
161  // Loop on the boundary regions
162  for (std::size_t i = 0; i < nBndRegions; ++i)
163  {
164  // Number of boundary regions related to region 'i'
165  std::size_t nBndEdges =
166  pFields[0]->GetBndCondExpansions()[i]->GetExpSize();
167 
168  if (pFields[0]->GetBndConditions()[i]->GetBoundaryConditionType() ==
170  {
171  continue;
172  }
173 
174  // Get value from interior
175  for (std::size_t e = 0; e < nBndEdges; ++e)
176  {
177  std::size_t nBndEdgePts = pFields[0]
178  ->GetBndCondExpansions()[i]
179  ->GetExp(e)
180  ->GetTotPoints();
181 
182  std::size_t id2 = pFields[0]->GetTrace()->GetPhys_Offset(
183  pFields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
184 
185  Vmath::Vcopy(nBndEdgePts, &Fwd[id2], 1, &Bwd[id2], 1);
186  }
187  }
188  // Get average of traces
189  Array<OneD, NekDouble> traceH{nTracePts, 1.0};
190  m_traceOneOverH = Array<OneD, NekDouble>{nTracePts, 1.0};
191  Vmath::Svtsvtp(nTracePts, 0.5, Fwd, 1, 0.5, Bwd, 1, traceH, 1);
192  // Multiply by coefficient = - C11 / h
193  m_session->LoadParameter("LDGNSc11", m_C11, 1.0);
194  Vmath::Sdiv(nTracePts, -m_C11, traceH, 1, m_traceOneOverH, 1);
195 }
196 
197 /**
198  * @brief Calculate weak DG Diffusion in the LDG form for the
199  * Navier-Stokes (NS) equations:
200  *
201  * \f$ \langle\psi, \hat{u}\cdot n\rangle
202  * - \langle\nabla\psi \cdot u\rangle
203  * \langle\phi, \hat{q}\cdot n\rangle -
204  * (\nabla \phi \cdot q) \rangle \f$
205  *
206  * The equations that need a diffusion operator are those related
207  * with the velocities and with the energy.
208  *
209  */
211  const std::size_t nConvectiveFields,
213  const Array<OneD, Array<OneD, NekDouble>> &inarray,
214  Array<OneD, Array<OneD, NekDouble>> &outarray,
215  const Array<OneD, Array<OneD, NekDouble>> &pFwd,
216  const Array<OneD, Array<OneD, NekDouble>> &pBwd)
217 {
218  std::size_t nCoeffs = fields[0]->GetNcoeffs();
219  Array<OneD, Array<OneD, NekDouble>> tmp2{nConvectiveFields};
220  for (std::size_t i = 0; i < nConvectiveFields; ++i)
221  {
222  tmp2[i] = Array<OneD, NekDouble>{nCoeffs, 0.0};
223  }
224  v_DiffuseCoeffs(nConvectiveFields, fields, inarray, tmp2, pFwd, pBwd);
225  for (std::size_t i = 0; i < nConvectiveFields; ++i)
226  {
227  fields[i]->BwdTrans(tmp2[i], outarray[i]);
228  }
229 }
230 
232  const std::size_t nConvectiveFields,
234  const Array<OneD, Array<OneD, NekDouble>> &inarray,
235  Array<OneD, Array<OneD, NekDouble>> &outarray,
236  const Array<OneD, Array<OneD, NekDouble>> &pFwd,
237  const Array<OneD, Array<OneD, NekDouble>> &pBwd)
238 {
239  std::size_t nDim = fields[0]->GetCoordim(0);
240  std::size_t nPts = fields[0]->GetTotPoints();
241  std::size_t nCoeffs = fields[0]->GetNcoeffs();
242  std::size_t nScalars = inarray.size();
243  std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
244 
245  TensorOfArray3D<NekDouble> derivativesO1{m_spaceDim};
246 
247  for (std::size_t j = 0; j < m_spaceDim; ++j)
248  {
249  derivativesO1[j] = Array<OneD, Array<OneD, NekDouble>>{nScalars};
250 
251  for (std::size_t i = 0; i < nScalars; ++i)
252  {
253  derivativesO1[j][i] = Array<OneD, NekDouble>{nPts, 0.0};
254  }
255  }
256 
257  DiffuseCalcDerivative(fields, inarray, derivativesO1, pFwd, pBwd);
258 
259  // Initialisation viscous tensor
261  Array<OneD, Array<OneD, NekDouble>> viscousFlux{nConvectiveFields};
262 
263  for (std::size_t j = 0; j < m_spaceDim; ++j)
264  {
266  for (std::size_t i = 0; i < nScalars + 1; ++i)
267  {
268  m_viscTensor[j][i] = Array<OneD, NekDouble>{nPts, 0.0};
269  }
270  }
271 
272  for (std::size_t i = 0; i < nConvectiveFields; ++i)
273  {
274  viscousFlux[i] = Array<OneD, NekDouble>{nTracePts, 0.0};
275  }
276 
277  DiffuseVolumeFlux(fields, inarray, derivativesO1, m_viscTensor);
278 
279  // Compute u from q_{\eta} and q_{\xi}
280  // Obtain numerical fluxes
281  DiffuseTraceFlux(fields, inarray, derivativesO1, m_viscTensor, viscousFlux,
282  pFwd, pBwd);
283 
284  Array<OneD, NekDouble> tmpOut{nCoeffs};
286 
287  for (std::size_t i = 0; i < nConvectiveFields; ++i)
288  {
289  // Temporary fix to call collection op
290  // we can reorder m_viscTensor later
291  for (std::size_t j = 0; j < nDim; ++j)
292  {
293  tmpIn[j] = m_viscTensor[j][i];
294  }
295  fields[i]->IProductWRTDerivBase(tmpIn, tmpOut);
296 
297  // Evaulate <\phi, \hat{F}\cdot n> - outarray[i]
298  Vmath::Neg(nCoeffs, tmpOut, 1);
299  fields[i]->AddTraceIntegral(viscousFlux[i], tmpOut);
300  fields[i]->SetPhysState(false);
301  fields[i]->MultiplyByElmtInvMass(tmpOut, outarray[i]);
302  }
303 }
304 
307  const Array<OneD, Array<OneD, NekDouble>> &inarray,
309  const Array<OneD, Array<OneD, NekDouble>> &pFwd,
310  const Array<OneD, Array<OneD, NekDouble>> &pBwd)
311 {
312  std::size_t nDim = fields[0]->GetCoordim(0);
313  std::size_t nCoeffs = fields[0]->GetNcoeffs();
314  std::size_t nScalars = inarray.size();
315  std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
316  std::size_t nConvectiveFields = fields.size();
317 
318  Array<OneD, NekDouble> tmp1{nCoeffs};
319  Array<OneD, Array<OneD, NekDouble>> tmp2{nConvectiveFields};
320  TensorOfArray3D<NekDouble> numericalFluxO1{m_spaceDim};
321 
322  for (std::size_t j = 0; j < m_spaceDim; ++j)
323  {
324  numericalFluxO1[j] = Array<OneD, Array<OneD, NekDouble>>{nScalars};
325 
326  for (std::size_t i = 0; i < nScalars; ++i)
327  {
328  numericalFluxO1[j][i] = Array<OneD, NekDouble>{nTracePts, 0.0};
329  }
330  }
331 
332  NumericalFluxO1(fields, inarray, numericalFluxO1, pFwd, pBwd);
333 
334  for (std::size_t j = 0; j < nDim; ++j)
335  {
336  for (std::size_t i = 0; i < nScalars; ++i)
337  {
338  fields[i]->IProductWRTDerivBase(j, inarray[i], tmp1);
339  Vmath::Neg(nCoeffs, tmp1, 1);
340  fields[i]->AddTraceIntegral(numericalFluxO1[j][i], tmp1);
341  fields[i]->SetPhysState(false);
342  fields[i]->MultiplyByElmtInvMass(tmp1, tmp1);
343  fields[i]->BwdTrans(tmp1, qfields[j][i]);
344  }
345  }
346  // For 3D Homogeneous 1D only take derivatives in 3rd direction
347  if (m_diffDim == 1)
348  {
349  for (std::size_t i = 0; i < nScalars; ++i)
350  {
351  qfields[2][i] = m_homoDerivs[i];
352  }
353  }
354 }
355 
358  const Array<OneD, Array<OneD, NekDouble>> &inarray,
360  Array<OneD, int> &nonZeroIndex)
361 {
362 
363  boost::ignore_unused(fields, nonZeroIndex);
364  m_fluxVectorNS(inarray, qfields, VolumeFlux);
365 }
366 
369  const Array<OneD, Array<OneD, NekDouble>> &inarray,
371  Array<OneD, Array<OneD, NekDouble>> &TraceFlux,
372  const Array<OneD, Array<OneD, NekDouble>> &pFwd,
373  const Array<OneD, Array<OneD, NekDouble>> &pBwd,
374  Array<OneD, int> &nonZeroIndex)
375 {
376  boost::ignore_unused(inarray, qfields, nonZeroIndex);
377  NumericalFluxO2(fields, VolumeFlux, TraceFlux, pFwd, pBwd);
378 }
379 
380 /**
381  * @brief Builds the numerical flux for the 1st order derivatives
382  *
383  */
386  const Array<OneD, Array<OneD, NekDouble>> &inarray,
387  TensorOfArray3D<NekDouble> &numericalFluxO1,
388  const Array<OneD, Array<OneD, NekDouble>> &pFwd,
389  const Array<OneD, Array<OneD, NekDouble>> &pBwd)
390 {
391  std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
392  std::size_t nScalars = inarray.size();
393 
394  // Upwind
395  Array<OneD, Array<OneD, NekDouble>> numflux{nScalars};
396  for (std::size_t i = 0; i < nScalars; ++i)
397  {
398  numflux[i] = {pFwd[i]};
399  }
400 
401  // Modify the values in case of boundary interfaces
402  if (fields[0]->GetBndCondExpansions().size())
403  {
404  ApplyBCsO1(fields, inarray, pFwd, pBwd, numflux);
405  }
406 
407  // Splitting the numerical flux into the dimensions
408  for (std::size_t j = 0; j < m_spaceDim; ++j)
409  {
410  for (std::size_t i = 0; i < nScalars; ++i)
411  {
412  Vmath::Vmul(nTracePts, m_traceNormals[j], 1, numflux[i], 1,
413  numericalFluxO1[j][i], 1);
414  }
415  }
416 }
417 
418 /**
419  * @brief Imposes appropriate bcs for the 1st order derivatives
420  *
421  */
424  const Array<OneD, Array<OneD, NekDouble>> &inarray,
425  const Array<OneD, Array<OneD, NekDouble>> &pFwd,
426  const Array<OneD, Array<OneD, NekDouble>> &pBwd,
428 {
429  boost::ignore_unused(pBwd);
430 
431  std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
432  std::size_t nScalars = inarray.size();
433 
434  Array<OneD, NekDouble> tmp1{nTracePts, 0.0};
435  Array<OneD, NekDouble> tmp2{nTracePts, 0.0};
436  Array<OneD, NekDouble> Tw{nTracePts, m_Twall};
437 
438  Array<OneD, Array<OneD, NekDouble>> scalarVariables{nScalars};
439 
440  // Extract internal values of the scalar variables for Neumann bcs
441  for (std::size_t i = 0; i < nScalars; ++i)
442  {
443  scalarVariables[i] = Array<OneD, NekDouble>{nTracePts, 0.0};
444  }
445 
446  // Compute boundary conditions for velocities
447  for (std::size_t i = 0; i < nScalars - 1; ++i)
448  {
449  // Note that cnt has to loop on nBndRegions and nBndEdges
450  // and has to be reset to zero per each equation
451  std::size_t cnt = 0;
452  std::size_t nBndRegions = fields[i + 1]->GetBndCondExpansions().size();
453  for (std::size_t j = 0; j < nBndRegions; ++j)
454  {
455  if (fields[i + 1]
456  ->GetBndConditions()[j]
457  ->GetBoundaryConditionType() == SpatialDomains::ePeriodic)
458  {
459  continue;
460  }
461 
462  std::size_t nBndEdges =
463  fields[i + 1]->GetBndCondExpansions()[j]->GetExpSize();
464  for (std::size_t e = 0; e < nBndEdges; ++e)
465  {
466  std::size_t nBndEdgePts = fields[i + 1]
467  ->GetBndCondExpansions()[j]
468  ->GetExp(e)
469  ->GetTotPoints();
470 
471  std::size_t id1 =
472  fields[i + 1]->GetBndCondExpansions()[j]->GetPhys_Offset(e);
473 
474  std::size_t id2 = fields[0]->GetTrace()->GetPhys_Offset(
475  fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
476  cnt++));
477 
478  if (boost::iequals(
479  fields[i]->GetBndConditions()[j]->GetUserDefined(),
480  "WallViscous") ||
481  boost::iequals(
482  fields[i]->GetBndConditions()[j]->GetUserDefined(),
483  "WallAdiabatic"))
484  {
485  // Reinforcing bcs for velocity in case of Wall bcs
486  Vmath::Zero(nBndEdgePts, &scalarVariables[i][id2], 1);
487  }
488  else if (boost::iequals(
489  fields[i]->GetBndConditions()[j]->GetUserDefined(),
490  "Wall") ||
491  boost::iequals(
492  fields[i]->GetBndConditions()[j]->GetUserDefined(),
493  "Symmetry"))
494  {
495  // Symmetry bc: normal velocity is zero
496  // get all velocities at once because we need u.n
497  if (i == 0)
498  {
499  // tmp1 = -(u.n)
500  Vmath::Zero(nBndEdgePts, tmp1, 1);
501  for (std::size_t k = 0; k < nScalars - 1; ++k)
502  {
503  Vmath::Vdiv(nBndEdgePts,
504  &(fields[k + 1]
505  ->GetBndCondExpansions()[j]
506  ->UpdatePhys())[id1],
507  1,
508  &(fields[0]
509  ->GetBndCondExpansions()[j]
510  ->UpdatePhys())[id1],
511  1, &scalarVariables[k][id2], 1);
512  Vmath::Vvtvp(nBndEdgePts, &m_traceNormals[k][id2],
513  1, &scalarVariables[k][id2], 1,
514  &tmp1[0], 1, &tmp1[0], 1);
515  }
516  Vmath::Smul(nBndEdgePts, -1.0, &tmp1[0], 1, &tmp1[0],
517  1);
518 
519  // u_i - (u.n)n_i
520  for (std::size_t k = 0; k < nScalars - 1; ++k)
521  {
522  Vmath::Vvtvp(nBndEdgePts, &tmp1[0], 1,
523  &m_traceNormals[k][id2], 1,
524  &scalarVariables[k][id2], 1,
525  &scalarVariables[k][id2], 1);
526  }
527  }
528  }
529  else if (fields[i]
530  ->GetBndConditions()[j]
531  ->GetBoundaryConditionType() ==
533  {
534  // Imposing velocity bcs if not Wall
535  Vmath::Vdiv(nBndEdgePts,
536  &(fields[i + 1]
537  ->GetBndCondExpansions()[j]
538  ->UpdatePhys())[id1],
539  1,
540  &(fields[0]
541  ->GetBndCondExpansions()[j]
542  ->UpdatePhys())[id1],
543  1, &scalarVariables[i][id2], 1);
544  }
545 
546  // For Dirichlet boundary condition: uflux = u_bcs
547  if (fields[i]
548  ->GetBndConditions()[j]
549  ->GetBoundaryConditionType() ==
551  {
552  Vmath::Vcopy(nBndEdgePts, &scalarVariables[i][id2], 1,
553  &fluxO1[i][id2], 1);
554  }
555 
556  // For Neumann boundary condition: uflux = u_+
557  else if ((fields[i]->GetBndConditions()[j])
558  ->GetBoundaryConditionType() ==
560  {
561  Vmath::Vcopy(nBndEdgePts, &pFwd[i][id2], 1, &fluxO1[i][id2],
562  1);
563  }
564 
565  // Building kinetic energy to be used for T bcs
566  Vmath::Vmul(nBndEdgePts, &scalarVariables[i][id2], 1,
567  &scalarVariables[i][id2], 1, &tmp1[id2], 1);
568 
569  Vmath::Smul(nBndEdgePts, 0.5, &tmp1[id2], 1, &tmp1[id2], 1);
570 
571  Vmath::Vadd(nBndEdgePts, &tmp2[id2], 1, &tmp1[id2], 1,
572  &tmp2[id2], 1);
573  }
574  }
575  }
576 
577  // Compute boundary conditions for temperature
578  std::size_t cnt = 0;
579  std::size_t nBndRegions = fields[nScalars]->GetBndCondExpansions().size();
580  for (std::size_t j = 0; j < nBndRegions; ++j)
581  {
582  if (fields[nScalars]
583  ->GetBndConditions()[j]
584  ->GetBoundaryConditionType() == SpatialDomains::ePeriodic)
585  {
586  continue;
587  }
588 
589  std::size_t nBndEdges =
590  fields[nScalars]->GetBndCondExpansions()[j]->GetExpSize();
591  for (std::size_t e = 0; e < nBndEdges; ++e)
592  {
593  std::size_t nBndEdgePts = fields[nScalars]
594  ->GetBndCondExpansions()[j]
595  ->GetExp(e)
596  ->GetTotPoints();
597 
598  std::size_t id1 =
599  fields[nScalars]->GetBndCondExpansions()[j]->GetPhys_Offset(e);
600 
601  std::size_t id2 = fields[0]->GetTrace()->GetPhys_Offset(
602  fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
603 
604  // Imposing Temperature Twall at the wall
605  if (boost::iequals(
606  fields[nScalars]->GetBndConditions()[j]->GetUserDefined(),
607  "WallViscous"))
608  {
609  Vmath::Vcopy(nBndEdgePts, &Tw[0], 1,
610  &scalarVariables[nScalars - 1][id2], 1);
611  }
612  // Imposing Temperature through condition on the Energy
613  // for no wall boundaries (e.g. farfield)
614  else if (fields[nScalars]
615  ->GetBndConditions()[j]
616  ->GetBoundaryConditionType() ==
618  {
619  // Use equation of state to evaluate temperature
620  NekDouble rho, ene;
621  for (std::size_t n = 0; n < nBndEdgePts; ++n)
622  {
623  rho = fields[0]
624  ->GetBndCondExpansions()[j]
625  ->GetPhys()[id1 + n];
626  ene = fields[nScalars]
627  ->GetBndCondExpansions()[j]
628  ->GetPhys()[id1 + n] /
629  rho -
630  tmp2[id2 + n];
631  scalarVariables[nScalars - 1][id2 + n] =
632  m_eos->GetTemperature(rho, ene);
633  }
634  }
635 
636  // For Dirichlet boundary condition: uflux = u_bcs
637  if (fields[nScalars]
638  ->GetBndConditions()[j]
639  ->GetBoundaryConditionType() ==
641  !boost::iequals(
642  fields[nScalars]->GetBndConditions()[j]->GetUserDefined(),
643  "WallAdiabatic"))
644  {
645  Vmath::Vcopy(nBndEdgePts, &scalarVariables[nScalars - 1][id2],
646  1, &fluxO1[nScalars - 1][id2], 1);
647  }
648 
649  // For Neumann boundary condition: uflux = u_+
650  else if (((fields[nScalars]->GetBndConditions()[j])
651  ->GetBoundaryConditionType() ==
653  boost::iequals(fields[nScalars]
654  ->GetBndConditions()[j]
655  ->GetUserDefined(),
656  "WallAdiabatic"))
657  {
658  Vmath::Vcopy(nBndEdgePts, &pFwd[nScalars - 1][id2], 1,
659  &fluxO1[nScalars - 1][id2], 1);
660  }
661  }
662  }
663 }
664 
665 /**
666  * @brief Build the numerical flux for the 2nd order derivatives
667  *
668  */
673  const Array<OneD, Array<OneD, NekDouble>> &uFwd,
674  const Array<OneD, Array<OneD, NekDouble>> &uBwd)
675 {
676  std::size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
677  std::size_t nVariables = fields.size();
678 
679  // Initialize penalty flux
680  Array<OneD, Array<OneD, NekDouble>> fluxPen{nVariables - 1};
681  for (std::size_t i = 0; i < nVariables - 1; ++i)
682  {
683  fluxPen[i] = Array<OneD, NekDouble>{nTracePts, 0.0};
684  }
685 
686  // Get penalty flux
687  m_fluxPenaltyNS(uFwd, uBwd, fluxPen);
688 
689  // Evaluate Riemann flux
690  // qflux = \hat{q} \cdot u = q \cdot n
691  // Notice: i = 1 (first row of the viscous tensor is zero)
692 
693  Array<OneD, NekDouble> qFwd{nTracePts};
694  Array<OneD, NekDouble> qBwd{nTracePts};
695  Array<OneD, NekDouble> qtemp{nTracePts};
696  Array<OneD, NekDouble> qfluxtemp{nTracePts, 0.0};
697  std::size_t nDim = fields[0]->GetCoordim(0);
698  for (std::size_t i = 1; i < nVariables; ++i)
699  {
700  qflux[i] = Array<OneD, NekDouble>{nTracePts, 0.0};
701  for (std::size_t j = 0; j < nDim; ++j)
702  {
703  // Compute qFwd and qBwd value of qfield in position 'ji'
704  fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
705 
706  // Downwind
707  Vmath::Vcopy(nTracePts, qBwd, 1, qfluxtemp, 1);
708 
709  // Multiply the Riemann flux by the trace normals
710  Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
711  qfluxtemp, 1);
712 
713  // Add penalty term
714  Vmath::Vvtvp(nTracePts, m_traceOneOverH, 1, fluxPen[i - 1], 1,
715  qfluxtemp, 1, qfluxtemp, 1);
716 
717  // Impose weak boundary condition with flux
718  if (fields[0]->GetBndCondExpansions().size())
719  {
720  ApplyBCsO2(fields, i, j, qfield[j][i], qFwd, qBwd, qfluxtemp);
721  }
722 
723  // Store the final flux into qflux
724  Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1, qflux[i], 1);
725  }
726  }
727 }
728 
729 /**
730  * @brief Imposes appropriate bcs for the 2nd order derivatives
731  *
732  */
735  const std::size_t var, const std::size_t dir,
736  const Array<OneD, const NekDouble> &qfield,
737  const Array<OneD, const NekDouble> &qFwd,
738  const Array<OneD, const NekDouble> &qBwd,
739  Array<OneD, NekDouble> &penaltyflux)
740 {
741  boost::ignore_unused(qfield, qBwd);
742 
743  std::size_t cnt = 0;
744  std::size_t nBndRegions = fields[var]->GetBndCondExpansions().size();
745  // Loop on the boundary regions to apply appropriate bcs
746  for (std::size_t i = 0; i < nBndRegions; ++i)
747  {
748  // Number of boundary regions related to region 'i'
749  std::size_t nBndEdges =
750  fields[var]->GetBndCondExpansions()[i]->GetExpSize();
751 
752  if (fields[var]->GetBndConditions()[i]->GetBoundaryConditionType() ==
754  {
755  continue;
756  }
757 
758  // Weakly impose bcs by modifying flux values
759  for (std::size_t e = 0; e < nBndEdges; ++e)
760  {
761  std::size_t nBndEdgePts = fields[var]
762  ->GetBndCondExpansions()[i]
763  ->GetExp(e)
764  ->GetTotPoints();
765 
766  std::size_t id2 = fields[0]->GetTrace()->GetPhys_Offset(
767  fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
768 
769  // In case of Dirichlet bcs:
770  // uflux = gD
771  if (fields[var]
772  ->GetBndConditions()[i]
773  ->GetBoundaryConditionType() ==
775  !boost::iequals(
776  fields[var]->GetBndConditions()[i]->GetUserDefined(),
777  "WallAdiabatic"))
778  {
779  Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
780  &qFwd[id2], 1, &penaltyflux[id2], 1);
781  }
782  // 3.4) In case of Neumann bcs:
783  // uflux = u+
784  else if ((fields[var]->GetBndConditions()[i])
785  ->GetBoundaryConditionType() ==
787  {
788  ASSERTL0(false, "Neumann bcs not implemented for LDGNS");
789 
790  /*
791  Vmath::Vmul(nBndEdgePts,
792  &m_traceNormals[dir][id2], 1,
793  &(fields[var]->
794  GetBndCondExpansions()[i]->
795  UpdatePhys())[id1], 1,
796  &penaltyflux[id2], 1);
797  */
798  }
799  else if (boost::iequals(
800  fields[var]->GetBndConditions()[i]->GetUserDefined(),
801  "WallAdiabatic"))
802  {
803  if ((var == m_spaceDim + 1))
804  {
805  Vmath::Zero(nBndEdgePts, &penaltyflux[id2], 1);
806  }
807  else
808  {
809 
810  Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
811  &qFwd[id2], 1, &penaltyflux[id2], 1);
812  }
813  }
814  }
815  }
816 }
817 
818 /**
819  * @brief Compute primary variables
820  *
821  */
824  const Array<OneD, Array<OneD, NekDouble>> &inarray,
825  Array<OneD, Array<OneD, NekDouble>> &primVar)
826 {
827  int nDim = fields[0]->GetCoordim(0);
828  int nPts = fields[0]->GetTotPoints();
829  for (int i = 0; i < nDim; ++i)
830  {
831  primVar[i] = Array<OneD, NekDouble>(nPts, 0.0);
832  primVar[i] = inarray[i];
833  }
834 }
835 } // end of namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
virtual void v_DiffuseCalcDerivative(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfields, const Array< OneD, Array< OneD, NekDouble >> &pFwd, const Array< OneD, Array< OneD, NekDouble >> &pBwd) override
Diffusion Flux, calculate the physical derivatives.
virtual void v_DiffuseCoeffs(const std::size_t nConvective, 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) override
Array< OneD, Array< OneD, NekDouble > > m_homoDerivs
EquationOfStateSharedPtr m_eos
Equation of system for computing temperature.
virtual void v_Diffuse(const std::size_t nConvective, 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) override
Calculate weak DG Diffusion in the LDG form for the Navier-Stokes (NS) equations:
LibUtilities::SessionReaderSharedPtr m_session
virtual void v_DiffuseTraceFlux(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfields, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, Array< OneD, NekDouble >> &TraceFlux, const Array< OneD, Array< OneD, NekDouble >> &pFwd, const Array< OneD, Array< OneD, NekDouble >> &pBwd, Array< OneD, int > &nonZeroIndex) override
Diffusion term Trace Flux.
Array< OneD, Array< OneD, NekDouble > > m_traceVel
void ApplyBCsO1(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, const Array< OneD, Array< OneD, NekDouble >> &pFwd, const Array< OneD, Array< OneD, NekDouble >> &pBwd, Array< OneD, Array< OneD, NekDouble >> &flux01)
Imposes appropriate bcs for the 1st order derivatives.
NekDouble m_C11
Penalty coefficient for LDGNS.
virtual void v_DiffuseVolumeFlux(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfields, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, int > &nonZeroIndex) override
Diffusion Volume Flux.
virtual void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
static std::string type
TensorOfArray3D< NekDouble > m_viscTensor
void NumericalFluxO1(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &numericalFluxO1, const Array< OneD, Array< OneD, NekDouble >> &pFwd, const Array< OneD, Array< OneD, NekDouble >> &pBwd)
Builds the numerical flux for the 1st order derivatives.
void ApplyBCsO2(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const std::size_t var, const std::size_t dir, const Array< OneD, const NekDouble > &qfield, const Array< OneD, const NekDouble > &qFwd, const Array< OneD, const NekDouble > &qBwd, Array< OneD, NekDouble > &penaltyflux)
Imposes appropriate bcs for the 2nd order derivatives.
void NumericalFluxO2(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, TensorOfArray3D< NekDouble > &qfield, Array< OneD, Array< OneD, NekDouble >> &qflux, const Array< OneD, Array< OneD, NekDouble >> &pFwd, const Array< OneD, Array< OneD, NekDouble >> &pBwd)
Build the numerical flux for the 2nd order derivatives.
Array< OneD, NekDouble > m_traceOneOverH
h scaling for penalty term
virtual void v_GetPrimVar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &primVar) override
Compute primary variables.
static DiffusionSharedPtr create(std::string diffType)
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
SpatialDomains::Geometry1DSharedPtr GetGeom1D() const
Definition: Expansion1D.h:108
SOLVER_UTILS_EXPORT void DiffuseVolumeFlux(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfields, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, int > &nonZeroIndex=NullInt1DArray)
Diffusion Volume FLux.
Definition: Diffusion.h:206
SOLVER_UTILS_EXPORT void DiffuseCalcDerivative(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfields, const Array< OneD, Array< OneD, NekDouble >> &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble >> &pBwd=NullNekDoubleArrayOfArray)
Definition: Diffusion.cpp:194
SOLVER_UTILS_EXPORT void DiffuseTraceFlux(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfields, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, Array< OneD, NekDouble >> &TraceFlux, const Array< OneD, Array< OneD, NekDouble >> &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble >> &pBwd=NullNekDoubleArrayOfArray, Array< OneD, int > &nonZeroIndex=NullInt1DArray)
Diffusion term Trace Flux.
Definition: Diffusion.h:217
DiffusionFluxVecCBNS m_fluxVectorNS
Definition: Diffusion.h:403
DiffusionFluxPenaltyNS m_fluxPenaltyNS
Definition: Diffusion.h:404
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:47
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:51
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:50
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:41
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
EquationOfStateFactory & GetEquationOfStateFactory()
Declaration of the equation of state factory singleton.
double NekDouble
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
svtvvtp (scalar times vector plus scalar times vector):
Definition: Vmath.cpp:751
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:209
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:574
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:359
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:248
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:324
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:284
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255