MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
fvstructuredsolver3drans.cpp
Go to the documentation of this file.
1// Copyright (C) 2024 The m-AIA AUTHORS
2//
3// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
4//
5// SPDX-License-Identifier: LGPL-3.0-only
6
7#include "fvstructuredsolver3drans.h"
8#include "globals.h"
9
10#include <cstdlib>
11#include "GRID/structuredpartition.h"
12#include "fvstructuredsolverwindowinfo.h"
13#if not defined(MAIA_MS_COMPILER)
14#include <unistd.h>
15#endif
16// temporaray
17#include <vector>
18
19using namespace std;
20
21
22FvStructuredSolver3DRans::FvStructuredSolver3DRans(FvStructuredSolver3D* solver)
23 : m_StructuredComm(solver->m_StructuredComm),
24 m_solverId(solver->m_solverId),
25 m_nCells(solver->m_nCells),
26 m_nPoints(solver->m_nPoints),
27 m_noCells(solver->m_noCells),
28 m_cells(solver->m_cells),
29 CV(solver->CV),
30 PV(solver->PV),
31 FQ(solver->FQ),
32 m_noGhostLayers(solver->m_noGhostLayers),
33 m_eps(solver->m_eps),
34 m_chi(solver->m_chi),
35 m_sutherlandConstant(solver->m_sutherlandConstant),
36 m_sutherlandPlusOne(solver->m_sutherlandPlusOne) {
37 m_solver = solver;
38
39 initFluxMethod();
40}
41
42FvStructuredSolver3DRans::~FvStructuredSolver3DRans() {}
43
44void FvStructuredSolver3DRans::initFluxMethod() {
45 // m_structuredBndryCndRans = new StructuredBndryCnd3DRans(m_solver, m_solver->m_noSpecies);
46 // set pointer to the right AUSM for the RANS equations;
47
60 m_ransMethod =
61 static_cast<RansMethod>(string2enum(Context::getSolverProperty<MString>("ransMethod", m_solverId, AT_)));
62 m_dsIsComputed = false;
63 switch(m_ransMethod) {
64 case RANS_SA_DV: {
65 mAlloc(m_cells->saFlux1, nDim, m_noCells, "m_cells->saFlux1", -999999.9, AT_);
66 mAlloc(m_cells->saFlux2, nDim, m_noCells, "m_cells->saFlux2", -999999.9, AT_);
67 mAlloc(m_cells->prodDest, m_noCells, "m_cells->prodDest", -999999.9, AT_);
68 compTurbVisc = &FvStructuredSolver3DRans::computeTurbViscosity_SA;
69 viscFluxMethod = &FvStructuredSolver3DRans::viscousFlux_SA;
70 switch(m_solver->CV->noVariables) {
71 case 5: {
72 reconstructSurfaceData = &FvStructuredSolver3DRans::Muscl_AusmDV<5>;
73 break;
74 }
75 case 6: {
76 reconstructSurfaceData = &FvStructuredSolver3DRans::Muscl_AusmDV<6>;
77 break;
78 }
79 case 7: {
80 reconstructSurfaceData = &FvStructuredSolver3DRans::Muscl_AusmDV<7>;
81 break;
82 }
83 default: {
84 stringstream errorMessage;
85 errorMessage << "Number of Variables " << m_solver->CV->noVariables
86 << " not implemented! in temlate Rans AUSM " << endl;
87 mTerm(1, AT_, errorMessage.str());
88 }
89 }
90 break;
91 }
92 case RANS_FS: {
93 viscFluxMethod = &FvStructuredSolver3DRans::viscousFlux_FS;
94 compTurbVisc = &FvStructuredSolver3DRans::computeTurbViscosity_FS;
95
96 if(m_solver->m_limiter) {
97 mAlloc(m_cells->ql, PV->noVariables, m_noCells, "m_cells->ql", F0, AT_);
98 mAlloc(m_cells->qr, PV->noVariables, m_noCells, "m_cells->qr", F0, AT_);
99 cout << "Using RANS with Fares-Schroeder model and limited AusmDV" << endl;
100 switch(m_solver->CV->noVariables) {
101 case 5: {
102 reconstructSurfaceData = &FvStructuredSolver3DRans::Muscl_AusmDV_Limited<5>;
103 break;
104 }
105 case 6: {
106 reconstructSurfaceData = &FvStructuredSolver3DRans::Muscl_AusmDV_Limited<6>;
107 break;
108 }
109 case 7: {
110 reconstructSurfaceData = &FvStructuredSolver3DRans::Muscl_AusmDV_Limited<7>;
111 break;
112 }
113 default: {
114 stringstream errorMessage;
115 errorMessage << "Number of Variables " << m_solver->CV->noVariables
116 << " not implemented! in temlate Rans AUSM " << endl;
117 mTerm(1, AT_, errorMessage.str());
118 }
119 }
120 } else {
121 cout << "Using RANS with Fares-Schroeder model and AusmDV" << endl;
122 switch(m_solver->CV->noVariables) {
123 case 5: {
124 reconstructSurfaceData = &FvStructuredSolver3DRans::Muscl_AusmDV<5>;
125 break;
126 }
127 case 6: {
128 reconstructSurfaceData = &FvStructuredSolver3DRans::Muscl_AusmDV<6>;
129 break;
130 }
131 case 7: {
132 reconstructSurfaceData = &FvStructuredSolver3DRans::Muscl_AusmDV<7>;
133 break;
134 }
135 default: {
136 stringstream errorMessage;
137 errorMessage << "Number of Variables " << m_solver->CV->noVariables
138 << " not implemented! in temlate Rans AUSM " << endl;
139 mTerm(1, AT_, errorMessage.str());
140 }
141 }
142 }
143 break;
144 }
145 default: {
146 mTerm(1, AT_, "RANS METHOD wsa not specified in properties");
147 break;
148 }
149 }
150}
151
152void FvStructuredSolver3DRans::Muscl(MInt) { (this->*reconstructSurfaceData)(); }
153
154
155void FvStructuredSolver3DRans::viscousFluxRANS() { (this->*viscFluxMethod)(); }
156
157
158void FvStructuredSolver3DRans::viscousFlux_SA() {
159 computeTurbViscosity_SA();
160
161 // call the standard LES viscous flux
162 m_solver->viscousFluxLES<>();
163
164 // OTHER variables required to calculate the laminar viscous fluxes
165 const MFloat rRe = F1 / m_solver->m_Re0;
166
167 MFloat* const RESTRICT u = ALIGNED_F(&m_cells->pvariables[PV->U][0]);
168 MFloat* const RESTRICT v = ALIGNED_F(&m_cells->pvariables[PV->V][0]);
169 MFloat* const RESTRICT w = ALIGNED_F(&m_cells->pvariables[PV->W][0]);
170 MFloat* const RESTRICT rho = ALIGNED_F(&m_cells->pvariables[PV->RHO][0]);
171 MFloat* const RESTRICT nuTilde = ALIGNED_F(&m_cells->pvariables[PV->RANS_VAR[0]][0]);
172 MFloat* const RESTRICT muLam = ALIGNED_F(&m_cells->fq[FQ->MU_L][0]);
173 MFloat* const RESTRICT T = ALIGNED_F(&m_cells->temperature[0]);
174
175 MFloat* const* const RESTRICT eflux = ALIGNED_F(m_cells->eFlux);
176 MFloat* const* const RESTRICT fflux = ALIGNED_F(m_cells->fFlux);
177 MFloat* const* const RESTRICT gflux = ALIGNED_F(m_cells->gFlux);
178 MFloat* const* const RESTRICT sa_1flux = ALIGNED_F(m_cells->saFlux1);
179 MFloat* const* const RESTRICT sa_2flux = ALIGNED_F(m_cells->saFlux2);
180
181 for(MInt k = m_noGhostLayers - 1; k < m_nCells[0] - m_noGhostLayers + 1; k++) {
182 for(MInt j = m_noGhostLayers - 1; j < m_nCells[1] - m_noGhostLayers + 1; j++) {
183 for(MInt i = m_noGhostLayers - 1; i < m_nCells[2] - m_noGhostLayers + 1; i++) {
184 // get the adjacent cells;
185 const MInt IJK = cellIndex(i, j, k);
186 const MInt IPJK = cellIndex((i + 1), j, k);
187 const MInt IPJPK = cellIndex((i + 1), (j + 1), k);
188 const MInt IJPK = cellIndex(i, (j + 1), k);
189 const MInt IJKP = cellIndex(i, j, (k + 1));
190 const MInt IPJKP = cellIndex((i + 1), j, (k + 1));
191 const MInt IPJPKP = cellIndex((i + 1), (j + 1), (k + 1));
192 const MInt IJPKP = cellIndex(i, (j + 1), (k + 1));
193
194 const MInt IMJK = cellIndex((i - 1), j, k);
195 const MInt IJMK = cellIndex(i, (j - 1), k);
196 const MInt IJKM = cellIndex(i, j, (k - 1));
197
198 const MFloat cornerMetrics[9] = {
199 m_cells->cornerMetrics[0][IJK], m_cells->cornerMetrics[1][IJK], m_cells->cornerMetrics[2][IJK],
200 m_cells->cornerMetrics[3][IJK], m_cells->cornerMetrics[4][IJK], m_cells->cornerMetrics[5][IJK],
201 m_cells->cornerMetrics[6][IJK], m_cells->cornerMetrics[7][IJK], m_cells->cornerMetrics[8][IJK]};
202
203
204 const MFloat dnutldxi = F1B4
205 * (nuTilde[IPJPKP] + nuTilde[IPJPK] + nuTilde[IPJKP] + nuTilde[IPJK] - nuTilde[IJPKP]
206 - nuTilde[IJPK] - nuTilde[IJKP] - nuTilde[IJK]);
207 const MFloat dnutldet = F1B4
208 * (nuTilde[IPJPKP] + nuTilde[IJPKP] + nuTilde[IPJPK] + nuTilde[IJPK] - nuTilde[IPJKP]
209 - nuTilde[IJKP] - nuTilde[IPJK] - nuTilde[IJK]);
210 const MFloat dnutldze = F1B4
211 * (nuTilde[IPJPKP] + nuTilde[IJPKP] + nuTilde[IPJKP] + nuTilde[IJKP] - nuTilde[IPJPK]
212 - nuTilde[IJPK] - nuTilde[IPJK] - nuTilde[IJK]);
213
214
215 const MFloat nutldAvg = F1B8
216 * (nuTilde[IPJPKP] + nuTilde[IJPKP] + nuTilde[IJPK] + nuTilde[IPJPK] + nuTilde[IPJKP]
217 + nuTilde[IJKP] + nuTilde[IJK] + nuTilde[IPJK]);
218
219 const MFloat nuLamAvg = F1B8
220 * (muLam[IPJPKP] / rho[IPJPKP] + muLam[IJPKP] / rho[IJPKP] + muLam[IJPK] / rho[IJPK]
221 + muLam[IPJPK] / rho[IPJPK] + muLam[IPJKP] / rho[IPJKP] + muLam[IJKP] / rho[IJKP]
222 + muLam[IJK] / rho[IJK] + muLam[IPJK] / rho[IPJK]);
223
224 const MFloat dnutldx = dnutldxi * cornerMetrics[xsd * 3 + xsd] + dnutldet * cornerMetrics[ysd * 3 + xsd]
225 + dnutldze * cornerMetrics[zsd * 3 + xsd];
226
227 const MFloat dnutldy = dnutldxi * cornerMetrics[xsd * 3 + ysd] + dnutldet * cornerMetrics[ysd * 3 + ysd]
228 + dnutldze * cornerMetrics[zsd * 3 + ysd];
229
230 const MFloat dnutldz = dnutldxi * cornerMetrics[xsd * 3 + zsd] + dnutldet * cornerMetrics[ysd * 3 + zsd]
231 + dnutldze * cornerMetrics[zsd * 3 + zsd];
232
233 const MFloat Frj = rRe / m_cells->cornerJac[IJK];
234
235 const MFloat sax1 = Frj * (nuLamAvg + (1.0 + RM_SA_DV::cb2) * nutldAvg)
236 * (dnutldx * cornerMetrics[xsd * 3 + xsd] + dnutldy * cornerMetrics[xsd * 3 + ysd]
237 + dnutldz * cornerMetrics[xsd * 3 + zsd]);
238
239 const MFloat say1 = Frj * (nuLamAvg + (1.0 + RM_SA_DV::cb2) * nutldAvg)
240 * (dnutldx * cornerMetrics[ysd * 3 + xsd] + dnutldy * cornerMetrics[ysd * 3 + ysd]
241 + dnutldz * cornerMetrics[ysd * 3 + zsd]);
242
243 const MFloat saz1 = Frj * (nuLamAvg + (1.0 + RM_SA_DV::cb2) * nutldAvg)
244 * (dnutldx * cornerMetrics[zsd * 3 + xsd] + dnutldy * cornerMetrics[zsd * 3 + ysd]
245 + dnutldz * cornerMetrics[zsd * 3 + zsd]);
246
247 const MFloat sax2 = -Frj * RM_SA_DV::cb2
248 * (dnutldx * cornerMetrics[xsd * 3 + xsd] + dnutldy * cornerMetrics[xsd * 3 + ysd]
249 + dnutldz * cornerMetrics[xsd * 3 + zsd]);
250 const MFloat say2 = -Frj * RM_SA_DV::cb2
251 * (dnutldx * cornerMetrics[ysd * 3 + xsd] + dnutldy * cornerMetrics[ysd * 3 + ysd]
252 + dnutldz * cornerMetrics[ysd * 3 + zsd]);
253
254 const MFloat saz2 = -Frj * RM_SA_DV::cb2
255 * (dnutldx * cornerMetrics[zsd * 3 + xsd] + dnutldy * cornerMetrics[zsd * 3 + ysd]
256 + dnutldz * cornerMetrics[zsd * 3 + zsd]);
257
258 // for dwdy
259 const MFloat dw1 = w[IPJK] - w[IMJK];
260 const MFloat dw2 = w[IJPK] - w[IJMK];
261 const MFloat dw3 = w[IJKP] - w[IJKM];
262
263 // for dvdz
264 const MFloat dv1 = v[IPJK] - v[IMJK];
265 const MFloat dv2 = v[IJPK] - v[IJMK];
266 const MFloat dv3 = v[IJKP] - v[IJKM];
267
268 // for dudz
269 const MFloat du1 = u[IPJK] - u[IMJK];
270 const MFloat du2 = u[IJPK] - u[IJMK];
271 const MFloat du3 = u[IJKP] - u[IJKM];
272
273 const MFloat vorti =
274 (m_cells->cellMetrics[xsd * 3 + ysd][IJK] * dw1) + (m_cells->cellMetrics[ysd * 3 + ysd][IJK] * dw2)
275 + (m_cells->cellMetrics[zsd * 3 + ysd][IJK] * dw3) -
276
277 (m_cells->cellMetrics[xsd * 3 + zsd][IJK] * dv1) - (m_cells->cellMetrics[ysd * 3 + zsd][IJK] * dv2)
278 - (m_cells->cellMetrics[zsd * 3 + zsd][IJK] * dv3);
279
280 const MFloat vortj =
281 (m_cells->cellMetrics[xsd * 3 + zsd][IJK] * du1) + (m_cells->cellMetrics[ysd * 3 + zsd][IJK] * du2)
282 + (m_cells->cellMetrics[zsd * 3 + zsd][IJK] * du3) -
283
284 (m_cells->cellMetrics[xsd * 3 + xsd][IJK] * dw1) - (m_cells->cellMetrics[ysd * 3 + xsd][IJK] * dw2)
285 - (m_cells->cellMetrics[zsd * 3 + xsd][IJK] * dw3);
286
287 const MFloat vortk =
288 (m_cells->cellMetrics[xsd * 3 + xsd][IJK] * dv1) + (m_cells->cellMetrics[ysd * 3 + xsd][IJK] * dv2)
289 + (m_cells->cellMetrics[zsd * 3 + xsd][IJK] * dv3) -
290
291 (m_cells->cellMetrics[xsd * 3 + ysd][IJK] * du1) - (m_cells->cellMetrics[ysd * 3 + ysd][IJK] * du2)
292 - (m_cells->cellMetrics[zsd * 3 + ysd][IJK] * du3);
293
294 MFloat s = (vorti * vorti) + (vortj * vortj) + (vortk * vortk);
295 s = F1B2 * sqrt(s) / m_cells->cellJac[IJK];
296
297 // assuming wall distance function
298 const MFloat distance = m_cells->fq[FQ->WALLDISTANCE][IJK];
299 const MFloat Fdist2 = 1.0 / (distance * distance);
300 const MFloat chi = nuTilde[IJK] * rho[IJK] / (SUTHERLANDLAW(T[IJK]) / rho[IJK]);
301 const MFloat chip3 = chi * chi * chi;
302 const MFloat Fv1 = chip3 / (chip3 + RM_SA_DV::cv1to3);
303 const MFloat Fv2 = F1 - (chi / (F1 + chi * Fv1));
304
305 const MFloat term = nuTilde[IJK] * Fdist2 * RM_SA_DV::Fkap2;
306 const MFloat stilde = s + term * Fv2 * rRe;
307 const MFloat r = min(10.0, rRe * term / stilde);
308
309 const MFloat g = r + RM_SA_DV::cw2 * (pow(r, 6) - r);
310 const MFloat Fwterm = (1 + RM_SA_DV::cw3to6) / (pow(g, 6) + RM_SA_DV::cw3to6);
311 const MFloat Fw = g * pow(Fwterm, (1.0 / 6.0));
312 const MFloat prodValue = rho[IJK] * RM_SA_DV::cb1 * (F1 - RM_SA_DV::Ft2) * stilde * nuTilde[IJK];
313 const MFloat destValue = rRe * rho[IJK] * (RM_SA_DV::cw1 * Fw - RM_SA_DV::cb1 * RM_SA_DV::Fkap2 * RM_SA_DV::Ft2)
314 * pow(nuTilde[IJK], 2.0) * Fdist2;
315
316 m_cells->prodDest[IJK] = (prodValue - destValue) * m_cells->cellJac[IJK];
317
318 eflux[0][IJK] = sax1; // diffusion of nutilde for every cell
319 eflux[1][IJK] = sax2; // diffusion of nutilde for every cell
320
321 fflux[0][IJK] = say1; // diffusion of nutilde for every cell
322 fflux[1][IJK] = say2; // diffusion of nutilde for every cell
323
324 gflux[0][IJK] = saz1; // diffusion of nutilde for every cell
325 gflux[1][IJK] = saz2; // diffusion of nutilde for every cell
326 }
327 }
328 }
329
330 for(MInt k = m_noGhostLayers; k < m_nCells[0] - m_noGhostLayers; k++) {
331 for(MInt j = m_noGhostLayers; j < m_nCells[1] - m_noGhostLayers; j++) {
332 for(MInt i = m_noGhostLayers - 1; i < m_nCells[2] - m_noGhostLayers; i++) {
333 const MInt IJK = cellIndex(i, j, k);
334 const MInt IJMK = cellIndex(i, (j - 1), k);
335 const MInt IJKM = cellIndex(i, j, (k - 1));
336 const MInt IJMKM = cellIndex(i, (j - 1), (k - 1));
337
338 sa_1flux[0][IJK] = F1B4 * (eflux[0][IJK] + eflux[0][IJKM] + eflux[0][IJMK] + eflux[0][IJMKM]);
339 sa_2flux[0][IJK] = F1B4 * (eflux[1][IJK] + eflux[1][IJKM] + eflux[1][IJMK] + eflux[1][IJMKM]);
340 }
341 }
342 }
343
344
345 for(MInt k = m_noGhostLayers; k < m_nCells[0] - m_noGhostLayers; k++) {
346 for(MInt j = m_noGhostLayers - 1; j < m_nCells[1] - m_noGhostLayers; j++) {
347 for(MInt i = m_noGhostLayers; i < m_nCells[2] - m_noGhostLayers; i++) {
348 const MInt IJK = cellIndex(i, j, k);
349 const MInt IMJK = cellIndex((i - 1), j, k);
350 const MInt IJKM = cellIndex(i, j, (k - 1));
351 const MInt IMJKM = cellIndex((i - 1), j, (k - 1));
352
353 sa_1flux[1][IJK] = F1B4 * (fflux[0][IJK] + fflux[0][IJKM] + fflux[0][IMJK] + fflux[0][IMJKM]);
354 sa_2flux[1][IJK] = F1B4 * (fflux[1][IJK] + fflux[1][IJKM] + fflux[1][IMJK] + fflux[1][IMJKM]);
355 }
356 }
357 }
358
359 for(MInt k = m_noGhostLayers - 1; k < m_nCells[0] - m_noGhostLayers; k++) {
360 for(MInt j = m_noGhostLayers; j < m_nCells[1] - m_noGhostLayers; j++) {
361 for(MInt i = m_noGhostLayers; i < m_nCells[2] - m_noGhostLayers; i++) {
362 const MInt IJK = cellIndex(i, j, k);
363 const MInt IMJK = cellIndex((i - 1), j, k);
364 const MInt IJMK = cellIndex(i, (j - 1), k);
365 const MInt IMJMK = cellIndex((i - 1), (j - 1), k);
366
367 sa_1flux[2][IJK] = F1B4 * (gflux[0][IJK] + gflux[0][IMJK] + gflux[0][IJMK] + gflux[0][IMJMK]);
368 sa_2flux[2][IJK] = F1B4 * (gflux[1][IJK] + gflux[1][IMJK] + gflux[1][IJMK] + gflux[1][IMJMK]);
369 }
370 }
371 }
372
373 // separate loop for adding the prodn nad destrn terms for tur kin viscosity transport variable
374 for(MInt k = m_noGhostLayers; k < m_nCells[0] - m_noGhostLayers; k++) {
375 for(MInt j = m_noGhostLayers; j < m_nCells[1] - m_noGhostLayers; j++) {
376 for(MInt i = m_noGhostLayers; i < m_nCells[2] - m_noGhostLayers; i++) {
377 const MInt IJK = cellIndex(i, j, k);
378 const MInt IMJK = cellIndex(i - 1, j, k);
379 const MInt IJMK = cellIndex(i, j - 1, k);
380 const MInt IJKM = cellIndex(i, j, k - 1);
381 const MFloat dissipation_term =
382 (((sa_1flux[0][IJK] - sa_1flux[0][IMJK]) + ((sa_2flux[0][IJK] - sa_2flux[0][IMJK]) * nuTilde[IJK]))
383 * rho[IJK] * RM_SA_DV::Fsigma
384 + ((sa_1flux[1][IJK] - sa_1flux[1][IJMK]) + ((sa_2flux[1][IJK] - sa_2flux[1][IJMK]) * nuTilde[IJK]))
385 * rho[IJK] * RM_SA_DV::Fsigma
386 + ((sa_1flux[2][IJK] - sa_1flux[2][IJKM]) + ((sa_2flux[2][IJK] - sa_2flux[2][IJKM]) * nuTilde[IJK]))
387 * rho[IJK] * RM_SA_DV::Fsigma);
388
389 m_cells->rightHandSide[CV->RANS_FIRST][IJK] += dissipation_term;
390 m_cells->rightHandSide[CV->RANS_FIRST][IJK] += m_cells->prodDest[IJK];
391 }
392 }
393 }
394}
395
396
397void FvStructuredSolver3DRans::viscousFlux_FS() {
398 const MFloat eps = 1e-16;
399 const MFloat epss = 1e-34;
400
401 computeTurbViscosity_FS();
402
403 // call the standard LES viscous flux
404 m_solver->viscousFluxLES<>();
405
406 // OTHER variables required to calculate the laminar viscous fluxes
407 const MInt noCells = m_noCells;
408 const MFloat rRe = F1 / m_solver->m_Re0;
409
410 MFloat* const RESTRICT u = ALIGNED_F(&m_cells->pvariables[PV->U][0]);
411 MFloat* const RESTRICT v = ALIGNED_F(&m_cells->pvariables[PV->V][0]);
412 MFloat* const RESTRICT w = ALIGNED_F(&m_cells->pvariables[PV->W][0]);
413 MFloat* const RESTRICT rho = ALIGNED_F(&m_cells->pvariables[PV->RHO][0]);
414 MFloat* const RESTRICT nuTilde = ALIGNED_F(&m_cells->pvariables[PV->RANS_VAR[0]][0]);
415 MFloat* const RESTRICT muLam = ALIGNED_F(&m_cells->fq[FQ->MU_L][0]);
416
417 MFloat* const* const RESTRICT eflux = ALIGNED_F(m_cells->eFlux);
418 MFloat* const* const RESTRICT fflux = ALIGNED_F(m_cells->fFlux);
419 MFloat* const* const RESTRICT gflux = ALIGNED_F(m_cells->gFlux);
420
421 MFloatScratchSpace uvwn(noCells, 4, nDim, AT_, "uvwn");
422
423 for(MInt k = m_noGhostLayers; k < m_nCells[0] - m_noGhostLayers; k++) {
424 for(MInt j = m_noGhostLayers; j < m_nCells[1] - m_noGhostLayers; j++) {
425 for(MInt i = m_noGhostLayers; i < m_nCells[2] - m_noGhostLayers; i++) {
426 // get the adjacent cells;
427 const MInt IJK = cellIndex(i, j, k);
428 const MInt IPJK = cellIndex((i + 1), j, k);
429 const MInt IJPK = cellIndex(i, (j + 1), k);
430 const MInt IJKP = cellIndex(i, j, (k + 1));
431
432 const MInt IMJK = cellIndex((i - 1), j, k);
433 const MInt IJMK = cellIndex(i, (j - 1), k);
434 const MInt IJKM = cellIndex(i, j, (k - 1));
435
436 const MFloat cellMetrics[9] = {
437 m_cells->cellMetrics[0][IJK], m_cells->cellMetrics[1][IJK], m_cells->cellMetrics[2][IJK],
438 m_cells->cellMetrics[3][IJK], m_cells->cellMetrics[4][IJK], m_cells->cellMetrics[5][IJK],
439 m_cells->cellMetrics[6][IJK], m_cells->cellMetrics[7][IJK], m_cells->cellMetrics[8][IJK]};
440
441 const MFloat fjac = F1B2 / m_cells->cellJac[IJK];
442
443 const MFloat dudxi = fjac * (u[IPJK] - u[IMJK]);
444 const MFloat dudet = fjac * (u[IJPK] - u[IJMK]);
445 const MFloat dudze = fjac * (u[IJKP] - u[IJKM]);
446
447 const MFloat dvdxi = fjac * (v[IPJK] - v[IMJK]);
448 const MFloat dvdet = fjac * (v[IJPK] - v[IJMK]);
449 const MFloat dvdze = fjac * (v[IJKP] - v[IJKM]);
450
451 const MFloat dwdxi = fjac * (w[IPJK] - w[IMJK]);
452 const MFloat dwdet = fjac * (w[IJPK] - w[IJMK]);
453 const MFloat dwdze = fjac * (w[IJKP] - w[IJKM]);
454
455 const MFloat dnutdxi = fjac * (nuTilde[IPJK] - nuTilde[IMJK]);
456 const MFloat dnutdet = fjac * (nuTilde[IJPK] - nuTilde[IJMK]);
457 const MFloat dnutdze = fjac * (nuTilde[IJKP] - nuTilde[IJKM]);
458
459
460 uvwn(IJK, 0, 0) = cellMetrics[0] * dudxi + cellMetrics[3] * dudet + cellMetrics[6] * dudze;
461 uvwn(IJK, 0, 1) = cellMetrics[1] * dudxi + cellMetrics[4] * dudet + cellMetrics[7] * dudze;
462 uvwn(IJK, 0, 2) = cellMetrics[2] * dudxi + cellMetrics[5] * dudet + cellMetrics[8] * dudze;
463
464 uvwn(IJK, 1, 0) = cellMetrics[0] * dvdxi + cellMetrics[3] * dvdet + cellMetrics[6] * dvdze;
465 uvwn(IJK, 1, 1) = cellMetrics[1] * dvdxi + cellMetrics[4] * dvdet + cellMetrics[7] * dvdze;
466 uvwn(IJK, 1, 2) = cellMetrics[2] * dvdxi + cellMetrics[5] * dvdet + cellMetrics[8] * dvdze;
467
468 uvwn(IJK, 2, 0) = cellMetrics[0] * dwdxi + cellMetrics[3] * dwdet + cellMetrics[6] * dwdze;
469 uvwn(IJK, 2, 1) = cellMetrics[1] * dwdxi + cellMetrics[4] * dwdet + cellMetrics[7] * dwdze;
470 uvwn(IJK, 2, 2) = cellMetrics[2] * dwdxi + cellMetrics[5] * dwdet + cellMetrics[8] * dwdze;
471
472 uvwn(IJK, 3, 0) = cellMetrics[0] * dnutdxi + cellMetrics[3] * dnutdet + cellMetrics[6] * dnutdze;
473 uvwn(IJK, 3, 1) = cellMetrics[1] * dnutdxi + cellMetrics[4] * dnutdet + cellMetrics[7] * dnutdze;
474 uvwn(IJK, 3, 2) = cellMetrics[2] * dnutdxi + cellMetrics[5] * dnutdet + cellMetrics[8] * dnutdze;
475 }
476 }
477 }
478
479
480 for(MInt var = 0; var < 4; var++) {
481 // i start
482 for(MInt k = m_noGhostLayers - 1; k < m_nCells[0] - m_noGhostLayers + 1; k++) {
483 for(MInt j = m_noGhostLayers - 1; j < m_nCells[1] - m_noGhostLayers + 1; j++) {
484 const MInt II1 = cellIndex(m_noGhostLayers - 1, j, k);
485 const MInt II2 = cellIndex(m_noGhostLayers, j, k);
486 const MInt II3 = cellIndex(m_noGhostLayers + 1, j, k);
487
488 uvwn(II1, var, 0) = 2.0 * uvwn(II2, var, 0) - uvwn(II3, var, 0);
489 uvwn(II1, var, 1) = 2.0 * uvwn(II2, var, 1) - uvwn(II3, var, 1);
490 uvwn(II1, var, 2) = 2.0 * uvwn(II2, var, 2) - uvwn(II3, var, 2);
491 }
492 }
493
494 // i end
495 for(MInt k = m_noGhostLayers - 1; k < m_nCells[0] - m_noGhostLayers + 1; k++) {
496 for(MInt j = m_noGhostLayers - 1; j < m_nCells[1] - m_noGhostLayers + 1; j++) {
497 const MInt II1 = cellIndex(m_nCells[2] - m_noGhostLayers + 1, j, k);
498 const MInt II2 = cellIndex(m_nCells[2] - m_noGhostLayers, j, k);
499 const MInt II3 = cellIndex(m_nCells[2] - m_noGhostLayers - 1, j, k);
500
501 uvwn(II1, var, 0) = 2.0 * uvwn(II2, var, 0) - uvwn(II3, var, 0);
502 uvwn(II1, var, 1) = 2.0 * uvwn(II2, var, 1) - uvwn(II3, var, 1);
503 uvwn(II1, var, 2) = 2.0 * uvwn(II2, var, 2) - uvwn(II3, var, 2);
504 }
505 }
506
507 // j start
508 for(MInt k = m_noGhostLayers - 1; k < m_nCells[0] - m_noGhostLayers + 1; k++) {
509 for(MInt i = m_noGhostLayers - 1; i < m_nCells[2] - m_noGhostLayers + 1; i++) {
510 const MInt II1 = cellIndex(i, m_noGhostLayers - 1, k);
511 const MInt II2 = cellIndex(i, m_noGhostLayers, k);
512 const MInt II3 = cellIndex(i, m_noGhostLayers + 1, k);
513
514 uvwn(II1, var, 0) = 2.0 * uvwn(II2, var, 0) - uvwn(II3, var, 0);
515 uvwn(II1, var, 1) = 2.0 * uvwn(II2, var, 1) - uvwn(II3, var, 1);
516 uvwn(II1, var, 2) = 2.0 * uvwn(II2, var, 2) - uvwn(II3, var, 2);
517 }
518 }
519
520 // j end
521 for(MInt k = m_noGhostLayers - 1; k < m_nCells[0] - m_noGhostLayers + 1; k++) {
522 for(MInt i = m_noGhostLayers - 1; i < m_nCells[2] - m_noGhostLayers + 1; i++) {
523 const MInt II1 = cellIndex(i, m_nCells[1] - m_noGhostLayers + 1, k);
524 const MInt II2 = cellIndex(i, m_nCells[1] - m_noGhostLayers, k);
525 const MInt II3 = cellIndex(i, m_nCells[1] - m_noGhostLayers - 1, k);
526
527 uvwn(II1, var, 0) = 2.0 * uvwn(II2, var, 0) - uvwn(II3, var, 0);
528 uvwn(II1, var, 1) = 2.0 * uvwn(II2, var, 1) - uvwn(II3, var, 1);
529 uvwn(II1, var, 2) = 2.0 * uvwn(II2, var, 2) - uvwn(II3, var, 2);
530 }
531 }
532
533
534 // k start
535 for(MInt j = m_noGhostLayers - 1; j < m_nCells[1] - m_noGhostLayers + 1; j++) {
536 for(MInt i = m_noGhostLayers - 1; i < m_nCells[2] - m_noGhostLayers + 1; i++) {
537 const MInt II1 = cellIndex(i, j, m_noGhostLayers - 1);
538 const MInt II2 = cellIndex(i, j, m_noGhostLayers);
539 const MInt II3 = cellIndex(i, j, m_noGhostLayers + 1);
540
541 uvwn(II1, var, 0) = 2.0 * uvwn(II2, var, 0) - uvwn(II3, var, 0);
542 uvwn(II1, var, 1) = 2.0 * uvwn(II2, var, 1) - uvwn(II3, var, 1);
543 uvwn(II1, var, 2) = 2.0 * uvwn(II2, var, 2) - uvwn(II3, var, 2);
544 }
545 }
546
547 // k end
548 for(MInt j = m_noGhostLayers - 1; j < m_nCells[1] - m_noGhostLayers + 1; j++) {
549 for(MInt i = m_noGhostLayers - 1; i < m_nCells[2] - m_noGhostLayers + 1; i++) {
550 const MInt II1 = cellIndex(i, j, m_nCells[0] - m_noGhostLayers + 1);
551 const MInt II2 = cellIndex(i, j, m_nCells[0] - m_noGhostLayers);
552 const MInt II3 = cellIndex(i, j, m_nCells[0] - m_noGhostLayers - 1);
553
554 uvwn(II1, var, 0) = 2.0 * uvwn(II2, var, 0) - uvwn(II3, var, 0);
555 uvwn(II1, var, 1) = 2.0 * uvwn(II2, var, 1) - uvwn(II3, var, 1);
556 uvwn(II1, var, 2) = 2.0 * uvwn(II2, var, 2) - uvwn(II3, var, 2);
557 }
558 }
559 }
560
561 MFloatScratchSpace omega(noCells, AT_, "omega");
562 MFloatScratchSpace dumm(noCells, AT_, "dumm");
563 omega.fill(0.0);
564 dumm.fill(0.0);
565
566 for(MInt j = 0; j < nDim; j++) {
567 for(MInt i = 0; i < nDim; i++) {
568 for(MInt kk = m_noGhostLayers - 1; kk < m_nCells[0] - m_noGhostLayers + 1; kk++) {
569 for(MInt jj = m_noGhostLayers - 1; jj < m_nCells[1] - m_noGhostLayers + 1; jj++) {
570 for(MInt ii = m_noGhostLayers - 1; ii < m_nCells[2] - m_noGhostLayers + 1; ii++) {
571 const MInt IJK = cellIndex(ii, jj, kk);
572 const MFloat Sij = F1B2 * (uvwn(IJK, i, j) + uvwn(IJK, j, i));
573 dumm(IJK) += Sij * Sij;
574 }
575 }
576 }
577 }
578 }
579
580 MFloat fac = 1.0 / sqrt(RM_FS::fabetcs);
581
582 for(MInt kk = m_noGhostLayers - 1; kk < m_nCells[0] - m_noGhostLayers + 1; kk++) {
583 for(MInt jj = m_noGhostLayers - 1; jj < m_nCells[1] - m_noGhostLayers + 1; jj++) {
584 for(MInt ii = m_noGhostLayers - 1; ii < m_nCells[2] - m_noGhostLayers + 1; ii++) {
585 const MInt IJK = cellIndex(ii, jj, kk);
586 omega(IJK) = mMax(eps, fac * sqrt(2.0 * dumm(IJK)));
587 }
588 }
589 }
590
591 dumm.fill(0.0);
592
593 for(MInt k = 0; k < nDim; k++) {
594 for(MInt j = 0; j < nDim; j++) {
595 for(MInt i = 0; i < nDim; i++) {
596 for(MInt kk = m_noGhostLayers; kk < m_nCells[0] - m_noGhostLayers; kk++) {
597 for(MInt jj = m_noGhostLayers; jj < m_nCells[1] - m_noGhostLayers; jj++) {
598 for(MInt ii = m_noGhostLayers; ii < m_nCells[2] - m_noGhostLayers; ii++) {
599 const MInt IJK = cellIndex(ii, jj, kk);
600 const MFloat Oij = 0.5 * (uvwn(IJK, i, j) - uvwn(IJK, j, i));
601 const MFloat Ojk = 0.5 * (uvwn(IJK, j, k) - uvwn(IJK, k, j));
602 const MFloat Ski = 0.5 * (uvwn(IJK, k, i) + uvwn(IJK, i, k));
603 dumm(IJK) += Oij * Ojk * Ski;
604 }
605 }
606 }
607 }
608 }
609 }
610
611 fac = 1.0 / pow(RM_FS::fabetcs, 3.0);
612
613 for(MInt k = m_noGhostLayers; k < m_nCells[0] - m_noGhostLayers; k++) {
614 for(MInt j = m_noGhostLayers; j < m_nCells[1] - m_noGhostLayers; j++) {
615 for(MInt i = m_noGhostLayers; i < m_nCells[2] - m_noGhostLayers; i++) {
616 const MInt IJK = cellIndex(i, j, k);
617 const MInt IPJK = cellIndex((i + 1), j, k);
618 const MInt IJPK = cellIndex(i, (j + 1), k);
619 const MInt IJKP = cellIndex(i, j, (k + 1));
620
621 const MInt IMJK = cellIndex((i - 1), j, k);
622 const MInt IJMK = cellIndex(i, (j - 1), k);
623 const MInt IJKM = cellIndex(i, j, (k - 1));
624
625 const MFloat fjac = 0.5 / m_cells->cellJac[IJK];
626
627 const MFloat cellMetrics[9] = {
628 m_cells->cellMetrics[0][IJK], m_cells->cellMetrics[1][IJK], m_cells->cellMetrics[2][IJK],
629 m_cells->cellMetrics[3][IJK], m_cells->cellMetrics[4][IJK], m_cells->cellMetrics[5][IJK],
630 m_cells->cellMetrics[6][IJK], m_cells->cellMetrics[7][IJK], m_cells->cellMetrics[8][IJK]};
631
632 const MFloat ddxi = fjac * (omega(IPJK) - omega(IMJK));
633 const MFloat ddeta = fjac * (omega(IJPK) - omega(IJMK));
634 const MFloat ddzeta = fjac * (omega(IJKP) - omega(IJKM));
635
636 const MFloat domdx = cellMetrics[0] * ddxi + cellMetrics[3] * ddeta + cellMetrics[6] * ddzeta;
637 const MFloat domdy = cellMetrics[1] * ddxi + cellMetrics[4] * ddeta + cellMetrics[7] * ddzeta;
638 const MFloat domdz = cellMetrics[2] * ddxi + cellMetrics[5] * ddeta + cellMetrics[8] * ddzeta;
639
640
641 const MFloat dudx = uvwn(IJK, 0, 0);
642 const MFloat dudy = uvwn(IJK, 0, 1);
643 const MFloat dudz = uvwn(IJK, 0, 2);
644
645 const MFloat dvdx = uvwn(IJK, 1, 0);
646 const MFloat dvdy = uvwn(IJK, 1, 1);
647 const MFloat dvdz = uvwn(IJK, 1, 2);
648
649 const MFloat dwdx = uvwn(IJK, 2, 0);
650 const MFloat dwdy = uvwn(IJK, 2, 1);
651 const MFloat dwdz = uvwn(IJK, 2, 2);
652
653 const MFloat dndx = uvwn(IJK, 3, 0);
654 const MFloat dndy = uvwn(IJK, 3, 1);
655 const MFloat dndz = uvwn(IJK, 3, 2);
656
657 const MFloat CDNOM = dndx * domdx + dndy * domdy + dndz * domdz;
658 const MFloat fpsik = 0.0;
659 const MFloat crdif = 0.0;
660
661 const MFloat betas = RM_FS::fabetcs * (1.0 + RM_FS::fapsik1 * fpsik) / (1.0 + RM_FS::fapsik2 * fpsik);
662 const MFloat beta = RM_FS::fabetc;
663
664 const MFloat P =
665 (dudy * (dudy + dvdx) + dudz * (dudz + dwdx) + dvdx * (dvdx + dudy) + dvdz * (dvdz + dwdy)
666 + dwdx * (dwdx + dudz) + dwdy * (dwdy + dvdz) + 2.0 * dudx * dudx + 2.0 * dvdy * dvdy + 2.0 * dwdz * dwdz);
667
668 const MFloat fv2t = 1.0;
669 const MFloat prod1 = (1.0 - RM_FS::faalpha * fv2t) * nuTilde[IJK] * rho[IJK] / omega(IJK) * P;
670 MFloat prod = prod1 - (1.0 - RM_FS::faalpha * fv2t) * F2B3 * nuTilde[IJK] * (dudx + dvdy + dwdz) * rho[IJK];
671 const MFloat dest = (betas - beta) * nuTilde[IJK] * omega(IJK) * rho[IJK] + rho[IJK] * crdif * 0.125;
672 const MFloat prodwall = mMin(
673 0.0, 2.0 * rRe * rho[IJK] / omega(IJK) * (muLam[IJK] / rho[IJK] + RM_FS::fasigma * nuTilde[IJK]) * CDNOM);
674
675 prod += prodwall;
676
677 m_cells->rightHandSide[CV->RANS_FIRST][IJK] += m_cells->cellJac[IJK] * (prod - dest);
678 }
679 }
680 }
681
682
683 for(MInt k = m_noGhostLayers - 1; k < m_nCells[0] - m_noGhostLayers + 1; k++) {
684 for(MInt j = m_noGhostLayers - 1; j < m_nCells[1] - m_noGhostLayers + 1; j++) {
685 for(MInt i = m_noGhostLayers - 1; i < m_nCells[2] - m_noGhostLayers + 1; i++) {
686 const MInt IJK = cellIndex(i, j, k);
687 const MInt IPJK = cellIndex((i + 1), j, k);
688 const MInt IPJPK = cellIndex((i + 1), (j + 1), k);
689 const MInt IJPK = cellIndex(i, (j + 1), k);
690 const MInt IJKP = cellIndex(i, j, (k + 1));
691 const MInt IPJKP = cellIndex((i + 1), j, (k + 1));
692 const MInt IPJPKP = cellIndex((i + 1), (j + 1), (k + 1));
693 const MInt IJPKP = cellIndex(i, (j + 1), (k + 1));
694
695 const MFloat cornerMetrics[9] = {
696 m_cells->cornerMetrics[0][IJK], m_cells->cornerMetrics[1][IJK], m_cells->cornerMetrics[2][IJK],
697 m_cells->cornerMetrics[3][IJK], m_cells->cornerMetrics[4][IJK], m_cells->cornerMetrics[5][IJK],
698 m_cells->cornerMetrics[6][IJK], m_cells->cornerMetrics[7][IJK], m_cells->cornerMetrics[8][IJK]};
699
700 const MFloat nutldAvg = F1B8
701 * (nuTilde[IPJPKP] + nuTilde[IJPKP] + nuTilde[IJPK] + nuTilde[IPJPK] + nuTilde[IPJKP]
702 + nuTilde[IJKP] + nuTilde[IJK] + nuTilde[IPJK]);
703
704 const MFloat nuLamAvg = F1B8
705 * (muLam[IPJPKP] / rho[IPJPKP] + muLam[IJPKP] / rho[IJPKP] + muLam[IJPK] / rho[IJPK]
706 + muLam[IPJPK] / rho[IPJPK] + muLam[IPJKP] / rho[IPJKP] + muLam[IJKP] / rho[IJKP]
707 + muLam[IJK] / rho[IJK] + muLam[IPJK] / rho[IPJK]);
708
709 const MFloat dnutldxi = F1B4
710 * (nuTilde[IPJPKP] + nuTilde[IPJPK] + nuTilde[IPJKP] + nuTilde[IPJK] - nuTilde[IJPKP]
711 - nuTilde[IJPK] - nuTilde[IJKP] - nuTilde[IJK]);
712 const MFloat dnutldet = F1B4
713 * (nuTilde[IPJPKP] + nuTilde[IJPKP] + nuTilde[IPJPK] + nuTilde[IJPK] - nuTilde[IPJKP]
714 - nuTilde[IJKP] - nuTilde[IPJK] - nuTilde[IJK]);
715 const MFloat dnutldze = F1B4
716 * (nuTilde[IPJPKP] + nuTilde[IJPKP] + nuTilde[IPJKP] + nuTilde[IJKP] - nuTilde[IPJPK]
717 - nuTilde[IJPK] - nuTilde[IPJK] - nuTilde[IJK]);
718
719 const MFloat dnutldx = (dnutldxi * cornerMetrics[xsd * 3 + xsd] + dnutldet * cornerMetrics[ysd * 3 + xsd]
720 + dnutldze * cornerMetrics[zsd * 3 + xsd]);
721
722 const MFloat dnutldy = (dnutldxi * cornerMetrics[xsd * 3 + ysd] + dnutldet * cornerMetrics[ysd * 3 + ysd]
723 + dnutldze * cornerMetrics[zsd * 3 + ysd]);
724
725 const MFloat dnutldz = (dnutldxi * cornerMetrics[xsd * 3 + zsd] + dnutldet * cornerMetrics[ysd * 3 + zsd]
726 + dnutldze * cornerMetrics[zsd * 3 + zsd]);
727
728 const MFloat Frj = rRe / mMax(fabs(m_cells->cornerJac[IJK]), epss);
729
730 eflux[6][IJK] = (Frj * (nuLamAvg + RM_FS::fasigma * nutldAvg)
731 * (dnutldx * cornerMetrics[xsd * 3 + xsd] + dnutldy * cornerMetrics[xsd * 3 + ysd]
732 + dnutldz * cornerMetrics[xsd * 3 + zsd]));
733
734 fflux[6][IJK] = (Frj * (nuLamAvg + RM_FS::fasigma * nutldAvg)
735 * (dnutldx * cornerMetrics[ysd * 3 + xsd] + dnutldy * cornerMetrics[ysd * 3 + ysd]
736 + dnutldz * cornerMetrics[ysd * 3 + zsd]));
737
738 gflux[6][IJK] = (Frj * (nuLamAvg + RM_FS::fasigma * nutldAvg)
739 * (dnutldx * cornerMetrics[zsd * 3 + xsd] + dnutldy * cornerMetrics[zsd * 3 + ysd]
740 + dnutldz * cornerMetrics[zsd * 3 + zsd]));
741 }
742 }
743 }
744
745 for(MInt k = m_noGhostLayers; k < m_nCells[0] - m_noGhostLayers; k++) {
746 for(MInt j = m_noGhostLayers; j < m_nCells[1] - m_noGhostLayers; j++) {
747 for(MInt i = m_noGhostLayers - 1; i < m_nCells[2] - m_noGhostLayers; i++) {
748 const MInt IJK = cellIndex(i, j, k);
749 const MInt IJMK = cellIndex(i, (j - 1), k);
750 const MInt IJKM = cellIndex(i, j, (k - 1));
751 const MInt IJMKM = cellIndex(i, (j - 1), (k - 1));
752
753 eflux[3][IJK] = F1B4 * (eflux[6][IJK] + eflux[6][IJKM] + eflux[6][IJMK] + eflux[6][IJMKM]);
754 }
755 }
756 }
757
758
759 for(MInt k = m_noGhostLayers; k < m_nCells[0] - m_noGhostLayers; k++) {
760 for(MInt j = m_noGhostLayers - 1; j < m_nCells[1] - m_noGhostLayers; j++) {
761 for(MInt i = m_noGhostLayers; i < m_nCells[2] - m_noGhostLayers; i++) {
762 const MInt IJK = cellIndex(i, j, k);
763 const MInt IMJK = cellIndex((i - 1), j, k);
764 const MInt IJKM = cellIndex(i, j, (k - 1));
765 const MInt IMJKM = cellIndex((i - 1), j, (k - 1));
766
767 fflux[3][IJK] = F1B4 * (fflux[6][IJK] + fflux[6][IJKM] + fflux[6][IMJK] + fflux[6][IMJKM]);
768 }
769 }
770 }
771
772 for(MInt k = m_noGhostLayers - 1; k < m_nCells[0] - m_noGhostLayers; k++) {
773 for(MInt j = m_noGhostLayers; j < m_nCells[1] - m_noGhostLayers; j++) {
774 for(MInt i = m_noGhostLayers; i < m_nCells[2] - m_noGhostLayers; i++) {
775 const MInt IJK = cellIndex(i, j, k);
776 const MInt IMJK = cellIndex((i - 1), j, k);
777 const MInt IJMK = cellIndex(i, (j - 1), k);
778 const MInt IMJMK = cellIndex((i - 1), (j - 1), k);
779
780 gflux[3][IJK] = F1B4 * (gflux[6][IJK] + gflux[6][IMJK] + gflux[6][IJMK] + gflux[6][IMJMK]);
781 }
782 }
783 }
784
785 // separate loop for adding the dissipation term for tur kin viscosity transport variable
786 for(MInt k = m_noGhostLayers; k < m_nCells[0] - m_noGhostLayers; k++) {
787 for(MInt j = m_noGhostLayers; j < m_nCells[1] - m_noGhostLayers; j++) {
788 for(MInt i = m_noGhostLayers; i < m_nCells[2] - m_noGhostLayers; i++) {
789 const MInt IJK = cellIndex(i, j, k);
790 const MInt IMJK = cellIndex(i - 1, j, k);
791 const MInt IJMK = cellIndex(i, j - 1, k);
792 const MInt IJKM = cellIndex(i, j, k - 1);
793 const MFloat dissipation_term =
794 ((eflux[3][IJK] - eflux[3][IMJK]) + (fflux[3][IJK] - fflux[3][IJMK]) + (gflux[3][IJK] - gflux[3][IJKM]));
795
796 m_cells->rightHandSide[CV->RANS_FIRST][IJK] += dissipation_term * rho[IJK];
797 }
798 }
799 }
800}
801
802
803void FvStructuredSolver3DRans::computeTurbViscosity() { (this->*compTurbVisc)(); }
804
805void FvStructuredSolver3DRans::computeTurbViscosity_SA() {
806 // OTHER variables required to calculate the laminar viscous fluxes
807 const MFloat transPos = m_solver->m_ransTransPos;
808 MFloat* const RESTRICT p = &m_cells->pvariables[PV->P][0];
809 MFloat* const RESTRICT rho = &m_cells->pvariables[PV->RHO][0];
810 MFloat* const RESTRICT nuTilde = &m_cells->pvariables[PV->RANS_VAR[0]][0];
811 MFloat* const RESTRICT T = &m_cells->temperature[0];
812
813 for(MInt i = 0; i < m_noCells; i++) {
814 if(m_cells->coordinates[0][i] <= transPos) {
815 nuTilde[i] = 0.0;
816 } else {
817 nuTilde[i] = mMin(nuTilde[i], 3000.0);
818 }
819 T[i] = m_solver->m_gamma * p[i] / rho[i];
820 // decode the kinematic turbulent viscosity from the turb dynamic visc arrays
821 const MFloat nuLaminar = SUTHERLANDLAW(T[i]) / rho[i];
822 const MFloat chi = nuTilde[i] / (nuLaminar);
823 const MFloat fv1 = pow(chi, 3) / (pow(chi, 3) + RM_SA_DV::cv1to3);
824 const MFloat nuTurb = fv1 * nuTilde[i];
825 m_cells->fq[FQ->NU_T][i] = nuTurb;
826 m_cells->fq[FQ->MU_T][i] = rho[i] * nuTurb;
827 }
828}
829
830void FvStructuredSolver3DRans::computeTurbViscosity_FS() {
831 // OTHER variables required to calculate the laminar viscous fluxes
832 const MFloat transPos = m_solver->m_ransTransPos;
833 MFloat* const RESTRICT p = &m_cells->pvariables[PV->P][0];
834 MFloat* const RESTRICT rho = &m_cells->pvariables[PV->RHO][0];
835 MFloat* const RESTRICT nuTilde = &m_cells->pvariables[PV->RANS_VAR[0]][0];
836 MFloat* const RESTRICT T = &m_cells->temperature[0];
837
838 for(MInt i = 0; i < m_noCells; i++) {
839 if(m_cells->coordinates[0][i] <= transPos) {
840 nuTilde[i] = 0.0;
841 } else {
842 nuTilde[i] = mMin(nuTilde[i], 1200.0);
843 }
844 T[i] = m_solver->m_gamma * p[i] / rho[i];
845 const MFloat nuLaminar = SUTHERLANDLAW(T[i]) / rho[i];
846 const MFloat chi = nuTilde[i] / nuLaminar;
847 const MFloat nuTurb = nuTilde[i] * pow(chi, 3.0) / (pow(chi, 3.0) + RM_FS::facv1to3);
848 m_cells->fq[FQ->NU_T][i] = nuTurb;
849 m_cells->fq[FQ->MU_T][i] = nuTurb * rho[i];
850 }
851}
852
853template <MInt noVars>
854void FvStructuredSolver3DRans::Muscl_AusmDV() {
855 TRACE();
856
857 // stencil identifier
858 const MInt IJK[3] = {1, m_nCells[2], m_nCells[1] * m_nCells[2]};
859
860 const MFloat* const RESTRICT x = ALIGNED_F(m_cells->coordinates[0]);
861 const MFloat* const RESTRICT y = ALIGNED_F(m_cells->coordinates[1]);
862 const MFloat* const RESTRICT z = ALIGNED_F(m_cells->coordinates[2]);
863 const MFloat* const* const RESTRICT vars = ALIGNED_F(m_cells->pvariables);
864 MFloat* const* const RESTRICT dss = ALIGNED_F(m_cells->dss);
865 MFloat* const* const RESTRICT flux = ALIGNED_F(m_cells->flux);
866 MFloat* const RESTRICT cellRhs = ALIGNED_MF(m_cells->rightHandSide[0]);
867
868 const MUint noCells = m_noCells;
869 const MFloat gamma = m_solver->m_gamma;
870 const MFloat gammaMinusOne = gamma - 1.0;
871 const MFloat FgammaMinusOne = F1 / gammaMinusOne;
872
873 const MInt noCellsI = m_nCells[2] - 2;
874 const MInt noCellsJ = m_nCells[1] - 2;
875 const MInt noCellsK = m_nCells[0] - 2;
876
877 const MInt noCellsIP1 = m_nCells[2] - 1;
878 const MInt noCellsJP1 = m_nCells[1] - 1;
879 const MInt noCellsKP1 = m_nCells[0] - 1;
880
881
882 if(!m_dsIsComputed) {
883 for(MInt dim = 0; dim < nDim; dim++) {
884#ifdef _OPENMP
885#pragma omp parallel for
886#endif
887 for(MInt k = 0; k < noCellsKP1; k++) {
888 for(MInt j = 0; j < noCellsJP1; j++) {
889 for(MInt i = 0; i < noCellsIP1; i++) {
890 const MInt I = cellIndex(i, j, k);
891 const MInt IP1 = I + IJK[dim];
892 dss[dim][I] = sqrt(POW2(x[IP1] - x[I]) + POW2(y[IP1] - y[I]) + POW2(z[IP1] - z[I]));
893 }
894 }
895 }
896 }
897
898 m_dsIsComputed = true;
899 }
900
901 for(MInt dim = 0; dim < nDim; dim++) {
902#ifdef _OPENMP
903#pragma omp parallel for
904#endif
905 for(MInt k = 1; k < noCellsK; k++) {
906 for(MInt j = 1; j < noCellsJ; j++) {
907#if defined(MAIA_INTEL_COMPILER)
908#pragma ivdep
909#pragma vector always
910#endif
911 for(MInt i = 1; i < noCellsI; i++) {
912 const MInt I = cellIndex(i, j, k);
913 const MInt IP1 = I + IJK[dim];
914 const MInt IM1 = I - IJK[dim];
915 const MInt IP2 = I + 2 * IJK[dim];
916
917 const MFloat DS = dss[dim][I];
918 const MFloat DSM1 = dss[dim][IM1];
919 const MFloat DSP1 = dss[dim][IP1];
920
921 const MFloat DSP = DS / POW2(DSP1 + DS);
922 const MFloat DSM = DS / POW2(DSM1 + DS);
923
924 // unrolled the loop so the compiler
925 // can optimize better
926 const MFloat DQU = vars[PV->U][IP1] - vars[PV->U][I];
927 const MFloat DQPU = vars[PV->U][IP2] - vars[PV->U][IP1];
928 const MFloat DQMU = vars[PV->U][I] - vars[PV->U][IM1];
929 MFloat UL = vars[PV->U][I] + DSM * (DSM1 * DQU + DS * DQMU);
930 MFloat UR = vars[PV->U][IP1] - DSP * (DS * DQPU + DSP1 * DQU);
931
932 const MFloat DQV = vars[PV->V][IP1] - vars[PV->V][I];
933 const MFloat DQPV = vars[PV->V][IP2] - vars[PV->V][IP1];
934 const MFloat DQMV = vars[PV->V][I] - vars[PV->V][IM1];
935 MFloat VL = vars[PV->V][I] + DSM * (DSM1 * DQV + DS * DQMV);
936 MFloat VR = vars[PV->V][IP1] - DSP * (DS * DQPV + DSP1 * DQV);
937
938 const MFloat DQW = vars[PV->W][IP1] - vars[PV->W][I];
939 const MFloat DQPW = vars[PV->W][IP2] - vars[PV->W][IP1];
940 const MFloat DQMW = vars[PV->W][I] - vars[PV->W][IM1];
941 MFloat WL = vars[PV->W][I] + DSM * (DSM1 * DQW + DS * DQMW);
942 MFloat WR = vars[PV->W][IP1] - DSP * (DS * DQPW + DSP1 * DQW);
943
944 const MFloat DQP = vars[PV->P][IP1] - vars[PV->P][I];
945 const MFloat DQPP = vars[PV->P][IP2] - vars[PV->P][IP1];
946 const MFloat DQMP = vars[PV->P][I] - vars[PV->P][IM1];
947 const MFloat PL = vars[PV->P][I] + DSM * (DSM1 * DQP + DS * DQMP);
948 const MFloat PR = vars[PV->P][IP1] - DSP * (DS * DQPP + DSP1 * DQP);
949
950 const MFloat DQRHO = vars[PV->RHO][IP1] - vars[PV->RHO][I];
951 const MFloat DQPRHO = vars[PV->RHO][IP2] - vars[PV->RHO][IP1];
952 const MFloat DQMRHO = vars[PV->RHO][I] - vars[PV->RHO][IM1];
953 const MFloat RHOL = vars[PV->RHO][I] + DSM * (DSM1 * DQRHO + DS * DQMRHO);
954 const MFloat RHOR = vars[PV->RHO][IP1] - DSP * (DS * DQPRHO + DSP1 * DQRHO);
955
956 const MFloat DQNUTILDE = vars[PV->RANS_VAR[0]][IP1] - vars[PV->RANS_VAR[0]][I];
957 const MFloat DQPNUTILDE = vars[PV->RANS_VAR[0]][IP2] - vars[PV->RANS_VAR[0]][IP1];
958 const MFloat DQMNUTILDE = vars[PV->RANS_VAR[0]][I] - vars[PV->RANS_VAR[0]][IM1];
959 const MFloat NUTILDEL = vars[PV->RANS_VAR[0]][I] + DSM * (DSM1 * DQNUTILDE + DS * DQMNUTILDE);
960 const MFloat NUTILDER = vars[PV->RANS_VAR[0]][IP1] - DSP * (DS * DQPNUTILDE + DSP1 * DQNUTILDE);
961 const MFloat surf0 = m_cells->surfaceMetrics[dim * 3 + 0][I];
962 const MFloat surf1 = m_cells->surfaceMetrics[dim * 3 + 1][I];
963 const MFloat surf2 = m_cells->surfaceMetrics[dim * 3 + 2][I];
964 const MFloat dxdtau = m_cells->dxt[dim][I];
965
966 const MFloat FRHOL = F1 / RHOL;
967 const MFloat FRHOR = F1 / RHOR;
968
969 const MFloat PLfRHOL = PL / RHOL;
970 const MFloat PRfRHOR = PR / RHOR;
971 const MFloat e0 = PLfRHOL * FgammaMinusOne + 0.5 * (POW2(UL) + POW2(VL) + POW2(WL)) + PLfRHOL;
972 const MFloat e1 = PRfRHOR * FgammaMinusOne + 0.5 * (POW2(UR) + POW2(VR) + POW2(WR)) + PRfRHOR;
973
974
975 // compute lenght of metric vector for normalization
976 const MFloat DGRAD = sqrt(POW2(surf0) + POW2(surf1) + POW2(surf2));
977 const MFloat FDGRAD = F1 / DGRAD;
978
979 // scale by metric length to get velocity in the new basis (get normalized basis vectors)
980 const MFloat UUL = ((UL * surf0 + VL * surf1 + WL * surf2) - dxdtau) * FDGRAD;
981
982
983 const MFloat UUR = ((UR * surf0 + VR * surf1 + WR * surf2) - dxdtau) * FDGRAD;
984
985 MFloat AL = FRHOL * PL;
986 MFloat AR = FRHOR * PR;
987
988 const MFloat FALR = 2.0 / (AL + AR);
989 const MFloat ALPHAL = AL * FALR;
990 const MFloat ALPHAR = AR * FALR;
991
992 AL = sqrt(gamma * AL);
993 AR = sqrt(gamma * AR);
994 AL = mMax(AL, AR);
995 AR = AL;
996
997 const MFloat XMAL = UUL / AL;
998 const MFloat XMAR = UUR / AR;
999
1000 AL = AL * DGRAD;
1001 AR = AR * DGRAD;
1002
1003 const MFloat RHOAL = AL * RHOL;
1004 const MFloat RHOAR = AR * RHOR;
1005
1006 const MFloat FDV = 0.3;
1007 const MFloat DXDXEZ = m_cells->coordinates[0][IP1] - m_cells->coordinates[0][I];
1008 const MFloat DYDXEZ = m_cells->coordinates[1][IP1] - m_cells->coordinates[1][I];
1009 const MFloat DZDXEZ = m_cells->coordinates[2][IP1] - m_cells->coordinates[2][I];
1010 MFloat SV = 2.0 * DGRAD / (m_cells->cellJac[I] + m_cells->cellJac[IP1]) * (FDV + (F1 - FDV) * getPSI(I, dim));
1011 const MFloat SV1 = F1 * SV * DXDXEZ;
1012 const MFloat SV2 = F1 * SV * DYDXEZ;
1013 const MFloat SV3 = F1 * SV * DZDXEZ;
1014
1015 const MFloat XMAL1 = mMin(F1, mMax(-F1, XMAL));
1016 const MFloat XMAR1 = mMin(F1, mMax(-F1, XMAR));
1017
1018 MFloat FXMA = F1B2 * (XMAL1 + fabs(XMAL1));
1019 const MFloat XMALP = ALPHAL * (F1B4 * POW2(XMAL1 + F1) - FXMA) + FXMA + (mMax(F1, XMAL) - F1);
1020 FXMA = F1B2 * (XMAR1 - fabs(XMAR1));
1021 const MFloat XMARM = ALPHAR * (-F1B4 * POW2(XMAR1 - F1) - FXMA) + FXMA + (mMin(-F1, XMAR) + F1);
1022
1023 const MFloat FLP = PL * ((F2 - XMAL1) * POW2(F1 + XMAL1));
1024 const MFloat FRP = PR * ((F2 + XMAR1) * POW2(F1 - XMAR1));
1025 const MFloat PLR = F1B4 * (FLP + FRP);
1026
1027 const MFloat RHOUL = XMALP * RHOAL;
1028 const MFloat RHOUR = XMARM * RHOAR;
1029 const MFloat RHOU = RHOUL + RHOUR;
1030 const MFloat RHOU2 = F1B2 * RHOU;
1031 const MFloat ARHOU2 = fabs(RHOU2);
1032
1033 const MFloat UUL2 = SV1 * UUL;
1034 const MFloat UUR2 = SV1 * UUR;
1035 UL = UL - UUL2;
1036 UR = UR - UUR2;
1037 const MFloat UUL3 = SV2 * UUL;
1038 const MFloat UUR3 = SV2 * UUR;
1039 VL = VL - UUL3;
1040 VR = VR - UUR3;
1041 const MFloat UUL4 = SV3 * UUL;
1042 const MFloat UUR4 = SV3 * UUR;
1043 WL = WL - UUL4;
1044 WR = WR - UUR4;
1045
1046 flux[CV->RHO_U][I] = RHOU2 * (UL + UR) + ARHOU2 * (UL - UR) + PLR * surf0 + RHOUL * UUL2 + RHOUR * UUR2;
1047 flux[CV->RHO_V][I] = RHOU2 * (VL + VR) + ARHOU2 * (VL - VR) + PLR * surf1 + RHOUL * UUL3 + RHOUR * UUR3;
1048 flux[CV->RHO_W][I] = RHOU2 * (WL + WR) + ARHOU2 * (WL - WR) + PLR * surf2 + RHOUL * UUL4 + RHOUR * UUR4;
1049 flux[CV->RHO_E][I] = RHOU2 * (e0 + e1) + ARHOU2 * (e0 - e1) + PLR * dxdtau;
1050 flux[CV->RHO][I] = RHOU;
1051 flux[CV->RANS_VAR[0]][I] = RHOU2 * (NUTILDEL + NUTILDER) + ARHOU2 * (NUTILDEL - NUTILDER);
1052 }
1053 }
1054 }
1055
1056 // FLUX BALANCE
1057 for(MUint v = 0; v < noVars; v++) {
1058#ifdef _OPENMP
1059#pragma omp parallel for
1060#endif
1061 for(MInt k = m_noGhostLayers; k < m_nCells[0] - m_noGhostLayers; k++) {
1062 for(MInt j = m_noGhostLayers; j < m_nCells[1] - m_noGhostLayers; j++) {
1063#if defined(MAIA_INTEL_COMPILER)
1064#pragma ivdep
1065#pragma vector always
1066#endif
1067 for(MInt i = m_noGhostLayers; i < m_nCells[2] - m_noGhostLayers; i++) {
1068 const MInt I = cellIndex(i, j, k);
1069 const MInt IM1 = I - IJK[dim];
1070 const MUint offset = v * noCells;
1071 MFloat* const RESTRICT rhs = ALIGNED_F(cellRhs + offset);
1072 rhs[I] += flux[v][IM1] - flux[v][I];
1073 }
1074 }
1075 }
1076 }
1077 }
1078}
1079
1080template void FvStructuredSolver3DRans::Muscl_AusmDV<5>();
1081template void FvStructuredSolver3DRans::Muscl_AusmDV<6>();
1082template void FvStructuredSolver3DRans::Muscl_AusmDV<7>();
1083
1084
1085template <MInt noVars>
1086void FvStructuredSolver3DRans::Muscl_AusmDV_Limited() {
1087 TRACE();
1088
1089 // stencil identifier
1090 const MInt IJK[3] = {1, m_nCells[2], m_nCells[1] * m_nCells[2]};
1091
1092 const MFloat* const RESTRICT x = ALIGNED_F(m_cells->coordinates[0]);
1093 const MFloat* const RESTRICT y = ALIGNED_F(m_cells->coordinates[1]);
1094 const MFloat* const RESTRICT z = ALIGNED_F(m_cells->coordinates[2]);
1095 const MFloat* const* const RESTRICT vars = ALIGNED_F(m_cells->pvariables);
1096 MFloat* const* const RESTRICT ql = ALIGNED_F(m_cells->ql);
1097 MFloat* const* const RESTRICT qr = ALIGNED_F(m_cells->qr);
1098 MFloat* const* const RESTRICT dss = ALIGNED_F(m_cells->dss);
1099
1100 MFloat* const* const RESTRICT flux = ALIGNED_F(m_cells->flux);
1101 MFloat* const RESTRICT cellRhs = ALIGNED_MF(m_cells->rightHandSide[0]);
1102
1103 const MFloat EPSLIM = 1e-14;
1104 const MFloat EPSS = 1e-34;
1105 const MUint noCells = m_noCells;
1106 const MFloat gamma = m_solver->m_gamma;
1107 const MFloat gammaMinusOne = gamma - 1.0;
1108 const MFloat FgammaMinusOne = F1 / gammaMinusOne;
1109
1110 const MInt noCellsI = m_nCells[2] - 2;
1111 const MInt noCellsJ = m_nCells[1] - 2;
1112 const MInt noCellsK = m_nCells[0] - 2;
1113
1114 const MInt noCellsIP1 = m_nCells[2] - 1;
1115 const MInt noCellsJP1 = m_nCells[1] - 1;
1116 const MInt noCellsKP1 = m_nCells[0] - 1;
1117
1118 if(!m_dsIsComputed) {
1119 for(MInt dim = 0; dim < nDim; dim++) {
1120#ifdef _OPENMP
1121#pragma omp parallel for
1122#endif
1123 for(MInt k = 0; k < noCellsKP1; k++) {
1124 for(MInt j = 0; j < noCellsJP1; j++) {
1125 for(MInt i = 0; i < noCellsIP1; i++) {
1126 const MInt I = cellIndex(i, j, k);
1127 const MInt IP1 = I + IJK[dim];
1128 dss[dim][I] = sqrt(POW2(x[IP1] - x[I]) + POW2(y[IP1] - y[I]) + POW2(z[IP1] - z[I]));
1129 }
1130 }
1131 }
1132 }
1133
1134 m_dsIsComputed = true;
1135 }
1136
1137 for(MInt dim = 0; dim < nDim; dim++) {
1138#ifdef _OPENMP
1139#pragma omp parallel for
1140#endif
1141 for(MInt k = 1; k < noCellsKP1; k++) {
1142 for(MInt j = 1; j < noCellsJP1; j++) {
1143 for(MInt i = 1; i < noCellsIP1; i++) {
1144 const MInt I = cellIndex(i, j, k);
1145 const MInt IP1 = I + IJK[dim];
1146 const MInt IM1 = I - IJK[dim];
1147
1148 const MFloat DS = dss[dim][I];
1149 const MFloat DSL = dss[dim][IM1];
1150 const MFloat DSR = dss[dim][I];
1151 const MFloat FDS = DS / POW2(DSL + DSR) * 2.0;
1152
1153 // unrolled the loop so the compiler
1154 // can optimize better
1155 const MFloat DQLU = DSL * (vars[PV->U][IP1] - vars[PV->U][I]);
1156 const MFloat DQRU = DSR * (vars[PV->U][I] - vars[PV->U][IM1]);
1157 const MFloat PHIU = mMax(F0, DQLU * DQRU / (POW2(DQLU) + POW2(DQRU) + EPSLIM * mMax(EPSS, DSL * DSR)));
1158 ql[PV->U][I] = vars[PV->U][I] + FDS * (DQLU + DQRU) * PHIU;
1159 qr[PV->U][IM1] = vars[PV->U][I] - FDS * (DQLU + DQRU) * PHIU;
1160
1161 const MFloat DQLV = DSL * (vars[PV->V][IP1] - vars[PV->V][I]);
1162 const MFloat DQRV = DSR * (vars[PV->V][I] - vars[PV->V][IM1]);
1163 const MFloat PHIV = mMax(F0, DQLV * DQRV / (POW2(DQLV) + POW2(DQRV) + EPSLIM * mMax(EPSS, DSL * DSR)));
1164 ql[PV->V][I] = vars[PV->V][I] + FDS * (DQLV + DQRV) * PHIV;
1165 qr[PV->V][IM1] = vars[PV->V][I] - FDS * (DQLV + DQRV) * PHIV;
1166
1167 const MFloat DQLW = DSL * (vars[PV->W][IP1] - vars[PV->W][I]);
1168 const MFloat DQRW = DSR * (vars[PV->W][I] - vars[PV->W][IM1]);
1169 const MFloat PHIW = mMax(F0, DQLW * DQRW / (POW2(DQLW) + POW2(DQRW) + EPSLIM * mMax(EPSS, DSL * DSR)));
1170 ql[PV->W][I] = vars[PV->W][I] + FDS * (DQLW + DQRW) * PHIW;
1171 qr[PV->W][IM1] = vars[PV->W][I] - FDS * (DQLW + DQRW) * PHIW;
1172
1173 const MFloat DQLP = DSL * (vars[PV->P][IP1] - vars[PV->P][I]);
1174 const MFloat DQRP = DSR * (vars[PV->P][I] - vars[PV->P][IM1]);
1175 const MFloat PHIP = mMax(F0, DQLP * DQRP / (POW2(DQLP) + POW2(DQRP) + EPSLIM * mMax(EPSS, DSL * DSR)));
1176 ql[PV->P][I] = vars[PV->P][I] + FDS * (DQLP + DQRP) * PHIP;
1177 qr[PV->P][IM1] = vars[PV->P][I] - FDS * (DQLP + DQRP) * PHIP;
1178
1179 const MFloat DQLRHO = DSL * (vars[PV->RHO][IP1] - vars[PV->RHO][I]);
1180 const MFloat DQRRHO = DSR * (vars[PV->RHO][I] - vars[PV->RHO][IM1]);
1181 const MFloat PHIRHO =
1182 mMax(F0, DQLRHO * DQRRHO / (POW2(DQLRHO) + POW2(DQRRHO) + EPSLIM * mMax(EPSS, DSL * DSR)));
1183 ql[PV->RHO][I] = vars[PV->RHO][I] + FDS * (DQLRHO + DQRRHO) * PHIRHO;
1184 qr[PV->RHO][IM1] = vars[PV->RHO][I] - FDS * (DQLRHO + DQRRHO) * PHIRHO;
1185
1186 const MFloat DQLNUTILDE = DSL * (vars[PV->RANS_VAR[0]][IP1] - vars[PV->RANS_VAR[0]][I]);
1187 const MFloat DQRNUTILDE = DSR * (vars[PV->RANS_VAR[0]][I] - vars[PV->RANS_VAR[0]][IM1]);
1188 const MFloat PHINUTILDE = mMax(
1189 F0, DQLNUTILDE * DQRNUTILDE / (POW2(DQLNUTILDE) + POW2(DQRNUTILDE) + EPSLIM * mMax(EPSS, DSL * DSR)));
1190 ql[PV->RANS_VAR[0]][I] = vars[PV->RANS_VAR[0]][I] + FDS * (DQLNUTILDE + DQRNUTILDE) * PHINUTILDE;
1191 qr[PV->RANS_VAR[0]][IM1] = vars[PV->RANS_VAR[0]][I] - FDS * (DQLNUTILDE + DQRNUTILDE) * PHINUTILDE;
1192 }
1193 }
1194 }
1195
1196
1197 for(MInt k = 1; k < noCellsK; k++) {
1198 for(MInt j = 1; j < noCellsJ; j++) {
1199 for(MInt i = 1; i < noCellsI; i++) {
1200 const MInt I = cellIndex(i, j, k);
1201 const MInt IP1 = I + IJK[dim];
1202
1203 MFloat UL = ql[PV->U][I];
1204 MFloat UR = qr[PV->U][I];
1205
1206 MFloat VL = ql[PV->V][I];
1207 MFloat VR = qr[PV->V][I];
1208
1209 MFloat WL = ql[PV->W][I];
1210 MFloat WR = qr[PV->W][I];
1211
1212 const MFloat PL = ql[PV->P][I];
1213 const MFloat PR = qr[PV->P][I];
1214
1215 const MFloat RHOL = ql[PV->RHO][I];
1216 const MFloat RHOR = qr[PV->RHO][I];
1217
1218 const MFloat NUTILDEL = ql[PV->RANS_VAR[0]][I];
1219 const MFloat NUTILDER = qr[PV->RANS_VAR[0]][I];
1220
1221 const MFloat surf0 = m_cells->surfaceMetrics[dim * 3 + 0][I];
1222 const MFloat surf1 = m_cells->surfaceMetrics[dim * 3 + 1][I];
1223 const MFloat surf2 = m_cells->surfaceMetrics[dim * 3 + 2][I];
1224 const MFloat dxdtau = m_cells->dxt[dim][I];
1225
1226 const MFloat FRHOL = F1 / RHOL;
1227 const MFloat FRHOR = F1 / RHOR;
1228
1229 const MFloat PLfRHOL = PL / RHOL;
1230 const MFloat PRfRHOR = PR / RHOR;
1231 const MFloat e0 = PLfRHOL * FgammaMinusOne + 0.5 * (POW2(UL) + POW2(VL) + POW2(WL)) + PLfRHOL;
1232 const MFloat e1 = PRfRHOR * FgammaMinusOne + 0.5 * (POW2(UR) + POW2(VR) + POW2(WR)) + PRfRHOR;
1233
1234
1235 // compute lenght of metric vector for normalization
1236 const MFloat DGRAD = sqrt(POW2(surf0) + POW2(surf1) + POW2(surf2));
1237 const MFloat FDGRAD = F1 / DGRAD;
1238
1239 // scale by metric length to get velocity in the new basis (get normalized basis vectors)
1240 const MFloat UUL = ((UL * surf0 + VL * surf1 + WL * surf2) - dxdtau) * FDGRAD;
1241
1242
1243 const MFloat UUR = ((UR * surf0 + VR * surf1 + WR * surf2) - dxdtau) * FDGRAD;
1244
1245 MFloat AL = FRHOL * PL;
1246 MFloat AR = FRHOR * PR;
1247
1248 const MFloat FALR = 2.0 / (AL + AR);
1249 const MFloat ALPHAL = AL * FALR;
1250 const MFloat ALPHAR = AR * FALR;
1251
1252 AL = sqrt(gamma * AL);
1253 AR = sqrt(gamma * AR);
1254 AL = mMax(AL, AR);
1255 AR = AL;
1256
1257 const MFloat XMAL = UUL / AL;
1258 const MFloat XMAR = UUR / AR;
1259
1260 AL = AL * DGRAD;
1261 AR = AR * DGRAD;
1262
1263 const MFloat RHOAL = AL * RHOL;
1264 const MFloat RHOAR = AR * RHOR;
1265
1266 const MFloat FDV = 0.3;
1267 const MFloat DXDXEZ = m_cells->coordinates[0][IP1] - m_cells->coordinates[0][I];
1268 const MFloat DYDXEZ = m_cells->coordinates[1][IP1] - m_cells->coordinates[1][I];
1269 const MFloat DZDXEZ = m_cells->coordinates[2][IP1] - m_cells->coordinates[2][I];
1270 MFloat SV = 2.0 * DGRAD / (m_cells->cellJac[I] + m_cells->cellJac[IP1]) * (FDV + (F1 - FDV) * getPSI(I, dim));
1271 const MFloat SV1 = F1 * SV * DXDXEZ;
1272 const MFloat SV2 = F1 * SV * DYDXEZ;
1273 const MFloat SV3 = F1 * SV * DZDXEZ;
1274
1275 const MFloat XMAL1 = mMin(F1, mMax(-F1, XMAL));
1276 const MFloat XMAR1 = mMin(F1, mMax(-F1, XMAR));
1277
1278 MFloat FXMA = F1B2 * (XMAL1 + fabs(XMAL1));
1279 const MFloat XMALP = ALPHAL * (F1B4 * POW2(XMAL1 + F1) - FXMA) + FXMA + (mMax(F1, XMAL) - F1);
1280 FXMA = F1B2 * (XMAR1 - fabs(XMAR1));
1281 const MFloat XMARM = ALPHAR * (-F1B4 * POW2(XMAR1 - F1) - FXMA) + FXMA + (mMin(-F1, XMAR) + F1);
1282
1283 const MFloat FLP = PL * ((F2 - XMAL1) * POW2(F1 + XMAL1));
1284 const MFloat FRP = PR * ((F2 + XMAR1) * POW2(F1 - XMAR1));
1285 const MFloat PLR = F1B4 * (FLP + FRP);
1286
1287 const MFloat RHOUL = XMALP * RHOAL;
1288 const MFloat RHOUR = XMARM * RHOAR;
1289 const MFloat RHOU = RHOUL + RHOUR;
1290 const MFloat RHOU2 = F1B2 * RHOU;
1291 const MFloat ARHOU2 = fabs(RHOU2);
1292
1293 const MFloat UUL2 = SV1 * UUL;
1294 const MFloat UUR2 = SV1 * UUR;
1295 UL = UL - UUL2;
1296 UR = UR - UUR2;
1297 const MFloat UUL3 = SV2 * UUL;
1298 const MFloat UUR3 = SV2 * UUR;
1299 VL = VL - UUL3;
1300 VR = VR - UUR3;
1301 const MFloat UUL4 = SV3 * UUL;
1302 const MFloat UUR4 = SV3 * UUR;
1303 WL = WL - UUL4;
1304 WR = WR - UUR4;
1305
1306 flux[CV->RHO_U][I] = RHOU2 * (UL + UR) + ARHOU2 * (UL - UR) + PLR * surf0 + RHOUL * UUL2 + RHOUR * UUR2;
1307 flux[CV->RHO_V][I] = RHOU2 * (VL + VR) + ARHOU2 * (VL - VR) + PLR * surf1 + RHOUL * UUL3 + RHOUR * UUR3;
1308 flux[CV->RHO_W][I] = RHOU2 * (WL + WR) + ARHOU2 * (WL - WR) + PLR * surf2 + RHOUL * UUL4 + RHOUR * UUR4;
1309 flux[CV->RHO_E][I] = RHOU2 * (e0 + e1) + ARHOU2 * (e0 - e1) + PLR * dxdtau;
1310 flux[CV->RHO][I] = RHOU;
1311 flux[CV->RANS_VAR[0]][I] = RHOU2 * (NUTILDEL + NUTILDER) + ARHOU2 * (NUTILDEL - NUTILDER);
1312 }
1313 }
1314 }
1315
1316 // FLUX BALANCE
1317 for(MUint v = 0; v < noVars; v++) {
1318#ifdef _OPENMP
1319#pragma omp parallel for
1320#endif
1321 for(MInt k = m_noGhostLayers; k < m_nCells[0] - m_noGhostLayers; k++) {
1322 for(MInt j = m_noGhostLayers; j < m_nCells[1] - m_noGhostLayers; j++) {
1323#if defined(MAIA_INTEL_COMPILER)
1324#pragma ivdep
1325#pragma vector always
1326#endif
1327 for(MInt i = m_noGhostLayers; i < m_nCells[2] - m_noGhostLayers; i++) {
1328 const MInt I = cellIndex(i, j, k);
1329 const MInt IM1 = I - IJK[dim];
1330 const MUint offset = v * noCells;
1331 MFloat* const RESTRICT rhs = ALIGNED_F(cellRhs + offset);
1332 rhs[I] += flux[v][IM1] - flux[v][I];
1333 }
1334 }
1335 }
1336 }
1337 }
1338}
1339
1340template void FvStructuredSolver3DRans::Muscl_AusmDV_Limited<5>();
1341template void FvStructuredSolver3DRans::Muscl_AusmDV_Limited<6>();
1342template void FvStructuredSolver3DRans::Muscl_AusmDV_Limited<7>();
1343
1344inline MInt FvStructuredSolver3DRans::cellIndex(MInt i, MInt j, MInt k) {
1345 return i + (j + k * m_nCells[1]) * m_nCells[2];
1346}
1347
1348inline MInt FvStructuredSolver3DRans::getCellIdfromCell(MInt origin, MInt incI, MInt incJ, MInt incK) {
1349 return origin + incI + incJ * m_nCells[2] + incK * m_nCells[2] * m_nCells[1];
1350}
1351
1352inline MFloat FvStructuredSolver3DRans::getPSI(MInt I, MInt dim) {
1353 const MFloat FK = 18.0;
1354 const MInt IJK[3] = {1, m_nCells[2], m_nCells[1] * m_nCells[2]};
1355 const MInt IP1 = I + IJK[dim];
1356 const MInt IM1 = I - IJK[dim];
1357 const MInt IP2 = I + 2 * IJK[dim];
1358
1359 const MFloat PIM2 = m_cells->pvariables[PV->P][IM1];
1360 const MFloat PIM1 = m_cells->pvariables[PV->P][I];
1361 const MFloat PIP2 = m_cells->pvariables[PV->P][IP2];
1362 const MFloat PIP1 = m_cells->pvariables[PV->P][IP1];
1363
1364 const MFloat PSI =
1365 mMin(F1,
1366 FK
1367 * mMax(mMax(fabs((PIM2 - PIM1) / mMin(PIM2, PIM1)), fabs((PIM1 - PIP1) / mMin(PIM1, PIP1))),
1368 fabs((PIP1 - PIP2) / mMin(PIP1, PIP2))));
1369 return PSI;
1370}
mAlloc
void mAlloc(T *&a, const MLong N, const MString &objectName, MString function)
allocates memory for one-dimensional array 'a' of size N
Definition: alloc.h:173
FvStructuredSolver3D
3D structured solver class
Definition: fvstructuredsolver3d.h:32
FvStructuredSolver3D::viscousFluxLES
void viscousFluxLES()
Viscous flux computation.
Definition: fvstructuredsolver3d.cpp:5020
FvStructuredSolver3DRans::m_solver
FvStructuredSolver3D * m_solver
Definition: fvstructuredsolver3drans.h:46
FvStructuredSolver3DRans::reconstructSurfaceData
void(FvStructuredSolver3DRans::* reconstructSurfaceData)()
Definition: fvstructuredsolver3drans.h:29
FvStructuredSolver3DRans::FQ
std::unique_ptr< StructuredFQVariables > & FQ
Definition: fvstructuredsolver3drans.h:56
FvStructuredSolver3DRans::FvStructuredSolver3DRans
FvStructuredSolver3DRans(FvStructuredSolver3D *solver)
Definition: fvstructuredsolver3drans.cpp:22
FvStructuredSolver3DRans::getCellIdfromCell
MInt getCellIdfromCell(MInt origin, MInt incI, MInt incJ, MInt incK)
Definition: fvstructuredsolver3drans.cpp:1348
FvStructuredSolver3DRans::cellIndex
MInt cellIndex(MInt i, MInt j, MInt k)
Definition: fvstructuredsolver3drans.cpp:1344
FvStructuredSolver3DRans::nDim
static constexpr const MInt nDim
Definition: fvstructuredsolver3drans.h:60
FvStructuredSolver3DRans::viscFluxMethod
void(FvStructuredSolver3DRans::* viscFluxMethod)()
Definition: fvstructuredsolver3drans.h:35
FvStructuredSolver3DRans::compTurbVisc
void(FvStructuredSolver3DRans::* compTurbVisc)()
Definition: fvstructuredsolver3drans.h:36
FvStructuredSolver3DRans::CV
std::unique_ptr< MConservativeVariables< 3 > > & CV
Definition: fvstructuredsolver3drans.h:54
FvStructuredSolver3DRans::PV
std::unique_ptr< MPrimitiveVariables< 3 > > & PV
Definition: fvstructuredsolver3drans.h:55
FvStructuredSolver::CV
std::unique_ptr< MConservativeVariables< nDim > > CV
Definition: fvstructuredsolver.h:332
ScratchSpace
This class is a ScratchSpace.
Definition: scratch.h:758
ScratchSpace::fill
void fill(T val)
fill the scratch with a given value
Definition: scratch.h:311
StructuredCell::cornerMetrics
MFloat ** cornerMetrics
Definition: fvstructuredcell.h:43
StructuredCell::surfaceMetrics
MFloat ** surfaceMetrics
Definition: fvstructuredcell.h:44
StructuredCell::rightHandSide
MFloat ** rightHandSide
Definition: fvstructuredcell.h:22
StructuredCell::coordinates
MFloat ** coordinates
Definition: fvstructuredcell.h:15
StructuredCell::temperature
MFloat * temperature
Definition: fvstructuredcell.h:20
StructuredCell::pvariables
MFloat ** pvariables
Definition: fvstructuredcell.h:18
StructuredCell::cellMetrics
MFloat ** cellMetrics
Definition: fvstructuredcell.h:42
string2enum
MInt string2enum(MString theString)
This global function translates strings in their corresponding enum values (integer values)....
Definition: enums.cpp:20
RansMethod
RansMethod
Definition: enums.h:51
RANS_SA_DV
@ RANS_SA_DV
Definition: enums.h:54
RANS_FS
@ RANS_FS
Definition: enums.h:55
mTerm
void mTerm(const MInt errorCode, const MString &location, const MString &message)
Definition: functions.cpp:29
POW2
constexpr Real POW2(const Real x)
Definition: functions.h:119
mMin
constexpr T mMin(const T &x, const T &y)
Definition: functions.h:90
mMax
constexpr T mMax(const T &x, const T &y)
Definition: functions.h:94
MInt
int32_t MInt
Definition: maiatypes.h:62
MUint
uint32_t MUint
Definition: maiatypes.h:63
MFloat
double MFloat
Definition: maiatypes.h:52
RANSModelConstants< RANS_FS >::faalpha
static constexpr MFloat faalpha
Definition: fvransmodelconstants.h:48
RANSModelConstants< RANS_FS >::fapsik2
static constexpr MFloat fapsik2
Definition: fvransmodelconstants.h:44
RANSModelConstants< RANS_FS >::fabetcs
static constexpr MFloat fabetcs
Definition: fvransmodelconstants.h:41
RANSModelConstants< RANS_FS >::fapsik1
static constexpr MFloat fapsik1
Definition: fvransmodelconstants.h:43
RANSModelConstants< RANS_FS >::fasigma
static constexpr MFloat fasigma
Definition: fvransmodelconstants.h:47
RANSModelConstants< RANS_FS >::facv1to3
static constexpr MFloat facv1to3
Definition: fvransmodelconstants.h:50
RANSModelConstants< RANS_FS >::fabetc
static constexpr MFloat fabetc
Definition: fvransmodelconstants.h:42
RANSModelConstants< RANS_SA_DV >::cw1
static constexpr MFloat cw1
Definition: fvransmodelconstants.h:28
RANSModelConstants< RANS_SA_DV >::Fsigma
static constexpr MFloat Fsigma
Definition: fvransmodelconstants.h:24
RANSModelConstants< RANS_SA_DV >::cv1to3
static constexpr MFloat cv1to3
Definition: fvransmodelconstants.h:31
RANSModelConstants< RANS_SA_DV >::cw3to6
static constexpr MFloat cw3to6
Definition: fvransmodelconstants.h:30
RANSModelConstants< RANS_SA_DV >::Ft2
static constexpr MFloat Ft2
Definition: fvransmodelconstants.h:32
RANSModelConstants< RANS_SA_DV >::cb2
static constexpr MFloat cb2
Definition: fvransmodelconstants.h:25
RANSModelConstants< RANS_SA_DV >::Fkap2
static constexpr MFloat Fkap2
Definition: fvransmodelconstants.h:26
RANSModelConstants< RANS_SA_DV >::cw2
static constexpr MFloat cw2
Definition: fvransmodelconstants.h:29
RANSModelConstants< RANS_SA_DV >::cb1
static constexpr MFloat cb1
Definition: fvransmodelconstants.h:27
y
define array structures
Definition: geometrycontexttypes.h:14