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