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