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 
35 #include <CompressibleFlowSolver/EquationSystems/NavierStokesCFE.h>
36 
37 using namespace std;
38 
39 namespace Nektar
40 {
41  string NavierStokesCFE::className =
42  SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
43  "NavierStokesCFE", NavierStokesCFE::create,
44  "NavierStokes equations in conservative variables.");
45 
46  NavierStokesCFE::NavierStokesCFE(
47  const LibUtilities::SessionReaderSharedPtr& pSession,
48  const SpatialDomains::MeshGraphSharedPtr& pGraph)
49  : UnsteadySystem(pSession, pGraph),
50  CompressibleFlowSystem(pSession, pGraph)
51  {
52  }
53 
54  NavierStokesCFE::~NavierStokesCFE()
55  {
56 
57  }
58 
59  /**
60  * @brief Initialization object for CompressibleFlowSystem class.
61  */
62  void NavierStokesCFE::v_InitObject()
63  {
64  CompressibleFlowSystem::v_InitObject();
65 
66  // rest of initialisation is in this routine so it can also be called
67  // in NavierStokesImplicitCFE initialisation
68  InitObject_Explicit();
69  }
70 
71  void NavierStokesCFE::InitObject_Explicit()
72  {
73  // Get gas constant from session file and compute Cp
74  NekDouble gasConstant;
75  m_session->LoadParameter ("GasConstant", gasConstant, 287.058);
76  m_Cp = m_gamma / (m_gamma - 1.0) * gasConstant;
77  m_Cv = m_Cp / m_gamma;
78 
79  m_session->LoadParameter ("Twall", m_Twall, 300.15);
80 
81  // Viscosity
82  int nPts = m_fields[0]->GetNpoints();
83  m_session->LoadSolverInfo("ViscosityType", m_ViscosityType, "Constant");
84  m_session->LoadParameter ("mu", m_muRef, 1.78e-05);
85  m_mu = Array<OneD, NekDouble>(nPts, m_muRef);
86  if (m_ViscosityType == "Variable")
87  {
88  m_is_mu_variable = true;
89  }
90 
91  // Thermal conductivity or Prandtl
92  if ( m_session->DefinesParameter("thermalConductivity"))
93  {
94  ASSERTL0( !m_session->DefinesParameter("Pr"),
95  "Cannot define both Pr and thermalConductivity.");
96 
97  m_session->LoadParameter ("thermalConductivity",
98  m_thermalConductivityRef);
99  m_Prandtl = m_Cp * m_muRef / m_thermalConductivityRef;
100  }
101  else
102  {
103  m_session->LoadParameter ("Pr", m_Prandtl, 0.72);
104  m_thermalConductivityRef = m_Cp * m_muRef / m_Prandtl;
105  }
106  m_thermalConductivity =
107  Array<OneD, NekDouble>(nPts, m_thermalConductivityRef);
108 
109  // Artificial viscosity parameter
110  m_session->LoadParameter("mu0", m_mu0, 1.0);
111 
112  // load smoothing tipe
113  m_session->LoadSolverInfo("Smoothing", m_smoothing, "Off");
114  if (m_smoothing == "C0")
115  {
116  m_C0ProjectExp = MemoryManager<MultiRegions::ContField>::
117  AllocateSharedPtr(m_session,m_graph,m_session->GetVariable(0));
118  }
119  // load physical sensor type
120  m_session->LoadSolverInfo("PhysicalSensorType", m_physicalSensorType,
121  "Off");
122 
123  string diffName;
124  m_session->LoadSolverInfo("DiffusionType", diffName, "LDGNS");
125 
126  m_diffusion = SolverUtils::GetDiffusionFactory()
127  .CreateInstance(diffName, diffName);
128  if ("InteriorPenalty" == diffName)
129  {
130  m_is_diffIP = true;
131  SetBoundaryConditionsBwdWeight();
132  }
133 
134  if ("LDGNS" == diffName||
135  "LDGNS3DHomogeneous1D" == diffName)
136  {
137  m_diffusion->SetFluxPenaltyNS(&NavierStokesCFE::
138  v_GetFluxPenalty, this);
139  }
140 
141  if (m_specHP_dealiasing)
142  {
143  m_diffusion->SetFluxVectorNS(
144  &NavierStokesCFE::v_GetViscousFluxVectorDeAlias,
145  this);
146  }
147  else
148  {
149  m_diffusion->SetFluxVectorNS(&NavierStokesCFE::
150  v_GetViscousFluxVector, this);
151  }
152 
153  m_diffusion->SetDiffusionFluxCons(
154  &NavierStokesCFE::GetViscousFluxVectorConservVar<false>, this);
155 
156  m_diffusion->SetDiffusionFluxConsTrace(
157  &NavierStokesCFE::GetViscousFluxVectorConservVar<true>, this);
158 
159  m_diffusion->SetSpecialBndTreat(
160  &NavierStokesCFE::SpecialBndTreat, this);
161 
162  m_diffusion->SetDiffusionSymmFluxCons(
163  &NavierStokesCFE::GetViscousSymmtrFluxConservVar, this);
164 
165  if (m_shockCaptureType != "Off")
166  {
167  m_diffusion->SetArtificialDiffusionVector(
168  &NavierStokesCFE::GetArtificialViscosity, this);
169  }
170 
171  m_diffusion->SetCalcViscosity(
172  &NavierStokesCFE::CalcViscosity, this);
173 
174  // Concluding initialisation of diffusion operator
175  m_diffusion->InitObject (m_session, m_fields);
176 
177  }
178 
179  void NavierStokesCFE::v_ExtraFldOutput(
180  std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
181  std::vector<std::string> &variables)
182  {
183  bool extraFields;
184  m_session->MatchSolverInfo("OutputExtraFields","True",
185  extraFields, true);
186  if (extraFields)
187  {
188  const int nPhys = m_fields[0]->GetNpoints();
189  const int nCoeffs = m_fields[0]->GetNcoeffs();
190  Array<OneD, Array<OneD, NekDouble> > tmp(m_fields.size());
191 
192  for (int i = 0; i < m_fields.size(); ++i)
193  {
194  tmp[i] = m_fields[i]->GetPhys();
195  }
196 
197  Array<OneD, Array<OneD, NekDouble> > velocity(m_spacedim);
198  Array<OneD, Array<OneD, NekDouble> > velFwd (m_spacedim);
199  for (int i = 0; i < m_spacedim; ++i)
200  {
201  velocity[i] = Array<OneD, NekDouble> (nPhys);
202  velFwd[i] = Array<OneD, NekDouble> (nCoeffs);
203  }
204 
205  Array<OneD, NekDouble> pressure(nPhys), temperature(nPhys);
206  Array<OneD, NekDouble> entropy(nPhys);
207  Array<OneD, NekDouble> soundspeed(nPhys), mach(nPhys);
208  Array<OneD, NekDouble> sensor(nPhys), SensorKappa(nPhys);
209 
210  m_varConv->GetVelocityVector(tmp, velocity);
211  m_varConv->GetPressure (tmp, pressure);
212  m_varConv->GetTemperature(tmp, temperature);
213  m_varConv->GetEntropy (tmp, entropy);
214  m_varConv->GetSoundSpeed(tmp, soundspeed);
215  m_varConv->GetMach (tmp, soundspeed, mach);
216 
217  int sensorOffset;
218  m_session->LoadParameter ("SensorOffset", sensorOffset, 1);
219  m_varConv->GetSensor (m_fields[0], tmp, sensor, SensorKappa,
220  sensorOffset);
221 
222  Array<OneD, NekDouble> pFwd(nCoeffs), TFwd(nCoeffs);
223  Array<OneD, NekDouble> sFwd(nCoeffs);
224  Array<OneD, NekDouble> aFwd(nCoeffs), mFwd(nCoeffs);
225  Array<OneD, NekDouble> sensFwd(nCoeffs);
226 
227  string velNames[3] = {"u", "v", "w"};
228  for (int i = 0; i < m_spacedim; ++i)
229  {
230  m_fields[0]->FwdTrans_IterPerExp(velocity[i], velFwd[i]);
231  variables.push_back(velNames[i]);
232  fieldcoeffs.push_back(velFwd[i]);
233  }
234 
235  m_fields[0]->FwdTrans_IterPerExp(pressure, pFwd);
236  m_fields[0]->FwdTrans_IterPerExp(temperature,TFwd);
237  m_fields[0]->FwdTrans_IterPerExp(entropy, sFwd);
238  m_fields[0]->FwdTrans_IterPerExp(soundspeed, aFwd);
239  m_fields[0]->FwdTrans_IterPerExp(mach, mFwd);
240  m_fields[0]->FwdTrans_IterPerExp(sensor, sensFwd);
241 
242  variables.push_back ("p");
243  variables.push_back ("T");
244  variables.push_back ("s");
245  variables.push_back ("a");
246  variables.push_back ("Mach");
247  variables.push_back ("Sensor");
248  fieldcoeffs.push_back(pFwd);
249  fieldcoeffs.push_back(TFwd);
250  fieldcoeffs.push_back(sFwd);
251  fieldcoeffs.push_back(aFwd);
252  fieldcoeffs.push_back(mFwd);
253  fieldcoeffs.push_back(sensFwd);
254 
255  if (m_artificialDiffusion)
256  {
257  // Get min h/p
258  // m_artificialDiffusion->m_hOverP = GetElmtMinHP();
259  // reuse pressure
260  Array<OneD, NekDouble> sensorFwd(nCoeffs);
261  m_artificialDiffusion->GetArtificialViscosity(tmp, pressure);
262  m_fields[0]->FwdTrans_IterPerExp(pressure, sensorFwd);
263 
264  variables.push_back ("ArtificialVisc");
265  fieldcoeffs.push_back(sensorFwd);
266 
267  }
268 
269  if (m_shockCaptureType == "Physical")
270  {
271  Array<OneD, NekDouble> muavFwd(nCoeffs);
272  m_fields[0]->FwdTrans_IterPerExp(m_muav, muavFwd);
273  variables.push_back ("ArtificialVisc");
274  fieldcoeffs.push_back(muavFwd);
275 
276  // Debug Ducros
277  // div square
278  Array<OneD, NekDouble> dv2Fwd(nCoeffs);
279  m_fields[0]->FwdTrans_IterPerExp(m_diffusion->m_divVelSquare,
280  dv2Fwd);
281  variables.push_back ("divVelSquare");
282  fieldcoeffs.push_back(dv2Fwd);
283  // curl square
284  Array<OneD, NekDouble> cv2Fwd(nCoeffs);
285  m_fields[0]->FwdTrans_IterPerExp(m_diffusion->m_curlVelSquare,
286  cv2Fwd);
287  variables.push_back ("curlVelSquare");
288  fieldcoeffs.push_back(cv2Fwd);
289  // Ducros
290  Array<OneD, NekDouble> duc(nPhys,1.0);
291  Ducros(duc);
292  Array<OneD, NekDouble> ducFwd(nCoeffs);
293  m_fields[0]->FwdTrans_IterPerExp(duc, ducFwd);
294  variables.push_back ("Ducros");
295  fieldcoeffs.push_back(ducFwd);
296 
297  }
298  }
299  }
300 
301  void NavierStokesCFE::v_DoDiffusion(
302  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
303  Array<OneD, Array<OneD, NekDouble>> &outarray,
304  const Array<OneD, const Array<OneD, NekDouble>> &pFwd,
305  const Array<OneD, const Array<OneD, NekDouble>> &pBwd)
306  {
307  size_t nvariables = inarray.size();
308  size_t npoints = GetNpoints();
309  size_t nTracePts = GetTraceTotPoints();
310 
311  // this should be preallocated
312  Array<OneD, Array<OneD, NekDouble> > outarrayDiff(nvariables);
313  for (int i = 0; i < nvariables; ++i)
314  {
315  outarrayDiff[i] = Array<OneD, NekDouble>(npoints, 0.0);
316  }
317 
318  // get artificial viscosity
319  if (m_shockCaptureType == "Physical" && m_CalcPhysicalAV)
320  {
321  GetPhysicalAV(inarray);
322  }
323 
324  if (m_is_diffIP)
325  {
326  if (m_bndEvaluateTime < 0.0)
327  {
328  NEKERROR(ErrorUtil::efatal, "m_bndEvaluateTime not setup");
329  }
330  m_diffusion->Diffuse(nvariables, m_fields, inarray, outarrayDiff,
331  m_bndEvaluateTime,
332  pFwd, pBwd);
333  for (int i = 0; i < nvariables; ++i)
334  {
335  Vmath::Vadd(npoints,
336  outarrayDiff[i], 1,
337  outarray[i], 1,
338  outarray[i], 1);
339  }
340  }
341  else
342  {
343  Array<OneD, Array<OneD, NekDouble> > inarrayDiff(nvariables-1);
344  Array<OneD, Array<OneD, NekDouble> > inFwd(nvariables-1);
345  Array<OneD, Array<OneD, NekDouble> > inBwd(nvariables-1);
346 
347 
348  for (int i = 0; i < nvariables - 1; ++i)
349  {
350  inarrayDiff[i] = Array<OneD, NekDouble>{npoints};
351  inFwd[i] = Array<OneD, NekDouble>{nTracePts};
352  inBwd[i] = Array<OneD, NekDouble>{nTracePts};
353  }
354 
355  // Extract pressure
356  // (use inarrayDiff[0] as a temporary storage for the pressure)
357  m_varConv->GetPressure(inarray, inarrayDiff[0]);
358 
359  // Extract temperature
360  m_varConv->GetTemperature(inarray, inarrayDiff[nvariables - 2]);
361 
362  // Extract velocities
363  m_varConv->GetVelocityVector(inarray, inarrayDiff);
364 
365  // Repeat calculation for trace space
366  if (pFwd == NullNekDoubleArrayOfArray ||
367  pBwd == NullNekDoubleArrayOfArray)
368  {
369  inFwd = NullNekDoubleArrayOfArray;
370  inBwd = NullNekDoubleArrayOfArray;
371  }
372  else
373  {
374  m_varConv->GetPressure(pFwd, inFwd[0]);
375  m_varConv->GetPressure(pBwd, inBwd[0]);
376 
377  m_varConv->GetTemperature(pFwd, inFwd[nvariables - 2]);
378  m_varConv->GetTemperature(pBwd, inBwd[nvariables - 2]);
379 
380  m_varConv->GetVelocityVector(pFwd, inFwd);
381  m_varConv->GetVelocityVector(pBwd, inBwd);
382  }
383 
384  // Diffusion term in physical rhs form
385  m_diffusion->Diffuse(nvariables, m_fields, inarrayDiff,
386  outarrayDiff, inFwd, inBwd);
387 
388  for (int i = 0; i < nvariables; ++i)
389  {
390  Vmath::Vadd(npoints,
391  outarrayDiff[i], 1,
392  outarray[i], 1,
393  outarray[i], 1);
394  }
395  }
396  }
397 
398 
399 
400  /**
401  * @brief Return the flux vector for the LDG diffusion problem.
402  * \todo Complete the viscous flux vector
403  */
404  void NavierStokesCFE::v_GetViscousFluxVector(
405  const Array<OneD, const Array<OneD, NekDouble>> &physfield,
406  TensorOfArray3D<NekDouble> &derivativesO1,
407  TensorOfArray3D<NekDouble> &viscousTensor)
408  {
409  // Auxiliary variables
410  size_t nScalar = physfield.size();
411  size_t nPts = physfield[0].size();
412  Array<OneD, NekDouble > divVel (nPts, 0.0);
413 
414  // Stokes hypothesis
415  const NekDouble lambda = -2.0/3.0;
416 
417  // Update viscosity and thermal conductivity
418  GetViscosityAndThermalCondFromTemp(physfield[nScalar-1], m_mu,
419  m_thermalConductivity);
420 
421  // Add artificial viscosity if wanted
422  if (m_shockCaptureType == "Physical")
423  {
424  // Apply Ducros sensor
425  if (m_physicalSensorType == "Ducros" && m_CalcPhysicalAV)
426  {
427  Ducros(m_muav);
428  }
429  // Apply approximate c0 smoothing
430  if (m_smoothing == "C0" && m_CalcPhysicalAV)
431  {
432  C0Smooth(m_muav);
433  }
434  Vmath::Vadd(nPts, m_mu, 1, m_muav, 1, m_mu, 1);
435  // Freeze AV for Implicit time stepping
436  if (m_explicitDiffusion == false)
437  {
438  m_CalcPhysicalAV = false;
439  }
440  }
441 
442  // Velocity divergence
443  for (int j = 0; j < m_spacedim; ++j)
444  {
445  Vmath::Vadd(nPts, divVel, 1, derivativesO1[j][j], 1,
446  divVel, 1);
447  }
448 
449  // Velocity divergence scaled by lambda * mu
450  Vmath::Smul(nPts, lambda, divVel, 1, divVel, 1);
451  Vmath::Vmul(nPts, m_mu, 1, divVel, 1, divVel, 1);
452 
453  // Viscous flux vector for the rho equation = 0
454  for (int i = 0; i < m_spacedim; ++i)
455  {
456  Vmath::Zero(nPts, viscousTensor[i][0], 1);
457  }
458 
459  // Viscous stress tensor (for the momentum equations)
460  for (int i = 0; i < m_spacedim; ++i)
461  {
462  for (int j = i; j < m_spacedim; ++j)
463  {
464  Vmath::Vadd(nPts, derivativesO1[i][j], 1,
465  derivativesO1[j][i], 1,
466  viscousTensor[i][j+1], 1);
467 
468  Vmath::Vmul(nPts, m_mu, 1,
469  viscousTensor[i][j+1], 1,
470  viscousTensor[i][j+1], 1);
471 
472  if (i == j)
473  {
474  // Add divergence term to diagonal
475  Vmath::Vadd(nPts, viscousTensor[i][j+1], 1,
476  divVel, 1,
477  viscousTensor[i][j+1], 1);
478  }
479  else
480  {
481  // Copy to make symmetric
482  Vmath::Vcopy(nPts, viscousTensor[i][j+1], 1,
483  viscousTensor[j][i+1], 1);
484  }
485  }
486  }
487 
488  // Terms for the energy equation
489  for (int i = 0; i < m_spacedim; ++i)
490  {
491  Vmath::Zero(nPts, viscousTensor[i][m_spacedim+1], 1);
492  // u_j * tau_ij
493  for (int j = 0; j < m_spacedim; ++j)
494  {
495  Vmath::Vvtvp(nPts, physfield[j], 1,
496  viscousTensor[i][j+1], 1,
497  viscousTensor[i][m_spacedim+1], 1,
498  viscousTensor[i][m_spacedim+1], 1);
499  }
500  // Add k*T_i
501  Vmath::Vvtvp(nPts, m_thermalConductivity, 1,
502  derivativesO1[i][m_spacedim], 1,
503  viscousTensor[i][m_spacedim+1], 1,
504  viscousTensor[i][m_spacedim+1], 1);
505  }
506  }
507 
508  /**
509  * @brief Return the flux vector for the LDG diffusion problem.
510  * \todo Complete the viscous flux vector
511  */
512  void NavierStokesCFE::v_GetViscousFluxVectorDeAlias(
513  const Array<OneD, const Array<OneD, NekDouble>> &physfield,
514  TensorOfArray3D<NekDouble> &derivativesO1,
515  TensorOfArray3D<NekDouble> &viscousTensor)
516  {
517  // Factor to rescale 1d points in dealiasing.
518  NekDouble OneDptscale = 2;
519  // Get number of points to dealias a cubic non-linearity
520  size_t nScalar = physfield.size();
521  int nPts = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
522  size_t nPts_orig = physfield[0].size();
523 
524  // Auxiliary variables
525  Array<OneD, NekDouble > divVel (nPts, 0.0);
526 
527  // Stokes hypothesis
528  const NekDouble lambda = -2.0/3.0;
529 
530  // Update viscosity and thermal conductivity
531  GetViscosityAndThermalCondFromTemp(physfield[nScalar-1], m_mu,
532  m_thermalConductivity);
533 
534  // Add artificial viscosity if wanted
535  if (m_shockCaptureType == "Physical")
536  {
537  // Apply Ducros sensor
538  if (m_physicalSensorType == "Ducros" && m_CalcPhysicalAV)
539  {
540  Ducros(m_muav);
541  }
542  // Apply approximate c0 smoothing
543  if (m_smoothing == "C0" && m_CalcPhysicalAV)
544  {
545  C0Smooth(m_muav);
546  }
547  Vmath::Vadd(nPts, m_mu, 1, m_muav, 1, m_mu, 1);
548  // Freeze AV for Implicit time stepping
549  if (m_explicitDiffusion == false)
550  {
551  m_CalcPhysicalAV = false;
552  }
553  }
554 
555  // Interpolate inputs and initialise interpolated output
556  Array<OneD, Array<OneD, NekDouble> > vel_interp(m_spacedim);
557  TensorOfArray3D<NekDouble> deriv_interp(m_spacedim);
558  TensorOfArray3D<NekDouble> out_interp(m_spacedim);
559  for (int i = 0; i < m_spacedim; ++i)
560  {
561  // Interpolate velocity
562  vel_interp[i] = Array<OneD, NekDouble> (nPts);
563  m_fields[0]->PhysInterp1DScaled(
564  OneDptscale, physfield[i], vel_interp[i]);
565 
566  // Interpolate derivatives
567  deriv_interp[i] = Array<OneD, Array<OneD, NekDouble> >
568  (m_spacedim+1);
569  for (int j = 0; j < m_spacedim+1; ++j)
570  {
571  deriv_interp[i][j] = Array<OneD, NekDouble> (nPts);
572  m_fields[0]->PhysInterp1DScaled(
573  OneDptscale, derivativesO1[i][j], deriv_interp[i][j]);
574  }
575 
576  // Output (start from j=1 since flux is zero for rho)
577  out_interp[i] = Array<OneD, Array<OneD, NekDouble> > (m_spacedim+2);
578  for (int j = 1; j < m_spacedim+2; ++j)
579  {
580  out_interp[i][j] = Array<OneD, NekDouble> (nPts);
581  }
582  }
583 
584  // Velocity divergence
585  for (int j = 0; j < m_spacedim; ++j)
586  {
587  Vmath::Vadd(nPts, divVel, 1, deriv_interp[j][j], 1,
588  divVel, 1);
589  }
590 
591  // Velocity divergence scaled by lambda * mu
592  Vmath::Smul(nPts, lambda, divVel, 1, divVel, 1);
593  Vmath::Vmul(nPts, m_mu, 1, divVel, 1, divVel, 1);
594 
595  // Viscous flux vector for the rho equation = 0 (no need to dealias)
596  for (int i = 0; i < m_spacedim; ++i)
597  {
598  Vmath::Zero(nPts_orig, viscousTensor[i][0], 1);
599  }
600 
601  // Viscous stress tensor (for the momentum equations)
602  for (int i = 0; i < m_spacedim; ++i)
603  {
604  for (int j = i; j < m_spacedim; ++j)
605  {
606  Vmath::Vadd(nPts, deriv_interp[i][j], 1,
607  deriv_interp[j][i], 1,
608  out_interp[i][j+1], 1);
609 
610  Vmath::Vmul(nPts, m_mu, 1,
611  out_interp[i][j+1], 1,
612  out_interp[i][j+1], 1);
613 
614  if (i == j)
615  {
616  // Add divergence term to diagonal
617  Vmath::Vadd(nPts, out_interp[i][j+1], 1,
618  divVel, 1,
619  out_interp[i][j+1], 1);
620  }
621  else
622  {
623  // Make symmetric
624  out_interp[j][i+1] = out_interp[i][j+1];
625  }
626  }
627  }
628 
629  // Terms for the energy equation
630  for (int i = 0; i < m_spacedim; ++i)
631  {
632  Vmath::Zero(nPts, out_interp[i][m_spacedim+1], 1);
633  // u_j * tau_ij
634  for (int j = 0; j < m_spacedim; ++j)
635  {
636  Vmath::Vvtvp(nPts, vel_interp[j], 1,
637  out_interp[i][j+1], 1,
638  out_interp[i][m_spacedim+1], 1,
639  out_interp[i][m_spacedim+1], 1);
640  }
641  // Add k*T_i
642  Vmath::Vvtvp(nPts, m_thermalConductivity, 1,
643  deriv_interp[i][m_spacedim], 1,
644  out_interp[i][m_spacedim+1], 1,
645  out_interp[i][m_spacedim+1], 1);
646  }
647 
648  // Project to original space
649  for (int i = 0; i < m_spacedim; ++i)
650  {
651  for (int j = 1; j < m_spacedim+2; ++j)
652  {
653  m_fields[0]->PhysGalerkinProjection1DScaled(
654  OneDptscale,
655  out_interp[i][j],
656  viscousTensor[i][j]);
657  }
658  }
659  }
660 
661 
662  /**
663  * @brief For very special treatment. For general boundaries it does nothing
664  * But for WallViscous and WallAdiabatic, special treatment is needed
665  * because they get the same Bwd value, special treatment is needed for
666  * boundary treatment of diffusion flux
667  * Note: This special treatment could be removed by seperating
668  * WallViscous and WallAdiabatic into two different classes.
669  *
670  */
671  void NavierStokesCFE::SpecialBndTreat(
672  Array<OneD, Array<OneD, NekDouble>> &consvar)
673  {
674  size_t nConvectiveFields = consvar.size();
675  int ndens = 0;
676  int nengy = nConvectiveFields - 1;
677 
678  Array<OneD, Array<OneD, NekDouble>> bndCons {nConvectiveFields};
679  Array<OneD, NekDouble> bndTotEngy;
680  Array<OneD, NekDouble> bndPressure;
681  Array<OneD, NekDouble> bndRho;
682  Array<OneD, NekDouble> bndIntEndy;
683  int nLengthArray = 0;
684 
685  int cnt = 0;
686  int nBndRegions = m_fields[nengy]->
687  GetBndCondExpansions().size();
688  for (int j = 0; j < nBndRegions; ++j)
689  {
690  if (m_fields[nengy]->GetBndConditions()[j]->
691  GetBoundaryConditionType() ==
692  SpatialDomains::ePeriodic)
693  {
694  continue;
695  }
696 
697  size_t nBndEdges = m_fields[nengy]->
698  GetBndCondExpansions()[j]->GetExpSize();
699  for (int e = 0; e < nBndEdges; ++e)
700  {
701  size_t nBndEdgePts = m_fields[nengy]->
702  GetBndCondExpansions()[j]->GetExp(e)->GetTotPoints();
703 
704  int id2 = m_fields[0]->GetTrace()->
705  GetPhys_Offset(m_fields[0]->GetTraceMap()->
706  GetBndCondIDToGlobalTraceID(cnt++));
707 
708  // Imposing Temperature Twall at the wall
709  if (boost::iequals(m_fields[nengy]->GetBndConditions()[j]->
710  GetUserDefined(), "WallViscous"))
711  {
712  if (nBndEdgePts != nLengthArray)
713  {
714  for (int i = 0; i < nConvectiveFields; ++i)
715  {
716  bndCons[i] = Array<OneD, NekDouble>
717  {nBndEdgePts, 0.0};
718  }
719  bndTotEngy = Array<OneD, NekDouble> {nBndEdgePts, 0.0};
720  bndPressure = Array<OneD, NekDouble> {nBndEdgePts, 0.0};
721  bndRho = Array<OneD, NekDouble> {nBndEdgePts, 0.0};
722  bndIntEndy = Array<OneD, NekDouble> {nBndEdgePts, 0.0};
723  nLengthArray = nBndEdgePts;
724  }
725  else
726  {
727  Vmath::Zero(nLengthArray, bndPressure, 1);
728  Vmath::Zero(nLengthArray, bndRho , 1);
729  Vmath::Zero(nLengthArray, bndIntEndy , 1);
730  }
731 
732  Array<OneD, NekDouble> tmp;
733 
734  for (int k = 0; k < nConvectiveFields; ++k)
735  {
736  Vmath::Vcopy(nBndEdgePts, tmp = consvar[k] + id2, 1,
737  bndCons[k], 1);
738  }
739 
740  m_varConv->GetPressure(bndCons, bndPressure);
741  Vmath::Fill(nLengthArray, m_Twall, bndTotEngy, 1);
742  m_varConv->GetRhoFromPT(bndPressure, bndTotEngy, bndRho);
743  m_varConv->GetEFromRhoP(bndRho, bndPressure, bndIntEndy);
744  m_varConv->GetDynamicEnergy(bndCons, bndTotEngy);
745 
746  Vmath::Vvtvp(nBndEdgePts, bndIntEndy, 1, bndCons[ndens], 1,
747  bndTotEngy, 1, bndTotEngy, 1);
748 
749  Vmath::Vcopy(nBndEdgePts,
750  bndTotEngy, 1,
751  tmp = consvar[nengy] + id2, 1);
752  }
753  }
754  }
755  }
756 
757  /**
758  * @brief Calculate and return the ArtificialViscosity for shock-capturing.
759  */
760  void NavierStokesCFE::GetArtificialViscosity(
761  const Array<OneD, Array<OneD, NekDouble>> &inarray,
762  Array<OneD, NekDouble> &muav)
763  {
764  m_artificialDiffusion->GetArtificialViscosity(inarray, muav);
765  }
766 
767  /**
768  * @brief Calculate and return the Symmetric flux in IP method.
769  */
770  void NavierStokesCFE::GetViscousSymmtrFluxConservVar(
771  const int nDim,
772  const Array<OneD, Array<OneD, NekDouble> > &inaverg,
773  const Array<OneD, Array<OneD, NekDouble > > &inarray,
774  TensorOfArray3D<NekDouble> &outarray,
775  Array< OneD, int > &nonZeroIndex,
776  const Array<OneD, Array<OneD, NekDouble> > &normal)
777  {
778  size_t nConvectiveFields = inarray.size();
779  size_t nPts = inaverg[nConvectiveFields - 1].size();
780  nonZeroIndex = Array<OneD, int>{nConvectiveFields - 1, 0};
781  for (int i = 0; i < nConvectiveFields - 1; ++i)
782  {
783  nonZeroIndex[i] = i + 1;
784  }
785 
786  std::vector<NekDouble> inAvgTmp(nConvectiveFields);
787  std::vector<NekDouble> inTmp(nConvectiveFields);
788  std::vector<NekDouble> outTmp(nConvectiveFields);
789  for (int d = 0; d < nDim; ++d)
790  {
791  for (int nderiv = 0; nderiv < nDim; ++nderiv)
792  {
793  for (size_t p = 0; p < nPts; ++p)
794  {
795  // rearrenge data
796  for (int f = 0; f < nConvectiveFields; ++f)
797  {
798  inAvgTmp[f] = inaverg[f][p];
799  inTmp[f] = inarray[f][p];
800  }
801 
802  // get temp
803  NekDouble temperature = m_varConv->GetTemperature(inTmp.data());
804  // get viscosity
805  NekDouble mu;
806  GetViscosityFromTempKernel(temperature, mu);
807 
808  GetViscousFluxBilinearFormKernel(nDim, d, nderiv,
809  inAvgTmp.data(), inTmp.data(), mu, outTmp.data());
810 
811  for (int f = 0; f < nConvectiveFields; ++f)
812  {
813  outarray[d][f][p] += normal[d][p] * outTmp[f];
814  }
815  }
816  }
817  }
818  }
819 
820  /**
821  * @brief Calculate the physical artificial viscosity
822  *
823  * @param physfield Input field.
824  */
825  void NavierStokesCFE::GetPhysicalAV(
826  const Array<OneD, const Array<OneD, NekDouble>> &physfield)
827  {
828  int nPts = physfield[0].size();
829  int nElements = m_fields[0]->GetExpSize();
830  Array <OneD, NekDouble > hOverP(nElements, 0.0);
831  hOverP = GetElmtMinHP();
832 
833  // Determine the maximum wavespeed
834  Array <OneD, NekDouble > Lambdas(nPts, 0.0);
835  Array <OneD, NekDouble > soundspeed(nPts, 0.0);
836  Array <OneD, NekDouble > absVelocity(nPts, 0.0);
837  m_varConv->GetSoundSpeed(physfield, soundspeed);
838  m_varConv->GetAbsoluteVelocity(physfield, absVelocity);
839 
840  Vmath::Vadd(nPts, absVelocity, 1, soundspeed, 1, Lambdas, 1);
841 
842  // Compute sensor based on rho
843  Array<OneD, NekDouble> Sensor(nPts, 0.0);
844  m_varConv->GetSensor(m_fields[0], physfield, Sensor, m_muav, 1);
845 
846  Array<OneD, NekDouble> tmp;
847  for (int e = 0; e < nElements; e++)
848  {
849  int physOffset = m_fields[0]->GetPhys_Offset(e);
850  int nElmtPoints = m_fields[0]->GetExp(e)->GetTotPoints();
851 
852  // Compute the maximum wave speed
853  NekDouble LambdaElmt = Vmath::Vmax(nElmtPoints, tmp = Lambdas
854  + physOffset, 1);
855 
856  // Compute average bounded density
857  NekDouble rhoAve = Vmath::Vsum(nElmtPoints, tmp = physfield[0]
858  + physOffset, 1);
859  rhoAve = rhoAve / nElmtPoints;
860  rhoAve = Smath::Smax(rhoAve , 1.0e-4, 1.0e+4);
861 
862  // Scale sensor by coeff, h/p, and density
863  LambdaElmt *= m_mu0 * hOverP[e] * rhoAve;
864  Vmath::Smul(nElmtPoints, LambdaElmt, tmp = m_muav + physOffset, 1,
865  tmp = m_muav + physOffset, 1);
866  }
867  }
868 
869  void NavierStokesCFE::CalcViscosity(
870  const Array<OneD, const Array<OneD, NekDouble>> &inaverg,
871  Array<OneD, NekDouble> &mu)
872  {
873  int nConvectiveFields = inaverg.size();
874  int nPts = inaverg[nConvectiveFields-1].size();
875 
876  if (m_ViscosityType == "Variable")
877  {
878  Array<OneD, NekDouble> tmp(nPts,0.0);
879  m_varConv->GetTemperature(inaverg,tmp);
880  m_varConv->GetDynamicViscosity(tmp, mu);
881  }
882  else
883  {
884  //mu may be on volume or trace
885  Vmath::Fill(nPts, m_mu[0], mu, 1);
886  }
887  }
888  /**
889  * @brief Make field C0.
890  *
891  * @param field Input Field
892  */
893  void NavierStokesCFE::C0Smooth( Array<OneD, NekDouble> &field )
894  {
895  int nCoeffs = m_C0ProjectExp->GetNcoeffs();
896  Array<OneD, NekDouble> muFwd(nCoeffs);
897  Array<OneD, NekDouble> weights(nCoeffs, 1.0);
898  // Assemble global expansion coefficients for viscosity
899  m_C0ProjectExp->FwdTrans_IterPerExp(field,
900  m_C0ProjectExp->UpdateCoeffs());
901  m_C0ProjectExp->Assemble();
902  Vmath::Vcopy(nCoeffs, m_C0ProjectExp->GetCoeffs(), 1, muFwd, 1);
903  // Global coefficients
904  Vmath::Vcopy(nCoeffs, weights, 1,
905  m_C0ProjectExp->UpdateCoeffs(), 1);
906  // This is the sign vector
907  m_C0ProjectExp->GlobalToLocal();
908  // Get weights
909  m_C0ProjectExp->Assemble();
910  // Divide
911  Vmath::Vdiv(nCoeffs, muFwd, 1, m_C0ProjectExp->GetCoeffs(), 1,
912  m_C0ProjectExp->UpdateCoeffs(), 1);
913  // Get local coefficients
914  m_C0ProjectExp->GlobalToLocal();
915  // Get C0 field
916  m_C0ProjectExp->BwdTrans_IterPerExp(
917  m_C0ProjectExp->GetCoeffs(), field);
918  }
919 
920  /**
921  * @brief Applied Ducros (anti-vorticity) sensor.
922  *
923  * @param field Input Field
924  */
925  void NavierStokesCFE::Ducros( Array<OneD, NekDouble> &field )
926  {
927  int nPts = m_fields[0]->GetTotPoints();
928  Array<OneD, NekDouble> denDuc(nPts, NekConstants::kNekZeroTol);
929  Array<OneD, NekDouble> ducros(nPts, 0.0);
930  Vmath::Vadd(nPts, denDuc, 1, m_diffusion->m_divVelSquare, 1,
931  denDuc, 1);
932  Vmath::Vadd(nPts, denDuc, 1, m_diffusion->m_curlVelSquare, 1,
933  denDuc, 1);
934  Vmath::Vdiv(nPts, m_diffusion->m_divVelSquare, 1, denDuc, 1,
935  ducros, 1);
936  // Average in cell
937  Array<OneD, NekDouble> tmp;
938  for (int e = 0; e < m_fields[0]->GetExpSize(); e++)
939  {
940  int nElmtPoints = m_fields[0]->GetExp(e)->GetTotPoints();
941  int physOffset = m_fields[0]->GetPhys_Offset(e);
942 
943  NekDouble eAve = Vmath::Vsum(nElmtPoints,
944  tmp = ducros + physOffset, 1);
945  eAve = eAve / nElmtPoints;
946  Vmath::Fill(nElmtPoints, eAve, tmp = ducros + physOffset, 1);
947  }
948  Vmath::Vmul(nPts, ducros, 1, field, 1, field, 1);
949  }
950 
951 
952  /**
953  * @brief Return the penalty vector for the LDGNS diffusion problem.
954  */
955  void NavierStokesCFE::v_GetFluxPenalty(
956  const Array<OneD, const Array<OneD, NekDouble>> &uFwd,
957  const Array<OneD, const Array<OneD, NekDouble>> &uBwd,
958  Array<OneD, Array<OneD, NekDouble>> &penaltyCoeff)
959  {
960  unsigned int nTracePts = uFwd[0].size();
961 
962  // Compute average temperature
963  unsigned int nVariables = uFwd.size();
964  Array<OneD, NekDouble> tAve{nTracePts, 0.0};
965  Vmath::Svtsvtp(nTracePts, 0.5, uFwd[nVariables-1], 1,
966  0.5, uBwd[nVariables-1], 1, tAve, 1);
967 
968  // Get average viscosity and thermal conductivity
969  Array<OneD, NekDouble> muAve{nTracePts, 0.0};
970  Array<OneD, NekDouble> tcAve{nTracePts, 0.0};
971 
972  GetViscosityAndThermalCondFromTemp(tAve, muAve, tcAve);
973 
974  // Compute penalty term
975  for (int i = 0; i < nVariables; ++i)
976  {
977  // Get jump of u variables
978  Vmath::Vsub(nTracePts, uFwd[i], 1, uBwd[i], 1, penaltyCoeff[i], 1);
979  // Multiply by variable coefficient = {coeff} ( u^+ - u^- )
980  if ( i < nVariables-1 )
981  {
982  Vmath::Vmul(nTracePts, muAve, 1, penaltyCoeff[i], 1,
983  penaltyCoeff[i], 1);
984  }
985  else
986  {
987  Vmath::Vmul(nTracePts, tcAve, 1, penaltyCoeff[i], 1,
988  penaltyCoeff[i], 1);
989  }
990  }
991  }
992 
993 
994  /**
995  * @brief Update viscosity
996  * todo: add artificial viscosity here
997  */
998  void NavierStokesCFE::GetViscosityAndThermalCondFromTemp(
999  const Array<OneD, NekDouble> &temperature,
1000  Array<OneD, NekDouble> &mu,
1001  Array<OneD, NekDouble> &thermalCond)
1002  {
1003  int nPts = temperature.size();
1004 
1005  for (size_t p = 0; p < nPts; ++p)
1006  {
1007  GetViscosityAndThermalCondFromTempKernel(temperature[p], mu[p],
1008  thermalCond[p]);
1009  }
1010  }
1011 
1012 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
Nektar::CompressibleFlowSystem::m_diffusion
SolverUtils::DiffusionSharedPtr m_diffusion
Definition: CompressibleFlowSystem.h:99
Nektar::CompressibleFlowSystem::GetElmtMinHP
Array< OneD, NekDouble > GetElmtMinHP(void)
Function to get estimate of min h/p factor per element.
Definition: CompressibleFlowSystem.cpp:1039
Nektar::CompressibleFlowSystem::m_varConv
VariableConverterSharedPtr m_varConv
Definition: CompressibleFlowSystem.h:121
Nektar::CompressibleFlowSystem::SetBoundaryConditionsBwdWeight
void SetBoundaryConditionsBwdWeight()
Set up a weight on physical boundaries for boundary condition applications.
Definition: CompressibleFlowSystem.cpp:426
Nektar::CompressibleFlowSystem::v_InitObject
virtual void v_InitObject()
Initialization object for CompressibleFlowSystem class.
Definition: CompressibleFlowSystem.cpp:56
Nektar::CompressibleFlowSystem::m_artificialDiffusion
ArtificialDiffusionSharedPtr m_artificialDiffusion
Definition: CompressibleFlowSystem.h:100
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:200
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:145
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:171
Nektar::NavierStokesCFE::CalcViscosity
void CalcViscosity(const Array< OneD, const Array< OneD, NekDouble >> &inaverg, Array< OneD, NekDouble > &mu)
Definition: NavierStokesCFE.cpp:869
Nektar::NavierStokesCFE::v_GetViscousFluxVector
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.
Definition: NavierStokesCFE.cpp:404
Nektar::NavierStokesCFE::SpecialBndTreat
void SpecialBndTreat(Array< OneD, Array< OneD, NekDouble > > &consvar)
For very special treatment. For general boundaries it does nothing But for WallViscous and WallAdiaba...
Definition: NavierStokesCFE.cpp:671
Nektar::NavierStokesCFE::GetViscousSymmtrFluxConservVar
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.
Definition: NavierStokesCFE.cpp:770
Nektar::NavierStokesCFE::m_thermalConductivity
Array< OneD, NekDouble > m_thermalConductivity
Definition: NavierStokesCFE.h:95
Nektar::NavierStokesCFE::GetPhysicalAV
void GetPhysicalAV(const Array< OneD, const Array< OneD, NekDouble >> &physfield)
Calculate the physical artificial viscosity.
Definition: NavierStokesCFE.cpp:825
Nektar::NavierStokesCFE::m_physicalSensorType
std::string m_physicalSensorType
Definition: NavierStokesCFE.h:83
Nektar::NavierStokesCFE::m_thermalConductivityRef
NekDouble m_thermalConductivityRef
Definition: NavierStokesCFE.h:93
Nektar::NavierStokesCFE::m_is_mu_variable
bool m_is_mu_variable
flag to switch between constant viscosity and Sutherland an enum could be added for more options
Definition: NavierStokesCFE.h:73
Nektar::NavierStokesCFE::v_ExtraFldOutput
virtual void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble >> &fieldcoeffs, std::vector< std::string > &variables)
Definition: NavierStokesCFE.cpp:179
Nektar::NavierStokesCFE::v_GetViscousFluxVectorDeAlias
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.
Definition: NavierStokesCFE.cpp:512
Nektar::NavierStokesCFE::m_C0ProjectExp
MultiRegions::ContFieldSharedPtr m_C0ProjectExp
Definition: NavierStokesCFE.h:85
Nektar::NavierStokesCFE::GetViscosityAndThermalCondFromTemp
void GetViscosityAndThermalCondFromTemp(const Array< OneD, NekDouble > &temperature, Array< OneD, NekDouble > &mu, Array< OneD, NekDouble > &thermalCond)
Update viscosity todo: add artificial viscosity here.
Definition: NavierStokesCFE.cpp:998
Nektar::NavierStokesCFE::C0Smooth
void C0Smooth(Array< OneD, NekDouble > &field)
Make field C0.
Definition: NavierStokesCFE.cpp:893
Nektar::NavierStokesCFE::v_DoDiffusion
virtual void v_DoDiffusion(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const Array< OneD, const Array< OneD, NekDouble >> &pFwd, const Array< OneD, const Array< OneD, NekDouble >> &pBwd)
Definition: NavierStokesCFE.cpp:301
Nektar::NavierStokesCFE::m_is_diffIP
bool m_is_diffIP
flag to switch between IP and LDG an enum could be added for more options
Definition: NavierStokesCFE.h:76
Nektar::NavierStokesCFE::Ducros
void Ducros(Array< OneD, NekDouble > &field)
Applied Ducros (anti-vorticity) sensor.
Definition: NavierStokesCFE.cpp:925
Nektar::NavierStokesCFE::m_mu
Array< OneD, NekDouble > m_mu
Definition: NavierStokesCFE.h:94
Nektar::NavierStokesCFE::GetViscousFluxBilinearFormKernel
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.
Definition: NavierStokesCFE.h:217
Nektar::NavierStokesCFE::GetViscosityFromTempKernel
void GetViscosityFromTempKernel(const T &temperature, T &mu)
Definition: NavierStokesCFE.h:189
Nektar::NavierStokesCFE::v_GetFluxPenalty
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.
Definition: NavierStokesCFE.cpp:955
Nektar::NavierStokesCFE::GetArtificialViscosity
void GetArtificialViscosity(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &muav)
Calculate and return the ArtificialViscosity for shock-capturing.
Definition: NavierStokesCFE.cpp:760
Nektar::NavierStokesCFE::GetViscosityAndThermalCondFromTempKernel
void GetViscosityAndThermalCondFromTempKernel(const T &temperature, T &mu, T &thermalCond)
Definition: NavierStokesCFE.h:175
Nektar::NavierStokesCFE::v_InitObject
virtual void v_InitObject()
Initialization object for CompressibleFlowSystem class.
Definition: NavierStokesCFE.cpp:62
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:376
Nektar::SolverUtils::EquationSystem::m_graph
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
Definition: EquationSystem.h:347
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:688
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:343
Nektar::SolverUtils::EquationSystem::m_specHP_dealiasing
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
Definition: EquationSystem.h:396
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:337
Nektar::SolverUtils::EquationSystem::GetNpoints
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:723
Nektar::SolverUtils::EquationSystem::GetPhys_Offset
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
Definition: EquationSystem.h:703
Nektar::SolverUtils::UnsteadySystem
Base class for unsteady solvers.
Definition: UnsteadySystem.h:49
Nektar::SolverUtils::UnsteadySystem::m_CalcPhysicalAV
bool m_CalcPhysicalAV
flag to update artificial viscosity
Definition: UnsteadySystem.h:144
Nektar::SolverUtils::UnsteadySystem::m_explicitDiffusion
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
Definition: UnsteadySystem.h:88
CG_Iterations.pressure
pressure
Definition: CG_Iterations.py:229
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:120
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:51
Nektar::SolverUtils::GetDiffusionFactory
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:41
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:87
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::NullNekDoubleArrayOfArray
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
Definition: SharedArray.hpp:821
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Smath::Smax
T Smax(const T a, const T b, const T k)
Return the soft max of between two scalars.
Definition: Smath.hpp:53
Vmath::Svtsvtp
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:691
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:192
Vmath::Vvtvp
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:513
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:846
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:322
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:257
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:436
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:892
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:372
Generated on Wed Nov 24 2021 10:33:27 for Nektar++ by doxygen 1.9.1