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