Nektar++
NavierStokesCFE.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File NavierStokesCFE.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: Navier Stokes equations in conservative variables
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
36 
37 using namespace std;
38 
39 namespace Nektar
40 {
41 string NavierStokesCFE::className =
43  "NavierStokesCFE", NavierStokesCFE::create,
44  "NavierStokes equations in conservative variables.");
45 
46 NavierStokesCFE::NavierStokesCFE(
49  : UnsteadySystem(pSession, pGraph), CompressibleFlowSystem(pSession, pGraph)
50 {
51 }
52 
54 {
55 }
56 
57 /**
58  * @brief Initialization object for CompressibleFlowSystem class.
59  */
60 void NavierStokesCFE::v_InitObject(bool DeclareFields)
61 {
63 
64  // rest of initialisation is in this routine so it can also be called
65  // in NavierStokesImplicitCFE initialisation
67 }
68 
70 {
71  // Get gas constant from session file and compute Cp
72  NekDouble gasConstant;
73  m_session->LoadParameter("GasConstant", gasConstant, 287.058);
74  m_Cp = m_gamma / (m_gamma - 1.0) * gasConstant;
75  m_Cv = m_Cp / m_gamma;
76 
77  m_session->LoadParameter("Twall", m_Twall, 300.15);
78 
79  // Viscosity
80  m_session->LoadSolverInfo("ViscosityType", m_ViscosityType, "Constant");
81  m_session->LoadParameter("mu", m_muRef, 1.78e-05);
82  if (m_ViscosityType == "Variable")
83  {
84  m_is_mu_variable = true;
85  }
86 
87  // Thermal conductivity or Prandtl
88  if (m_session->DefinesParameter("thermalConductivity"))
89  {
90  ASSERTL0(!m_session->DefinesParameter("Pr"),
91  "Cannot define both Pr and thermalConductivity.");
92 
93  m_session->LoadParameter("thermalConductivity",
96  }
97  else
98  {
99  m_session->LoadParameter("Pr", m_Prandtl, 0.72);
101  }
102 
103  if (m_shockCaptureType == "Physical")
104  {
105  m_is_shockCaptPhys = true;
106  }
107 
108  string diffName;
109  m_session->LoadSolverInfo("DiffusionType", diffName, "LDGNS");
110 
111  m_diffusion =
112  SolverUtils::GetDiffusionFactory().CreateInstance(diffName, diffName);
113  if ("InteriorPenalty" == diffName)
114  {
115  m_is_diffIP = true;
117  }
118 
119  if ("LDGNS" == diffName || "LDGNS3DHomogeneous1D" == diffName)
120  {
121  m_diffusion->SetFluxPenaltyNS(&NavierStokesCFE::v_GetFluxPenalty, this);
122  }
123 
125  {
126  m_diffusion->SetFluxVectorNS(
128  }
129  else
130  {
132  this);
133  }
134 
135  m_diffusion->SetDiffusionFluxCons(
136  &NavierStokesCFE::GetViscousFluxVectorConservVar<false>, this);
137 
138  m_diffusion->SetDiffusionFluxConsTrace(
139  &NavierStokesCFE::GetViscousFluxVectorConservVar<true>, this);
140 
141  m_diffusion->SetSpecialBndTreat(&NavierStokesCFE::SpecialBndTreat, this);
142 
143  m_diffusion->SetDiffusionSymmFluxCons(
145 
146  // Concluding initialisation of diffusion operator
147  m_diffusion->InitObject(m_session, m_fields);
148 }
149 
151  const Array<OneD, Array<OneD, NekDouble>> &inarray,
152  Array<OneD, Array<OneD, NekDouble>> &outarray,
153  const Array<OneD, Array<OneD, NekDouble>> &pFwd,
154  const Array<OneD, Array<OneD, NekDouble>> &pBwd)
155 {
156  size_t nvariables = inarray.size();
157  size_t npoints = GetNpoints();
158  size_t nTracePts = GetTraceTotPoints();
159 
160  // this should be preallocated
161  Array<OneD, Array<OneD, NekDouble>> outarrayDiff(nvariables);
162  for (int i = 0; i < nvariables; ++i)
163  {
164  outarrayDiff[i] = Array<OneD, NekDouble>(npoints, 0.0);
165  }
166 
167  // Set artificial viscosity based on NS viscous tensor
168  if (m_is_shockCaptPhys)
169  {
170  Array<OneD, NekDouble> div(npoints), curlSquare(npoints);
171  GetDivCurlSquared(m_fields, inarray, div, curlSquare, pFwd, pBwd);
172 
173  // Set volume and trace artificial viscosity
174  m_varConv->SetAv(m_fields, inarray, div, curlSquare);
175  }
176 
177  if (m_is_diffIP)
178  {
179  if (m_bndEvaluateTime < 0.0)
180  {
181  NEKERROR(ErrorUtil::efatal, "m_bndEvaluateTime not setup");
182  }
183  m_diffusion->Diffuse(nvariables, m_fields, inarray, outarrayDiff,
184  m_bndEvaluateTime, pFwd, pBwd);
185  for (int i = 0; i < nvariables; ++i)
186  {
187  Vmath::Vadd(npoints, outarrayDiff[i], 1, outarray[i], 1,
188  outarray[i], 1);
189  }
190  }
191  else
192  {
193  // Get primitive variables [u,v,w,T]
194  Array<OneD, Array<OneD, NekDouble>> inarrayDiff(nvariables - 1);
195  Array<OneD, Array<OneD, NekDouble>> inFwd(nvariables - 1);
196  Array<OneD, Array<OneD, NekDouble>> inBwd(nvariables - 1);
197 
198  for (int i = 0; i < nvariables - 1; ++i)
199  {
200  inarrayDiff[i] = Array<OneD, NekDouble>{npoints};
201  inFwd[i] = Array<OneD, NekDouble>{nTracePts};
202  inBwd[i] = Array<OneD, NekDouble>{nTracePts};
203  }
204 
205  // Extract pressure
206  // (use inarrayDiff[0] as a temporary storage for the pressure)
207  m_varConv->GetPressure(inarray, inarrayDiff[0]);
208 
209  // Extract temperature
210  m_varConv->GetTemperature(inarray, inarrayDiff[nvariables - 2]);
211 
212  // Extract velocities
213  m_varConv->GetVelocityVector(inarray, inarrayDiff);
214 
215  // Repeat calculation for trace space
216  if (pFwd == NullNekDoubleArrayOfArray ||
218  {
221  }
222  else
223  {
224  m_varConv->GetPressure(pFwd, inFwd[0]);
225  m_varConv->GetPressure(pBwd, inBwd[0]);
226 
227  m_varConv->GetTemperature(pFwd, inFwd[nvariables - 2]);
228  m_varConv->GetTemperature(pBwd, inBwd[nvariables - 2]);
229 
230  m_varConv->GetVelocityVector(pFwd, inFwd);
231  m_varConv->GetVelocityVector(pBwd, inBwd);
232  }
233 
234  // Diffusion term in physical rhs form
235  m_diffusion->Diffuse(nvariables, m_fields, inarrayDiff, outarrayDiff,
236  inFwd, inBwd);
237 
238  for (int i = 0; i < nvariables; ++i)
239  {
240  Vmath::Vadd(npoints, outarrayDiff[i], 1, outarray[i], 1,
241  outarray[i], 1);
242  }
243  }
244 }
245 
246 /**
247  * @brief Return the flux vector for the LDG diffusion problem.
248  * \todo Complete the viscous flux vector
249  */
251  const Array<OneD, const Array<OneD, NekDouble>> &physfield,
252  TensorOfArray3D<NekDouble> &derivativesO1,
253  TensorOfArray3D<NekDouble> &viscousTensor)
254 {
255  // Auxiliary variables
256  size_t nScalar = physfield.size();
257  size_t nPts = physfield[0].size();
258  Array<OneD, NekDouble> divVel(nPts, 0.0);
259 
260  // Stokes hypothesis
261  const NekDouble lambda = -2.0 / 3.0;
262 
263  // Update viscosity and thermal conductivity
264  Array<OneD, NekDouble> mu(nPts, 0.0);
265  Array<OneD, NekDouble> thermalConductivity(nPts, 0.0);
266  GetViscosityAndThermalCondFromTemp(physfield[nScalar - 1], mu,
267  thermalConductivity);
268 
269  // Velocity divergence
270  for (int j = 0; j < m_spacedim; ++j)
271  {
272  Vmath::Vadd(nPts, divVel, 1, derivativesO1[j][j], 1, divVel, 1);
273  }
274 
275  // Velocity divergence scaled by lambda * mu
276  Vmath::Smul(nPts, lambda, divVel, 1, divVel, 1);
277  Vmath::Vmul(nPts, mu, 1, divVel, 1, divVel, 1);
278 
279  // Viscous flux vector for the rho equation = 0
280  for (int i = 0; i < m_spacedim; ++i)
281  {
282  Vmath::Zero(nPts, viscousTensor[i][0], 1);
283  }
284 
285  // Viscous stress tensor (for the momentum equations)
286  for (int i = 0; i < m_spacedim; ++i)
287  {
288  for (int j = i; j < m_spacedim; ++j)
289  {
290  Vmath::Vadd(nPts, derivativesO1[i][j], 1, derivativesO1[j][i], 1,
291  viscousTensor[i][j + 1], 1);
292 
293  Vmath::Vmul(nPts, mu, 1, viscousTensor[i][j + 1], 1,
294  viscousTensor[i][j + 1], 1);
295 
296  if (i == j)
297  {
298  // Add divergence term to diagonal
299  Vmath::Vadd(nPts, viscousTensor[i][j + 1], 1, divVel, 1,
300  viscousTensor[i][j + 1], 1);
301  }
302  else
303  {
304  // Copy to make symmetric
305  Vmath::Vcopy(nPts, viscousTensor[i][j + 1], 1,
306  viscousTensor[j][i + 1], 1);
307  }
308  }
309  }
310 
311  // Terms for the energy equation
312  for (int i = 0; i < m_spacedim; ++i)
313  {
314  Vmath::Zero(nPts, viscousTensor[i][m_spacedim + 1], 1);
315  // u_j * tau_ij
316  for (int j = 0; j < m_spacedim; ++j)
317  {
318  Vmath::Vvtvp(nPts, physfield[j], 1, viscousTensor[i][j + 1], 1,
319  viscousTensor[i][m_spacedim + 1], 1,
320  viscousTensor[i][m_spacedim + 1], 1);
321  }
322  // Add k*T_i
323  Vmath::Vvtvp(nPts, thermalConductivity, 1, derivativesO1[i][m_spacedim],
324  1, viscousTensor[i][m_spacedim + 1], 1,
325  viscousTensor[i][m_spacedim + 1], 1);
326  }
327 }
328 
329 /**
330  * @brief Return the flux vector for the LDG diffusion problem.
331  * \todo Complete the viscous flux vector
332  */
334  const Array<OneD, const Array<OneD, NekDouble>> &physfield,
335  TensorOfArray3D<NekDouble> &derivativesO1,
336  TensorOfArray3D<NekDouble> &viscousTensor)
337 {
338  // Factor to rescale 1d points in dealiasing.
339  NekDouble OneDptscale = 2;
340  // Get number of points to dealias a cubic non-linearity
341  size_t nScalar = physfield.size();
342  int nPts = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
343  size_t nPts_orig = physfield[0].size();
344 
345  // Auxiliary variables
346  Array<OneD, NekDouble> divVel(nPts, 0.0);
347 
348  // Stokes hypothesis
349  const NekDouble lambda = -2.0 / 3.0;
350 
351  // Update viscosity and thermal conductivity
352  Array<OneD, NekDouble> mu(nPts, 0.0);
353  Array<OneD, NekDouble> thermalConductivity(nPts, 0.0);
354  GetViscosityAndThermalCondFromTemp(physfield[nScalar - 1], mu,
355  thermalConductivity);
356 
357  // Interpolate inputs and initialise interpolated output
361  for (int i = 0; i < m_spacedim; ++i)
362  {
363  // Interpolate velocity
364  vel_interp[i] = Array<OneD, NekDouble>(nPts);
365  m_fields[0]->PhysInterp1DScaled(OneDptscale, physfield[i],
366  vel_interp[i]);
367 
368  // Interpolate derivatives
369  deriv_interp[i] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim + 1);
370  for (int j = 0; j < m_spacedim + 1; ++j)
371  {
372  deriv_interp[i][j] = Array<OneD, NekDouble>(nPts);
373  m_fields[0]->PhysInterp1DScaled(OneDptscale, derivativesO1[i][j],
374  deriv_interp[i][j]);
375  }
376 
377  // Output (start from j=1 since flux is zero for rho)
378  out_interp[i] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim + 2);
379  for (int j = 1; j < m_spacedim + 2; ++j)
380  {
381  out_interp[i][j] = Array<OneD, NekDouble>(nPts);
382  }
383  }
384 
385  // Velocity divergence
386  for (int j = 0; j < m_spacedim; ++j)
387  {
388  Vmath::Vadd(nPts, divVel, 1, deriv_interp[j][j], 1, divVel, 1);
389  }
390 
391  // Velocity divergence scaled by lambda * mu
392  Vmath::Smul(nPts, lambda, divVel, 1, divVel, 1);
393  Vmath::Vmul(nPts, mu, 1, divVel, 1, divVel, 1);
394 
395  // Viscous flux vector for the rho equation = 0 (no need to dealias)
396  for (int i = 0; i < m_spacedim; ++i)
397  {
398  Vmath::Zero(nPts_orig, viscousTensor[i][0], 1);
399  }
400 
401  // Viscous stress tensor (for the momentum equations)
402  for (int i = 0; i < m_spacedim; ++i)
403  {
404  for (int j = i; j < m_spacedim; ++j)
405  {
406  Vmath::Vadd(nPts, deriv_interp[i][j], 1, deriv_interp[j][i], 1,
407  out_interp[i][j + 1], 1);
408 
409  Vmath::Vmul(nPts, mu, 1, out_interp[i][j + 1], 1,
410  out_interp[i][j + 1], 1);
411 
412  if (i == j)
413  {
414  // Add divergence term to diagonal
415  Vmath::Vadd(nPts, out_interp[i][j + 1], 1, divVel, 1,
416  out_interp[i][j + 1], 1);
417  }
418  else
419  {
420  // Make symmetric
421  out_interp[j][i + 1] = out_interp[i][j + 1];
422  }
423  }
424  }
425 
426  // Terms for the energy equation
427  for (int i = 0; i < m_spacedim; ++i)
428  {
429  Vmath::Zero(nPts, out_interp[i][m_spacedim + 1], 1);
430  // u_j * tau_ij
431  for (int j = 0; j < m_spacedim; ++j)
432  {
433  Vmath::Vvtvp(nPts, vel_interp[j], 1, out_interp[i][j + 1], 1,
434  out_interp[i][m_spacedim + 1], 1,
435  out_interp[i][m_spacedim + 1], 1);
436  }
437  // Add k*T_i
438  Vmath::Vvtvp(nPts, thermalConductivity, 1, deriv_interp[i][m_spacedim],
439  1, out_interp[i][m_spacedim + 1], 1,
440  out_interp[i][m_spacedim + 1], 1);
441  }
442 
443  // Project to original space
444  for (int i = 0; i < m_spacedim; ++i)
445  {
446  for (int j = 1; j < m_spacedim + 2; ++j)
447  {
448  m_fields[0]->PhysGalerkinProjection1DScaled(
449  OneDptscale, out_interp[i][j], viscousTensor[i][j]);
450  }
451  }
452 }
453 
454 /**
455  * @brief For very special treatment. For general boundaries it does nothing
456  * But for WallViscous and WallAdiabatic, special treatment is needed
457  * because they get the same Bwd value, special treatment is needed for
458  * boundary treatment of diffusion flux
459  * Note: This special treatment could be removed by seperating
460  * WallViscous and WallAdiabatic into two different classes.
461  *
462  */
464  Array<OneD, Array<OneD, NekDouble>> &consvar)
465 {
466  size_t nConvectiveFields = consvar.size();
467  int ndens = 0;
468  int nengy = nConvectiveFields - 1;
469 
470  Array<OneD, Array<OneD, NekDouble>> bndCons{nConvectiveFields};
471  Array<OneD, NekDouble> bndTotEngy;
472  Array<OneD, NekDouble> bndPressure;
473  Array<OneD, NekDouble> bndRho;
474  Array<OneD, NekDouble> bndIntEndy;
475  int nLengthArray = 0;
476 
477  int cnt = 0;
478  int nBndRegions = m_fields[nengy]->GetBndCondExpansions().size();
479  for (int j = 0; j < nBndRegions; ++j)
480  {
481  if (m_fields[nengy]
482  ->GetBndConditions()[j]
483  ->GetBoundaryConditionType() == SpatialDomains::ePeriodic)
484  {
485  continue;
486  }
487 
488  size_t nBndEdges =
489  m_fields[nengy]->GetBndCondExpansions()[j]->GetExpSize();
490  for (int e = 0; e < nBndEdges; ++e)
491  {
492  size_t nBndEdgePts = m_fields[nengy]
493  ->GetBndCondExpansions()[j]
494  ->GetExp(e)
495  ->GetTotPoints();
496 
497  int id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
498  m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
499 
500  // Imposing Temperature Twall at the wall
501  if (boost::iequals(
502  m_fields[nengy]->GetBndConditions()[j]->GetUserDefined(),
503  "WallViscous"))
504  {
505  if (nBndEdgePts != nLengthArray)
506  {
507  for (int i = 0; i < nConvectiveFields; ++i)
508  {
509  bndCons[i] = Array<OneD, NekDouble>{nBndEdgePts, 0.0};
510  }
511  bndTotEngy = Array<OneD, NekDouble>{nBndEdgePts, 0.0};
512  bndPressure = Array<OneD, NekDouble>{nBndEdgePts, 0.0};
513  bndRho = Array<OneD, NekDouble>{nBndEdgePts, 0.0};
514  bndIntEndy = Array<OneD, NekDouble>{nBndEdgePts, 0.0};
515  nLengthArray = nBndEdgePts;
516  }
517  else
518  {
519  Vmath::Zero(nLengthArray, bndPressure, 1);
520  Vmath::Zero(nLengthArray, bndRho, 1);
521  Vmath::Zero(nLengthArray, bndIntEndy, 1);
522  }
523 
525 
526  for (int k = 0; k < nConvectiveFields; ++k)
527  {
528  Vmath::Vcopy(nBndEdgePts, tmp = consvar[k] + id2, 1,
529  bndCons[k], 1);
530  }
531 
532  m_varConv->GetPressure(bndCons, bndPressure);
533  Vmath::Fill(nLengthArray, m_Twall, bndTotEngy, 1);
534  m_varConv->GetRhoFromPT(bndPressure, bndTotEngy, bndRho);
535  m_varConv->GetEFromRhoP(bndRho, bndPressure, bndIntEndy);
536  m_varConv->GetDynamicEnergy(bndCons, bndTotEngy);
537 
538  Vmath::Vvtvp(nBndEdgePts, bndIntEndy, 1, bndCons[ndens], 1,
539  bndTotEngy, 1, bndTotEngy, 1);
540 
541  Vmath::Vcopy(nBndEdgePts, bndTotEngy, 1,
542  tmp = consvar[nengy] + id2, 1);
543  }
544  }
545  }
546 }
547 
548 /**
549  *
550  * @brief Calculate and return the Symmetric flux in IP method.
551  */
553  const int nDim, const Array<OneD, Array<OneD, NekDouble>> &inaverg,
554  const Array<OneD, Array<OneD, NekDouble>> &inarray,
555  TensorOfArray3D<NekDouble> &outarray, Array<OneD, int> &nonZeroIndex,
556  const Array<OneD, Array<OneD, NekDouble>> &normal)
557 {
558  size_t nConvectiveFields = inarray.size();
559  size_t nPts = inaverg[nConvectiveFields - 1].size();
560  nonZeroIndex = Array<OneD, int>{nConvectiveFields - 1, 0};
561  for (int i = 0; i < nConvectiveFields - 1; ++i)
562  {
563  nonZeroIndex[i] = i + 1;
564  }
565 
566  Array<OneD, NekDouble> mu(nPts, 0.0);
567  Array<OneD, NekDouble> thermalConductivity(nPts, 0.0);
568  Array<OneD, NekDouble> temperature(nPts, 0.0);
569  m_varConv->GetTemperature(inaverg, temperature);
570  GetViscosityAndThermalCondFromTemp(temperature, mu, thermalConductivity);
571 
572  std::vector<NekDouble> inAvgTmp(nConvectiveFields);
573  std::vector<NekDouble> inTmp(nConvectiveFields);
574  std::vector<NekDouble> outTmp(nConvectiveFields);
575  for (int d = 0; d < nDim; ++d)
576  {
577  for (int nderiv = 0; nderiv < nDim; ++nderiv)
578  {
579  for (size_t p = 0; p < nPts; ++p)
580  {
581  // rearrenge data
582  for (int f = 0; f < nConvectiveFields; ++f)
583  {
584  inAvgTmp[f] = inaverg[f][p];
585  inTmp[f] = inarray[f][p];
586  }
587 
588  GetViscousFluxBilinearFormKernel(nDim, d, nderiv,
589  inAvgTmp.data(), inTmp.data(),
590  mu[p], outTmp.data());
591 
592  for (int f = 0; f < nConvectiveFields; ++f)
593  {
594  outarray[d][f][p] += normal[nderiv][p] * outTmp[f];
595  }
596  }
597  }
598  }
599 }
600 
601 /**
602  * @brief Return the penalty vector for the LDGNS diffusion problem.
603  */
605  const Array<OneD, const Array<OneD, NekDouble>> &uFwd,
606  const Array<OneD, const Array<OneD, NekDouble>> &uBwd,
607  Array<OneD, Array<OneD, NekDouble>> &penaltyCoeff)
608 {
609  unsigned int nTracePts = uFwd[0].size();
610 
611  // Compute average temperature
612  unsigned int nVariables = uFwd.size();
613  Array<OneD, NekDouble> tAve{nTracePts, 0.0};
614  Vmath::Svtsvtp(nTracePts, 0.5, uFwd[nVariables - 1], 1, 0.5,
615  uBwd[nVariables - 1], 1, tAve, 1);
616 
617  // Get average viscosity and thermal conductivity
618  Array<OneD, NekDouble> muAve{nTracePts, 0.0};
619  Array<OneD, NekDouble> tcAve{nTracePts, 0.0};
620 
621  GetViscosityAndThermalCondFromTemp(tAve, muAve, tcAve);
622 
623  // Compute penalty term
624  for (int i = 0; i < nVariables; ++i)
625  {
626  // Get jump of u variables
627  Vmath::Vsub(nTracePts, uFwd[i], 1, uBwd[i], 1, penaltyCoeff[i], 1);
628  // Multiply by variable coefficient = {coeff} ( u^+ - u^- )
629  if (i < nVariables - 1)
630  {
631  Vmath::Vmul(nTracePts, muAve, 1, penaltyCoeff[i], 1,
632  penaltyCoeff[i], 1);
633  }
634  else
635  {
636  Vmath::Vmul(nTracePts, tcAve, 1, penaltyCoeff[i], 1,
637  penaltyCoeff[i], 1);
638  }
639  }
640 }
641 
642 /**
643  * @brief Update viscosity
644  * todo: add artificial viscosity here
645  */
647  const Array<OneD, NekDouble> &temperature, Array<OneD, NekDouble> &mu,
648  Array<OneD, NekDouble> &thermalCond)
649 {
650  auto nPts = temperature.size();
651 
652  for (size_t p = 0; p < nPts; ++p)
653  {
654  GetViscosityAndThermalCondFromTempKernel(temperature[p], mu[p],
655  thermalCond[p]);
656  }
657 
658  // Add artificial viscosity if wanted
659  // move this above and add in kernel
660  if (m_is_shockCaptPhys)
661  {
662  auto nTracePts = m_fields[0]->GetTrace()->GetTotPoints();
663  if (nPts != nTracePts)
664  {
665  Vmath::Vadd(nPts, mu, 1, m_varConv->GetAv(), 1, mu, 1);
666  }
667  else
668  {
669  Vmath::Vadd(nPts, mu, 1, m_varConv->GetAvTrace(), 1, mu, 1);
670  }
671 
672  // Update thermal conductivity
673  NekDouble tRa = m_Cp / m_Prandtl;
674  Vmath::Smul(nPts, tRa, mu, 1, thermalCond, 1);
675  }
676 }
677 
678 /**
679  * @brief Get divergence and curl squared
680  *
681  * @param input
682  * fields -> expansion list pointer
683  * cnsVar -> conservative variables
684  * cnsVarFwd -> forward trace of conservative variables
685  * cnsVarBwd -> backward trace of conservative variables
686  * @paran output
687  * divSquare -> divergence
688  * curlSquare -> curl squared
689  *
690  */
693  const Array<OneD, Array<OneD, NekDouble>> &cnsVar,
695  const Array<OneD, Array<OneD, NekDouble>> &cnsVarFwd,
696  const Array<OneD, Array<OneD, NekDouble>> &cnsVarBwd)
697 {
698  auto nDim = fields[0]->GetCoordim(0);
699  auto nVar = cnsVar.size();
700  auto nPts = cnsVar[0].size();
701  auto nPtsTrc = cnsVarFwd[0].size();
702 
703  // These should be allocated once
704  Array<OneD, Array<OneD, NekDouble>> primVar(nVar - 1), primVarFwd(nVar - 1),
705  primVarBwd(nVar - 1);
706 
707  for (unsigned short d = 0; d < nVar - 2; ++d)
708  {
709  primVar[d] = Array<OneD, NekDouble>(nPts, 0.0);
710  primVarFwd[d] = Array<OneD, NekDouble>(nPtsTrc, 0.0);
711  primVarBwd[d] = Array<OneD, NekDouble>(nPtsTrc, 0.0);
712  }
713  auto ergLoc = nVar - 2;
714  primVar[ergLoc] = Array<OneD, NekDouble>(nPts, 0.0);
715  primVarFwd[ergLoc] = Array<OneD, NekDouble>(nPtsTrc, 0.0);
716  primVarBwd[ergLoc] = Array<OneD, NekDouble>(nPtsTrc, 0.0);
717 
718  // Get primitive variables [u,v,w,T=0]
719  // Possibly should be changed to [rho,u,v,w,T] to make IP and LDGNS more
720  // consistent with each other
721  for (unsigned short d = 0; d < nVar - 2; ++d)
722  {
723  // Volume
724  for (size_t p = 0; p < nPts; ++p)
725  {
726  primVar[d][p] = cnsVar[d + 1][p] / cnsVar[0][p];
727  }
728  // Trace
729  for (size_t p = 0; p < nPtsTrc; ++p)
730  {
731  primVarFwd[d][p] = cnsVarFwd[d + 1][p] / cnsVarFwd[0][p];
732  primVarBwd[d][p] = cnsVarBwd[d + 1][p] / cnsVarBwd[0][p];
733  }
734  }
735 
736  // this should be allocated once
738  for (unsigned short j = 0; j < nDim; ++j)
739  {
740  primVarDer[j] = Array<OneD, Array<OneD, NekDouble>>(nVar - 1);
741  for (unsigned short i = 0; i < nVar - 1; ++i)
742  {
743  primVarDer[j][i] = Array<OneD, NekDouble>(nPts, 0.0);
744  }
745  }
746 
747  // Get derivative tensor
748  m_diffusion->DiffuseCalcDerivative(fields, primVar, primVarDer, primVarFwd,
749  primVarBwd);
750 
751  // Get div curl squared
752  GetDivCurlFromDvelT(primVarDer, div, curlSquare);
753 }
754 
755 /**
756  * @brief Get divergence and curl from velocity derivative tensor
757  *
758  */
761  Array<OneD, NekDouble> &curlSquare)
762 {
763  auto nDim = pVarDer.size();
764  auto nPts = div.size();
765 
766  // div velocity
767  for (size_t p = 0; p < nPts; ++p)
768  {
769  NekDouble divTmp = 0;
770  for (unsigned short j = 0; j < nDim; ++j)
771  {
772  divTmp += pVarDer[j][j][p];
773  }
774  div[p] = divTmp;
775  }
776 
777  // |curl velocity| ** 2
778  if (nDim > 2)
779  {
780  for (size_t p = 0; p < nPts; ++p)
781  {
782  // curl[0] 3/2 - 2/3
783  NekDouble curl032 = pVarDer[2][1][p]; // load 1x
784  NekDouble curl023 = pVarDer[1][2][p]; // load 1x
785  NekDouble curl0 = curl032 - curl023;
786  // square curl[0]
787  NekDouble curl0sqr = curl0 * curl0;
788 
789  // curl[1] 3/1 - 1/3
790  NekDouble curl131 = pVarDer[2][0][p]; // load 1x
791  NekDouble curl113 = pVarDer[0][2][p]; // load 1x
792  NekDouble curl1 = curl131 - curl113;
793  // square curl[1]
794  NekDouble curl1sqr = curl1 * curl1;
795 
796  // curl[2] 1/2 - 2/1
797  NekDouble curl212 = pVarDer[0][1][p]; // load 1x
798  NekDouble curl221 = pVarDer[1][0][p]; // load 1x
799  NekDouble curl2 = curl212 - curl221;
800  // square curl[2]
801  NekDouble curl2sqr = curl2 * curl2;
802 
803  // reduce
804  curl0sqr += curl1sqr + curl2sqr;
805  // store
806  curlSquare[p] = curl0sqr; // store 1x
807  }
808  }
809  else if (nDim > 1)
810  {
811  for (size_t p = 0; p < nPts; ++p)
812  {
813  // curl[2] 1/2
814  NekDouble c212 = pVarDer[0][1][p]; // load 1x
815  // curl[2] 2/1
816  NekDouble c221 = pVarDer[1][0][p]; // load 1x
817  // curl[2] 1/2 - 2/1
818  NekDouble curl = c212 - c221;
819  // square curl[2]
820  curlSquare[p] = curl * curl; // store 1x
821  }
822  }
823  else
824  {
825  Vmath::Fill(nPts, 0.0, curlSquare, 1);
826  }
827 }
828 
830  std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
831  std::vector<std::string> &variables)
832 {
833  bool extraFields;
834  m_session->MatchSolverInfo("OutputExtraFields", "True", extraFields, true);
835  if (extraFields)
836  {
837  const int nPhys = m_fields[0]->GetNpoints();
838  const int nCoeffs = m_fields[0]->GetNcoeffs();
840 
841  for (int i = 0; i < m_fields.size(); ++i)
842  {
843  cnsVar[i] = m_fields[i]->GetPhys();
844  }
845 
848  for (int i = 0; i < m_spacedim; ++i)
849  {
850  velocity[i] = Array<OneD, NekDouble>(nPhys);
851  velFwd[i] = Array<OneD, NekDouble>(nCoeffs);
852  }
853 
854  Array<OneD, NekDouble> pressure(nPhys), temperature(nPhys);
855  Array<OneD, NekDouble> entropy(nPhys);
856  Array<OneD, NekDouble> soundspeed(nPhys), mach(nPhys);
857  Array<OneD, NekDouble> sensor(nPhys), SensorKappa(nPhys);
858 
859  m_varConv->GetVelocityVector(cnsVar, velocity);
860  m_varConv->GetPressure(cnsVar, pressure);
861  m_varConv->GetTemperature(cnsVar, temperature);
862  m_varConv->GetEntropy(cnsVar, entropy);
863  m_varConv->GetSoundSpeed(cnsVar, soundspeed);
864  m_varConv->GetMach(cnsVar, soundspeed, mach);
865 
866  int sensorOffset;
867  m_session->LoadParameter("SensorOffset", sensorOffset, 1);
868  m_varConv->GetSensor(m_fields[0], cnsVar, sensor, SensorKappa,
869  sensorOffset);
870 
871  Array<OneD, NekDouble> pFwd(nCoeffs), TFwd(nCoeffs);
872  Array<OneD, NekDouble> sFwd(nCoeffs);
873  Array<OneD, NekDouble> aFwd(nCoeffs), mFwd(nCoeffs);
874  Array<OneD, NekDouble> sensFwd(nCoeffs);
875 
876  string velNames[3] = {"u", "v", "w"};
877  for (int i = 0; i < m_spacedim; ++i)
878  {
879  m_fields[0]->FwdTransLocalElmt(velocity[i], velFwd[i]);
880  variables.push_back(velNames[i]);
881  fieldcoeffs.push_back(velFwd[i]);
882  }
883 
884  m_fields[0]->FwdTransLocalElmt(pressure, pFwd);
885  m_fields[0]->FwdTransLocalElmt(temperature, TFwd);
886  m_fields[0]->FwdTransLocalElmt(entropy, sFwd);
887  m_fields[0]->FwdTransLocalElmt(soundspeed, aFwd);
888  m_fields[0]->FwdTransLocalElmt(mach, mFwd);
889  m_fields[0]->FwdTransLocalElmt(sensor, sensFwd);
890 
891  variables.push_back("p");
892  variables.push_back("T");
893  variables.push_back("s");
894  variables.push_back("a");
895  variables.push_back("Mach");
896  variables.push_back("Sensor");
897  fieldcoeffs.push_back(pFwd);
898  fieldcoeffs.push_back(TFwd);
899  fieldcoeffs.push_back(sFwd);
900  fieldcoeffs.push_back(aFwd);
901  fieldcoeffs.push_back(mFwd);
902  fieldcoeffs.push_back(sensFwd);
903 
905  {
906  // reuse pressure
907  Array<OneD, NekDouble> sensorFwd(nCoeffs);
908  m_artificialDiffusion->GetArtificialViscosity(cnsVar, pressure);
909  m_fields[0]->FwdTransLocalElmt(pressure, sensorFwd);
910 
911  variables.push_back("ArtificialVisc");
912  fieldcoeffs.push_back(sensorFwd);
913  }
914 
915  if (m_is_shockCaptPhys)
916  {
917 
918  Array<OneD, Array<OneD, NekDouble>> cnsVarFwd(m_fields.size()),
919  cnsVarBwd(m_fields.size());
920 
921  for (int i = 0; i < m_fields.size(); ++i)
922  {
923  cnsVarFwd[i] = Array<OneD, NekDouble>(GetTraceTotPoints());
924  cnsVarBwd[i] = Array<OneD, NekDouble>(GetTraceTotPoints());
925  m_fields[i]->GetFwdBwdTracePhys(cnsVar[i], cnsVarFwd[i],
926  cnsVarBwd[i]);
927  }
928 
929  Array<OneD, NekDouble> div(nPhys), curlSquare(nPhys);
930  GetDivCurlSquared(m_fields, cnsVar, div, curlSquare, cnsVarFwd,
931  cnsVarBwd);
932 
933  Array<OneD, NekDouble> divFwd(nCoeffs, 0.0);
934  m_fields[0]->FwdTransLocalElmt(div, divFwd);
935  variables.push_back("div");
936  fieldcoeffs.push_back(divFwd);
937 
938  Array<OneD, NekDouble> curlFwd(nCoeffs, 0.0);
939  m_fields[0]->FwdTransLocalElmt(curlSquare, curlFwd);
940  variables.push_back("curl^2");
941  fieldcoeffs.push_back(curlFwd);
942 
943  m_varConv->SetAv(m_fields, cnsVar, div, curlSquare);
944 
945  Array<OneD, NekDouble> muavFwd(nCoeffs);
946  m_fields[0]->FwdTransLocalElmt(m_varConv->GetAv(), muavFwd);
947  variables.push_back("ArtificialVisc");
948  fieldcoeffs.push_back(muavFwd);
949  }
950  }
951 }
952 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
virtual void v_InitObject(bool DeclareFields=true)
Initialization object for CompressibleFlowSystem class.
SolverUtils::DiffusionSharedPtr m_diffusion
VariableConverterSharedPtr m_varConv
void SetBoundaryConditionsBwdWeight()
Set up a weight on physical boundaries for boundary condition applications.
ArtificialDiffusionSharedPtr m_artificialDiffusion
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
virtual void v_GetViscousFluxVector(const Array< OneD, const Array< OneD, NekDouble >> &physfield, TensorOfArray3D< NekDouble > &derivatives, TensorOfArray3D< NekDouble > &viscousTensor)
Return the flux vector for the LDG diffusion problem.
void GetDivCurlFromDvelT(const TensorOfArray3D< NekDouble > &pVarDer, Array< OneD, NekDouble > &div, Array< OneD, NekDouble > &curlSquare)
Get divergence and curl from velocity derivative tensor.
virtual void v_InitObject(bool DeclareField=true) override
Initialization object for CompressibleFlowSystem class.
void SpecialBndTreat(Array< OneD, Array< OneD, NekDouble >> &consvar)
For very special treatment. For general boundaries it does nothing But for WallViscous and WallAdiaba...
void GetDivCurlSquared(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble >> &cnsVar, Array< OneD, NekDouble > &div, Array< OneD, NekDouble > &curlSquare, const Array< OneD, Array< OneD, NekDouble >> &cnsVarFwd, const Array< OneD, Array< OneD, NekDouble >> &cnsVarBwd)
Get divergence and curl squared.
virtual void v_DoDiffusion(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
Apply artificial diffusion (Laplacian operator)
NekDouble m_thermalConductivityRef
bool m_is_mu_variable
flag to switch between constant viscosity and Sutherland an enum could be added for more options
bool m_is_shockCaptPhys
flag for shock capturing switch on/off an enum could be added for more options
virtual void v_GetViscousFluxVectorDeAlias(const Array< OneD, const Array< OneD, NekDouble >> &physfield, TensorOfArray3D< NekDouble > &derivatives, TensorOfArray3D< NekDouble > &viscousTensor)
Return the flux vector for the LDG diffusion problem.
void GetViscosityAndThermalCondFromTemp(const Array< OneD, NekDouble > &temperature, Array< OneD, NekDouble > &mu, Array< OneD, NekDouble > &thermalCond)
Update viscosity todo: add artificial viscosity here.
bool m_is_diffIP
flag to switch between IP and LDG an enum could be added for more options
void GetViscousSymmtrFluxConservVar(const int nSpaceDim, const Array< OneD, Array< OneD, NekDouble >> &inaverg, const Array< OneD, Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &outarray, Array< OneD, int > &nonZeroIndex, const Array< OneD, Array< OneD, NekDouble >> &normals)
Calculate and return the Symmetric flux in IP method.
virtual void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables) override
void GetViscousFluxBilinearFormKernel(const unsigned short nDim, const unsigned short FluxDirection, const unsigned short DerivDirection, const T *inaverg, const T *injumpp, const T &mu, T *outarray)
Calculate diffusion flux using the Jacobian form.
virtual void v_GetFluxPenalty(const Array< OneD, const Array< OneD, NekDouble >> &uFwd, const Array< OneD, const Array< OneD, NekDouble >> &uBwd, Array< OneD, Array< OneD, NekDouble >> &penaltyCoeff)
Return the penalty vector for the LDGNS diffusion problem.
void GetViscosityAndThermalCondFromTempKernel(const T &temperature, T &mu, T &thermalCond)
int m_spacedim
Spatial dimension (>= expansion dim).
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetNpoints()
Base class for unsteady solvers.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:41
EquationSystemFactory & GetEquationSystemFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
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)
vvtvvtp (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 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 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
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:419