Nektar++
NavierStokesImplicitCFE.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File NavierStokesImplicitCFE.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/NavierStokesImplicitCFE.h>
36 
37 using namespace std;
38 
39 namespace Nektar
40 {
41 string NavierStokesImplicitCFE::className =
42  SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
43  "NavierStokesImplicitCFE", NavierStokesImplicitCFE::create,
44  "NavierStokes equations in conservative variables.");
45 
46 NavierStokesImplicitCFE::NavierStokesImplicitCFE(
47  const LibUtilities::SessionReaderSharedPtr &pSession,
48  const SpatialDomains::MeshGraphSharedPtr &pGraph)
49  : UnsteadySystem(pSession, pGraph),
50  CompressibleFlowSystem(pSession, pGraph),
51  NavierStokesCFE(pSession, pGraph), CFSImplicit(pSession, pGraph)
52 {
53 }
54 
55 NavierStokesImplicitCFE::~NavierStokesImplicitCFE()
56 {
57 }
58 
59 /**
60  * @brief Initialization object for CompressibleFlowSystem class.
61  */
62 void NavierStokesImplicitCFE::v_InitObject(bool DeclareFields)
63 {
64  CFSImplicit::v_InitObject(DeclareFields);
65 
66  NavierStokesCFE::InitObject_Explicit();
67 
68  m_GetdFlux_dDeriv_Array = Array<OneD, GetdFlux_dDeriv>(m_spacedim);
69  switch (m_spacedim)
70  {
71  case 2:
72  /* code */
73  m_GetdFlux_dDeriv_Array[0] =
74  std::bind(&NavierStokesImplicitCFE::GetdFlux_dQx_2D, this,
75  std::placeholders::_1, std::placeholders::_2,
76  std::placeholders::_3, std::placeholders::_4);
77 
78  m_GetdFlux_dDeriv_Array[1] =
79  std::bind(&NavierStokesImplicitCFE::GetdFlux_dQy_2D, this,
80  std::placeholders::_1, std::placeholders::_2,
81  std::placeholders::_3, std::placeholders::_4);
82  break;
83  case 3:
84  /* code */
85  m_GetdFlux_dDeriv_Array[0] =
86  std::bind(&NavierStokesImplicitCFE::GetdFlux_dQx_3D, this,
87  std::placeholders::_1, std::placeholders::_2,
88  std::placeholders::_3, std::placeholders::_4);
89 
90  m_GetdFlux_dDeriv_Array[1] =
91  std::bind(&NavierStokesImplicitCFE::GetdFlux_dQy_3D, this,
92  std::placeholders::_1, std::placeholders::_2,
93  std::placeholders::_3, std::placeholders::_4);
94  m_GetdFlux_dDeriv_Array[2] =
95  std::bind(&NavierStokesImplicitCFE::GetdFlux_dQz_3D, this,
96  std::placeholders::_1, std::placeholders::_2,
97  std::placeholders::_3, std::placeholders::_4);
98 
99  break;
100 
101  default:
102 
103  break;
104  }
105 }
106 
107 void NavierStokesImplicitCFE::v_DoDiffusionCoeff(
108  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
109  Array<OneD, Array<OneD, NekDouble>> &outarray,
110  const Array<OneD, const Array<OneD, NekDouble>> &pFwd,
111  const Array<OneD, const Array<OneD, NekDouble>> &pBwd)
112 {
113  size_t nvariables = inarray.size();
114  size_t npoints = GetNpoints();
115  size_t ncoeffs = GetNcoeffs();
116  size_t nTracePts = GetTraceTotPoints();
117 
118  Array<OneD, Array<OneD, NekDouble>> outarrayDiff{nvariables};
119  for (int i = 0; i < nvariables; ++i)
120  {
121  outarrayDiff[i] = Array<OneD, NekDouble>{ncoeffs, 0.0};
122  }
123 
124  if (m_is_diffIP)
125  {
126  if (m_bndEvaluateTime < 0.0)
127  {
128  NEKERROR(ErrorUtil::efatal, "m_bndEvaluateTime not setup");
129  }
130  m_diffusion->DiffuseCoeffs(nvariables, m_fields, inarray, outarrayDiff,
131  m_bndEvaluateTime, pFwd, pBwd);
132  for (int i = 0; i < nvariables; ++i)
133  {
134  Vmath::Vadd(ncoeffs, outarrayDiff[i], 1, outarray[i], 1,
135  outarray[i], 1);
136  }
137  }
138  else
139  {
140  Array<OneD, Array<OneD, NekDouble>> inarrayDiff{nvariables - 1};
141  Array<OneD, Array<OneD, NekDouble>> inFwd{nvariables - 1};
142  Array<OneD, Array<OneD, NekDouble>> inBwd{nvariables - 1};
143 
144  for (int i = 0; i < nvariables - 1; ++i)
145  {
146  inarrayDiff[i] = Array<OneD, NekDouble>{npoints};
147  inFwd[i] = Array<OneD, NekDouble>{nTracePts};
148  inBwd[i] = Array<OneD, NekDouble>{nTracePts};
149  }
150 
151  // Extract pressure
152  // (use inarrayDiff[0] as a temporary storage for the pressure)
153  m_varConv->GetPressure(inarray, inarrayDiff[0]);
154 
155  // Extract temperature
156  m_varConv->GetTemperature(inarray, inarrayDiff[nvariables - 2]);
157 
158  // Extract velocities
159  m_varConv->GetVelocityVector(inarray, inarrayDiff);
160 
161  // Repeat calculation for trace space
162  if (pFwd == NullNekDoubleArrayOfArray ||
163  pBwd == NullNekDoubleArrayOfArray)
164  {
165  inFwd = NullNekDoubleArrayOfArray;
166  inBwd = NullNekDoubleArrayOfArray;
167  }
168  else
169  {
170  m_varConv->GetPressure(pFwd, inFwd[0]);
171  m_varConv->GetPressure(pBwd, inBwd[0]);
172 
173  m_varConv->GetTemperature(pFwd, inFwd[nvariables - 2]);
174  m_varConv->GetTemperature(pBwd, inBwd[nvariables - 2]);
175 
176  m_varConv->GetVelocityVector(pFwd, inFwd);
177  m_varConv->GetVelocityVector(pBwd, inBwd);
178  }
179 
180  // Diffusion term in physical rhs form
181  m_diffusion->Diffuse(nvariables, m_fields, inarrayDiff, outarrayDiff,
182  inFwd, inBwd);
183 
184  for (int i = 0; i < nvariables; ++i)
185  {
186  Vmath::Vadd(npoints, outarrayDiff[i], 1, outarray[i], 1,
187  outarray[i], 1);
188  }
189 
190  // Laplacian operator based artificial viscosity
191  if (m_artificialDiffusion)
192  {
193  m_artificialDiffusion->DoArtificialDiffusionCoeff(inarray,
194  outarray);
195  }
196  }
197 }
198 
199 /**
200  * @brief return part of viscous Jacobian:
201  * \todo flux derived with Qx=[drho_dx,drhou_dx,drhov_dx,drhoE_dx]
202  * Input:
203  * normals:Point normals
204  * U=[rho,rhou,rhov,rhoE]
205  * Output: 2D 3*4 Matrix (flux with rho is zero)
206  */
207 void NavierStokesImplicitCFE::GetdFlux_dQx_2D(
208  const Array<OneD, NekDouble> &normals, const NekDouble &mu,
209  const Array<OneD, NekDouble> &U, DNekMatSharedPtr &OutputMatrix)
210 {
211  NekDouble nx = normals[0];
212  NekDouble ny = normals[1];
213  NekDouble rho = U[0];
214  NekDouble orho = 1.0 / rho;
215  NekDouble u = U[1] * orho;
216  NekDouble v = U[2] * orho;
217  NekDouble E = U[3] * orho;
218  NekDouble q2 = u * u + v * v;
219  NekDouble gamma = m_gamma;
220  // q_x=-kappa*dT_dx;
221  NekDouble Pr = m_Prandtl;
222  NekDouble oPr = 1.0 / Pr;
223  // To notice, here is positive, which is consistent with
224  //"SYMMETRIC INTERIOR PENALTY DG METHODS FOR THE COMPRESSIBLE
225  // NAVIER-STOKES EQUATIONS"
226  // But opposite to "I Do like CFD"
227  NekDouble tmp = mu * orho;
228  NekDouble tmp2 = gamma * oPr;
229  NekDouble OneThird, TwoThird, FourThird;
230  OneThird = 1.0 / 3.0;
231  TwoThird = 2.0 * OneThird;
232  FourThird = 4.0 * OneThird;
233 
234  Array<OneD, NekDouble> tmpArray;
235  tmpArray = OutputMatrix->GetPtr();
236  int nrow = OutputMatrix->GetRows();
237 
238  tmpArray[0 + 0 * nrow] = tmp * (-FourThird * u * nx - v * ny);
239  tmpArray[0 + 1 * nrow] = tmp * (FourThird * nx);
240  tmpArray[0 + 2 * nrow] = tmp * ny;
241  tmpArray[0 + 3 * nrow] = 0.0;
242  tmpArray[1 + 0 * nrow] = tmp * (-v * nx + TwoThird * u * ny);
243  tmpArray[1 + 1 * nrow] = tmp * (-TwoThird * ny);
244  tmpArray[1 + 2 * nrow] = tmp * nx;
245  tmpArray[1 + 3 * nrow] = 0.0;
246  tmpArray[2 + 0 * nrow] =
247  (FourThird * u * u + v * v + tmp2 * (E - q2)) * nx +
248  OneThird * u * v * ny;
249  tmpArray[2 + 0 * nrow] = -tmp * (*OutputMatrix)(2, 0);
250  tmpArray[2 + 1 * nrow] = (FourThird - tmp2) * u * nx - TwoThird * v * ny;
251  tmpArray[2 + 1 * nrow] = tmp * (*OutputMatrix)(2, 1);
252  tmpArray[2 + 2 * nrow] = (1 - tmp2) * v * nx + u * ny;
253  tmpArray[2 + 2 * nrow] = tmp * (*OutputMatrix)(2, 2);
254  tmpArray[2 + 3 * nrow] = tmp * tmp2 * nx;
255 }
256 
257 /**
258  * @brief return part of viscous Jacobian:
259  * \todo flux derived with Qx=[drho_dy,drhou_dy,drhov_dy,drhoE_dy]
260  * Input:
261  * normals:Point normals
262  * U=[rho,rhou,rhov,rhoE]
263  * Output: 2D 3*4 Matrix (flux with rho is zero)
264  */
265 void NavierStokesImplicitCFE::GetdFlux_dQy_2D(
266  const Array<OneD, NekDouble> &normals, const NekDouble &mu,
267  const Array<OneD, NekDouble> &U, DNekMatSharedPtr &OutputMatrix)
268 {
269  NekDouble nx = normals[0];
270  NekDouble ny = normals[1];
271  NekDouble rho = U[0];
272  NekDouble orho = 1.0 / rho;
273  NekDouble u = U[1] * orho;
274  NekDouble v = U[2] * orho;
275  NekDouble E = U[3] * orho;
276  NekDouble q2 = u * u + v * v;
277  NekDouble gamma = m_gamma;
278  // q_x=-kappa*dT_dx;
279  NekDouble Pr = m_Prandtl;
280  NekDouble oPr = 1.0 / Pr;
281  // To notice, here is positive, which is consistent with
282  //"SYMMETRIC INTERIOR PENALTY DG METHODS FOR THE COMPRESSIBLE
283  // NAVIER-STOKES EQUATIONS"
284  // But opposite to "I Do like CFD"
285  NekDouble tmp = mu * orho;
286  NekDouble tmp2 = gamma * oPr;
287  NekDouble OneThird, TwoThird, FourThird;
288  OneThird = 1.0 / 3.0;
289  TwoThird = 2.0 * OneThird;
290  FourThird = 4.0 * OneThird;
291 
292  Array<OneD, NekDouble> tmpArray;
293  tmpArray = OutputMatrix->GetPtr();
294  int nrow = OutputMatrix->GetRows();
295 
296  tmpArray[0 + 0 * nrow] = tmp * (TwoThird * v * nx - u * ny);
297  tmpArray[0 + 1 * nrow] = tmp * ny;
298  tmpArray[0 + 2 * nrow] = tmp * (-TwoThird) * nx;
299  tmpArray[0 + 3 * nrow] = 0.0;
300  tmpArray[1 + 0 * nrow] = tmp * (-u * nx - FourThird * v * ny);
301  tmpArray[1 + 1 * nrow] = tmp * nx;
302  tmpArray[1 + 2 * nrow] = tmp * (FourThird * ny);
303  tmpArray[1 + 3 * nrow] = 0.0;
304  tmpArray[2 + 0 * nrow] = OneThird * u * v * nx +
305  (FourThird * v * v + u * u + tmp2 * (E - q2)) * ny;
306  tmpArray[2 + 0 * nrow] = -tmp * (*OutputMatrix)(2, 0);
307  tmpArray[2 + 1 * nrow] = (1 - tmp2) * u * ny + v * nx;
308  tmpArray[2 + 1 * nrow] = tmp * (*OutputMatrix)(2, 1);
309  tmpArray[2 + 2 * nrow] = (FourThird - tmp2) * v * ny - TwoThird * u * nx;
310  tmpArray[2 + 2 * nrow] = tmp * (*OutputMatrix)(2, 2);
311  tmpArray[2 + 3 * nrow] = tmp * tmp2 * ny;
312 }
313 
314 /**
315  * @brief return part of viscous Jacobian derived with
316  * Qx=[drho_dx,drhou_dx,drhov_dx,drhow_dx,drhoE_dx]
317  * Input:
318  * normals:Point normals
319  * U=[rho,rhou,rhov,rhow,rhoE]
320  * dir: means whether derive with
321  * Qx=[drho_dx,drhou_dx,drhov_dx,drhow_dx,drhoE_dx]
322  * Output: 3D 4*5 Matrix (flux about rho is zero)
323  * OutputMatrix(dir=0)= dF_dQx;
324  */
325 void NavierStokesImplicitCFE::GetdFlux_dQx_3D(
326  const Array<OneD, NekDouble> &normals, const NekDouble &mu,
327  const Array<OneD, NekDouble> &U, DNekMatSharedPtr &OutputMatrix)
328 {
329  NekDouble nx = normals[0];
330  NekDouble ny = normals[1];
331  NekDouble nz = normals[2];
332  NekDouble rho = U[0];
333  NekDouble orho = 1.0 / rho;
334  NekDouble u = U[1] * orho;
335  NekDouble v = U[2] * orho;
336  NekDouble w = U[3] * orho;
337  NekDouble E = U[4] * orho;
338  NekDouble q2 = u * u + v * v + w * w;
339  NekDouble gamma = m_gamma;
340  // q_x=-kappa*dT_dx;
341  NekDouble Pr = m_Prandtl;
342  NekDouble oPr = 1.0 / Pr;
343  // To notice, here is positive, which is consistent with
344  //"SYMMETRIC INTERIOR PENALTY DG METHODS FOR THE COMPRESSIBLE
345  // NAVIER-STOKES EQUATIONS"
346  // But opposite to "I do like CFD"
347  NekDouble tmp = mu * orho;
348  NekDouble tmpx = tmp * nx;
349  NekDouble tmpy = tmp * ny;
350  NekDouble tmpz = tmp * nz;
351  NekDouble tmp2 = gamma * oPr;
352  NekDouble OneThird, TwoThird, FourThird;
353  OneThird = 1.0 / 3.0;
354  TwoThird = 2.0 * OneThird;
355  FourThird = 4.0 * OneThird;
356 
357  Array<OneD, NekDouble> tmpArray;
358  tmpArray = OutputMatrix->GetPtr();
359  int nrow = OutputMatrix->GetRows();
360 
361  tmpArray[0 + 0 * nrow] =
362  tmpx * (-FourThird * u) + tmpy * (-v) + tmpz * (-w);
363  tmpArray[0 + 1 * nrow] = tmpx * FourThird;
364  tmpArray[0 + 2 * nrow] = tmpy;
365  tmpArray[0 + 3 * nrow] = tmpz;
366  tmpArray[0 + 4 * nrow] = 0.0;
367  tmpArray[1 + 0 * nrow] = tmpx * (-v) + tmpy * (TwoThird * u);
368  tmpArray[1 + 1 * nrow] = tmpy * (-TwoThird);
369  tmpArray[1 + 2 * nrow] = tmpx;
370  tmpArray[1 + 3 * nrow] = 0.0;
371  tmpArray[1 + 4 * nrow] = 0.0;
372  tmpArray[2 + 0 * nrow] = tmpx * (-w) + tmpz * (TwoThird * u);
373  tmpArray[2 + 1 * nrow] = tmpz * (-TwoThird);
374  tmpArray[2 + 2 * nrow] = 0.0;
375  tmpArray[2 + 3 * nrow] = tmpx;
376  tmpArray[2 + 4 * nrow] = 0.0;
377  tmpArray[3 + 0 * nrow] =
378  -tmpx * (FourThird * u * u + v * v + w * w + tmp2 * (E - q2)) +
379  tmpy * (-OneThird * u * v) + tmpz * (-OneThird * u * w);
380  tmpArray[3 + 1 * nrow] = tmpx * (FourThird - tmp2) * u +
381  tmpy * (-TwoThird * v) + tmpz * (-TwoThird * w);
382  tmpArray[3 + 2 * nrow] = tmpx * (1.0 - tmp2) * v + tmpy * u;
383  tmpArray[3 + 3 * nrow] = tmpx * (1.0 - tmp2) * w + tmpz * u;
384  tmpArray[3 + 4 * nrow] = tmpx * tmp2;
385 }
386 
387 /**
388  * @brief return part of viscous Jacobian derived with
389  * Qy=[drho_dy,drhou_dy,drhov_dy,drhow_dy,drhoE_dy]
390  * Input:
391  * normals:Point normals
392  * U=[rho,rhou,rhov,rhow,rhoE]
393  * dir: means whether derive with
394  * Qy=[drho_dy,drhou_dy,drhov_dy,drhow_dy,drhoE_dy]
395  * Output: 3D 4*5 Matrix (flux about rho is zero)
396  * OutputMatrix(dir=1)= dF_dQy;
397  */
398 void NavierStokesImplicitCFE::GetdFlux_dQy_3D(
399  const Array<OneD, NekDouble> &normals, const NekDouble &mu,
400  const Array<OneD, NekDouble> &U, DNekMatSharedPtr &OutputMatrix)
401 {
402  NekDouble nx = normals[0];
403  NekDouble ny = normals[1];
404  NekDouble nz = normals[2];
405  NekDouble rho = U[0];
406  NekDouble orho = 1.0 / rho;
407  NekDouble u = U[1] * orho;
408  NekDouble v = U[2] * orho;
409  NekDouble w = U[3] * orho;
410  NekDouble E = U[4] * orho;
411  NekDouble q2 = u * u + v * v + w * w;
412  NekDouble gamma = m_gamma;
413  // q_x=-kappa*dT_dx;
414  NekDouble Pr = m_Prandtl;
415  NekDouble oPr = 1.0 / Pr;
416  // To notice, here is positive, which is consistent with
417  //"SYMMETRIC INTERIOR PENALTY DG METHODS FOR THE COMPRESSIBLE
418  // NAVIER-STOKES EQUATIONS"
419  // But opposite to "I do like CFD"
420  NekDouble tmp = mu * orho;
421  NekDouble tmpx = tmp * nx;
422  NekDouble tmpy = tmp * ny;
423  NekDouble tmpz = tmp * nz;
424  NekDouble tmp2 = gamma * oPr;
425  NekDouble OneThird, TwoThird, FourThird;
426  OneThird = 1.0 / 3.0;
427  TwoThird = 2.0 * OneThird;
428  FourThird = 4.0 * OneThird;
429 
430  Array<OneD, NekDouble> tmpArray;
431  tmpArray = OutputMatrix->GetPtr();
432  int nrow = OutputMatrix->GetRows();
433 
434  tmpArray[0 + 0 * nrow] = tmpx * (TwoThird * v) + tmpy * (-u);
435  tmpArray[0 + 1 * nrow] = tmpy;
436  tmpArray[0 + 2 * nrow] = tmpx * (-TwoThird);
437  tmpArray[0 + 3 * nrow] = 0.0;
438  tmpArray[0 + 4 * nrow] = 0.0;
439  tmpArray[1 + 0 * nrow] =
440  tmpx * (-u) + tmpy * (-FourThird * v) + tmpz * (-w);
441  tmpArray[1 + 1 * nrow] = tmpx;
442  tmpArray[1 + 2 * nrow] = tmpy * FourThird;
443  tmpArray[1 + 3 * nrow] = tmpz;
444  tmpArray[1 + 4 * nrow] = 0.0;
445  tmpArray[2 + 0 * nrow] = tmpy * (-w) + tmpz * (TwoThird * v);
446  tmpArray[2 + 1 * nrow] = 0.0;
447  tmpArray[2 + 2 * nrow] = tmpz * (-TwoThird);
448  tmpArray[2 + 3 * nrow] = tmpy;
449  tmpArray[2 + 4 * nrow] = 0.0;
450  tmpArray[3 + 0 * nrow] =
451  tmpx * (-OneThird * u * v) -
452  tmpy * (u * u + FourThird * v * v + w * w + tmp2 * (E - q2)) +
453  tmpz * (-OneThird * v * w);
454  tmpArray[3 + 1 * nrow] = tmpx * v + tmpy * (1 - tmp2) * u;
455  tmpArray[3 + 2 * nrow] = tmpx * (-TwoThird * u) +
456  tmpy * (FourThird - tmp2) * v +
457  tmpz * (-TwoThird * w);
458  tmpArray[3 + 3 * nrow] = tmpy * (1 - tmp2) * w + tmpz * v;
459  tmpArray[3 + 4 * nrow] = tmpy * tmp2;
460 }
461 
462 /**
463  * @brief return part of viscous Jacobian derived with
464  * Qz=[drho_dz,drhou_dz,drhov_dz,drhow_dz,drhoE_dz]
465  * Input:
466  * normals:Point normals
467  * U=[rho,rhou,rhov,rhow,rhoE]
468  * dir: means whether derive with
469  * Qz=[drho_dz,drhou_dz,drhov_dz,drhow_dz,drhoE_dz]
470  * Output: 3D 4*5 Matrix (flux about rho is zero)
471  * OutputMatrix(dir=2)= dF_dQz;
472  */
473 void NavierStokesImplicitCFE::GetdFlux_dQz_3D(
474  const Array<OneD, NekDouble> &normals, const NekDouble &mu,
475  const Array<OneD, NekDouble> &U, DNekMatSharedPtr &OutputMatrix)
476 {
477  NekDouble nx = normals[0];
478  NekDouble ny = normals[1];
479  NekDouble nz = normals[2];
480  NekDouble rho = U[0];
481  NekDouble orho = 1.0 / rho;
482  NekDouble u = U[1] * orho;
483  NekDouble v = U[2] * orho;
484  NekDouble w = U[3] * orho;
485  NekDouble E = U[4] * orho;
486  NekDouble q2 = u * u + v * v + w * w;
487  NekDouble gamma = m_gamma;
488  // q_x=-kappa*dT_dx;
489  NekDouble Pr = m_Prandtl;
490  NekDouble oPr = 1.0 / Pr;
491  // To notice, here is positive, which is consistent with
492  //"SYMMETRIC INTERIOR PENALTY DG METHODS FOR THE COMPRESSIBLE
493  // NAVIER-STOKES EQUATIONS"
494  // But opposite to "I do like CFD"
495  NekDouble tmp = mu * orho;
496  NekDouble tmpx = tmp * nx;
497  NekDouble tmpy = tmp * ny;
498  NekDouble tmpz = tmp * nz;
499  NekDouble tmp2 = gamma * oPr;
500  NekDouble OneThird, TwoThird, FourThird;
501  OneThird = 1.0 / 3.0;
502  TwoThird = 2.0 * OneThird;
503  FourThird = 4.0 * OneThird;
504 
505  Array<OneD, NekDouble> tmpArray;
506  tmpArray = OutputMatrix->GetPtr();
507  int nrow = OutputMatrix->GetRows();
508 
509  tmpArray[0 + 0 * nrow] = tmpx * (TwoThird * w) + tmpz * (-u);
510  tmpArray[0 + 1 * nrow] = tmpz;
511  tmpArray[0 + 2 * nrow] = 0.0;
512  tmpArray[0 + 3 * nrow] = tmpx * (-TwoThird);
513  tmpArray[0 + 4 * nrow] = 0.0;
514  tmpArray[1 + 0 * nrow] = tmpy * (TwoThird * w) + tmpz * (-v);
515  tmpArray[1 + 1 * nrow] = 0.0;
516  tmpArray[1 + 2 * nrow] = tmpz;
517  tmpArray[1 + 3 * nrow] = tmpy * (-TwoThird);
518  tmpArray[1 + 4 * nrow] = 0.0;
519  tmpArray[2 + 0 * nrow] =
520  tmpx * (-u) + tmpy * (-v) + tmpz * (-FourThird * w);
521  tmpArray[2 + 1 * nrow] = tmpx;
522  tmpArray[2 + 2 * nrow] = tmpy;
523  tmpArray[2 + 3 * nrow] = tmpz * FourThird;
524  tmpArray[2 + 4 * nrow] = 0.0;
525  tmpArray[3 + 0 * nrow] =
526  tmpx * (-OneThird * u * w) + tmpy * (-OneThird * v * w) -
527  tmpz * (u * u + v * v + FourThird * w * w + tmp2 * (E - q2));
528  tmpArray[3 + 1 * nrow] = tmpx * w + tmpz * (1 - tmp2) * u;
529  tmpArray[3 + 2 * nrow] = tmpy * w + tmpz * (1 - tmp2) * v;
530  tmpArray[3 + 3 * nrow] = tmpx * (-TwoThird * u) + tmpy * (-TwoThird * v) +
531  tmpz * (FourThird - tmp2) * w;
532  tmpArray[3 + 4 * nrow] = tmpz * tmp2;
533 }
534 
535 /**
536  * @brief return part of viscous Jacobian
537  * Input:
538  * normals:Point normals
539  * mu: dynamicviscosity
540  * dmu_dT: mu's derivative with T using Sutherland's law
541  * U=[rho,rhou,rhov,rhoE]
542  * Output: 3*4 Matrix (the flux about rho is zero)
543  * OutputMatrix dFLux_dU, the matrix sign is consistent with SIPG
544  */
545 void NavierStokesImplicitCFE::GetdFlux_dU_2D(
546  const Array<OneD, NekDouble> &normals, const NekDouble mu,
547  const NekDouble dmu_dT, const Array<OneD, NekDouble> &U,
548  const Array<OneD, const Array<OneD, NekDouble>> &qfield,
549  DNekMatSharedPtr &OutputMatrix)
550 {
551  Array<OneD, NekDouble> tmpArray;
552  tmpArray = OutputMatrix->GetPtr();
553  int nrow = OutputMatrix->GetRows();
554 
555  NekDouble nx = normals[0];
556  NekDouble ny = normals[1];
557  NekDouble U1 = U[0];
558  NekDouble U2 = U[1];
559  NekDouble U3 = U[2];
560  NekDouble U4 = U[3];
561  NekDouble dU1_dx = qfield[0][0];
562  NekDouble dU2_dx = qfield[0][1];
563  NekDouble dU3_dx = qfield[0][2];
564  NekDouble dU4_dx = qfield[0][3];
565  NekDouble dU1_dy = qfield[1][0];
566  NekDouble dU2_dy = qfield[1][1];
567  NekDouble dU3_dy = qfield[1][2];
568  NekDouble dU4_dy = qfield[1][3];
569  NekDouble gamma = m_gamma;
570  NekDouble Cv = m_Cv;
571  NekDouble Pr = m_Prandtl;
572  NekDouble oPr = 1.0 / Pr;
573 
574  NekDouble orho1, orho2, orho3, orho4;
575  NekDouble oCv = 1. / Cv;
576  orho1 = 1.0 / U1;
577  orho2 = orho1 * orho1;
578  orho3 = orho1 * orho2;
579  orho4 = orho2 * orho2;
580 
581  // Assume Fn=mu*Sn
582  // Sn=Sx*nx+Sy*ny
583 
584  NekDouble TwoThrid = 2. / 3.;
585  NekDouble FourThird = 2.0 * TwoThrid;
586  NekDouble u = U2 * orho1;
587  NekDouble v = U3 * orho1;
588  NekDouble du_dx = orho1 * (dU2_dx - u * dU1_dx);
589  NekDouble dv_dx = orho1 * (dU3_dx - v * dU1_dx);
590  NekDouble du_dy = orho1 * (dU2_dy - u * dU1_dy);
591  NekDouble dv_dy = orho1 * (dU3_dy - v * dU1_dy);
592  NekDouble s12 = FourThird * du_dx - TwoThrid * dv_dy;
593  NekDouble s13 = du_dy + dv_dx;
594  NekDouble s22 = s13;
595  NekDouble s23 = FourThird * dv_dy - TwoThrid * du_dx;
596  NekDouble snx = s12 * nx + s22 * ny;
597  NekDouble sny = s13 * nx + s23 * ny;
598  NekDouble snv = snx * u + sny * v;
599  NekDouble qx = -gamma * mu * oPr *
600  (orho1 * dU4_dx - U[3] * orho2 * dU1_dx -
601  u * (orho1 * dU2_dx - U[1] * orho2 * dU1_dx) -
602  v * (orho1 * dU3_dx - U[2] * orho2 * dU1_dx));
603  NekDouble qy = -gamma * mu * oPr *
604  (orho1 * dU4_dy - U[3] * orho2 * dU1_dy -
605  u * (orho1 * dU2_dy - U[1] * orho2 * dU1_dy) -
606  v * (orho1 * dU3_dy - U[2] * orho2 * dU1_dy));
607  NekDouble qn = qx * nx + qy * ny;
608 
609  // Term1 mu's derivative with U: dmu_dU*Sn
610  Array<OneD, NekDouble> tmp(3, 0.0);
611  tmp[0] = snx;
612  tmp[1] = sny;
613  tmp[2] = snv - qn / mu;
614  Array<OneD, NekDouble> dT_dU(4, 0.0);
615  dT_dU[0] = oCv * (-orho2 * U4 + orho3 * U2 * U2 + orho3 * U3 * U3);
616  dT_dU[1] = -oCv * orho2 * U2;
617  dT_dU[2] = -oCv * orho2 * U3;
618  dT_dU[3] = oCv * orho1;
619  for (int i = 0; i < 3; i++)
620  {
621  for (int j = 0; j < 4; j++)
622  {
623  tmpArray[i + j * nrow] = dmu_dT * dT_dU[j] * tmp[i];
624  }
625  }
626 
627  // Term 2 +mu*dSn_dU
628  NekDouble du_dx_dU1, du_dx_dU2;
629  NekDouble du_dy_dU1, du_dy_dU2;
630  NekDouble dv_dx_dU1, dv_dx_dU3;
631  NekDouble dv_dy_dU1, dv_dy_dU3;
632  NekDouble ds12_dU1, ds12_dU2, ds12_dU3;
633  NekDouble ds13_dU1, ds13_dU2, ds13_dU3;
634  NekDouble ds22_dU1, ds22_dU2, ds22_dU3;
635  NekDouble ds23_dU1, ds23_dU2, ds23_dU3;
636  NekDouble dsnx_dU1, dsnx_dU2, dsnx_dU3;
637  NekDouble dsny_dU1, dsny_dU2, dsny_dU3;
638  NekDouble dsnv_dU1, dsnv_dU2, dsnv_dU3;
639 
640  du_dx_dU1 = -orho2 * dU2_dx + 2 * orho3 * U2 * dU1_dx;
641  du_dx_dU2 = -orho2 * dU1_dx;
642  du_dy_dU1 = -orho2 * dU2_dy + 2 * orho3 * U2 * dU1_dy;
643  du_dy_dU2 = -orho2 * dU1_dy;
644  dv_dx_dU1 = -orho2 * dU3_dx + 2 * orho3 * U3 * dU1_dx;
645  dv_dx_dU3 = du_dx_dU2;
646  dv_dy_dU1 = -orho2 * dU3_dy + 2 * orho3 * U3 * dU1_dy;
647  dv_dy_dU3 = du_dy_dU2;
648  ds12_dU1 = FourThird * du_dx_dU1 - TwoThrid * dv_dy_dU1;
649  ds12_dU2 = FourThird * du_dx_dU2;
650  ds12_dU3 = -TwoThrid * dv_dy_dU3;
651  ds13_dU1 = du_dy_dU1 + dv_dx_dU1;
652  ds13_dU2 = du_dy_dU2;
653  ds13_dU3 = dv_dx_dU3;
654  ds22_dU1 = ds13_dU1;
655  ds22_dU2 = ds13_dU2;
656  ds22_dU3 = ds13_dU3;
657  ds23_dU1 = FourThird * dv_dy_dU1 - TwoThrid * du_dx_dU1;
658  ds23_dU2 = -TwoThrid * du_dx_dU2;
659  ds23_dU3 = FourThird * dv_dy_dU3;
660  dsnx_dU1 = ds12_dU1 * nx + ds22_dU1 * ny;
661  dsnx_dU2 = ds12_dU2 * nx + ds22_dU2 * ny;
662  dsnx_dU3 = ds12_dU3 * nx + ds22_dU3 * ny;
663  dsny_dU1 = ds13_dU1 * nx + ds23_dU1 * ny;
664  dsny_dU2 = ds13_dU2 * nx + ds23_dU2 * ny;
665  dsny_dU3 = ds13_dU3 * nx + ds23_dU3 * ny;
666  dsnv_dU1 =
667  u * dsnx_dU1 + v * dsny_dU1 - orho2 * U2 * snx - orho2 * U3 * sny;
668  dsnv_dU2 = u * dsnx_dU2 + v * dsny_dU2 + orho1 * snx;
669  dsnv_dU3 = u * dsnx_dU3 + v * dsny_dU3 + orho1 * sny;
670  tmpArray[0 + 0 * nrow] = tmpArray[0 + 0 * nrow] + mu * dsnx_dU1;
671  tmpArray[0 + 1 * nrow] = tmpArray[0 + 1 * nrow] + mu * dsnx_dU2;
672  tmpArray[0 + 2 * nrow] = tmpArray[0 + 2 * nrow] + mu * dsnx_dU3;
673  tmpArray[1 + 0 * nrow] = tmpArray[1 + 0 * nrow] + mu * dsny_dU1;
674  tmpArray[1 + 1 * nrow] = tmpArray[1 + 1 * nrow] + mu * dsny_dU2;
675  tmpArray[1 + 2 * nrow] = tmpArray[1 + 2 * nrow] + mu * dsny_dU3;
676  tmpArray[2 + 0 * nrow] = tmpArray[2 + 0 * nrow] + mu * dsnv_dU1;
677  tmpArray[2 + 1 * nrow] = tmpArray[2 + 1 * nrow] + mu * dsnv_dU2;
678  tmpArray[2 + 2 * nrow] = tmpArray[2 + 2 * nrow] + mu * dsnv_dU3;
679 
680  // Consider +qn's effect (does not include mu's effect)
681  NekDouble dqx_dU1, dqx_dU2, dqx_dU3, dqx_dU4;
682  NekDouble dqy_dU1, dqy_dU2, dqy_dU3, dqy_dU4;
683  NekDouble tmpx = -nx * mu * gamma * oPr;
684  dqx_dU1 = tmpx * (-orho2 * dU4_dx + 2 * orho3 * U4 * dU1_dx +
685  2 * orho3 * U2 * dU2_dx - 3 * orho4 * U2 * U2 * dU1_dx +
686  2 * orho3 * U3 * dU3_dx - 3 * orho4 * U3 * U3 * dU1_dx);
687  dqx_dU2 = tmpx * (-orho2 * dU2_dx + 2 * orho3 * U2 * dU1_dx);
688  dqx_dU3 = tmpx * (-orho2 * dU3_dx + 2 * orho3 * U3 * dU1_dx);
689  dqx_dU4 = -tmpx * orho2 * dU1_dx;
690  NekDouble tmpy = -ny * mu * gamma * oPr;
691  dqy_dU1 = tmpy * (-orho2 * dU4_dy + 2 * orho3 * U4 * dU1_dy +
692  2 * orho3 * U2 * dU2_dy - 3 * orho4 * U2 * U2 * dU1_dy +
693  2 * orho3 * U3 * dU3_dy - 3 * orho4 * U3 * U3 * dU1_dy);
694  dqy_dU2 = tmpy * (-orho2 * dU2_dy + 2 * orho3 * U2 * dU1_dy);
695  dqy_dU3 = tmpy * (-orho2 * dU3_dy + 2 * orho3 * U3 * dU1_dy);
696  dqy_dU4 = -tmpy * orho2 * dU1_dy;
697  tmpArray[2 + 0 * nrow] = tmpArray[2 + 0 * nrow] - dqx_dU1 - dqy_dU1;
698  tmpArray[2 + 1 * nrow] = tmpArray[2 + 1 * nrow] - dqx_dU2 - dqy_dU2;
699  tmpArray[2 + 2 * nrow] = tmpArray[2 + 2 * nrow] - dqx_dU3 - dqy_dU3;
700  tmpArray[2 + 3 * nrow] = tmpArray[2 + 3 * nrow] - dqx_dU4 - dqy_dU4;
701 }
702 
703 /**
704  * @brief return part of viscous Jacobian
705  * Input:
706  * normals:Point normals
707  * mu: dynamicviscosity
708  * dmu_dT: mu's derivative with T using Sutherland's law
709  * U=[rho,rhou,rhov,rhow,rhoE]
710  * Output: 4*5 Matrix (the flux about rho is zero)
711  * OutputMatrix dFLux_dU, the matrix sign is consistent with SIPG
712  */
713 void NavierStokesImplicitCFE::GetdFlux_dU_3D(
714  const Array<OneD, NekDouble> &normals, const NekDouble mu,
715  const NekDouble dmu_dT, const Array<OneD, NekDouble> &U,
716  const Array<OneD, const Array<OneD, NekDouble>> &qfield,
717  DNekMatSharedPtr &OutputMatrix)
718 {
719  Array<OneD, NekDouble> tmpArray;
720  tmpArray = OutputMatrix->GetPtr();
721  int nrow = OutputMatrix->GetRows();
722 
723  NekDouble nx = normals[0];
724  NekDouble ny = normals[1];
725  NekDouble nz = normals[2];
726  NekDouble U1 = U[0];
727  NekDouble U2 = U[1];
728  NekDouble U3 = U[2];
729  NekDouble U4 = U[3];
730  NekDouble U5 = U[4];
731  NekDouble dU1_dx = qfield[0][0];
732  NekDouble dU2_dx = qfield[0][1];
733  NekDouble dU3_dx = qfield[0][2];
734  NekDouble dU4_dx = qfield[0][3];
735  NekDouble dU5_dx = qfield[0][4];
736  NekDouble dU1_dy = qfield[1][0];
737  NekDouble dU2_dy = qfield[1][1];
738  NekDouble dU3_dy = qfield[1][2];
739  NekDouble dU4_dy = qfield[1][3];
740  NekDouble dU5_dy = qfield[1][4];
741  NekDouble dU1_dz = qfield[2][0];
742  NekDouble dU2_dz = qfield[2][1];
743  NekDouble dU3_dz = qfield[2][2];
744  NekDouble dU4_dz = qfield[2][3];
745  NekDouble dU5_dz = qfield[2][4];
746  NekDouble gamma = m_gamma;
747  NekDouble Cv = m_Cv;
748  NekDouble Pr = m_Prandtl;
749  NekDouble oPr = 1.0 / Pr;
750 
751  NekDouble orho1, orho2, orho3, orho4;
752  NekDouble oCv = 1. / Cv;
753  orho1 = 1.0 / U1;
754  orho2 = orho1 * orho1;
755  orho3 = orho1 * orho2;
756  orho4 = orho2 * orho2;
757 
758  // Assume Fn=mu*Sn
759  // Sn=Sx*nx+Sy*ny+Sz*nz
760  NekDouble TwoThrid = 2. / 3.;
761  NekDouble FourThird = 2.0 * TwoThrid;
762  NekDouble tmp2 = gamma * mu * oPr;
763  NekDouble u = U2 * orho1;
764  NekDouble v = U3 * orho1;
765  NekDouble w = U4 * orho1;
766  NekDouble du_dx = orho1 * (dU2_dx - u * dU1_dx);
767  NekDouble dv_dx = orho1 * (dU3_dx - v * dU1_dx);
768  NekDouble dw_dx = orho1 * (dU4_dx - w * dU1_dx);
769  NekDouble du_dy = orho1 * (dU2_dy - u * dU1_dy);
770  NekDouble dv_dy = orho1 * (dU3_dy - v * dU1_dy);
771  NekDouble dw_dy = orho1 * (dU4_dy - w * dU1_dy);
772  NekDouble du_dz = orho1 * (dU2_dz - u * dU1_dz);
773  NekDouble dv_dz = orho1 * (dU3_dz - v * dU1_dz);
774  NekDouble dw_dz = orho1 * (dU4_dz - w * dU1_dz);
775  NekDouble s12 = FourThird * du_dx - TwoThrid * dv_dy - TwoThrid * dw_dz;
776  NekDouble s13 = du_dy + dv_dx;
777  NekDouble s14 = dw_dx + du_dz;
778  NekDouble s22 = s13;
779  NekDouble s23 = FourThird * dv_dy - TwoThrid * du_dx - TwoThrid * dw_dz;
780  NekDouble s24 = dv_dz + dw_dy;
781  NekDouble s32 = s14;
782  NekDouble s33 = s24;
783  NekDouble s34 = FourThird * dw_dz - TwoThrid * du_dx - TwoThrid * dv_dy;
784  NekDouble snx = s12 * nx + s22 * ny + s32 * nz;
785  NekDouble sny = s13 * nx + s23 * ny + s33 * nz;
786  NekDouble snz = s14 * nz + s24 * ny + s34 * nz;
787  NekDouble snv = snx * u + sny * v + snz * w;
788  NekDouble qx = -tmp2 * (orho1 * dU5_dx - U5 * orho2 * dU1_dx -
789  u * (orho1 * dU2_dx - U2 * orho2 * dU1_dx) -
790  v * (orho1 * dU3_dx - U3 * orho2 * dU1_dx) -
791  w * (orho1 * dU4_dx - U4 * orho2 * dU1_dx));
792  NekDouble qy = -tmp2 * (orho1 * dU5_dy - U5 * orho2 * dU1_dy -
793  u * (orho1 * dU2_dy - U2 * orho2 * dU1_dy) -
794  v * (orho1 * dU3_dy - U3 * orho2 * dU1_dy) -
795  w * (orho1 * dU4_dy - U4 * orho2 * dU1_dy));
796  NekDouble qz = -tmp2 * (orho1 * dU5_dz - U5 * orho2 * dU1_dz -
797  u * (orho1 * dU2_dz - U2 * orho2 * dU1_dz) -
798  v * (orho1 * dU3_dz - U3 * orho2 * dU1_dz) -
799  w * (orho1 * dU4_dz - U4 * orho2 * dU1_dz));
800  NekDouble qn = qx * nx + qy * ny + qz * nz;
801 
802  // Term1 mu's derivative with U: dmu_dU*Sn
803  Array<OneD, NekDouble> tmp(4, 0.0);
804  tmp[0] = snx;
805  tmp[1] = sny;
806  tmp[2] = snz;
807  tmp[3] = snv - qn / mu;
808  Array<OneD, NekDouble> dT_dU(5, 0.0);
809  dT_dU[0] = oCv * (-orho2 * U5 + orho3 * U2 * U2 + orho3 * U3 * U3 +
810  orho3 * U4 * U4);
811  dT_dU[1] = -oCv * orho2 * U2;
812  dT_dU[2] = -oCv * orho2 * U3;
813  dT_dU[3] = -oCv * orho2 * U4;
814  dT_dU[4] = oCv * orho1;
815  for (int i = 0; i < 4; i++)
816  {
817  for (int j = 0; j < 5; j++)
818  {
819  tmpArray[i + j * nrow] = dmu_dT * dT_dU[j] * tmp[i];
820  }
821  }
822 
823  // Term 2 +mu*dSn_dU
824  NekDouble du_dx_dU1, du_dx_dU2;
825  NekDouble du_dy_dU1, du_dy_dU2;
826  NekDouble du_dz_dU1, du_dz_dU2;
827  NekDouble dv_dx_dU1, dv_dx_dU3;
828  NekDouble dv_dy_dU1, dv_dy_dU3;
829  NekDouble dv_dz_dU1, dv_dz_dU3;
830  NekDouble dw_dx_dU1, dw_dx_dU4;
831  NekDouble dw_dy_dU1, dw_dy_dU4;
832  NekDouble dw_dz_dU1, dw_dz_dU4;
833  NekDouble ds12_dU1, ds12_dU2, ds12_dU3, ds12_dU4;
834  NekDouble ds13_dU1, ds13_dU2, ds13_dU3;
835  NekDouble ds14_dU1, ds14_dU2, ds14_dU4;
836  NekDouble ds22_dU1, ds22_dU2, ds22_dU3;
837  NekDouble ds23_dU1, ds23_dU2, ds23_dU3, ds23_dU4;
838  NekDouble ds24_dU1, ds24_dU3, ds24_dU4;
839  NekDouble ds32_dU1, ds32_dU2, ds32_dU4;
840  NekDouble ds33_dU1, ds33_dU3, ds33_dU4;
841  NekDouble ds34_dU1, ds34_dU2, ds34_dU3, ds34_dU4;
842  NekDouble dsnx_dU1, dsnx_dU2, dsnx_dU3, dsnx_dU4;
843  NekDouble dsny_dU1, dsny_dU2, dsny_dU3, dsny_dU4;
844  NekDouble dsnz_dU1, dsnz_dU2, dsnz_dU3, dsnz_dU4;
845  NekDouble dsnv_dU1, dsnv_dU2, dsnv_dU3, dsnv_dU4;
846 
847  du_dx_dU1 = -orho2 * dU2_dx + 2 * orho3 * U2 * dU1_dx;
848  du_dx_dU2 = -orho2 * dU1_dx;
849  du_dy_dU1 = -orho2 * dU2_dy + 2 * orho3 * U2 * dU1_dy;
850  du_dy_dU2 = -orho2 * dU1_dy;
851  du_dz_dU1 = -orho2 * dU2_dz + 2 * orho3 * U2 * dU1_dz;
852  du_dz_dU2 = -orho2 * dU1_dz;
853  dv_dx_dU1 = -orho2 * dU3_dx + 2 * orho3 * U3 * dU1_dx;
854  dv_dx_dU3 = -orho2 * dU1_dx;
855  dv_dy_dU1 = -orho2 * dU3_dy + 2 * orho3 * U3 * dU1_dy;
856  dv_dy_dU3 = -orho2 * dU1_dy;
857  dv_dz_dU1 = -orho2 * dU3_dz + 2 * orho3 * U3 * dU1_dz;
858  dv_dz_dU3 = -orho2 * dU1_dz;
859  dw_dx_dU1 = -orho2 * dU4_dx + 2 * orho3 * U4 * dU1_dx;
860  dw_dx_dU4 = -orho2 * dU1_dx;
861  dw_dy_dU1 = -orho2 * dU4_dy + 2 * orho3 * U4 * dU1_dy;
862  dw_dy_dU4 = -orho2 * dU1_dy;
863  dw_dz_dU1 = -orho2 * dU4_dz + 2 * orho3 * U4 * dU1_dz;
864  dw_dz_dU4 = -orho2 * dU1_dz;
865  ds12_dU1 =
866  FourThird * du_dx_dU1 - TwoThrid * dv_dy_dU1 - TwoThrid * dw_dz_dU1;
867  ds12_dU2 = FourThird * du_dx_dU2;
868  ds12_dU3 = -TwoThrid * dv_dy_dU3;
869  ds12_dU4 = -TwoThrid * dw_dz_dU4;
870  ds13_dU1 = du_dy_dU1 + dv_dx_dU1;
871  ds13_dU2 = du_dy_dU2;
872  ds13_dU3 = dv_dx_dU3;
873  ds14_dU1 = dw_dx_dU1 + du_dz_dU1;
874  ds14_dU2 = du_dz_dU2;
875  ds14_dU4 = dw_dx_dU4;
876  ds22_dU1 = du_dy_dU1 + dv_dx_dU1;
877  ds22_dU2 = du_dy_dU2;
878  ds22_dU3 = dv_dx_dU3;
879  ds23_dU1 =
880  FourThird * dv_dy_dU1 - TwoThrid * du_dx_dU1 - TwoThrid * dw_dz_dU1;
881  ds23_dU2 = -TwoThrid * du_dx_dU2;
882  ds23_dU3 = FourThird * dv_dy_dU3;
883  ds23_dU4 = -TwoThrid * dw_dz_dU4;
884  ds24_dU1 = dv_dz_dU1 + dw_dy_dU1;
885  ds24_dU3 = dv_dz_dU3;
886  ds24_dU4 = dw_dy_dU4;
887  ds32_dU1 = dw_dx_dU1 + du_dz_dU1;
888  ds32_dU2 = du_dz_dU2;
889  ds32_dU4 = dw_dx_dU4;
890  ds33_dU1 = dv_dz_dU1 + dw_dy_dU1;
891  ds33_dU3 = dv_dz_dU3;
892  ds33_dU4 = dw_dy_dU4;
893  ds34_dU1 =
894  FourThird * dw_dz_dU1 - TwoThrid * du_dx_dU1 - TwoThrid * dv_dy_dU1;
895  ds34_dU2 = -TwoThrid * du_dx_dU2;
896  ds34_dU3 = -TwoThrid * dv_dy_dU3;
897  ds34_dU4 = FourThird * dw_dz_dU4;
898  dsnx_dU1 = ds12_dU1 * nx + ds22_dU1 * ny + ds32_dU1 * nz;
899  dsnx_dU2 = ds12_dU2 * nx + ds22_dU2 * ny + ds32_dU2 * nz;
900  dsnx_dU3 = ds12_dU3 * nx + ds22_dU3 * ny;
901  dsnx_dU4 = ds12_dU4 * nx + ds32_dU4 * nz;
902  dsny_dU1 = ds13_dU1 * nx + ds23_dU1 * ny + ds33_dU1 * nz;
903  dsny_dU2 = ds13_dU2 * nx + ds23_dU2 * ny;
904  dsny_dU3 = ds13_dU3 * nx + ds23_dU3 * ny + ds33_dU3 * nz;
905  dsny_dU4 = ds23_dU4 * ny + ds33_dU4 * nz;
906  dsnz_dU1 = ds14_dU1 * nx + ds24_dU1 * ny + ds34_dU1 * nz;
907  dsnz_dU2 = ds14_dU2 * nx + ds34_dU2 * nz;
908  dsnz_dU3 = ds24_dU3 * ny + ds34_dU3 * nz;
909  //? why there is value if 2D
910  dsnz_dU4 = ds14_dU4 * nx + ds24_dU4 * ny + ds34_dU4 * nz;
911  dsnv_dU1 = u * dsnx_dU1 + v * dsny_dU1 + w * dsnz_dU1 - orho2 * U2 * snx -
912  orho2 * U3 * sny - orho2 * U4 * snz;
913  dsnv_dU2 = u * dsnx_dU2 + v * dsny_dU2 + w * dsnz_dU2 + orho1 * snx;
914  dsnv_dU3 = u * dsnx_dU3 + v * dsny_dU3 + w * dsnz_dU3 + orho1 * sny;
915  dsnv_dU4 = u * dsnx_dU4 + v * dsny_dU4 + w * dsnz_dU4 + orho1 * snz;
916  tmpArray[0 + 0 * nrow] = tmpArray[0 + 0 * nrow] + mu * dsnx_dU1;
917  tmpArray[0 + 1 * nrow] = tmpArray[0 + 1 * nrow] + mu * dsnx_dU2;
918  tmpArray[0 + 2 * nrow] = tmpArray[0 + 2 * nrow] + mu * dsnx_dU3;
919  tmpArray[0 + 3 * nrow] = tmpArray[0 + 3 * nrow] + mu * dsnx_dU4;
920  tmpArray[1 + 0 * nrow] = tmpArray[1 + 0 * nrow] + mu * dsny_dU1;
921  tmpArray[1 + 1 * nrow] = tmpArray[1 + 1 * nrow] + mu * dsny_dU2;
922  tmpArray[1 + 2 * nrow] = tmpArray[1 + 2 * nrow] + mu * dsny_dU3;
923  tmpArray[1 + 3 * nrow] = tmpArray[1 + 3 * nrow] + mu * dsny_dU4;
924  tmpArray[2 + 0 * nrow] = tmpArray[2 + 0 * nrow] + mu * dsnz_dU1;
925  tmpArray[2 + 1 * nrow] = tmpArray[2 + 1 * nrow] + mu * dsnz_dU2;
926  tmpArray[2 + 2 * nrow] = tmpArray[2 + 2 * nrow] + mu * dsnz_dU3;
927  tmpArray[2 + 3 * nrow] = tmpArray[2 + 3 * nrow] + mu * dsnz_dU4;
928  tmpArray[3 + 0 * nrow] = tmpArray[3 + 0 * nrow] + mu * dsnv_dU1;
929  tmpArray[3 + 1 * nrow] = tmpArray[3 + 1 * nrow] + mu * dsnv_dU2;
930  tmpArray[3 + 2 * nrow] = tmpArray[3 + 2 * nrow] + mu * dsnv_dU3;
931  tmpArray[3 + 3 * nrow] = tmpArray[3 + 3 * nrow] + mu * dsnv_dU4;
932 
933  // Consider heat flux qn's effect (does not include mu's effect)
934  NekDouble dqx_dU1, dqx_dU2, dqx_dU3, dqx_dU4, dqx_dU5;
935  NekDouble dqy_dU1, dqy_dU2, dqy_dU3, dqy_dU4, dqy_dU5;
936  NekDouble dqz_dU1, dqz_dU2, dqz_dU3, dqz_dU4, dqz_dU5;
937  NekDouble tmpx = -nx * tmp2;
938  dqx_dU1 = tmpx * (-orho2 * dU5_dx + 2 * orho3 * U5 * dU1_dx +
939  2 * orho3 * U2 * dU2_dx - 3 * orho4 * U2 * U2 * dU1_dx +
940  2 * orho3 * U3 * dU3_dx - 3 * orho4 * U3 * U3 * dU1_dx +
941  2 * orho3 * U4 * dU4_dx - 3 * orho4 * U4 * U4 * dU1_dx);
942  dqx_dU2 = tmpx * (-orho2 * dU2_dx + 2 * orho3 * U2 * dU1_dx);
943  dqx_dU3 = tmpx * (-orho2 * dU3_dx + 2 * orho3 * U3 * dU1_dx);
944  dqx_dU4 = tmpx * (-orho2 * dU4_dx + 2 * orho3 * U4 * dU1_dx);
945  dqx_dU5 = -tmpx * orho2 * dU1_dx;
946  NekDouble tmpy = -ny * tmp2;
947  dqy_dU1 = tmpy * (-orho2 * dU5_dy + 2 * orho3 * U5 * dU1_dy +
948  2 * orho3 * U2 * dU2_dy - 3 * orho4 * U2 * U2 * dU1_dy +
949  2 * orho3 * U3 * dU3_dy - 3 * orho4 * U3 * U3 * dU1_dy +
950  2 * orho3 * U4 * dU4_dy - 3 * orho4 * U4 * U4 * dU1_dy);
951  dqy_dU2 = tmpy * (-orho2 * dU2_dy + 2 * orho3 * U2 * dU1_dy);
952  dqy_dU3 = tmpy * (-orho2 * dU3_dy + 2 * orho3 * U3 * dU1_dy);
953  dqy_dU4 = tmpy * (-orho2 * dU4_dy + 2 * orho3 * U4 * dU1_dy);
954  dqy_dU5 = -tmpy * orho2 * dU1_dy;
955  NekDouble tmpz = -nz * tmp2;
956  dqz_dU1 = tmpz * (-orho2 * dU5_dz + 2 * orho3 * U5 * dU1_dz +
957  2 * orho3 * U2 * dU2_dz - 3 * orho4 * U2 * U2 * dU1_dz +
958  2 * orho3 * U3 * dU3_dz - 3 * orho4 * U3 * U3 * dU1_dz +
959  2 * orho3 * U4 * dU4_dz - 3 * orho4 * U4 * U4 * dU1_dz);
960  dqz_dU2 = tmpz * (-orho2 * dU2_dz + 2 * orho3 * U2 * dU1_dz);
961  dqz_dU3 = tmpz * (-orho2 * dU3_dz + 2 * orho3 * U3 * dU1_dz);
962  dqz_dU4 = tmpz * (-orho2 * dU4_dz + 2 * orho3 * U4 * dU1_dz);
963  dqz_dU5 = -tmpz * orho2 * dU1_dz;
964  tmpArray[3 + 0 * nrow] =
965  tmpArray[3 + 0 * nrow] - dqx_dU1 - dqy_dU1 - dqz_dU1;
966  tmpArray[3 + 1 * nrow] =
967  tmpArray[3 + 1 * nrow] - dqx_dU2 - dqy_dU2 - dqz_dU2;
968  tmpArray[3 + 2 * nrow] =
969  tmpArray[3 + 2 * nrow] - dqx_dU3 - dqy_dU3 - dqz_dU3;
970  tmpArray[3 + 3 * nrow] =
971  tmpArray[3 + 3 * nrow] - dqx_dU4 - dqy_dU4 - dqz_dU4;
972  tmpArray[3 + 4 * nrow] =
973  tmpArray[3 + 4 * nrow] - dqx_dU5 - dqy_dU5 - dqz_dU5;
974 }
975 
976 void NavierStokesImplicitCFE::v_MinusDiffusionFluxJacPoint(
977  const int nConvectiveFields, const int nElmtPnt,
978  const Array<OneD, const Array<OneD, NekDouble>> &locVars,
979  const TensorOfArray3D<NekDouble> &locDerv,
980  const Array<OneD, NekDouble> &locmu, const Array<OneD, NekDouble> &locDmuDT,
981  const Array<OneD, NekDouble> &normals, DNekMatSharedPtr &wspMat,
982  Array<OneD, Array<OneD, NekDouble>> &PntJacArray)
983 {
984  int nSpaceDim = m_graph->GetSpaceDimension();
985 
986  NekDouble pointmu = 0.0;
987  NekDouble pointDmuDT = 0.0;
988  Array<OneD, NekDouble> pointVar(nConvectiveFields, 0.0);
989  Array<OneD, Array<OneD, NekDouble>> pointDerv(nSpaceDim);
990  for (int j = 0; j < nSpaceDim; j++)
991  {
992  pointDerv[j] = Array<OneD, NekDouble>(nConvectiveFields, 0.0);
993  }
994 
995  Array<OneD, NekDouble> wspMatData = wspMat->GetPtr();
996  Array<OneD, NekDouble> tmp1;
997  Array<OneD, NekDouble> tmp2;
998  Array<OneD, NekDouble> tmp3;
999 
1000  for (int npnt = 0; npnt < nElmtPnt; npnt++)
1001  {
1002  for (int j = 0; j < nConvectiveFields; j++)
1003  {
1004  pointVar[j] = locVars[j][npnt];
1005  }
1006  for (int j = 0; j < nSpaceDim; j++)
1007  {
1008  for (int k = 0; k < nConvectiveFields; k++)
1009  {
1010  pointDerv[j][k] = locDerv[j][k][npnt];
1011  }
1012  }
1013 
1014  pointmu = locmu[npnt];
1015  pointDmuDT = locDmuDT[npnt];
1016 
1017  v_GetDiffusionFluxJacPoint(pointVar, pointDerv, pointmu, pointDmuDT,
1018  normals, wspMat);
1019  for (int j = 0; j < nConvectiveFields; j++)
1020  {
1021  int noffset = j * nConvectiveFields;
1022 
1023  Vmath::Vsub(nConvectiveFields - 1,
1024  tmp1 = PntJacArray[npnt] + (noffset + 1), 1,
1025  tmp2 = wspMatData + (noffset - j), 1,
1026  tmp3 = PntJacArray[npnt] + (noffset + 1), 1);
1027  }
1028  }
1029 }
1030 
1031 void NavierStokesImplicitCFE::v_GetDiffusionFluxJacPoint(
1032  const Array<OneD, NekDouble> &conservVar,
1033  const Array<OneD, const Array<OneD, NekDouble>> &conseDeriv,
1034  const NekDouble mu, const NekDouble DmuDT,
1035  const Array<OneD, NekDouble> &normals, DNekMatSharedPtr &fluxJac)
1036 {
1037  switch (m_spacedim)
1038  {
1039  case 2:
1040  GetdFlux_dU_2D(normals, mu, DmuDT, conservVar, conseDeriv, fluxJac);
1041  break;
1042 
1043  case 3:
1044  GetdFlux_dU_3D(normals, mu, DmuDT, conservVar, conseDeriv, fluxJac);
1045  break;
1046 
1047  default:
1048  NEKERROR(ErrorUtil::efatal, "v_GetDiffusionFluxJacPoint not coded");
1049  break;
1050  }
1051 }
1052 
1053 void NavierStokesImplicitCFE::v_GetFluxDerivJacDirctn(
1054  const MultiRegions::ExpListSharedPtr &explist,
1055  const Array<OneD, const Array<OneD, NekDouble>> &normals,
1056  const int nDervDir,
1057  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1058  TensorOfArray5D<NekDouble> &ElmtJacArray, const int nFluxDir)
1059 {
1060  int nConvectiveFields = inarray.size();
1061  std::shared_ptr<LocalRegions::ExpansionVector> expvect = explist->GetExp();
1062  int nTotElmt = (*expvect).size();
1063  int nPts = explist->GetTotPoints();
1064  int nSpaceDim = m_graph->GetSpaceDimension();
1065 
1066  // Auxiliary variables
1067  Array<OneD, NekDouble> mu(nPts, 0.0);
1068  Array<OneD, NekDouble> thermalConductivity(nPts, 0.0);
1069  Array<OneD, NekDouble> temperature(nPts, 0.0);
1070  m_varConv->GetTemperature(inarray, temperature);
1071  GetViscosityAndThermalCondFromTemp(temperature, mu, thermalConductivity);
1072 
1073  NekDouble pointmu = 0.0;
1074  Array<OneD, NekDouble> locmu;
1075  Array<OneD, NekDouble> pointVar(nConvectiveFields, 0.0);
1076  Array<OneD, Array<OneD, NekDouble>> locVars(nConvectiveFields);
1077  Array<OneD, NekDouble> pointnormals(nSpaceDim, 0.0);
1078  Array<OneD, Array<OneD, NekDouble>> locnormal(nSpaceDim);
1079 
1080  DNekMatSharedPtr PointFJac = MemoryManager<DNekMat>::AllocateSharedPtr(
1081  nConvectiveFields - 1, nConvectiveFields);
1082  Array<OneD, NekDouble> PointFJac_data = PointFJac->GetPtr();
1083 
1084  for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
1085  {
1086  int nElmtPnt = (*expvect)[nelmt]->GetTotPoints();
1087  int noffest = explist->GetPhys_Offset(nelmt);
1088 
1089  for (int j = 0; j < nConvectiveFields; j++)
1090  {
1091  locVars[j] = inarray[j] + noffest;
1092  }
1093 
1094  for (int j = 0; j < nSpaceDim; j++)
1095  {
1096  locnormal[j] = normals[j] + noffest;
1097  }
1098 
1099  locmu = mu + noffest;
1100  for (int npnt = 0; npnt < nElmtPnt; npnt++)
1101  {
1102  for (int j = 0; j < nConvectiveFields; j++)
1103  {
1104  pointVar[j] = locVars[j][npnt];
1105  }
1106  for (int j = 0; j < nSpaceDim; j++)
1107  {
1108  pointnormals[j] = locnormal[j][npnt];
1109  }
1110 
1111  pointmu = locmu[npnt];
1112 
1113  m_GetdFlux_dDeriv_Array[nDervDir](pointnormals, pointmu, pointVar,
1114  PointFJac);
1115  for (int j = 0; j < nConvectiveFields; j++)
1116  {
1117  ElmtJacArray[0][j][nFluxDir][nelmt][npnt] = 0.0;
1118  }
1119  for (int j = 0; j < nConvectiveFields; j++)
1120  {
1121  int noffset = j * (nConvectiveFields - 1);
1122  for (int i = 0; i < nConvectiveFields - 1; i++)
1123  {
1124  ElmtJacArray[i + 1][j][nFluxDir][nelmt][npnt] =
1125  PointFJac_data[noffset + i];
1126  }
1127  }
1128  }
1129  }
1130 }
1131 
1132 void NavierStokesImplicitCFE::v_GetFluxDerivJacDirctnElmt(
1133  const int nConvectiveFields, const int nElmtPnt, const int nDervDir,
1134  const Array<OneD, const Array<OneD, NekDouble>> &locVars,
1135  const Array<OneD, NekDouble> &locmu,
1136  const Array<OneD, const Array<OneD, NekDouble>> &locnormal,
1137  DNekMatSharedPtr &wspMat, Array<OneD, Array<OneD, NekDouble>> &PntJacArray)
1138 {
1139  int nSpaceDim = m_graph->GetSpaceDimension();
1140 
1141  NekDouble pointmu = 0.0;
1142  Array<OneD, NekDouble> pointVar(nConvectiveFields, 0.0);
1143  Array<OneD, NekDouble> pointnormals(nSpaceDim, 0.0);
1144 
1145  Array<OneD, NekDouble> wspMatData = wspMat->GetPtr();
1146 
1147  Array<OneD, NekDouble> tmp1;
1148  Array<OneD, NekDouble> tmp2;
1149 
1150  for (int npnt = 0; npnt < nElmtPnt; npnt++)
1151  {
1152  for (int j = 0; j < nConvectiveFields; j++)
1153  {
1154  pointVar[j] = locVars[j][npnt];
1155  }
1156  for (int j = 0; j < nSpaceDim; j++)
1157  {
1158  pointnormals[j] = locnormal[j][npnt];
1159  }
1160 
1161  pointmu = locmu[npnt];
1162 
1163  m_GetdFlux_dDeriv_Array[nDervDir](pointnormals, pointmu, pointVar,
1164  wspMat);
1165  Vmath::Zero(nConvectiveFields, PntJacArray[npnt], nConvectiveFields);
1166  for (int j = 0; j < nConvectiveFields; j++)
1167  {
1168  int noffset = j * (nConvectiveFields - 1);
1169  Vmath::Vcopy((nConvectiveFields - 1), tmp1 = wspMatData + noffset,
1170  1, tmp2 = PntJacArray[npnt] + noffset + j + 1, 1);
1171  }
1172  }
1173 }
1174 
1175 void NavierStokesImplicitCFE::v_GetFluxDerivJacDirctn(
1176  const MultiRegions::ExpListSharedPtr &explist,
1177  const Array<OneD, const Array<OneD, NekDouble>> &normals,
1178  const int nDervDir,
1179  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1180  Array<OneD, Array<OneD, DNekMatSharedPtr>> &ElmtJac)
1181 {
1182  int nConvectiveFields = inarray.size();
1183  std::shared_ptr<LocalRegions::ExpansionVector> expvect = explist->GetExp();
1184  int nTotElmt = (*expvect).size();
1185  int nPts = explist->GetTotPoints();
1186  int nSpaceDim = m_graph->GetSpaceDimension();
1187 
1188  // Debug
1189  if (!ElmtJac.size())
1190  {
1191  ElmtJac = Array<OneD, Array<OneD, DNekMatSharedPtr>>(nTotElmt);
1192  for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
1193  {
1194  int nElmtPnt = (*expvect)[nelmt]->GetTotPoints();
1195  ElmtJac[nelmt] = Array<OneD, DNekMatSharedPtr>(nElmtPnt);
1196  for (int npnt = 0; npnt < nElmtPnt; npnt++)
1197  {
1198  ElmtJac[nelmt][npnt] =
1199  MemoryManager<DNekMat>::AllocateSharedPtr(
1200  nConvectiveFields, nConvectiveFields);
1201  }
1202  }
1203  }
1204 
1205  // Auxiliary variables
1206  Array<OneD, NekDouble> mu(nPts, 0.0);
1207  Array<OneD, NekDouble> thermalConductivity(nPts, 0.0);
1208  Array<OneD, NekDouble> temperature(nPts, 0.0);
1209  m_varConv->GetTemperature(inarray, temperature);
1210  GetViscosityAndThermalCondFromTemp(temperature, mu, thermalConductivity);
1211 
1212  // What about thermal conductivity?
1213 
1214  NekDouble pointmu = 0.0;
1215  Array<OneD, NekDouble> locmu;
1216  Array<OneD, NekDouble> pointVar(nConvectiveFields, 0.0);
1217  Array<OneD, Array<OneD, NekDouble>> locVars(nConvectiveFields);
1218  Array<OneD, NekDouble> pointnormals(nSpaceDim, 0.0);
1219  Array<OneD, Array<OneD, NekDouble>> locnormal(nSpaceDim);
1220 
1221  DNekMatSharedPtr PointFJac = MemoryManager<DNekMat>::AllocateSharedPtr(
1222  nConvectiveFields - 1, nConvectiveFields);
1223  Array<OneD, NekDouble> tmpMatinnData, tmpMatoutData;
1224  Array<OneD, NekDouble> tmp1, tmp2;
1225 
1226  for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
1227  {
1228  int nElmtPnt = (*expvect)[nelmt]->GetTotPoints();
1229  int noffest = explist->GetPhys_Offset(nelmt);
1230 
1231  for (int j = 0; j < nConvectiveFields; j++)
1232  {
1233  locVars[j] = inarray[j] + noffest;
1234  }
1235 
1236  for (int j = 0; j < nSpaceDim; j++)
1237  {
1238  locnormal[j] = normals[j] + noffest;
1239  }
1240 
1241  locmu = mu + noffest;
1242  for (int npnt = 0; npnt < nElmtPnt; npnt++)
1243  {
1244  for (int j = 0; j < nConvectiveFields; j++)
1245  {
1246  pointVar[j] = locVars[j][npnt];
1247  }
1248  for (int j = 0; j < nSpaceDim; j++)
1249  {
1250  pointnormals[j] = locnormal[j][npnt];
1251  }
1252 
1253  pointmu = locmu[npnt];
1254 
1255  m_GetdFlux_dDeriv_Array[nDervDir](pointnormals, pointmu, pointVar,
1256  PointFJac);
1257  tmpMatinnData = PointFJac->GetPtr();
1258  tmpMatoutData = ElmtJac[nelmt][npnt]->GetPtr();
1259 
1260  Vmath::Fill(nConvectiveFields, 0.0, tmpMatoutData,
1261  nConvectiveFields);
1262  for (int j = 0; j < nConvectiveFields; j++)
1263  {
1264  Vmath::Vcopy(
1265  nConvectiveFields - 1,
1266  tmp1 = tmpMatinnData + (j * (nConvectiveFields - 1)), 1,
1267  tmp2 = tmpMatoutData + (1 + j * nConvectiveFields), 1);
1268  }
1269  }
1270  }
1271 }
1272 
1273 void NavierStokesImplicitCFE::v_CalcPhysDeriv(
1274  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1275  TensorOfArray3D<NekDouble> &qfield)
1276 {
1277  int nConvectiveFields = m_fields.size();
1278  int npoints = GetTotPoints();
1279  const Array<OneD, Array<OneD, NekDouble>> pFwd;
1280  const Array<OneD, Array<OneD, NekDouble>> pBwd;
1281  if (!qfield.size())
1282  {
1283  qfield = TensorOfArray3D<NekDouble>(m_spacedim);
1284  for (int i = 0; i < m_spacedim; i++)
1285  {
1286  qfield[i] = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
1287  for (int j = 0; j < nConvectiveFields; j++)
1288  {
1289  qfield[i][j] = Array<OneD, NekDouble>(npoints, 0.0);
1290  }
1291  }
1292  }
1293  m_diffusion->DiffuseCalcDerivative(m_fields, inarray, qfield, pFwd, pBwd);
1294 }
1295 
1296 void NavierStokesImplicitCFE::v_CalcMuDmuDT(
1297  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1298  Array<OneD, NekDouble> &mu, Array<OneD, NekDouble> &DmuDT)
1299 {
1300  int nPts = mu.size();
1301 
1302  Array<OneD, NekDouble> thermalConductivity(nPts, 0.0);
1303  Array<OneD, NekDouble> temperature(nPts, 0.0);
1304  m_varConv->GetTemperature(inarray, temperature);
1305  GetViscosityAndThermalCondFromTemp(temperature, mu, thermalConductivity);
1306 
1307  if (m_ViscosityType == "Variable")
1308  {
1309  if (DmuDT.size() > 0)
1310  {
1311  m_varConv->GetDmuDT(temperature, mu, DmuDT);
1312  }
1313  }
1314  else
1315  {
1316  if (DmuDT.size() > 0)
1317  {
1318  Vmath::Zero(nPts, DmuDT, 1);
1319  }
1320  }
1321 }
1322 } // namespace Nektar
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::CFSImplicit::v_InitObject
virtual void v_InitObject(bool DeclareFields=true)
Initialization object for CFSImplicit class.
Definition: CompressibleFlowSystemImplicit.cpp:54
Nektar::CompressibleFlowSystem::m_diffusion
SolverUtils::DiffusionSharedPtr m_diffusion
Definition: CompressibleFlowSystem.h:94
Nektar::CompressibleFlowSystem::m_varConv
VariableConverterSharedPtr m_varConv
Definition: CompressibleFlowSystem.h:116
Nektar::CompressibleFlowSystem::m_artificialDiffusion
ArtificialDiffusionSharedPtr m_artificialDiffusion
Definition: CompressibleFlowSystem.h:95
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:168
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:646
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::NavierStokesImplicitCFE::v_GetDiffusionFluxJacPoint
virtual void v_GetDiffusionFluxJacPoint(const Array< OneD, NekDouble > &conservVar, const Array< OneD, const Array< OneD, NekDouble >> &conseDeriv, const NekDouble mu, const NekDouble DmuDT, const Array< OneD, NekDouble > &normals, DNekMatSharedPtr &fluxJac)
Definition: NavierStokesImplicitCFE.cpp:1031
Nektar::NavierStokesImplicitCFE::m_GetdFlux_dDeriv_Array
Array< OneD, GetdFlux_dDeriv > m_GetdFlux_dDeriv_Array
Definition: NavierStokesImplicitCFE.h:78
Nektar::NavierStokesImplicitCFE::v_CalcMuDmuDT
virtual void v_CalcMuDmuDT(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, NekDouble > &mu, Array< OneD, NekDouble > &DmuDT) override
Definition: NavierStokesImplicitCFE.cpp:1296
Nektar::NavierStokesImplicitCFE::GetdFlux_dQx_3D
void GetdFlux_dQx_3D(const Array< OneD, NekDouble > &normals, const NekDouble &mu, const Array< OneD, NekDouble > &U, DNekMatSharedPtr &OutputMatrix)
return part of viscous Jacobian derived with Qx=[drho_dx,drhou_dx,drhov_dx,drhow_dx,...
Definition: NavierStokesImplicitCFE.cpp:325
Nektar::NavierStokesImplicitCFE::v_DoDiffusionCoeff
virtual void v_DoDiffusionCoeff(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) override
Definition: NavierStokesImplicitCFE.cpp:107
Nektar::NavierStokesImplicitCFE::GetdFlux_dU_2D
void GetdFlux_dU_2D(const Array< OneD, NekDouble > &normals, const NekDouble mu, const NekDouble dmu_dT, const Array< OneD, NekDouble > &U, const Array< OneD, const Array< OneD, NekDouble >> &qfield, DNekMatSharedPtr &OutputMatrix)
return part of viscous Jacobian Input: normals:Point normals mu: dynamicviscosity dmu_dT: mu's deriva...
Definition: NavierStokesImplicitCFE.cpp:545
Nektar::NavierStokesImplicitCFE::v_GetFluxDerivJacDirctnElmt
virtual void v_GetFluxDerivJacDirctnElmt(const int nConvectiveFields, const int nElmtPnt, const int nDervDir, const Array< OneD, const Array< OneD, NekDouble >> &locVars, const Array< OneD, NekDouble > &locmu, const Array< OneD, const Array< OneD, NekDouble >> &locnormal, DNekMatSharedPtr &wspMat, Array< OneD, Array< OneD, NekDouble >> &PntJacArray) override
Definition: NavierStokesImplicitCFE.cpp:1132
Nektar::NavierStokesImplicitCFE::GetdFlux_dQx_2D
void GetdFlux_dQx_2D(const Array< OneD, NekDouble > &normals, const NekDouble &mu, const Array< OneD, NekDouble > &U, DNekMatSharedPtr &OutputMatrix)
return part of viscous Jacobian:
Definition: NavierStokesImplicitCFE.cpp:207
Nektar::NavierStokesImplicitCFE::v_GetFluxDerivJacDirctn
virtual void v_GetFluxDerivJacDirctn(const MultiRegions::ExpListSharedPtr &explist, const Array< OneD, const Array< OneD, NekDouble >> &normals, const int nDervDir, const Array< OneD, const Array< OneD, NekDouble >> &inarray, TensorOfArray5D< NekDouble > &ElmtJacArray, const int nFluxDir) override
Definition: NavierStokesImplicitCFE.cpp:1053
Nektar::NavierStokesImplicitCFE::v_InitObject
virtual void v_InitObject(bool DeclareFields=true) override
Initialization object for CompressibleFlowSystem class.
Definition: NavierStokesImplicitCFE.cpp:62
Nektar::NavierStokesImplicitCFE::GetdFlux_dQy_3D
void GetdFlux_dQy_3D(const Array< OneD, NekDouble > &normals, const NekDouble &mu, const Array< OneD, NekDouble > &U, DNekMatSharedPtr &OutputMatrix)
return part of viscous Jacobian derived with Qy=[drho_dy,drhou_dy,drhov_dy,drhow_dy,...
Definition: NavierStokesImplicitCFE.cpp:398
Nektar::NavierStokesImplicitCFE::v_MinusDiffusionFluxJacPoint
virtual void v_MinusDiffusionFluxJacPoint(const int nConvectiveFields, const int nElmtPnt, const Array< OneD, const Array< OneD, NekDouble >> &locVars, const TensorOfArray3D< NekDouble > &locDerv, const Array< OneD, NekDouble > &locmu, const Array< OneD, NekDouble > &locDmuDT, const Array< OneD, NekDouble > &normals, DNekMatSharedPtr &wspMat, Array< OneD, Array< OneD, NekDouble >> &PntJacArray) override
Definition: NavierStokesImplicitCFE.cpp:976
Nektar::NavierStokesImplicitCFE::v_CalcPhysDeriv
virtual void v_CalcPhysDeriv(const Array< OneD, const Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfield) override
Definition: NavierStokesImplicitCFE.cpp:1273
Nektar::NavierStokesImplicitCFE::GetdFlux_dQz_3D
void GetdFlux_dQz_3D(const Array< OneD, NekDouble > &normals, const NekDouble &mu, const Array< OneD, NekDouble > &U, DNekMatSharedPtr &OutputMatrix)
return part of viscous Jacobian derived with Qz=[drho_dz,drhou_dz,drhov_dz,drhow_dz,...
Definition: NavierStokesImplicitCFE.cpp:473
Nektar::NavierStokesImplicitCFE::GetdFlux_dQy_2D
void GetdFlux_dQy_2D(const Array< OneD, NekDouble > &normals, const NekDouble &mu, const Array< OneD, NekDouble > &U, DNekMatSharedPtr &OutputMatrix)
return part of viscous Jacobian:
Definition: NavierStokesImplicitCFE.cpp:265
Nektar::NavierStokesImplicitCFE::GetdFlux_dU_3D
void GetdFlux_dU_3D(const Array< OneD, NekDouble > &normals, const NekDouble mu, const NekDouble dmu_dT, const Array< OneD, NekDouble > &U, const Array< OneD, const Array< OneD, NekDouble >> &qfield, DNekMatSharedPtr &OutputMatrix)
return part of viscous Jacobian Input: normals:Point normals mu: dynamicviscosity dmu_dT: mu's deriva...
Definition: NavierStokesImplicitCFE.cpp:713
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:368
Nektar::SolverUtils::EquationSystem::m_graph
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
Definition: EquationSystem.h:339
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:695
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:335
Nektar::SolverUtils::EquationSystem::GetNpoints
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:730
Nektar::SolverUtils::EquationSystem::GetNcoeffs
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:662
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:720
Nektar::SolverUtils::UnsteadySystem
Base class for unsteady solvers.
Definition: UnsteadySystem.h:49
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:115
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:55
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:85
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
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:884
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
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:359
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
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::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
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:419
Generated on Tue Aug 23 2022 07:27:14 for Nektar++ by doxygen 1.9.1