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