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