MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
fvstructuredsolver2drans.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 "fvstructuredsolver2drans.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
21class FvStructuredSolver2D;
22
23FvStructuredSolver2DRans::FvStructuredSolver2DRans(FvStructuredSolver2D* solver)
24 : m_StructuredComm(solver->m_StructuredComm),
25 m_solverId(solver->m_solverId),
26 m_nCells(solver->m_nCells),
27 m_nPoints(solver->m_nPoints),
28 m_noCells(solver->m_noCells),
29 m_cells(solver->m_cells),
30 CV(solver->CV),
31 PV(solver->PV),
32 FQ(solver->FQ),
33 m_noGhostLayers(solver->m_noGhostLayers),
34 m_eps(solver->m_eps),
35 m_chi(solver->m_chi),
36 m_sutherlandConstant(solver->m_sutherlandConstant),
37 m_sutherlandPlusOne(solver->m_sutherlandPlusOne) {
38 m_solver = solver;
39 // set all the methods for the 2d RANS code here
40
41 initFluxMethod();
42}
43
44FvStructuredSolver2DRans::~FvStructuredSolver2DRans() {}
45
46void FvStructuredSolver2DRans::initFluxMethod() {
47 // set pointer to the right AUSM for the RANS equations;
48
61 m_ransMethod =
62 static_cast<RansMethod>(string2enum(Context::getSolverProperty<MString>("ransMethod", m_solverId, AT_)));
63
64 m_dsIsComputed = false;
65 switch(m_ransMethod) {
66 case RANS_SA:
67 case RANS_SA_DV: {
68 mAlloc(m_cells->saFlux1, nDim, m_noCells, "m_cells->saFlux1", -999999.9, AT_);
69 mAlloc(m_cells->saFlux2, nDim, m_noCells, "m_cells->saFlux2", -999999.9, AT_);
70 mAlloc(m_cells->prodDest, m_noCells, "m_cells->prodDest", -999999.9, AT_);
71 compTurbVisc = &FvStructuredSolver2DRans::computeTurbViscosity_SA;
72 viscFluxMethod = &FvStructuredSolver2DRans::viscousFlux_SA;
73 switch(m_solver->CV->noVariables) {
74 case 5: {
75 reconstructSurfaceData = &FvStructuredSolver2DRans::Muscl_AusmDV<5, 1>;
76 break;
77 }
78 case 6: {
79 reconstructSurfaceData = &FvStructuredSolver2DRans::Muscl_AusmDV<6, 1>;
80 break;
81 }
82 case 7: {
83 reconstructSurfaceData = &FvStructuredSolver2DRans::Muscl_AusmDV<7, 1>;
84 break;
85 }
86 default: {
87 stringstream errorMessage;
88 errorMessage << "Number of Variables " << m_solver->CV->noVariables
89 << " not implemented! in temlate Rans AUSM " << endl;
90 mTerm(1, AT_, errorMessage.str());
91 }
92 }
93 break;
94 }
95 case RANS_KEPSILON: {
96 setAndAllocate_KEPSILON();
97
98 compTurbVisc = &FvStructuredSolver2DRans::computeTurbViscosity_KEPSILON;
99 viscFluxMethod = &FvStructuredSolver2DRans::viscousFlux_KEPSILON2;
100
101 // TODO_SS labels:FV
102 MBool AusmDV = true;
103 AusmDV = Context::getSolverProperty<MBool>("AusmDV", m_solverId, AT_, &AusmDV);
104
105 if(m_solver->m_limiter) {
106 mAlloc(m_cells->ql, PV->noVariables, m_noCells, "m_cells->ql", F0, AT_);
107 mAlloc(m_cells->qr, PV->noVariables, m_noCells, "m_cells->qr", F0, AT_);
108 m_log << "Using RANS with K-Epsilon model and limited AusmDV" << endl;
109 switch(m_solver->CV->noVariables) {
110 case 6: {
111 if(m_solver->m_rans2eq_mode == "production")
112 reconstructSurfaceData = (AusmDV) ? &FvStructuredSolver2DRans::Muscl_AusmDV_Limited<6, 2, true>
113 : &FvStructuredSolver2DRans::Muscl_Ausm_Limited<6, 2, true>;
114 else
115 reconstructSurfaceData = (AusmDV) ? &FvStructuredSolver2DRans::Muscl_AusmDV_Limited<6, 2>
116 : &FvStructuredSolver2DRans::Muscl_Ausm_Limited<6, 2>;
117 break;
118 }
119 case 7: {
120 if(m_solver->m_rans2eq_mode == "production")
121 reconstructSurfaceData = (AusmDV) ? &FvStructuredSolver2DRans::Muscl_AusmDV_Limited<7, 2, true>
122 : &FvStructuredSolver2DRans::Muscl_Ausm_Limited<7, 2, true>;
123 else
124 reconstructSurfaceData = (AusmDV) ? &FvStructuredSolver2DRans::Muscl_AusmDV_Limited<7, 2>
125 : &FvStructuredSolver2DRans::Muscl_Ausm_Limited<7, 2>;
126 break;
127 }
128 default: {
129 stringstream errorMessage;
130 errorMessage << "Number of Variables " << m_solver->CV->noVariables
131 << " not implemented! in temlate Rans AUSM " << endl;
132 mTerm(1, AT_, errorMessage.str());
133 }
134 }
135 } else {
136 if(!AusmDV) mTerm(1, "Not implemented yet!");
137 m_log << "Using RANS with K-Epsilon model and AusmDV" << endl;
138 switch(m_solver->CV->noVariables) {
139 case 6: {
140 if(m_solver->m_rans2eq_mode == "production")
141 reconstructSurfaceData = &FvStructuredSolver2DRans::Muscl_AusmDV<6, 2, true>;
142 else
143 reconstructSurfaceData = &FvStructuredSolver2DRans::Muscl_AusmDV<6, 2>;
144 break;
145 }
146 case 7: {
147 if(m_solver->m_rans2eq_mode == "production")
148 reconstructSurfaceData = &FvStructuredSolver2DRans::Muscl_AusmDV<7, 2, true>;
149 else
150 reconstructSurfaceData = &FvStructuredSolver2DRans::Muscl_AusmDV<7, 2>;
151 break;
152 }
153 default: {
154 stringstream errorMessage;
155 errorMessage << "Number of Variables " << m_solver->CV->noVariables
156 << " not implemented! in temlate Rans AUSM " << endl;
157 mTerm(1, AT_, errorMessage.str());
158 }
159 }
160 }
161 break;
162 }
163 default: {
164 mTerm(1, AT_, "RANS METHOD wsa not specified in properties");
165 break;
166 }
167 }
168}
169
170void FvStructuredSolver2DRans::Muscl(MInt) { (this->*reconstructSurfaceData)(); }
171
172template <MInt noVars, MInt noRANS, MBool rans2eq_production_mode>
173void FvStructuredSolver2DRans::Muscl_AusmDV() {
174 TRACE();
175
176 // stencil identifier
177 const MInt IJ[2] = {1, m_nCells[1]};
178
179 const MFloat* const RESTRICT x = ALIGNED_F(m_cells->coordinates[0]);
180 const MFloat* const RESTRICT y = ALIGNED_F(m_cells->coordinates[1]);
181 const MFloat* const* const RESTRICT vars = ALIGNED_F(m_cells->pvariables);
182 MFloat* const* const RESTRICT dss = ALIGNED_F(m_cells->dss);
183 MFloat* const* const RESTRICT flux = ALIGNED_F(m_cells->flux);
184 MFloat* const RESTRICT cellRhs = ALIGNED_MF(m_cells->rightHandSide[0]);
185
186 const MUint noCells = m_noCells;
187 const MFloat gamma = m_solver->m_gamma;
188 const MFloat gammaMinusOne = gamma - 1.0;
189 const MFloat FgammaMinusOne = F1 / gammaMinusOne;
190
191 const MInt noCellsI = m_nCells[1] - 2;
192 const MInt noCellsJ = m_nCells[0] - 2;
193
194 const MInt noCellsIP1 = m_nCells[1] - 1;
195 const MInt noCellsJP1 = m_nCells[0] - 1;
196
197 if(!m_dsIsComputed) {
198 for(MInt dim = 0; dim < nDim; dim++) {
199#ifdef _OPENMP
200#pragma omp parallel for
201#endif
202 for(MInt j = 0; j < noCellsJP1; j++) {
203 for(MInt i = 0; i < noCellsIP1; i++) {
204 const MInt I = cellIndex(i, j);
205 const MInt IP1 = I + IJ[dim];
206 dss[dim][I] = sqrt(POW2(x[IP1] - x[I]) + POW2(y[IP1] - y[I]));
207 }
208 }
209 }
210
211 m_dsIsComputed = true;
212 }
213
214
215 for(MInt dim = 0; dim < nDim; dim++) {
216 for(MInt j = 1; j < noCellsJ; j++) {
217 for(MInt i = 1; i < noCellsI; i++) {
218 const MInt I = cellIndex(i, j);
219 const MInt IP1 = I + IJ[dim];
220 const MInt IM1 = I - IJ[dim];
221 const MInt IP2 = I + 2 * IJ[dim];
222
223 const MFloat DS = dss[dim][I];
224 const MFloat DSM1 = dss[dim][IM1];
225 const MFloat DSP1 = dss[dim][IP1];
226
227 const MFloat DSP = DS / POW2(DSP1 + DS);
228 const MFloat DSM = DS / POW2(DSM1 + DS);
229
230 // unrolled the loop so the compiler
231 // can optimize better
232 const MFloat DQU = vars[PV->U][IP1] - vars[PV->U][I];
233 const MFloat DQPU = vars[PV->U][IP2] - vars[PV->U][IP1];
234 const MFloat DQMU = vars[PV->U][I] - vars[PV->U][IM1];
235 MFloat UL = vars[PV->U][I] + DSM * (DSM1 * DQU + DS * DQMU);
236 MFloat UR = vars[PV->U][IP1] - DSP * (DS * DQPU + DSP1 * DQU);
237
238 const MFloat DQV = vars[PV->V][IP1] - vars[PV->V][I];
239 const MFloat DQPV = vars[PV->V][IP2] - vars[PV->V][IP1];
240 const MFloat DQMV = vars[PV->V][I] - vars[PV->V][IM1];
241 MFloat VL = vars[PV->V][I] + DSM * (DSM1 * DQV + DS * DQMV);
242 MFloat VR = vars[PV->V][IP1] - DSP * (DS * DQPV + DSP1 * DQV);
243
244 const MFloat DQP = vars[PV->P][IP1] - vars[PV->P][I];
245 const MFloat DQPP = vars[PV->P][IP2] - vars[PV->P][IP1];
246 const MFloat DQMP = vars[PV->P][I] - vars[PV->P][IM1];
247 const MFloat PL = vars[PV->P][I] + DSM * (DSM1 * DQP + DS * DQMP);
248 const MFloat PR = vars[PV->P][IP1] - DSP * (DS * DQPP + DSP1 * DQP);
249
250 const MFloat DQRHO = vars[PV->RHO][IP1] - vars[PV->RHO][I];
251 const MFloat DQPRHO = vars[PV->RHO][IP2] - vars[PV->RHO][IP1];
252 const MFloat DQMRHO = vars[PV->RHO][I] - vars[PV->RHO][IM1];
253 const MFloat RHOL = vars[PV->RHO][I] + DSM * (DSM1 * DQRHO + DS * DQMRHO);
254 const MFloat RHOR = vars[PV->RHO][IP1] - DSP * (DS * DQPRHO + DSP1 * DQRHO);
255
256 // in C++17 and later use "if constexpr"
257 // for SA-model: nutilde
258 // for k-epsilon-model: k and epsilon
259 MFloat RANSL[noRANS];
260 MFloat RANSR[noRANS];
261 for(MInt v = 0; v < noRANS; ++v) {
262 const MFloat DQRANSVAR = vars[PV->RANS_VAR[v]][IP1] - vars[PV->RANS_VAR[v]][I];
263 const MFloat DQPRANSVAR = vars[PV->RANS_VAR[v]][IP2] - vars[PV->RANS_VAR[v]][IP1];
264 const MFloat DQMRANSVAR = vars[PV->RANS_VAR[v]][I] - vars[PV->RANS_VAR[v]][IM1];
265 RANSL[v] = vars[PV->RANS_VAR[v]][I] + DSM * (DSM1 * DQRANSVAR + DS * DQMRANSVAR);
266 RANSR[v] = vars[PV->RANS_VAR[v]][IP1] - DSP * (DS * DQPRANSVAR + DSP1 * DQRANSVAR);
267 }
268 // const MFloat DQNUTILDE = vars[PV->RANS_VAR[0]][IP1]-vars[PV->RANS_VAR[0]][I];
269 // const MFloat DQPNUTILDE = vars[PV->RANS_VAR[0]][IP2]-vars[PV->RANS_VAR[0]][IP1];
270 // const MFloat DQMNUTILDE = vars[PV->RANS_VAR[0]][I]-vars[PV->RANS_VAR[0]][IM1];
271 // const MFloat NUTILDEL = vars[PV->RANS_VAR[0]][I]+DSM*(DSM1*DQNUTILDE+DS*DQMNUTILDE);
272 // const MFloat NUTILDER = vars[PV->RANS_VAR[0]][IP1]-DSP*(DS*DQPNUTILDE+DSP1*DQNUTILDE);
273
274 const MFloat surf0 = m_cells->surfaceMetrics[dim * 2 + 0][I];
275 const MFloat surf1 = m_cells->surfaceMetrics[dim * 2 + 1][I];
276
277 const MFloat FRHOL = F1 / RHOL;
278 const MFloat FRHOR = F1 / RHOR;
279
280 const MFloat PLfRHOL = PL / RHOL;
281 const MFloat PRfRHOR = PR / RHOR;
282 const MFloat e0 = (rans2eq_production_mode)
283 ? PLfRHOL * FgammaMinusOne + 0.5 * (POW2(UL) + POW2(VL)) + RANSL[0] + PLfRHOL
284 : PLfRHOL * FgammaMinusOne + 0.5 * (POW2(UL) + POW2(VL)) + PLfRHOL;
285 const MFloat e1 = (rans2eq_production_mode)
286 ? PRfRHOR * FgammaMinusOne + 0.5 * (POW2(UR) + POW2(VR)) + RANSR[0] + PRfRHOR
287 : PRfRHOR * FgammaMinusOne + 0.5 * (POW2(UR) + POW2(VR)) + PRfRHOR;
288
289
290 // compute lenght of metric vector for normalization
291 const MFloat DGRAD = sqrt(POW2(surf0) + POW2(surf1));
292 const MFloat FDGRAD = F1 / DGRAD;
293
294 // scale by metric length to get velocity in the new basis (get normalized basis vectors)
295 const MFloat UUL = ((UL * surf0 + VL * surf1)) * FDGRAD;
296 const MFloat UUR = ((UR * surf0 + VR * surf1)) * FDGRAD;
297
298 MFloat AL = FRHOL * PL;
299 MFloat AR = FRHOR * PR;
300
301 const MFloat FALR = 2.0 / (AL + AR);
302 const MFloat ALPHAL = AL * FALR;
303 const MFloat ALPHAR = AR * FALR;
304
305 AL = sqrt(gamma * AL);
306 AR = sqrt(gamma * AR);
307 AL = mMax(AL, AR);
308 AR = AL;
309
310 const MFloat XMAL = UUL / AL;
311 const MFloat XMAR = UUR / AR;
312
313 AL = AL * DGRAD;
314 AR = AR * DGRAD;
315
316 const MFloat RHOAL = AL * RHOL;
317 const MFloat RHOAR = AR * RHOR;
318
319 const MFloat FDV = 0.3;
320 const MFloat DXDXEZ = m_cells->coordinates[0][IP1] - m_cells->coordinates[0][I];
321 const MFloat DYDXEZ = m_cells->coordinates[1][IP1] - m_cells->coordinates[1][I];
322 MFloat SV = 2.0 * DGRAD / (m_cells->cellJac[I] + m_cells->cellJac[IP1]) * (FDV + (F1 - FDV) * getPSI(I, dim));
323 const MFloat SV1 = F1 * SV * DXDXEZ;
324 const MFloat SV2 = F1 * SV * DYDXEZ;
325
326 const MFloat XMAL1 = mMin(F1, mMax(-F1, XMAL));
327 const MFloat XMAR1 = mMin(F1, mMax(-F1, XMAR));
328
329 MFloat FXMA = F1B2 * (XMAL1 + fabs(XMAL1));
330 const MFloat XMALP = ALPHAL * (F1B4 * POW2(XMAL1 + F1) - FXMA) + FXMA + (mMax(F1, XMAL) - F1);
331 FXMA = F1B2 * (XMAR1 - fabs(XMAR1));
332 const MFloat XMARM = ALPHAR * (-F1B4 * POW2(XMAR1 - F1) - FXMA) + FXMA + (mMin(-F1, XMAR) + F1);
333
334 const MFloat FLP = PL * ((F2 - XMAL1) * POW2(F1 + XMAL1));
335 const MFloat FRP = PR * ((F2 + XMAR1) * POW2(F1 - XMAR1));
336 const MFloat PLR = F1B4 * (FLP + FRP);
337
338 const MFloat RHOUL = XMALP * RHOAL;
339 const MFloat RHOUR = XMARM * RHOAR;
340 const MFloat RHOU = RHOUL + RHOUR;
341 const MFloat RHOU2 = F1B2 * RHOU;
342 const MFloat ARHOU2 = fabs(RHOU2);
343
344 const MFloat UUL2 = SV1 * UUL;
345 const MFloat UUR2 = SV1 * UUR;
346 UL = UL - UUL2;
347 UR = UR - UUR2;
348 const MFloat UUL3 = SV2 * UUL;
349 const MFloat UUR3 = SV2 * UUR;
350 VL = VL - UUL3;
351 VR = VR - UUR3;
352
353 flux[CV->RHO_U][I] = RHOU2 * (UL + UR) + ARHOU2 * (UL - UR) + PLR * surf0 + RHOUL * UUL2 + RHOUR * UUR2;
354 flux[CV->RHO_V][I] = RHOU2 * (VL + VR) + ARHOU2 * (VL - VR) + PLR * surf1 + RHOUL * UUL3 + RHOUR * UUR3;
355 flux[CV->RHO_E][I] = RHOU2 * (e0 + e1) + ARHOU2 * (e0 - e1);
356 flux[CV->RHO][I] = RHOU;
357 for(MInt v = 0; v < noRANS; ++v) {
358 flux[CV->RANS_VAR[v]][I] = RHOU2 * (RANSL[v] + RANSR[v]) + ARHOU2 * (RANSL[v] - RANSR[v]);
359 // flux[I+noCells*CV->RANS_VAR[0]] = RHOU2*(NUTILDEL+NUTILDER) + ARHOU2*(NUTILDEL-NUTILDER);
360 }
361 }
362 }
363
364 // FLUX BALANCE
365 for(MUint v = 0; v < noVars; v++) {
366 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; j++) {
367 for(MInt i = m_noGhostLayers; i < m_nCells[1] - m_noGhostLayers; i++) {
368 const MInt I = cellIndex(i, j);
369 const MInt IM1 = I - IJ[dim];
370 const MUint offset = v * noCells;
371 MFloat* const RESTRICT rhs = ALIGNED_F(cellRhs + offset);
372 rhs[I] += flux[v][IM1] - flux[v][I];
373 }
374 }
375 }
376 }
377}
378template void FvStructuredSolver2DRans::Muscl_AusmDV<5, 1>();
379template void FvStructuredSolver2DRans::Muscl_AusmDV<6, 1>();
380template void FvStructuredSolver2DRans::Muscl_AusmDV<7, 1>();
381template void FvStructuredSolver2DRans::Muscl_AusmDV<6, 2>();
382template void FvStructuredSolver2DRans::Muscl_AusmDV<7, 2>();
383
384
385template <MInt noVars, MInt noRANS, MBool rans2eq_production_mode>
386void FvStructuredSolver2DRans::Muscl_AusmDV_Limited() {
387 TRACE();
388
389 // stencil identifier
390 const MInt IJ[2] = {1, m_nCells[1]};
391
392 const MFloat* const RESTRICT x = ALIGNED_F(m_cells->coordinates[0]);
393 const MFloat* const RESTRICT y = ALIGNED_F(m_cells->coordinates[1]);
394 const MFloat* const* const RESTRICT vars = ALIGNED_F(m_cells->pvariables);
395 MFloat* const* const RESTRICT ql = ALIGNED_F(m_cells->ql);
396 MFloat* const* const RESTRICT qr = ALIGNED_F(m_cells->qr);
397 MFloat* const* const RESTRICT dss = ALIGNED_F(m_cells->dss);
398
399 MFloat* const* const RESTRICT flux = ALIGNED_F(m_cells->flux);
400 MFloat* const RESTRICT cellRhs = ALIGNED_MF(m_cells->rightHandSide[0]);
401
402 const MFloat EPSLIM = 1e-14;
403 const MFloat EPSS = 1e-34;
404 const MUint noCells = m_noCells;
405 const MFloat gamma = m_solver->m_gamma;
406 const MFloat gammaMinusOne = gamma - 1.0;
407 const MFloat FgammaMinusOne = F1 / gammaMinusOne;
408
409 const MInt noCellsI = m_nCells[1] - 2;
410 const MInt noCellsJ = m_nCells[0] - 2;
411
412 const MInt noCellsIP1 = m_nCells[1] - 1;
413 const MInt noCellsJP1 = m_nCells[0] - 1;
414
415 if(!m_dsIsComputed) {
416 for(MInt dim = 0; dim < nDim; dim++) {
417#ifdef _OPENMP
418#pragma omp parallel for
419#endif
420 for(MInt j = 0; j < noCellsJP1; j++) {
421 for(MInt i = 0; i < noCellsIP1; i++) {
422 const MInt I = cellIndex(i, j);
423 const MInt IP1 = I + IJ[dim];
424 dss[dim][I] = sqrt(POW2(x[IP1] - x[I]) + POW2(y[IP1] - y[I]));
425 }
426 }
427 }
428
429 m_dsIsComputed = true;
430 }
431
432 for(MInt dim = 0; dim < nDim; dim++) {
433 for(MInt j = 1; j < noCellsJP1; j++) {
434 for(MInt i = 1; i < noCellsIP1; i++) {
435 const MInt I = cellIndex(i, j);
436 const MInt IP1 = I + IJ[dim];
437 const MInt IM1 = I - IJ[dim];
438
439 const MFloat DS = dss[dim][I];
440 const MFloat DSL = dss[dim][IM1];
441 const MFloat DSR = dss[dim][I];
442 const MFloat FDS = DS / POW2(DSL + DSR) * 2.0;
443
444 // unrolled the loop so the compiler
445 // can optimize better
446 const MFloat DQLU = DSL * (vars[PV->U][IP1] - vars[PV->U][I]);
447 const MFloat DQRU = DSR * (vars[PV->U][I] - vars[PV->U][IM1]);
448 const MFloat PHIU = mMax(F0, DQLU * DQRU / (POW2(DQLU) + POW2(DQRU) + EPSLIM * mMax(EPSS, DSL * DSR)));
449 ql[PV->U][I] = vars[PV->U][I] + FDS * (DQLU + DQRU) * PHIU;
450 qr[PV->U][IM1] = vars[PV->U][I] - FDS * (DQLU + DQRU) * PHIU;
451
452 const MFloat DQLV = DSL * (vars[PV->V][IP1] - vars[PV->V][I]);
453 const MFloat DQRV = DSR * (vars[PV->V][I] - vars[PV->V][IM1]);
454 const MFloat PHIV = mMax(F0, DQLV * DQRV / (POW2(DQLV) + POW2(DQRV) + EPSLIM * mMax(EPSS, DSL * DSR)));
455 ql[PV->V][I] = vars[PV->V][I] + FDS * (DQLV + DQRV) * PHIV;
456 qr[PV->V][IM1] = vars[PV->V][I] - FDS * (DQLV + DQRV) * PHIV;
457
458 const MFloat DQLP = DSL * (vars[PV->P][IP1] - vars[PV->P][I]);
459 const MFloat DQRP = DSR * (vars[PV->P][I] - vars[PV->P][IM1]);
460 const MFloat PHIP = mMax(F0, DQLP * DQRP / (POW2(DQLP) + POW2(DQRP) + EPSLIM * mMax(EPSS, DSL * DSR)));
461 ql[PV->P][I] = vars[PV->P][I] + FDS * (DQLP + DQRP) * PHIP;
462 qr[PV->P][IM1] = vars[PV->P][I] - FDS * (DQLP + DQRP) * PHIP;
463
464 const MFloat DQLRHO = DSL * (vars[PV->RHO][IP1] - vars[PV->RHO][I]);
465 const MFloat DQRRHO = DSR * (vars[PV->RHO][I] - vars[PV->RHO][IM1]);
466 const MFloat PHIRHO =
467 mMax(F0, DQLRHO * DQRRHO / (POW2(DQLRHO) + POW2(DQRRHO) + EPSLIM * mMax(EPSS, DSL * DSR)));
468 ql[PV->RHO][I] = vars[PV->RHO][I] + FDS * (DQLRHO + DQRRHO) * PHIRHO;
469 qr[PV->RHO][IM1] = vars[PV->RHO][I] - FDS * (DQLRHO + DQRRHO) * PHIRHO;
470
471 for(MInt v = 0; v < noRANS; ++v) {
472 const MFloat DQLRANSVAR = DSL * (vars[PV->RANS_VAR[v]][IP1] - vars[PV->RANS_VAR[v]][I]);
473 const MFloat DQRRANSVAR = DSR * (vars[PV->RANS_VAR[v]][I] - vars[PV->RANS_VAR[v]][IM1]);
474 const MFloat PHIRANSVAR = mMax(
475 F0, DQLRANSVAR * DQRRANSVAR / (POW2(DQLRANSVAR) + POW2(DQRRANSVAR) + EPSLIM * mMax(EPSS, DSL * DSR)));
476 ql[PV->RANS_VAR[v]][I] = vars[PV->RANS_VAR[v]][I] + FDS * (DQLRANSVAR + DQRRANSVAR) * PHIRANSVAR;
477 qr[PV->RANS_VAR[v]][IM1] = vars[PV->RANS_VAR[v]][I] - FDS * (DQLRANSVAR + DQRRANSVAR) * PHIRANSVAR;
478 }
479 }
480 }
481
482
483 for(MInt j = 1; j < noCellsJ; j++) {
484 for(MInt i = 1; i < noCellsI; i++) {
485 const MInt I = cellIndex(i, j);
486 const MInt IP1 = I + IJ[dim];
487
488 MFloat UL = ql[PV->U][I];
489 MFloat UR = qr[PV->U][I];
490
491 MFloat VL = ql[PV->V][I];
492 MFloat VR = qr[PV->V][I];
493
494 const MFloat PL = ql[PV->P][I];
495 const MFloat PR = qr[PV->P][I];
496
497 const MFloat RHOL = ql[PV->RHO][I];
498 const MFloat RHOR = qr[PV->RHO][I];
499
500 MFloat RANSL[noRANS];
501 MFloat RANSR[noRANS];
502 for(MInt v = 0; v < noRANS; ++v) {
503 RANSL[v] = ql[PV->RANS_VAR[v]][I];
504 RANSR[v] = qr[PV->RANS_VAR[v]][I];
505 }
506
507 const MFloat surf0 = m_cells->surfaceMetrics[dim * 2 + 0][I];
508 const MFloat surf1 = m_cells->surfaceMetrics[dim * 2 + 1][I];
509
510 const MFloat FRHOL = F1 / RHOL;
511 const MFloat FRHOR = F1 / RHOR;
512
513 const MFloat PLfRHOL = PL / RHOL;
514 const MFloat PRfRHOR = PR / RHOR;
515 const MFloat e0 = (rans2eq_production_mode)
516 ? PLfRHOL * FgammaMinusOne + 0.5 * (POW2(UL) + POW2(VL)) + RANSL[0] + PLfRHOL
517 : PLfRHOL * FgammaMinusOne + 0.5 * (POW2(UL) + POW2(VL)) + PLfRHOL;
518 const MFloat e1 = (rans2eq_production_mode)
519 ? PRfRHOR * FgammaMinusOne + 0.5 * (POW2(UR) + POW2(VR)) + RANSR[0] + PRfRHOR
520 : PRfRHOR * FgammaMinusOne + 0.5 * (POW2(UR) + POW2(VR)) + PRfRHOR;
521
522
523 // compute lenght of metric vector for normalization
524 const MFloat DGRAD = sqrt(POW2(surf0) + POW2(surf1));
525 const MFloat FDGRAD = F1 / DGRAD;
526
527 // scale by metric length to get velocity in the new basis (get normalized basis vectors)
528 const MFloat UUL = ((UL * surf0 + VL * surf1)) * FDGRAD;
529 const MFloat UUR = ((UR * surf0 + VR * surf1)) * FDGRAD;
530
531 MFloat AL = FRHOL * PL;
532 MFloat AR = FRHOR * PR;
533
534 const MFloat FALR = 2.0 / (AL + AR);
535 const MFloat ALPHAL = AL * FALR;
536 const MFloat ALPHAR = AR * FALR;
537
538 AL = sqrt(gamma * AL);
539 AR = sqrt(gamma * AR);
540 AL = mMax(AL, AR);
541 AR = AL;
542
543 const MFloat XMAL = UUL / AL;
544 const MFloat XMAR = UUR / AR;
545
546 AL = AL * DGRAD;
547 AR = AR * DGRAD;
548
549 const MFloat RHOAL = AL * RHOL;
550 const MFloat RHOAR = AR * RHOR;
551
552 const MFloat FDV = 0.3;
553 const MFloat DXDXEZ = m_cells->coordinates[0][IP1] - m_cells->coordinates[0][I];
554 const MFloat DYDXEZ = m_cells->coordinates[1][IP1] - m_cells->coordinates[1][I];
555 MFloat SV = 2.0 * DGRAD / (m_cells->cellJac[I] + m_cells->cellJac[IP1]) * (FDV + (F1 - FDV) * getPSI(I, dim));
556 const MFloat SV1 = F1 * SV * DXDXEZ;
557 const MFloat SV2 = F1 * SV * DYDXEZ;
558
559 const MFloat XMAL1 = mMin(F1, mMax(-F1, XMAL));
560 const MFloat XMAR1 = mMin(F1, mMax(-F1, XMAR));
561
562 MFloat FXMA = F1B2 * (XMAL1 + fabs(XMAL1));
563 const MFloat XMALP = ALPHAL * (F1B4 * POW2(XMAL1 + F1) - FXMA) + FXMA + (mMax(F1, XMAL) - F1);
564 FXMA = F1B2 * (XMAR1 - fabs(XMAR1));
565 const MFloat XMARM = ALPHAR * (-F1B4 * POW2(XMAR1 - F1) - FXMA) + FXMA + (mMin(-F1, XMAR) + F1);
566
567 const MFloat FLP = PL * ((F2 - XMAL1) * POW2(F1 + XMAL1));
568 const MFloat FRP = PR * ((F2 + XMAR1) * POW2(F1 - XMAR1));
569 const MFloat PLR = F1B4 * (FLP + FRP);
570
571 const MFloat RHOUL = XMALP * RHOAL;
572 const MFloat RHOUR = XMARM * RHOAR;
573 const MFloat RHOU = RHOUL + RHOUR;
574 const MFloat RHOU2 = F1B2 * RHOU;
575 const MFloat ARHOU2 = fabs(RHOU2);
576
577 const MFloat UUL2 = SV1 * UUL;
578 const MFloat UUR2 = SV1 * UUR;
579 UL = UL - UUL2;
580 UR = UR - UUR2;
581 const MFloat UUL3 = SV2 * UUL;
582 const MFloat UUR3 = SV2 * UUR;
583 VL = VL - UUL3;
584 VR = VR - UUR3;
585
586 flux[CV->RHO_U][I] = RHOU2 * (UL + UR) + ARHOU2 * (UL - UR) + PLR * surf0 + RHOUL * UUL2 + RHOUR * UUR2;
587 flux[CV->RHO_V][I] = RHOU2 * (VL + VR) + ARHOU2 * (VL - VR) + PLR * surf1 + RHOUL * UUL3 + RHOUR * UUR3;
588 flux[CV->RHO_E][I] = RHOU2 * (e0 + e1) + ARHOU2 * (e0 - e1);
589 flux[CV->RHO][I] = RHOU;
590 for(MInt v = 0; v < noRANS; ++v) {
591 flux[CV->RANS_VAR[v]][I] = RHOU2 * (RANSL[v] + RANSR[v]) + ARHOU2 * (RANSL[v] - RANSR[v]);
592 }
593 }
594 }
595
596 // FLUX BALANCE
597 for(MUint v = 0; v < noVars; v++) {
598 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; j++) {
599 for(MInt i = m_noGhostLayers; i < m_nCells[1] - m_noGhostLayers; i++) {
600 const MInt I = cellIndex(i, j);
601 const MInt IM1 = I - IJ[dim];
602 const MUint offset = v * noCells;
603 MFloat* const RESTRICT rhs = ALIGNED_F(cellRhs + offset);
604 rhs[I] += flux[v][IM1] - flux[v][I];
605 }
606 }
607 }
608 }
609}
610template void FvStructuredSolver2DRans::Muscl_AusmDV_Limited<6, 2>();
611template void FvStructuredSolver2DRans::Muscl_AusmDV_Limited<7, 2>();
612
613
614template <MInt noVars, MInt noRANS, MBool rans2eq_production_mode>
615void FvStructuredSolver2DRans::Muscl_Ausm_Limited() {
616 TRACE();
617
618 // stencil identifier
619 const MInt IJ[2] = {1, m_nCells[1]};
620
621 const MFloat* const RESTRICT x = ALIGNED_F(m_cells->coordinates[0]);
622 const MFloat* const RESTRICT y = ALIGNED_F(m_cells->coordinates[1]);
623 const MFloat* const* const RESTRICT vars = ALIGNED_F(m_cells->pvariables);
624 MFloat* const* const RESTRICT ql = ALIGNED_F(m_cells->ql);
625 MFloat* const* const RESTRICT qr = ALIGNED_F(m_cells->qr);
626 MFloat* const* const RESTRICT dss = ALIGNED_F(m_cells->dss);
627
628 MFloat* const* const RESTRICT flux = ALIGNED_F(m_cells->flux);
629 MFloat* const RESTRICT cellRhs = ALIGNED_MF(m_cells->rightHandSide[0]);
630
631 const MFloat EPSLIM = 1e-14;
632 const MFloat EPSS = 1e-34;
633 const MUint noCells = m_noCells;
634 const MFloat gamma = m_solver->m_gamma;
635 const MFloat gammaMinusOne = gamma - 1.0;
636 const MFloat FgammaMinusOne = F1 / gammaMinusOne;
637
638 const MInt noCellsI = m_nCells[1] - 2;
639 const MInt noCellsJ = m_nCells[0] - 2;
640
641 const MInt noCellsIP1 = m_nCells[1] - 1;
642 const MInt noCellsJP1 = m_nCells[0] - 1;
643
644 if(!m_dsIsComputed) {
645 for(MInt dim = 0; dim < nDim; dim++) {
646#ifdef _OPENMP
647#pragma omp parallel for
648#endif
649 for(MInt j = 0; j < noCellsJP1; j++) {
650 for(MInt i = 0; i < noCellsIP1; i++) {
651 const MInt I = cellIndex(i, j);
652 const MInt IP1 = I + IJ[dim];
653 dss[dim][I] = sqrt(POW2(x[IP1] - x[I]) + POW2(y[IP1] - y[I]));
654 }
655 }
656 }
657
658 m_dsIsComputed = true;
659 }
660
661 for(MInt dim = 0; dim < nDim; dim++) {
662 for(MInt j = 1; j < noCellsJP1; j++) {
663 for(MInt i = 1; i < noCellsIP1; i++) {
664 const MInt I = cellIndex(i, j);
665 const MInt IP1 = I + IJ[dim];
666 const MInt IM1 = I - IJ[dim];
667
668 const MFloat DS = dss[dim][I];
669 const MFloat DSL = dss[dim][IM1];
670 const MFloat DSR = dss[dim][I];
671 const MFloat FDS = DS / POW2(DSL + DSR) * 2.0;
672
673 // unrolled the loop so the compiler
674 // can optimize better
675 const MFloat DQLU = DSL * (vars[PV->U][IP1] - vars[PV->U][I]);
676 const MFloat DQRU = DSR * (vars[PV->U][I] - vars[PV->U][IM1]);
677 const MFloat PHIU = mMax(F0, DQLU * DQRU / (POW2(DQLU) + POW2(DQRU) + EPSLIM * mMax(EPSS, DSL * DSR)));
678 ql[PV->U][I] = vars[PV->U][I] + FDS * (DQLU + DQRU) * PHIU;
679 qr[PV->U][IM1] = vars[PV->U][I] - FDS * (DQLU + DQRU) * PHIU;
680
681 const MFloat DQLV = DSL * (vars[PV->V][IP1] - vars[PV->V][I]);
682 const MFloat DQRV = DSR * (vars[PV->V][I] - vars[PV->V][IM1]);
683 const MFloat PHIV = mMax(F0, DQLV * DQRV / (POW2(DQLV) + POW2(DQRV) + EPSLIM * mMax(EPSS, DSL * DSR)));
684 ql[PV->V][I] = vars[PV->V][I] + FDS * (DQLV + DQRV) * PHIV;
685 qr[PV->V][IM1] = vars[PV->V][I] - FDS * (DQLV + DQRV) * PHIV;
686
687 const MFloat DQLP = DSL * (vars[PV->P][IP1] - vars[PV->P][I]);
688 const MFloat DQRP = DSR * (vars[PV->P][I] - vars[PV->P][IM1]);
689 const MFloat PHIP = mMax(F0, DQLP * DQRP / (POW2(DQLP) + POW2(DQRP) + EPSLIM * mMax(EPSS, DSL * DSR)));
690 ql[PV->P][I] = vars[PV->P][I] + FDS * (DQLP + DQRP) * PHIP;
691 qr[PV->P][IM1] = vars[PV->P][I] - FDS * (DQLP + DQRP) * PHIP;
692
693 const MFloat DQLRHO = DSL * (vars[PV->RHO][IP1] - vars[PV->RHO][I]);
694 const MFloat DQRRHO = DSR * (vars[PV->RHO][I] - vars[PV->RHO][IM1]);
695 const MFloat PHIRHO =
696 mMax(F0, DQLRHO * DQRRHO / (POW2(DQLRHO) + POW2(DQRRHO) + EPSLIM * mMax(EPSS, DSL * DSR)));
697 ql[PV->RHO][I] = vars[PV->RHO][I] + FDS * (DQLRHO + DQRRHO) * PHIRHO;
698 qr[PV->RHO][IM1] = vars[PV->RHO][I] - FDS * (DQLRHO + DQRRHO) * PHIRHO;
699
700 for(MInt v = 0; v < noRANS; ++v) {
701 const MFloat DQLRANSVAR = DSL * (vars[PV->RANS_VAR[v]][IP1] - vars[PV->RANS_VAR[v]][I]);
702 const MFloat DQRRANSVAR = DSR * (vars[PV->RANS_VAR[v]][I] - vars[PV->RANS_VAR[v]][IM1]);
703 const MFloat PHIRANSVAR = mMax(
704 F0, DQLRANSVAR * DQRRANSVAR / (POW2(DQLRANSVAR) + POW2(DQRRANSVAR) + EPSLIM * mMax(EPSS, DSL * DSR)));
705 ql[PV->RANS_VAR[v]][I] = vars[PV->RANS_VAR[v]][I] + FDS * (DQLRANSVAR + DQRRANSVAR) * PHIRANSVAR;
706 qr[PV->RANS_VAR[v]][IM1] = vars[PV->RANS_VAR[v]][I] - FDS * (DQLRANSVAR + DQRRANSVAR) * PHIRANSVAR;
707 }
708 }
709 }
710
711
712 for(MInt j = 1; j < noCellsJ; j++) {
713 for(MInt i = 1; i < noCellsI; i++) {
714 const MInt I = cellIndex(i, j);
715
716 MFloat UL = ql[PV->U][I];
717 MFloat UR = qr[PV->U][I];
718
719 MFloat VL = ql[PV->V][I];
720 MFloat VR = qr[PV->V][I];
721
722 const MFloat PL = ql[PV->P][I];
723 const MFloat PR = qr[PV->P][I];
724
725 const MFloat RHOL = ql[PV->RHO][I];
726 const MFloat RHOR = qr[PV->RHO][I];
727
728 MFloat RANSL[noRANS];
729 MFloat RANSR[noRANS];
730 for(MInt v = 0; v < noRANS; ++v) {
731 RANSL[v] = ql[PV->RANS_VAR[v]][I];
732 RANSR[v] = qr[PV->RANS_VAR[v]][I];
733 }
734
735 const MFloat surf0 = m_cells->surfaceMetrics[dim * 2 + 0][I];
736 const MFloat surf1 = m_cells->surfaceMetrics[dim * 2 + 1][I];
737
738 // compute lenght of metric vector for normalization
739 const MFloat DGRAD = sqrt(POW2(surf0) + POW2(surf1));
740 const MFloat FDGRAD = F1 / DGRAD;
741
742 // scale by metric length to get velocity in the new basis (get normalized basis vectors)
743 const MFloat UUL = ((UL * surf0 + VL * surf1)) * FDGRAD;
744
745
746 const MFloat UUR = ((UR * surf0 + VR * surf1)) * FDGRAD;
747
748 // speed of sound
749 const MFloat AL = sqrt(gamma * mMax(m_eps, PL / mMax(m_eps, RHOL)));
750 const MFloat AR = sqrt(gamma * mMax(m_eps, PR / mMax(m_eps, RHOR)));
751
752 const MFloat MAL = UUL / AL;
753 const MFloat MAR = UUR / AR;
754
755 const MFloat MALR = F1B2 * (MAL + MAR);
756 const MFloat PLR = F1B2 * (PL + PR);
757
758 const MFloat RHO_AL = RHOL * AL;
759 const MFloat RHO_AR = RHOR * AR;
760
761 const MFloat PLfRHOL = PL / RHOL;
762 const MFloat PRfRHOR = PR / RHOR;
763 const MFloat e0 = (rans2eq_production_mode)
764 ? PLfRHOL * FgammaMinusOne + 0.5 * (POW2(UL) + POW2(VL)) + RANSL[0] + PLfRHOL
765 : PLfRHOL * FgammaMinusOne + 0.5 * (POW2(UL) + POW2(VL)) + PLfRHOL;
766 const MFloat e1 = (rans2eq_production_mode)
767 ? PRfRHOR * FgammaMinusOne + 0.5 * (POW2(UR) + POW2(VR)) + RANSR[0] + PRfRHOR
768 : PRfRHOR * FgammaMinusOne + 0.5 * (POW2(UR) + POW2(VR)) + PRfRHOR;
769
770 const MFloat RHOU = F1B2 * (MALR * (RHO_AL + RHO_AR) + fabs(MALR) * (RHO_AL - RHO_AR)) * DGRAD;
771 const MFloat RHOU2 = F1B2 * RHOU;
772 const MFloat ARHOU2 = fabs(RHOU2);
773
774 flux[CV->RHO_U][I] = RHOU2 * (UL + UR) + ARHOU2 * (UL - UR) + PLR * surf0;
775 flux[CV->RHO_V][I] = RHOU2 * (VL + VR) + ARHOU2 * (VL - VR) + PLR * surf1;
776 flux[CV->RHO_E][I] = RHOU2 * (e0 + e1) + ARHOU2 * (e0 - e1);
777 flux[CV->RHO][I] = RHOU;
778 for(MInt v = 0; v < noRANS; ++v) {
779 flux[CV->RANS_VAR[v]][I] = RHOU2 * (RANSL[v] + RANSR[v]) + ARHOU2 * (RANSL[v] - RANSR[v]);
780 }
781 }
782 }
783
784 // FLUX BALANCE
785 for(MUint v = 0; v < noVars; v++) {
786 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; j++) {
787 for(MInt i = m_noGhostLayers; i < m_nCells[1] - m_noGhostLayers; i++) {
788 const MInt I = cellIndex(i, j);
789 const MInt IM1 = I - IJ[dim];
790 const MUint offset = v * noCells;
791 MFloat* const RESTRICT rhs = ALIGNED_F(cellRhs + offset);
792 rhs[I] += flux[v][IM1] - flux[v][I];
793 }
794 }
795 }
796 }
797}
798template void FvStructuredSolver2DRans::Muscl_Ausm_Limited<6, 2>();
799template void FvStructuredSolver2DRans::Muscl_Ausm_Limited<7, 2>();
800
801
802void FvStructuredSolver2DRans::viscousFluxRANS() { (this->*viscFluxMethod)(); }
803
804
805void FvStructuredSolver2DRans::viscousFlux_SA() {
806 computeTurbViscosity_SA();
807
808 // call the standard LES viscous flux
809 m_solver->viscousFluxLES<>();
810
811 // OTHER variables required to calculate the laminar viscous fluxes
812 const MFloat rRe = F1 / m_solver->m_Re0;
813
814 MFloat* const* const RESTRICT eflux = ALIGNED_F(m_cells->eFlux);
815 MFloat* const* const RESTRICT fflux = ALIGNED_F(m_cells->fFlux);
816 MFloat* const* const RESTRICT sa_1flux = ALIGNED_F(m_cells->saFlux1);
817 MFloat* const* const RESTRICT sa_2flux = ALIGNED_F(m_cells->saFlux2);
818
819 MFloat* const RESTRICT u = ALIGNED_F(m_cells->pvariables[PV->U]);
820 MFloat* const RESTRICT v = ALIGNED_F(m_cells->pvariables[PV->V]);
821 MFloat* const RESTRICT rho = ALIGNED_F(m_cells->pvariables[PV->RHO]);
822 MFloat* const RESTRICT T = ALIGNED_F(m_cells->temperature);
823 MFloat* const RESTRICT nuTilde = ALIGNED_F(m_cells->pvariables[PV->RANS_VAR[0]]);
824 MFloat* const RESTRICT muLam = ALIGNED_F(m_cells->fq[FQ->MU_L]);
825
826 for(MInt j = m_noGhostLayers - 1; j < m_nCells[0] - m_noGhostLayers + 1; j++) {
827 for(MInt i = m_noGhostLayers - 1; i < m_nCells[1] - m_noGhostLayers + 1; i++) {
828 // get the adjacent cells;
829
830 const MInt IJ = cellIndex(i, j);
831 const MInt IPJ = cellIndex((i + 1), j);
832 const MInt IPJP = cellIndex((i + 1), (j + 1));
833 const MInt IJP = cellIndex(i, (j + 1));
834 const MInt IMJ = cellIndex((i - 1), j);
835 const MInt IJM = cellIndex(i, (j - 1));
836
837 const MFloat cornerMetrics[9] = {m_cells->cornerMetrics[0][IJ], m_cells->cornerMetrics[1][IJ],
838 m_cells->cornerMetrics[2][IJ], m_cells->cornerMetrics[3][IJ]};
839
840 const MFloat dnutdxi = F1B2 * (nuTilde[IPJP] + nuTilde[IPJ] - nuTilde[IJP] - nuTilde[IJ]);
841 const MFloat dnutdet = F1B2 * (nuTilde[IPJP] + nuTilde[IJP] - nuTilde[IPJ] - nuTilde[IJ]);
842
843 const MFloat nutldAvg = F1B4 * (nuTilde[IJP] + nuTilde[IPJP] + nuTilde[IJ] + nuTilde[IPJ]);
844
845 const MFloat nuLamAvg =
846 F1B4 * (muLam[IJP] / rho[IJP] + muLam[IPJP] / rho[IPJP] + muLam[IJ] / rho[IJ] + muLam[IPJ] / rho[IPJ]);
847
848
849 const MFloat dnutldx =
850 dnutdxi * m_cells->cornerMetrics[xsd * 2 + xsd][IJ] + dnutdet * m_cells->cornerMetrics[ysd * 2 + xsd][IJ];
851
852 const MFloat dnutldy =
853 dnutdxi * m_cells->cornerMetrics[xsd * 2 + ysd][IJ] + dnutdet * m_cells->cornerMetrics[ysd * 2 + ysd][IJ];
854
855
856 const MFloat Frj = rRe / m_cells->cornerJac[IJ];
857
858 const MFloat sax1 = Frj * (nuLamAvg + (1.0 + RM_SA_DV::cb2) * nutldAvg)
859 * (dnutldx * cornerMetrics[xsd * 2 + xsd] + dnutldy * cornerMetrics[xsd * 2 + ysd]);
860 const MFloat sax2 =
861 -Frj * RM_SA_DV::cb2 * (dnutldx * cornerMetrics[xsd * 2 + xsd] + dnutldy * cornerMetrics[xsd * 2 + ysd]);
862 const MFloat say1 = Frj * (nuLamAvg + (1.0 + RM_SA_DV::cb2) * nutldAvg)
863 * (dnutldx * cornerMetrics[ysd * 2 + xsd] + dnutldy * cornerMetrics[ysd * 2 + ysd]);
864 const MFloat say2 =
865 -Frj * RM_SA_DV::cb2 * (dnutldx * cornerMetrics[ysd * 2 + xsd] + dnutldy * cornerMetrics[ysd * 2 + ysd]);
866
867
868 // compute vorticity
869 const MFloat du1 = u[IPJ] - u[IMJ];
870 const MFloat du2 = u[IJP] - u[IJM];
871 const MFloat dv1 = v[IPJ] - v[IMJ];
872 const MFloat dv2 = v[IJP] - v[IJM];
873 const MFloat vortk =
874 (m_cells->cellMetrics[xsd * 2 + xsd][IJ] * dv1) + (m_cells->cellMetrics[ysd * 2 + xsd][IJ] * dv2)
875 - (m_cells->cellMetrics[xsd * 2 + ysd][IJ] * du1) - (m_cells->cellMetrics[ysd * 2 + ysd][IJ] * du2);
876 const MFloat s = F1B2 * fabs(vortk) / m_cells->cellJac[IJ];
877
878
879 const MFloat distance = m_cells->fq[FQ->WALLDISTANCE][IJ];
880 const MFloat Fdist2 = 1.0 / (distance * distance);
881 const MFloat chi = nuTilde[IJ] * rho[IJ] / (SUTHERLANDLAW(T[IJ]));
882 const MFloat chip3 = chi * chi * chi;
883 const MFloat Fv1 = chip3 / (chip3 + RM_SA_DV::cv1to3);
884 const MFloat Fv2 = F1 - (chi / (F1 + chi * Fv1));
885
886 const MFloat term = nuTilde[IJ] * Fdist2 * RM_SA_DV::Fkap2;
887 const MFloat stilde = s + term * Fv2 * rRe;
888 const MFloat r = min(10.0, rRe * term / stilde);
889
890 const MFloat g = r + RM_SA_DV::cw2 * (pow(r, 6) - r);
891 const MFloat Fwterm = (1 + RM_SA_DV::cw3to6) / (pow(g, 6) + RM_SA_DV::cw3to6);
892 const MFloat Fw = g * pow(Fwterm, (1.0 / 6.0));
893 const MFloat prodValue = rho[IJ] * RM_SA_DV::cb1 * (F1 - RM_SA_DV::Ft2) * stilde * nuTilde[IJ];
894 const MFloat destValue = rRe * rho[IJ] * (RM_SA_DV::cw1 * Fw - RM_SA_DV::cb1 * RM_SA_DV::Fkap2 * RM_SA_DV::Ft2)
895 * pow(nuTilde[IJ], 2.0) * Fdist2;
896
897 m_cells->prodDest[IJ] = (prodValue - destValue) * m_cells->cellJac[IJ];
898
899 eflux[0][IJ] = sax1; // diffusion of nutilde for every cell
900 eflux[1][IJ] = sax2; // diffusion of nutilde for every cell
901
902 fflux[0][IJ] = say1; // diffusion of nutilde for every cell
903 fflux[1][IJ] = say2; // diffusion of nutilde for every cell
904 }
905 }
906
910 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; ++j) {
911 for(MInt i = m_noGhostLayers - 1; i < m_nCells[1] - m_noGhostLayers; ++i) {
912 const MInt IJ = cellIndex(i, j);
913 const MInt IJM = cellIndex(i, (j - 1));
914
915 sa_1flux[0][IJ] = F1B2 * (eflux[0][IJ] + eflux[0][IJM]);
916 sa_2flux[0][IJ] = F1B2 * (eflux[1][IJ] + eflux[1][IJM]);
917 }
918 }
919
920 for(MInt j = m_noGhostLayers - 1; j < m_nCells[0] - m_noGhostLayers; ++j) {
921 for(MInt i = m_noGhostLayers; i < m_nCells[1] - m_noGhostLayers; ++i) {
922 const MInt IJ = cellIndex(i, j);
923 const MInt IMJ = cellIndex((i - 1), j);
924
925 sa_1flux[1][IJ] = F1B2 * (fflux[0][IJ] + fflux[0][IMJ]);
926 sa_2flux[1][IJ] = F1B2 * (fflux[1][IJ] + fflux[1][IMJ]);
927 }
928 }
929
930
931 // separate loop for adding the prodn nad destrn terms for tur kin viscosity transport variable
932 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; j++) {
933 for(MInt i = m_noGhostLayers; i < m_nCells[1] - m_noGhostLayers; i++) {
934 const MInt IJ = cellIndex(i, j);
935 const MInt IMJ = cellIndex((i - 1), j);
936 const MInt IJM = cellIndex(i, (j - 1));
937
938 const MFloat dissipation_term =
939 (((sa_1flux[0][IJ] - sa_1flux[0][IMJ]) + ((sa_2flux[0][IJ] - sa_2flux[0][IMJ]) * nuTilde[IJ])) * rho[IJ]
940 * RM_SA_DV::Fsigma
941 + ((sa_1flux[1][IJ] - sa_1flux[1][IJM]) + ((sa_2flux[1][IJ] - sa_2flux[1][IJM]) * nuTilde[IJ])) * rho[IJ]
942 * RM_SA_DV::Fsigma);
943
944 m_cells->rightHandSide[CV->RANS_VAR[0]][IJ] += dissipation_term;
945 m_cells->rightHandSide[CV->RANS_VAR[0]][IJ] += m_cells->prodDest[IJ];
946 }
947 }
948}
949
950void FvStructuredSolver2DRans::viscousFlux_KEPSILON() {
951 ASSERT(!m_solver->m_hasSingularity, "Get your ass out of this fucntion!!!");
952
953 // compute friction velocity at wall
954 m_solver->m_structuredBndryCnd->computeFrictionCoef();
955 // communicate wall properties to all cells
956 m_solver->m_structuredBndryCnd->distributeWallAndFPProperties();
957
958 computeTurbViscosity_KEPSILON();
959
960 // call the standard LES viscous flux
961 if(m_solver->m_viscCompact)
962 m_solver->viscousFluxLESCompact<true>();
963 else
964 m_solver->viscousFluxLES<true>();
965
966 // OTHER variables required to calculate the laminar viscous fluxes
967 static constexpr MFloat minMFloat =
968 1e-16; // std::min(std::numeric_limits<MFloat>::min(), 1.0/std::numeric_limits<MFloat>::max());
969 const MFloat rRe0 = F1 / m_solver->m_Re0;
970 const MFloat gamma = m_solver->m_gamma;
971 const MFloat fac = (m_solver->m_rans2eq_mode == "production") ? 1.0 : 0.0;
972 const MFloat fac_nonDim = m_solver->m_keps_nonDimType ? 1.0 : PV->UInfinity;
973
974 MFloat* const* const RESTRICT eflux = ALIGNED_F(m_cells->eFlux);
975 MFloat* const* const RESTRICT fflux = ALIGNED_F(m_cells->fFlux);
976 MFloat* const* const RESTRICT sa_1flux = ALIGNED_F(m_cells->saFlux1);
977 MFloat* const* const RESTRICT sa_2flux = ALIGNED_F(m_cells->saFlux2);
978
979 MFloat* const RESTRICT u = ALIGNED_F(m_cells->pvariables[PV->U]);
980 MFloat* const RESTRICT v = ALIGNED_F(m_cells->pvariables[PV->V]);
981 MFloat* const RESTRICT rho = ALIGNED_F(m_cells->pvariables[PV->RHO]);
982 MFloat* const RESTRICT p = ALIGNED_F(m_cells->pvariables[PV->P]);
983 // MFloat* const RESTRICT T = ALIGNED_F(m_cells->temperature);
984 MFloat* const RESTRICT TKE = ALIGNED_F(m_cells->pvariables[PV->RANS_VAR[0]]);
985 MFloat* const RESTRICT eps = ALIGNED_F(m_cells->pvariables[PV->RANS_VAR[1]]);
986 MFloat* const RESTRICT muLam = ALIGNED_F(m_cells->fq[FQ->MU_L]);
987 MFloat* const RESTRICT muTurb = ALIGNED_F(m_cells->fq[FQ->MU_T]);
988 const MFloat* const RESTRICT utau = &m_cells->fq[FQ->UTAU][0];
989 const MFloat* const RESTRICT wallDist = &m_cells->fq[FQ->WALLDISTANCE][0];
990
991 for(MInt j = m_noGhostLayers - 1; j < m_nCells[0] - m_noGhostLayers; j++) {
992 for(MInt i = m_noGhostLayers - 1; i < m_nCells[1] - m_noGhostLayers; i++) {
993 // get the adjacent cells;
994
995 const MInt IJ = cellIndex(i, j);
996 const MInt IPJ = cellIndex((i + 1), j);
997 const MInt IPJP = cellIndex((i + 1), (j + 1));
998 const MInt IJP = cellIndex(i, (j + 1));
999
1000 const MFloat cornerMetrics[nDim * nDim] = {m_cells->cornerMetrics[0][IJ], m_cells->cornerMetrics[1][IJ],
1001 m_cells->cornerMetrics[2][IJ], m_cells->cornerMetrics[3][IJ]};
1002
1003 // TODO_SS labels:FV if wall is along eta=constant, then at wall dTKEdxi and depsdxi must be zero! Check!!!
1004
1005 const MFloat dTKEdxi = F1B2 * (TKE[IPJP] + TKE[IPJ] - TKE[IJP] - TKE[IJ]);
1006 const MFloat dTKEdet = F1B2 * (TKE[IPJP] + TKE[IJP] - TKE[IPJ] - TKE[IJ]);
1007
1008 const MFloat dTKEdx = dTKEdxi * cornerMetrics[xsd * 2 + xsd] + dTKEdet * cornerMetrics[ysd * 2 + xsd];
1009
1010 const MFloat dTKEdy = dTKEdxi * cornerMetrics[xsd * 2 + ysd] + dTKEdet * cornerMetrics[ysd * 2 + ysd];
1011
1012
1013 const MFloat depsdxi = F1B2 * (eps[IPJP] + eps[IPJ] - eps[IJP] - eps[IJ]);
1014 const MFloat depsdet = F1B2 * (eps[IPJP] + eps[IJP] - eps[IPJ] - eps[IJ]);
1015
1016 const MFloat depsdx = depsdxi * cornerMetrics[xsd * 2 + xsd] + depsdet * cornerMetrics[ysd * 2 + xsd];
1017
1018 const MFloat depsdy = depsdxi * cornerMetrics[xsd * 2 + ysd] + depsdet * cornerMetrics[ysd * 2 + ysd];
1019
1020 // TODO_SS labels:FV muTurbAvg at wall surface must be zero! Check!
1021 const MFloat muLamAvg = F1B4 * (muLam[IJP] + muLam[IPJP] + muLam[IJ] + muLam[IPJ]);
1022 const MFloat muTurbAvg = F1B4 * (muTurb[IJP] + muTurb[IPJP] + muTurb[IJ] + muTurb[IPJ]);
1023
1024 const MFloat Frj = rRe0 / m_cells->cornerJac[IJ];
1025 const MFloat mu_k = muLamAvg + muTurbAvg * RM_KEPS::rsigma_k;
1026 const MFloat mu_eps = muLamAvg + muTurbAvg * RM_KEPS::rsigma_eps;
1027
1028 const MFloat sax1 = Frj * mu_k * (dTKEdx * cornerMetrics[xsd * 2 + xsd] + dTKEdy * cornerMetrics[xsd * 2 + ysd]);
1029 const MFloat sax2 =
1030 Frj * mu_eps * (depsdx * cornerMetrics[xsd * 2 + xsd] + depsdy * cornerMetrics[xsd * 2 + ysd]);
1031 const MFloat say1 = Frj * mu_k * (dTKEdx * cornerMetrics[ysd * 2 + xsd] + dTKEdy * cornerMetrics[ysd * 2 + ysd]);
1032 const MFloat say2 =
1033 Frj * mu_eps * (depsdx * cornerMetrics[ysd * 2 + xsd] + depsdy * cornerMetrics[ysd * 2 + ysd]);
1034
1035 eflux[0][IJ] = sax1; // diffusion of nutilde for every cell
1036 eflux[1][IJ] = sax2; // diffusion of nutilde for every cell
1037
1038 fflux[0][IJ] = say1; // diffusion of nutilde for every cell
1039 fflux[1][IJ] = say2; // diffusion of nutilde for every cell
1040 }
1041 }
1042
1046 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; ++j) {
1047 for(MInt i = m_noGhostLayers - 1; i < m_nCells[1] - m_noGhostLayers; ++i) {
1048 const MInt IJ = cellIndex(i, j);
1049 const MInt IJM = cellIndex(i, (j - 1));
1050
1051 sa_1flux[0][IJ] = F1B2 * (eflux[0][IJ] + eflux[0][IJM]);
1052 sa_2flux[0][IJ] = F1B2 * (eflux[1][IJ] + eflux[1][IJM]);
1053 }
1054 }
1055
1056 for(MInt j = m_noGhostLayers - 1; j < m_nCells[0] - m_noGhostLayers; ++j) {
1057 for(MInt i = m_noGhostLayers; i < m_nCells[1] - m_noGhostLayers; ++i) {
1058 const MInt IJ = cellIndex(i, j);
1059 const MInt IMJ = cellIndex((i - 1), j);
1060
1061 sa_1flux[1][IJ] = F1B2 * (fflux[0][IJ] + fflux[0][IMJ]);
1062 sa_2flux[1][IJ] = F1B2 * (fflux[1][IJ] + fflux[1][IMJ]);
1063 }
1064 }
1065
1066
1067 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; j++) {
1068 for(MInt i = m_noGhostLayers; i < m_nCells[1] - m_noGhostLayers; i++) {
1069 const MInt IJ = cellIndex(i, j);
1070 const MInt IMJ = cellIndex((i - 1), j);
1071 const MInt IJM = cellIndex(i, (j - 1));
1072 const MInt IPJ = cellIndex((i + 1), j);
1073 const MInt IJP = cellIndex(i, (j + 1));
1074
1078 const MFloat diffusion_TKE = (sa_1flux[0][IJ] - sa_1flux[0][IMJ]) + (sa_1flux[1][IJ] - sa_1flux[1][IJM]);
1079
1080 const MFloat diffusion_eps = (sa_2flux[0][IJ] - sa_2flux[0][IMJ]) + (sa_2flux[1][IJ] - sa_2flux[1][IJM]);
1081
1085 const MFloat invCellJac = 1.0 / m_cells->cellJac[IJ];
1086
1087 // compute
1088 const MFloat dudxi = 0.5 * (u[IPJ] - u[IMJ]); // TODO_SS labels:FV,totest check the 0.5; it should be correct
1089 const MFloat dudet = 0.5 * (u[IJP] - u[IJM]);
1090 const MFloat dvdxi = 0.5 * (v[IPJ] - v[IMJ]);
1091 const MFloat dvdet = 0.5 * (v[IJP] - v[IJM]);
1092
1093 const MFloat dudx =
1094 dudxi * m_cells->cellMetrics[xsd * 2 + xsd][IJ] + dudet * m_cells->cellMetrics[ysd * 2 + xsd][IJ];
1095 const MFloat dudy =
1096 dudxi * m_cells->cellMetrics[xsd * 2 + ysd][IJ] + dudet * m_cells->cellMetrics[ysd * 2 + ysd][IJ];
1097 const MFloat dvdx =
1098 dvdxi * m_cells->cellMetrics[xsd * 2 + xsd][IJ] + dvdet * m_cells->cellMetrics[ysd * 2 + xsd][IJ];
1099 const MFloat dvdy =
1100 dvdxi * m_cells->cellMetrics[xsd * 2 + ysd][IJ] + dvdet * m_cells->cellMetrics[ysd * 2 + ysd][IJ];
1101
1102 // Production P = rho*nu_t/Re * ( 2*dudx^2 + 2*dvdy^2 + (dudy+dvdx)^2 - (2/3)*(dudx+dvdy)^2)
1103 // - 2/3*rho*k * ( dudx + dvdy )
1104 const MFloat P = ((rRe0 / POW2(fac_nonDim)) * invCellJac * muTurb[IJ]
1105 * (2.0 * POW2(dudx) + 2.0 * POW2(dvdy) + POW2(dudy + dvdx) - fac * F2B3 * POW2(dudx + dvdy))
1106 - fac * F2B3 * rho[IJ] * TKE[IJ] * (dudx + dvdy))
1107 * invCellJac;
1108
1109 // Ret is only required for inner cells, for which eps & TKE must be positive
1110 const MFloat Ret =
1111 fac_nonDim * m_solver->m_Re0 * rho[IJ] * POW2(TKE[IJ]) / muLam[IJ] / std::max(eps[IJ], minMFloat);
1112 const MFloat f2 = 1 - 0.222222222222222 /*=0.4/1.8*/ * exp(-POW2(Ret) * 0.027777777777777777 /*=1/36*/);
1113
1114 const MFloat Mt2 = 1.5 * TKE[IJ] * rho[IJ] / (gamma * p[IJ]);
1115 const MFloat Fdist2 = 1.0 / POW2(wallDist[IJ]);
1116 const MFloat k_diss1 = fac_nonDim * rho[IJ] * eps[IJ] * (1 + Mt2);
1117 const MFloat k_diss2 = 2.0 * rRe0 * muLam[IJ] * TKE[IJ] * Fdist2;
1118
1119 const MFloat tau = m_cells->turbTimeScale[IJ]; // eps[IJ]/std::max(TKE[IJ], minMFloat); // time scale?!
1120 const MFloat eps_prod = RM_KEPS::C1 * tau * P;
1121 const MFloat yp = utau[IJ] * wallDist[IJ]; // utau[IJ]*wallDist[IJ]*rho[IJ]/muLam[IJ]; // TODO_SS labels:FV
1122 const MFloat exp_wall = exp(-RM_KEPS::C4 * m_solver->m_Re0 * yp);
1123
1124 // limiter
1125 MFloat f2_ = f2;
1126 MFloat f3 = 1.0;
1127 if(m_solver->m_solutionAnomaly && m_cells->isAnomCell[IJ]) {
1128 const MFloat invRhoEps = 1.0 / (rho[IJ] * eps[IJ]);
1129 const MFloat RHS1 = invRhoEps / tau * diffusion_eps / m_cells->cellJac[IJ]
1130 - invRhoEps * diffusion_TKE / m_cells->cellJac[IJ] + (RM_KEPS::C1 - 1.0) * P * invRhoEps
1131 + fac_nonDim * (1 + Mt2);
1132 const MFloat f2_limit = (RHS1 + k_diss2 * invRhoEps * (1 - exp_wall)) / (fac_nonDim * RM_KEPS::C2);
1133 f2_ = std::min(1.0, std::max(f2, f2_limit));
1134 const MFloat diff = std::max(0.0, f2_limit - 1.0);
1135 f3 = std::max(0.0, 1.0 - fac_nonDim * RM_KEPS::C2 * diff / (k_diss2 * invRhoEps * (1 - exp_wall)));
1136 }
1137
1138 const MFloat k_diss = k_diss1 + f3 * k_diss2;
1139 const MFloat eps_diss = tau * (fac_nonDim * RM_KEPS::C2 * f2_ * rho[IJ] * eps[IJ] + f3 * k_diss2 * exp_wall);
1140
1141 const MFloat prodDest_TKE = (P - k_diss) * m_cells->cellJac[IJ];
1142 const MFloat prodDest_eps = (eps_prod - eps_diss) * m_cells->cellJac[IJ];
1143
1147 m_cells->rightHandSide[CV->RANS_VAR[0]][IJ] += diffusion_TKE;
1148 m_cells->rightHandSide[CV->RANS_VAR[0]][IJ] += prodDest_TKE;
1149 m_cells->rightHandSide[CV->RANS_VAR[1]][IJ] += diffusion_eps;
1150 m_cells->rightHandSide[CV->RANS_VAR[1]][IJ] += prodDest_eps;
1151 }
1152 }
1153}
1154
1155void FvStructuredSolver2DRans::viscousFlux_KEPSILON2() {
1156 // compute friction velocity at wall
1157 m_solver->m_structuredBndryCnd->computeFrictionCoef();
1158 // communicate wall properties to all cells
1159 m_solver->m_structuredBndryCnd->distributeWallAndFPProperties();
1160
1161 computeTurbViscosity_KEPSILON();
1162
1163 // OTHER variables required to calculate the laminar viscous fluxes
1164 static constexpr MFloat minMFloat =
1165 1e-16; // std::min(std::numeric_limits<MFloat>::min(), 1.0/std::numeric_limits<MFloat>::max());
1166 const MFloat rRe0 = F1 / m_solver->m_Re0;
1167 const MFloat gamma = m_solver->m_gamma;
1168 const MFloat fac = (m_solver->m_rans2eq_mode == "production") ? 1.0 : 0.0;
1169 const MFloat fac_nonDim = m_solver->m_keps_nonDimType ? 1.0 : PV->UInfinity;
1170 const MFloat transPos = m_solver->m_ransTransPos;
1171
1172 MFloat* const* const RESTRICT eflux = ALIGNED_F(m_cells->eFlux);
1173 MFloat* const* const RESTRICT fflux = ALIGNED_F(m_cells->fFlux);
1174 MFloat* const* const RESTRICT sa_1flux = ALIGNED_F(m_cells->saFlux1);
1175 MFloat* const* const RESTRICT sa_2flux = ALIGNED_F(m_cells->saFlux2);
1176
1177 MFloat* const RESTRICT u = ALIGNED_F(m_cells->pvariables[PV->U]);
1178 MFloat* const RESTRICT v = ALIGNED_F(m_cells->pvariables[PV->V]);
1179 MFloat* const RESTRICT rho = ALIGNED_F(m_cells->pvariables[PV->RHO]);
1180 MFloat* const RESTRICT p = ALIGNED_F(m_cells->pvariables[PV->P]);
1181 // MFloat* const RESTRICT T = ALIGNED_F(m_cells->temperature);
1182 MFloat* const RESTRICT TKE = ALIGNED_F(m_cells->pvariables[PV->RANS_VAR[0]]);
1183 MFloat* const RESTRICT eps = ALIGNED_F(m_cells->pvariables[PV->RANS_VAR[1]]);
1184 MFloat* const RESTRICT muLam = ALIGNED_F(m_cells->fq[FQ->MU_L]);
1185 MFloat* const RESTRICT muTurb = ALIGNED_F(m_cells->fq[FQ->MU_T]);
1186 const MFloat* const RESTRICT utau = &m_cells->fq[FQ->UTAU][0];
1187 const MFloat* const RESTRICT wallDist = &m_cells->fq[FQ->WALLDISTANCE][0];
1188
1189 for(MInt j = m_noGhostLayers - 1; j < m_nCells[0] - m_noGhostLayers; j++) {
1190 for(MInt i = m_noGhostLayers - 1; i < m_nCells[1] - m_noGhostLayers; i++) {
1191 // get the adjacent cells;
1192
1193 const MInt IJ = cellIndex(i, j);
1194 const MInt IPJ = cellIndex((i + 1), j);
1195 const MInt IPJP = cellIndex((i + 1), (j + 1));
1196 const MInt IJP = cellIndex(i, (j + 1));
1197
1198 const MFloat cornerMetrics[nDim * nDim] = {m_cells->cornerMetrics[0][IJ], m_cells->cornerMetrics[1][IJ],
1199 m_cells->cornerMetrics[2][IJ], m_cells->cornerMetrics[3][IJ]};
1200
1201 // TODO_SS labels:FV,totest if wall is along eta=constant, then at wall dTKEdxi and depsdxi must be zero! Check!!!
1202
1203 const MFloat dTKEdxi = F1B2 * (TKE[IPJP] + TKE[IPJ] - TKE[IJP] - TKE[IJ]);
1204 const MFloat dTKEdet = F1B2 * (TKE[IPJP] + TKE[IJP] - TKE[IPJ] - TKE[IJ]);
1205
1206 const MFloat depsdxi = F1B2 * (eps[IPJP] + eps[IPJ] - eps[IJP] - eps[IJ]);
1207 const MFloat depsdet = F1B2 * (eps[IPJP] + eps[IJP] - eps[IPJ] - eps[IJ]);
1208
1209 const MFloat metricTerms = cornerMetrics[xsd * 2 + xsd] * cornerMetrics[ysd * 2 + xsd]
1210 + cornerMetrics[xsd * 2 + ysd] * cornerMetrics[ysd * 2 + ysd];
1211
1212 const MFloat invCornerJac = F1 / POW2(m_cells->cornerJac[IJ]);
1213
1214 const MFloat sax1 = invCornerJac * dTKEdet * metricTerms;
1215 const MFloat sax2 = invCornerJac * depsdet * metricTerms;
1216 const MFloat say1 = invCornerJac * dTKEdxi * metricTerms;
1217 const MFloat say2 = invCornerJac * depsdxi * metricTerms;
1218
1219 eflux[0][IJ] = sax1; // diffusion of nutilde for every cell
1220 eflux[1][IJ] = sax2; // diffusion of nutilde for every cell
1221
1222 fflux[0][IJ] = say1; // diffusion of nutilde for every cell
1223 fflux[1][IJ] = say2; // diffusion of nutilde for every cell
1224
1228 if(m_P_keps) {
1229 const MFloat dudxi = F1B2 * (u[IPJP] + u[IPJ] - u[IJP] - u[IJ]);
1230 const MFloat dudet = F1B2 * (u[IPJP] + u[IJP] - u[IPJ] - u[IJ]);
1231
1232 const MFloat dvdxi = F1B2 * (v[IPJP] + v[IPJ] - v[IJP] - v[IJ]);
1233 const MFloat dvdet = F1B2 * (v[IPJP] + v[IJP] - v[IPJ] - v[IJ]);
1234
1235 const MFloat S11 = dudxi * cornerMetrics[xsd * 2 + xsd] + dudet * cornerMetrics[ysd * 2 + xsd];
1236 const MFloat S12 = 0.5
1237 * (dudxi * cornerMetrics[xsd * 2 + ysd] + dudet * cornerMetrics[ysd * 2 + ysd]
1238 + dvdxi * cornerMetrics[xsd * 2 + xsd] + dvdet * cornerMetrics[ysd * 2 + xsd]);
1239 const MFloat S22 = dvdxi * cornerMetrics[xsd * 2 + ysd] + dvdet * cornerMetrics[ysd * 2 + ysd];
1240 const MFloat S = S11 * S11 + 2 * S12 * S12 + S22 * S22;
1241 // TODO_SS labels:FV,totest muTurbAvg at wall surface should be zero! Check!
1242 const MFloat muTurbAvg = F1B4 * (muTurb[IJP] + muTurb[IPJP] + muTurb[IJ] + muTurb[IPJ]);
1243
1244 // if ((m_solver->domainId()==7 && IJ==83) || (m_solver->domainId()==8 && IJ==203))
1245 // muTurbAvg = 0;
1246
1247 const MFloat mueOverRe = rRe0 * invCornerJac * muTurbAvg / POW2(fac_nonDim);
1248 // The simplified (incompressible) version reads: 2*rho*nu_t/Re * S_ij*S_ij
1249 m_cells->P_keps[IJ] = 2.0 * mueOverRe * S;
1250 }
1251 }
1252 }
1253
1254
1255 // diffusive flux correction for the singular points
1256 // m_hasSingularity=0 means no singular points in this solver, otherwise do flux correction
1257 if(m_solver->m_hasSingularity > 0) {
1258 diffusiveFluxCorrection();
1259 }
1260
1264 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; ++j) {
1265 for(MInt i = m_noGhostLayers - 1; i < m_nCells[1] - m_noGhostLayers; ++i) {
1266 const MInt IJ = cellIndex(i, j);
1267 const MInt IJM = cellIndex(i, (j - 1));
1268 const MInt IPJ = cellIndex(i + 1, j);
1269 const MFloat surf0 = m_cells->surfaceMetrics[0][IJ];
1270 const MFloat surf1 = m_cells->surfaceMetrics[1][IJ];
1271
1272 const MFloat dTKEdxi = TKE[IPJ] - TKE[IJ];
1273 const MFloat depsdxi = eps[IPJ] - eps[IJ];
1274
1275 const MFloat muLamAvg_xi = F1B2 * (muLam[IJ] + muLam[IPJ]);
1276 const MFloat muTurbAvg_xi = F1B2 * (muTurb[IJ] + muTurb[IPJ]);
1277
1278 // TODO_SS labels:FV here surface jac (m_cells->surfJac) verwenden!!!
1279 const MFloat surfJac =
1280 0.5
1281 * (m_cells->cornerJac[IJ] + m_cells->cornerJac[IJM]); // TODO_SS labels:FV,totest Is this the right way to go?
1282 const MFloat Frj = rRe0 * surfJac;
1283 const MFloat mu_k_xi = muLamAvg_xi + muTurbAvg_xi * RM_KEPS::rsigma_k;
1284 const MFloat mu_eps_xi = muLamAvg_xi + muTurbAvg_xi * RM_KEPS::rsigma_eps;
1285
1286 const MFloat temp_TKE = dTKEdxi * (POW2(surf0) + POW2(surf1)) * POW2(1.0 / surfJac);
1287 const MFloat temp_eps = depsdxi * (POW2(surf0) + POW2(surf1)) * POW2(1.0 / surfJac);
1288
1289 MFloat limiterVisc1 = 1.0;
1290 MFloat limiterVisc2 = 1.0;
1291 if(m_solver->m_limiterVisc) {
1292 // TODO_SS labels:FV nochmal richtig anschauen
1293 const MFloat mu_ref = F1B2 * (m_cells->fq[FQ->LIMITERVISC][IJ] + m_cells->fq[FQ->LIMITERVISC][IPJ]);
1294 limiterVisc1 = std::min(1.0, mu_ref / std::abs(rRe0 * mu_k_xi));
1295 limiterVisc2 = std::min(1.0, mu_ref / std::abs(rRe0 * mu_eps_xi));
1296 }
1297
1298 sa_1flux[0][IJ] = Frj * mu_k_xi * (temp_TKE + F1B2 * (eflux[0][IJ] + eflux[0][IJM])) * limiterVisc1;
1299 sa_2flux[0][IJ] = Frj * mu_eps_xi * (temp_eps + F1B2 * (eflux[1][IJ] + eflux[1][IJM])) * limiterVisc2;
1300 }
1301 }
1302
1303 for(MInt j = m_noGhostLayers - 1; j < m_nCells[0] - m_noGhostLayers; ++j) {
1304 for(MInt i = m_noGhostLayers; i < m_nCells[1] - m_noGhostLayers; ++i) {
1305 const MInt IJ = cellIndex(i, j);
1306 const MInt IMJ = cellIndex((i - 1), j);
1307 const MInt IJP = cellIndex(i, j + 1);
1308 const MFloat surf0 = m_cells->surfaceMetrics[2 + 0][IJ];
1309 const MFloat surf1 = m_cells->surfaceMetrics[2 + 1][IJ];
1310
1311 const MFloat dTKEdet = TKE[IJP] - TKE[IJ];
1312 const MFloat depsdet = eps[IJP] - eps[IJ];
1313
1314 const MFloat muLamAvg_eta = F1B2 * (muLam[IJ] + muLam[IJP]);
1315 const MFloat muTurbAvg_eta = F1B2 * (muTurb[IJ] + muTurb[IJP]);
1316
1317 // const MFloat limit = 1e12;
1318 // MFloat mu = muLamAvg_eta+muTurbAvg_eta;
1319 // m_cells->fq[FQ->VAR1][IJ] = 1.0;
1320 // if (m_cells->coordinates[0][IJ]>0.0) {
1321 // const MInt IJM = cellIndex(i,j-1);
1322 // const MFloat mu_ref =
1323 // limit*POW2(0.5*(m_cells->coordinates[1][IJP]-m_cells->coordinates[1][IJM]))/m_solver->m_timeRef;
1324 // m_cells->fq[FQ->VAR1][IJ] = std::min(mu_ref/mu, 1.0);
1325 // }
1326
1327 // TODO_SS labels:FV here surface jac (m_cells->surfJac) verwenden!!!
1328 const MFloat surfJac =
1329 0.5
1330 * (m_cells->cornerJac[IJ] + m_cells->cornerJac[IMJ]); // TODO_SS labels:FV,totest Is this the right way to go?
1331 const MFloat Frj = rRe0 * surfJac;
1332 const MFloat mu_k_eta = muLamAvg_eta + muTurbAvg_eta * RM_KEPS::rsigma_k;
1333 const MFloat mu_eps_eta = muLamAvg_eta + muTurbAvg_eta * RM_KEPS::rsigma_eps;
1334
1335 const MFloat temp_TKE = dTKEdet * (POW2(surf0) + POW2(surf1)) * POW2(1.0 / surfJac);
1336 const MFloat temp_eps = depsdet * (POW2(surf0) + POW2(surf1)) * POW2(1.0 / surfJac);
1337
1338 MFloat limiterVisc1 = 1.0;
1339 MFloat limiterVisc2 = 1.0;
1340 if(m_solver->m_limiterVisc) {
1341 // TODO_SS labels:FV nochmal richtig anschauen
1342 const MFloat mu_ref = F1B2 * (m_cells->fq[FQ->LIMITERVISC][IJ] + m_cells->fq[FQ->LIMITERVISC][IJP]);
1343 limiterVisc1 = std::min(1.0, mu_ref / std::abs(rRe0 * mu_k_eta));
1344 limiterVisc2 = std::min(1.0, mu_ref / std::abs(rRe0 * mu_eps_eta));
1345 }
1346
1347 sa_1flux[1][IJ] = Frj * mu_k_eta * (temp_TKE + F1B2 * (fflux[0][IJ] + fflux[0][IMJ])) * limiterVisc1;
1348 sa_2flux[1][IJ] = Frj * mu_eps_eta * (temp_eps + F1B2 * (fflux[1][IJ] + fflux[1][IMJ])) * limiterVisc2;
1349 }
1350 }
1351
1352
1353 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; j++) {
1354 for(MInt i = m_noGhostLayers; i < m_nCells[1] - m_noGhostLayers; i++) {
1355 const MInt IJ = cellIndex(i, j);
1356 const MInt IMJ = cellIndex((i - 1), j);
1357 const MInt IJM = cellIndex(i, (j - 1));
1358 const MInt IPJ = cellIndex((i + 1), j);
1359 const MInt IJP = cellIndex(i, (j + 1));
1360
1361 // Actualy skipping should not be necessary, but somehow for localTimeStepping, the code might react
1362 // to the non-meaningful values even though these are not used
1363 if(m_cells->coordinates[0][IJ] <= transPos) continue;
1364
1368 MFloat limiterVisc1 = 1.0;
1369 MFloat limiterVisc2 = 1.0;
1370 /* if (m_solver->m_limiterVisc) {
1371 //TODO_SS labels:FV nochmal richtig anschauen
1372 const MFloat mu_ref = m_cells->fq[FQ->LIMITERVISC][IJ];
1373 limiterVisc1 = std::min(1.0, mu_ref/std::abs(rRe0*(muLam[IJ]+muTurb[IJ]*RM_KEPS::rsigma_k)));
1374 limiterVisc2 = std::min(1.0, mu_ref/std::abs(rRe0*(muLam[IJ]+muTurb[IJ]*RM_KEPS::rsigma_eps)));
1375 }*/
1376
1377 const MFloat diffusion_TKE =
1378 (sa_1flux[0][IJ] - sa_1flux[0][IMJ]) + (sa_1flux[1][IJ] - sa_1flux[1][IJM]) * limiterVisc1;
1379
1380 const MFloat diffusion_eps =
1381 (sa_2flux[0][IJ] - sa_2flux[0][IMJ]) + (sa_2flux[1][IJ] - sa_2flux[1][IJM]) * limiterVisc2;
1382
1386 const MFloat invCellJac = 1.0 / m_cells->cellJac[IJ];
1387
1388 MFloat P;
1389 if(m_P_keps) {
1390 const MInt IMJM = cellIndex(i - 1, j - 1);
1391 P = F1B4 * (m_cells->P_keps[IJ] + m_cells->P_keps[IMJ] + m_cells->P_keps[IJM] + m_cells->P_keps[IMJM]);
1392 } else {
1393 // compute
1394 const MFloat dudxi = 0.5 * (u[IPJ] - u[IMJ]); // TODO_SS labels:FV,totest check the 0.5; it should be correct
1395 const MFloat dudet = 0.5 * (u[IJP] - u[IJM]);
1396 const MFloat dvdxi = 0.5 * (v[IPJ] - v[IMJ]);
1397 const MFloat dvdet = 0.5 * (v[IJP] - v[IJM]);
1398
1399 const MFloat dudx =
1400 dudxi * m_cells->cellMetrics[xsd * 2 + xsd][IJ] + dudet * m_cells->cellMetrics[ysd * 2 + xsd][IJ];
1401 const MFloat dudy =
1402 dudxi * m_cells->cellMetrics[xsd * 2 + ysd][IJ] + dudet * m_cells->cellMetrics[ysd * 2 + ysd][IJ];
1403 const MFloat dvdx =
1404 dvdxi * m_cells->cellMetrics[xsd * 2 + xsd][IJ] + dvdet * m_cells->cellMetrics[ysd * 2 + xsd][IJ];
1405 const MFloat dvdy =
1406 dvdxi * m_cells->cellMetrics[xsd * 2 + ysd][IJ] + dvdet * m_cells->cellMetrics[ysd * 2 + ysd][IJ];
1407
1408 // Production P = rho*nu_t/Re * ( 2*dudx^2 + 2*dvdy^2 + (dudy+dvdx)^2 - (2/3)*(dudx+dvdy)^2)
1409 // - 2/3*rho*k * ( dudx + dvdy )
1410 P = ((rRe0 / POW2(fac_nonDim)) * invCellJac * muTurb[IJ]
1411 * (2.0 * POW2(dudx) + 2.0 * POW2(dvdy) + POW2(dudy + dvdx) - fac * F2B3 * POW2(dudx + dvdy))
1412 - fac * F2B3 * rho[IJ] * TKE[IJ] * (dudx + dvdy))
1413 * invCellJac;
1414 }
1415
1416 // Ret is only required for inner cells, for which eps & TKE must be positive
1417 const MFloat Ret =
1418 fac_nonDim * m_solver->m_Re0 * rho[IJ] * POW2(TKE[IJ]) / muLam[IJ] / std::max(eps[IJ], minMFloat);
1419 const MFloat f2 = 1 - 0.222222222222222 /*=0.4/1.8*/ * exp(-POW2(Ret) * 0.027777777777777777 /*=1/36*/);
1420
1421 const MFloat Mt2 = 1.5 * TKE[IJ] * rho[IJ] / (gamma * p[IJ]);
1422 const MFloat Fdist2 = 1.0 / POW2(wallDist[IJ]);
1423 const MFloat k_diss1 = fac_nonDim * rho[IJ] * eps[IJ] * (1 + Mt2);
1424 const MFloat k_diss2 = 2.0 * rRe0 * muLam[IJ] * TKE[IJ] * Fdist2;
1425
1426 const MFloat tau = m_cells->turbTimeScale[IJ]; // eps[IJ]/std::max(TKE[IJ], minMFloat); // time scale?!
1427 const MFloat eps_prod = RM_KEPS::C1 * tau * P;
1428 const MFloat yp = utau[IJ] * wallDist[IJ]; // utau[IJ]*wallDist[IJ]*rho[IJ]/muLam[IJ]; // TODO_SS labels:FV
1429 const MFloat exp_wall = exp(-RM_KEPS::C4 * m_solver->m_Re0 * yp);
1430
1431 // limiter
1432 MFloat f2_ = f2;
1433 MFloat f3 = 1.0;
1434 if(m_solver->m_solutionAnomaly && m_cells->isAnomCell[IJ]) {
1435 const MFloat invRhoEps = 1.0 / std::max(rho[IJ] * eps[IJ], minMFloat);
1436 // The additional porous terms which might be part of the rhs are omitted
1437 const MFloat RHS1 = invRhoEps / tau * diffusion_eps / m_cells->cellJac[IJ]
1438 - invRhoEps * diffusion_TKE / m_cells->cellJac[IJ] + (RM_KEPS::C1 - 1.0) * P * invRhoEps
1439 + fac_nonDim * (1 + Mt2);
1440 const MFloat f2_limit = (RHS1 + k_diss2 * invRhoEps * (1 - exp_wall)) / (fac_nonDim * RM_KEPS::C2);
1441 f2_ = std::min(1.0, std::max(f2, f2_limit));
1442 const MFloat diff = std::max(0.0, f2_limit - 1.0);
1443 f3 = std::max(0.0, 1.0 - fac_nonDim * RM_KEPS::C2 * diff / (k_diss2 * invRhoEps * (1 - exp_wall)));
1444 }
1445
1446 const MFloat k_diss = k_diss1 + f3 * k_diss2;
1447 const MFloat eps_diss = tau * (fac_nonDim * RM_KEPS::C2 * f2_ * rho[IJ] * eps[IJ] + f3 * k_diss2 * exp_wall);
1448
1449 const MFloat prodDest_TKE = (P - k_diss) * m_cells->cellJac[IJ];
1450 const MFloat prodDest_eps = (eps_prod - eps_diss) * m_cells->cellJac[IJ];
1451
1455 // m_cells->fq[FQ->VAR1][IJ] = 2.0*POW2(dudx);//m_cells->rightHandSide[CV->RANS_VAR[0]][IJ];
1456 // m_cells->fq[FQ->VAR2][IJ] = dudy;//diffusion_TKE;
1457 // m_cells->fq[FQ->VAR3][IJ] = prodDest_TKE;
1458 m_cells->rightHandSide[CV->RANS_VAR[0]][IJ] += diffusion_TKE;
1459 m_cells->rightHandSide[CV->RANS_VAR[0]][IJ] += prodDest_TKE;
1460 m_cells->rightHandSide[CV->RANS_VAR[1]][IJ] += diffusion_eps;
1461 m_cells->rightHandSide[CV->RANS_VAR[1]][IJ] += prodDest_eps;
1462 // m_cells->fq[FQ->VAR4][IJ] = dvdx;//m_cells->rightHandSide[CV->RANS_VAR[0]][IJ];
1463 // m_cells->fq[FQ->VAR1][IJ] = rho[IJ]*v[IJ];
1464 // m_cells->fq[FQ->VAR2][IJ] = v[IJ];
1465 }
1466 }
1467
1468 // call the standard LES viscous flux
1469 if(m_solver->m_viscCompact)
1470 m_solver->viscousFluxLESCompact<true>();
1471 else
1472 m_solver->viscousFluxLES<true>();
1473}
1474
1475
1476void FvStructuredSolver2DRans::computeTurbViscosity() { (this->*compTurbVisc)(); }
1477
1478void FvStructuredSolver2DRans::computeTurbViscosity_SA() {
1479 // OTHER variables required to calculate the laminar viscous fluxes
1480 const MFloat transPos = m_solver->m_ransTransPos;
1481 MFloat* const RESTRICT p = &m_cells->pvariables[PV->P][0];
1482 MFloat* const RESTRICT rho = &m_cells->pvariables[PV->RHO][0];
1483 MFloat* const RESTRICT nuTilde = &m_cells->pvariables[PV->RANS_VAR[0]][0];
1484 MFloat* const RESTRICT T = &m_cells->temperature[0];
1485
1486 for(MInt i = 0; i < m_noCells; i++) {
1487 if(m_cells->coordinates[0][i] <= transPos) {
1488 nuTilde[i] = 0.0;
1489 } else {
1490 nuTilde[i] = mMin(mMax(nuTilde[i], 0.0), 3000.0);
1491 }
1492 T[i] = m_solver->m_gamma * p[i] / rho[i];
1493 // decode the kinematic turbulent viscosity from the turb dynamic visc arrays
1494 const MFloat nuLaminar = SUTHERLANDLAW(T[i]) / rho[i];
1495 const MFloat chi = nuTilde[i] / (nuLaminar);
1496 const MFloat fv1 = pow(chi, 3) / (pow(chi, 3) + RM_SA_DV::cv1to3);
1497 const MFloat nuTurb = fv1 * nuTilde[i];
1498 m_cells->fq[FQ->NU_T][i] = nuTurb;
1499 m_cells->fq[FQ->MU_T][i] = rho[i] * nuTurb;
1500 }
1501}
1502
1503void FvStructuredSolver2DRans::computeTurbViscosity_KEPSILON() {
1504 // nu_t is also required in the ghost cells at boundaries because of e.g. viscousFluxLES
1505 // OTHER variables required to calculate the laminar viscous fluxes
1506 static constexpr MFloat minMFloat =
1507 1e-16; // std::min(std::numeric_limits<MFloat>::min(), 1.0/std::numeric_limits<MFloat>::max());
1508 const MFloat transPos = m_solver->m_ransTransPos;
1509 const MFloat Re0 = m_solver->m_Re0;
1510 const MFloat fac_nonDim = m_solver->m_keps_nonDimType ? 1.0 : PV->UInfinity;
1511 static const MFloat sqrtF2B3 = sqrt(2.0) / 3.0;
1512 const MFloat c_wd_eff = m_solver->m_c_wd;
1513 MFloat* const RESTRICT p = &m_cells->pvariables[PV->P][0];
1514 MFloat* const RESTRICT rho = &m_cells->pvariables[PV->RHO][0];
1515 MFloat* const RESTRICT TKE = &m_cells->pvariables[PV->RANS_VAR[0]][0];
1516 MFloat* const RESTRICT eps = &m_cells->pvariables[PV->RANS_VAR[1]][0];
1517 MFloat* const RESTRICT T = &m_cells->temperature[0];
1518 MFloat* const RESTRICT u = &m_cells->pvariables[PV->U][0];
1519 MFloat* const RESTRICT v = &m_cells->pvariables[PV->V][0];
1520 const MFloat* const RESTRICT utau = &m_cells->fq[FQ->UTAU][0];
1521 const MFloat* const RESTRICT utau2 =
1522 (m_solver->m_blockType != "porous") ? &m_cells->fq[FQ->UTAU][0] : &m_cells->fq[FQ->UTAU2][0];
1523 const MFloat* const RESTRICT wallDist = &m_cells->fq[FQ->WALLDISTANCE][0];
1524 MFloat* const RESTRICT muLam = &m_cells->fq[FQ->MU_L][0];
1525 const MFloat* const RESTRICT por = &m_cells->fq[FQ->POROSITY][0];
1526 const MFloat* const RESTRICT Da = &m_cells->fq[FQ->DARCY][0];
1527
1528 for(MInt jj = m_noGhostLayers - 1; jj < m_nCells[0] - m_noGhostLayers + 1; jj++) {
1529 for(MInt ii = m_noGhostLayers - 1; ii < m_nCells[1] - m_noGhostLayers + 1; ii++) {
1530 const MInt i = cellIndex(ii, jj);
1531 if(m_cells->coordinates[0][i] <= transPos) {
1532 TKE[i] = PV->ransInfinity[0];
1533 eps[i] = PV->ransInfinity[1];
1534 } else {
1535 // TODO_SS labels:FV maybe else not necessary
1536 }
1537 T[i] = m_solver->m_gamma * p[i] / rho[i];
1538 // decode the kinematic turbulent viscosity from the turb dynamic visc arrays
1539 muLam[i] = SUTHERLANDLAW(T[i]); // TODO_SS labels:FV put this maybe somewhere else
1540 // const MFloat nuLaminar = muLam[i]/rho[i];
1541 const MFloat yp = utau[i] * wallDist[i];
1542 // TODO_SS labels:FV,totest Using Da & por in halo cells at fluid-porous interfaces might be not correct
1543 // TODO_SS labels:FV,toenhance In the following think if it makes sense to take the minmum of yp or wallDistance
1544 // const MFloat yp_ = (m_solver->m_blockType=="porous" && c_wd_eff*sqrt(Da[i]/por[i])<wallDist[i]) ?
1545 // utau2[i]*c_wd_eff*sqrt(Da[i]/por[i]) : yp;
1546 const MFloat yp_ = (m_solver->m_blockType == "porous" && utau2[i] * c_wd_eff * sqrt(Da[i] / por[i]) < yp)
1547 ? utau2[i] * c_wd_eff * sqrt(Da[i] / por[i])
1548 : yp;
1549 const MFloat f_mu = 1.0 - exp(-RM_KEPS::C3 * Re0 * yp_);
1550 const MFloat sgn = eps[i] < 0 ? -1.0 : 1.0;
1551 const MFloat eps_i = sgn * mMax(fabs(eps[i]), minMFloat); // eps & TKE can be negative in BC treatment
1552
1553 //
1554 // Note: actually the computation of the derivatives here, needs to be in the same way as in the
1555 // computation of the production term
1556 // TODO_SS labels:FV determine neighbor cells by i+i, etc.
1557 const MInt IPJ = getCellIdfromCell(i, 1, 0);
1558 const MInt IMJ = getCellIdfromCell(i, -1, 0);
1559 const MInt IJP = getCellIdfromCell(i, 0, 1);
1560 const MInt IJM = getCellIdfromCell(i, 0, -1);
1561 const MFloat dudxi = 0.5 * (u[IPJ] - u[IMJ]); // TODO_SS labels:FV,totest check the 0.5; it should be correct
1562 const MFloat dudet = 0.5 * (u[IJP] - u[IJM]);
1563 const MFloat dvdxi = 0.5 * (v[IPJ] - v[IMJ]);
1564 const MFloat dvdet = 0.5 * (v[IJP] - v[IJM]);
1565
1566 const MFloat dudx =
1567 dudxi * m_cells->cellMetrics[xsd * 2 + xsd][i] + dudet * m_cells->cellMetrics[ysd * 2 + xsd][i];
1568 const MFloat dudy =
1569 dudxi * m_cells->cellMetrics[xsd * 2 + ysd][i] + dudet * m_cells->cellMetrics[ysd * 2 + ysd][i];
1570 const MFloat dvdx =
1571 dvdxi * m_cells->cellMetrics[xsd * 2 + xsd][i] + dvdet * m_cells->cellMetrics[ysd * 2 + xsd][i];
1572 const MFloat dvdy =
1573 dvdxi * m_cells->cellMetrics[xsd * 2 + ysd][i] + dvdet * m_cells->cellMetrics[ysd * 2 + ysd][i];
1574
1575 const MFloat invCellJac = 1.0 / m_cells->cellJac[i];
1576 const MFloat S11 = F2B3 * dudx - F1B3 * dvdy;
1577 const MFloat S12 = 0.5 * (dudy + dvdx);
1578 const MFloat S22 = F2B3 * dvdy - F1B3 * dudx;
1579 const MFloat S = sqrt(S11 * S11 + 2 * S12 * S12 + S22 * S22) * invCellJac;
1580 const MFloat temp = sqrtF2B3 / std::max(S, minMFloat);
1581 const MFloat nuTurbMax = temp * fabs(TKE[i]) * Re0 * POW2(fac_nonDim);
1582
1583 const MFloat nuTurb_true = fabs(Re0 * RM_KEPS::C_mu * f_mu * POW2(TKE[i]) / eps_i * fac_nonDim);
1584 // m_cells->fq[FQ->VAR2][i] = nuTurbMax;
1585 // if (nuTurbMax<nuTurb_true)
1586 // m_cells->fq[FQ->VAR1][i] = true;
1587 // else
1588 // m_cells->fq[FQ->VAR1][i] = false;
1589 const MFloat nuTurb = sgn * std::min(nuTurbMax, nuTurb_true);
1590 m_cells->fq[FQ->NU_T][i] = nuTurb;
1591 m_cells->fq[FQ->MU_T][i] = rho[i] * nuTurb;
1592 // turbTimeScale is always positive, because, if TKE negative in boundary ghost cell, then also eps
1593 m_cells->turbTimeScale[i] =
1594 RM_KEPS::C_mu * f_mu * fabs(TKE[i]) * Re0 * fac_nonDim / std::max(fabs(nuTurb), minMFloat);
1595
1596 //
1597 /* MFloat tau1 = F4B3 * dudx - F2B3 * dvdy;
1598 MFloat tau2 = dudy + dvdx;
1599 const MFloat mueOverRe = 1.0/(Re0*m_cells->cellJac[i])*nuTurb;
1600 m_cells->fq[FQ->VAR3][i] = -mueOverRe*tau1+F2B3*TKE[i];
1601 m_cells->fq[FQ->VAR4][i] = -mueOverRe*tau2;*/
1602 }
1603 }
1604}
1605
1606
1607void FvStructuredSolver2DRans::setAndAllocate_KEPSILON() {
1608 mAlloc(m_cells->saFlux1, nDim, m_noCells, "m_cells->saFlux1", -999999.9, AT_);
1609 mAlloc(m_cells->saFlux2, nDim, m_noCells, "m_cells->saFlux2", -999999.9, AT_);
1610 mAlloc(m_cells->turbTimeScale, m_noCells, "m_cells->turbTimeScale", -999999.9, AT_);
1611
1612 // Modified production term computation required temporary buffer
1613 m_P_keps = Context::getSolverProperty<MBool>("P_keps", m_solverId, AT_, &m_P_keps);
1614 if(m_P_keps) {
1615 if(m_solver->m_hasSingularity) mTerm(1, "P_keps=true for grid with singularities not implemented yet!");
1616 mAlloc(m_cells->P_keps, m_noCells, "m_cells->P_keps", -999999.9, AT_);
1617 }
1618
1619 //
1620 m_solver->m_solutionAnomaly = false;
1621 m_solver->m_solutionAnomaly =
1622 Context::getSolverProperty<MBool>("solutionAnomaly", m_solverId, AT_, &m_solver->m_solutionAnomaly);
1623 if(m_solver->m_solutionAnomaly) {
1624 // TODO_SS labels:FV smoothly blend the anomaly area from the rest
1625 // By default if solution anomaly treatment is enabled all cells will be taken into account
1626 mAlloc(m_cells->isAnomCell, m_noCells, "m_cells->isAnomCell", true, AT_);
1627 if(Context::propertyExists("anomPoint", m_solverId) + Context::propertyExists("anomRadius") == 1)
1628 mTerm(1, "Either you specify 'anomPoint' and 'anomRadius' or neither of them!");
1629 if(Context::propertyExists("anomPoint", m_solverId)) {
1630 for(MInt i = 0; i < m_noCells; ++i) {
1631 m_cells->isAnomCell[i] = false;
1632 }
1633 const MInt noAnomPoints = Context::propertyLength("anomPoint", m_solverId);
1634 const MInt noAnomRadius = Context::propertyLength("anomRadius", m_solverId);
1635 if(noAnomPoints % nDim != 0 || noAnomPoints != nDim * noAnomRadius)
1636 mTerm(1, "IQ>10 mandatory for usage of k-epsilon model!!!");
1637 for(MInt n = 0; n < noAnomRadius; ++n) {
1638 const MFloat rPOW2 = POW2(Context::getSolverProperty<MFloat>("anomRadius", m_solverId, AT_, n));
1639 const MFloat c[nDim] = {Context::getSolverProperty<MFloat>("anomPoint", m_solverId, AT_, n * nDim),
1640 Context::getSolverProperty<MFloat>("anomPoint", m_solverId, AT_, n * nDim + 1)};
1641 for(MInt i = 0; i < m_noCells; ++i) {
1642 if(POW2(m_cells->coordinates[0][i] - c[0]) + POW2(m_cells->coordinates[1][i] - c[1]) < rPOW2)
1643 m_cells->isAnomCell[i] = true;
1644 }
1645 }
1646 }
1647 } // if(m_solver->m_solutionAnomaly)
1648}
1649
1650
1651void FvStructuredSolver2DRans::diffusiveFluxCorrection() {
1652 MFloat* const* const RESTRICT eflux = ALIGNED_F(m_cells->eFlux);
1653 MFloat* const* const RESTRICT fflux = ALIGNED_F(m_cells->fFlux);
1654 MFloat* const RESTRICT TKE_ = ALIGNED_F(m_cells->pvariables[PV->RANS_VAR[0]]);
1655 MFloat* const RESTRICT EPS_ = ALIGNED_F(m_cells->pvariables[PV->RANS_VAR[1]]);
1656
1657 const auto& singularity = m_solver->m_singularity;
1658
1659 // MInt dim = 0;
1660 MInt start[nDim], end[nDim], nghbr[10];
1661 MInt len1[nDim];
1662 // MInt totalCells;
1663
1664 for(MInt i = 0; i < m_solver->m_hasSingularity; ++i) {
1665 // only correct for bc 6000 not for bc 4000-5000
1666 if(singularity[i].BC == -6000) {
1667 // totalCells=1;
1668 for(MInt j = 0; j < nDim; j++) {
1669 len1[j] = singularity[i].end[j] - singularity[i].start[j];
1670 // if(len1[j]!=0) totalCells*=len1[j];
1671 }
1672
1673 for(MInt n = 0; n < nDim; ++n) {
1674 if(singularity[i].end[n] - singularity[i].start[n] > 1) {
1675 mTerm(1, "In 2D not possible!");
1676 // dim=n;
1677 // start[n]=singularity[i].start[n]+1;
1678 start[n] = singularity[i].start[n] + 1;
1679 end[n] = singularity[i].end[n] - 1;
1680 } else {
1681 start[n] = singularity[i].start[n];
1682 end[n] = singularity[i].end[n];
1683 }
1684 }
1685
1686 MFloat TKE[10], EPS[10];
1687 MFloat TKEcorner, EPScorner;
1688 for(MInt jj = start[1]; jj < end[1]; ++jj) {
1689 for(MInt ii = start[0]; ii < end[0]; ++ii) {
1690 MInt count = 0;
1691 // MInt temp[nDim]{};
1692
1693 const MInt IJ_ = cellIndex(ii + singularity[i].Viscous[0], jj + singularity[i].Viscous[1]);
1694 // const MFloat* const surf = m_cells->surfaceMetrics[IJ_];
1695 // dxidx, dxidy, detadx, detady
1696 const MFloat surfaceMetricsS[nDim * nDim] = {
1697 m_cells->surfaceMetricsSingularity[0][i], m_cells->surfaceMetricsSingularity[1][i],
1698 m_cells->surfaceMetricsSingularity[2][i], m_cells->surfaceMetricsSingularity[3][i]};
1699
1700 // temp[dim]=1;
1701 const MInt IJ = cellIndex(ii, jj);
1702 nghbr[count++] = IJ;
1703 // nghbr[count++]=cellIndex(ii+temp[0],jj+temp[1],kk+temp[2]);
1704
1705 const MInt IPMJ_ = getCellIdfromCell(IJ, singularity[i].Viscous[0], 0);
1706 const MInt IJPM_ = getCellIdfromCell(IJ, 0, singularity[i].Viscous[1]);
1707
1708 const MFloat surfaceMetrics[nDim * nDim] = {
1709 m_cells->surfaceMetrics[0][IPMJ_], m_cells->surfaceMetrics[1][IPMJ_], m_cells->surfaceMetrics[2][IJPM_],
1710 m_cells->surfaceMetrics[3][IJPM_]};
1711
1712 for(MInt m = 0; m < singularity[i].Nstar - 1; ++m) {
1713 MInt* change = singularity[i].displacement[m];
1714 nghbr[count++] = cellIndex(ii + change[0], jj + change[1]);
1715 // nghbr[count++]=cellIndex(ii+temp[0]+change[0],jj+temp[1]+change[1],kk+temp[2]+change[2]);
1716 }
1717
1718 if(count != singularity[i].Nstar) {
1719 cout << "what the hell! it is wrong!!!" << endl;
1720 }
1721
1722 for(MInt m = 0; m < singularity[i].Nstar; ++m) {
1723 TKE[m] = TKE_[nghbr[m]];
1724 EPS[m] = EPS_[nghbr[m]];
1725 }
1726
1727 TKEcorner = F0;
1728 EPScorner = F0;
1729
1730 MInt id2 = ii - start[0] + (jj - start[1]) * len1[0];
1731
1732 for(MInt n = 0; n < count; n++) {
1733 MInt ID = id2 * count + n;
1734 TKEcorner += singularity[i].ReconstructionConstants[0][ID] * TKE[n];
1735 EPScorner += singularity[i].ReconstructionConstants[0][ID] * EPS[n];
1736 }
1737
1738 const MInt sign_xi = 2 * singularity[i].Viscous[0] + 1;
1739 const MInt sign_eta = 2 * singularity[i].Viscous[1] + 1;
1740
1741 const MInt IPMJ = getCellIdfromCell(IJ, sign_xi, 0);
1742 const MInt IJPM = getCellIdfromCell(IJ, 0, sign_eta);
1743
1744 const MFloat dTKEdet = sign_eta * 0.5 * (TKEcorner - 0.5 * (TKE[0] + TKE_[IPMJ]));
1745 const MFloat dEPSdet = sign_eta * 0.5 * (EPScorner - 0.5 * (EPS[0] + EPS_[IPMJ]));
1746
1747 const MFloat dTKEdxi = sign_xi * 0.5 * (TKEcorner - 0.5 * (TKE[0] + TKE_[IJPM]));
1748 const MFloat dEPSdxi = sign_xi * 0.5 * (EPScorner - 0.5 * (EPS[0] + EPS_[IJPM]));
1749
1750 // const MFloat metricTerms = cornerMetrics[ xsd * 2 + xsd ]*cornerMetrics[ ysd * 2 + xsd ]
1751 // + cornerMetrics[ xsd * 2 + ysd ]*cornerMetrics[ ysd * 2 + ysd ];
1752
1753 const MFloat metricTerms_xi = surfaceMetrics[xsd * 2 + xsd] * surfaceMetricsS[ysd * 2 + xsd]
1754 + surfaceMetrics[xsd * 2 + ysd] * surfaceMetricsS[ysd * 2 + ysd];
1755 const MFloat metricTerms_eta = surfaceMetricsS[xsd * 2 + xsd] * surfaceMetrics[ysd * 2 + xsd]
1756 + surfaceMetricsS[xsd * 2 + ysd] * surfaceMetrics[ysd * 2 + ysd];
1757
1758 const MFloat invSurfJac_xi = F1 / POW2(m_cells->surfJac[IPMJ_]);
1759 const MFloat invSurfJac_eta = F1 / POW2(m_cells->surfJac[m_noCells + IJPM_]);
1760
1761 const MFloat sax1 = invSurfJac_xi * dTKEdet * metricTerms_xi;
1762 const MFloat sax2 = invSurfJac_xi * dEPSdet * metricTerms_xi;
1763 const MFloat say1 = invSurfJac_eta * dTKEdxi * metricTerms_eta;
1764 const MFloat say2 = invSurfJac_eta * dEPSdxi * metricTerms_eta;
1765
1766 eflux[0][IJ_] = sax1; // diffusion of nutilde for every cell
1767 eflux[1][IJ_] = sax2; // diffusion of nutilde for every cell
1768
1769 fflux[0][IJ_] = say1; // diffusion of nutilde for every cell
1770 fflux[1][IJ_] = say2; // diffusion of nutilde for every cell
1771
1772 // cout << globalTimeStep << "(" << m_solver->m_RKStep << ") dom=" << m_solver->domainId() <<
1773 // setprecision(10) << " x|y=" << m_cells->coordinates[0][IJ] << "|" << m_cells->coordinates[1][IJ]
1774 // << " eflux=" << eflux[ 0*noCells+IJ_ ] << "|" << eflux[ 1*noCells+IJ_ ] << " fflux=" << fflux[
1775 // 0*noCells+IJ_ ] << "|" << fflux[ 1*noCells+IJ_ ] << " metricTerms_xi=" << metricTerms_xi << "
1776 // metricTerms_eta=" << metricTerms_eta << " " << m_cells->surfaceMetricsSingularity[i][0] << "|"
1777 // << m_cells->surfaceMetricsSingularity[i][1] << "|" << m_cells->surfaceMetricsSingularity[i][2]
1778 // << "|" << m_cells->surfaceMetricsSingularity[i][3]
1779 // << " " << surfaceMetrics[0] << "|" << surfaceMetrics[1] << "|" << surfaceMetrics[2] << " " <<
1780 // surfaceMetrics[3] << endl;
1781 }
1782 }
1783 }
1784 }
1785}
1786
1787
1788inline MInt FvStructuredSolver2DRans::cellIndex(MInt i, MInt j) { return i + j * m_nCells[1]; }
1789
1790inline MInt FvStructuredSolver2DRans::getCellIdfromCell(MInt origin, MInt incI, MInt incJ) {
1791 return origin + incI + incJ * m_nCells[1];
1792}
1793
1794inline MFloat FvStructuredSolver2DRans::getPSI(MInt I, MInt dim) {
1795 const MFloat FK = 18.0;
1796 const MInt IJK[2] = {1, m_nCells[1]};
1797 const MInt IP1 = I + IJK[dim];
1798 const MInt IM1 = I - IJK[dim];
1799 const MInt IP2 = I + 2 * IJK[dim];
1800
1801 const MFloat PIM2 = m_cells->pvariables[PV->P][IM1];
1802 const MFloat PIM1 = m_cells->pvariables[PV->P][I];
1803 const MFloat PIP2 = m_cells->pvariables[PV->P][IP2];
1804 const MFloat PIP1 = m_cells->pvariables[PV->P][IP1];
1805
1806 const MFloat PSI =
1807 mMin(F1,
1808 FK
1809 * mMax(mMax(fabs((PIM2 - PIM1) / mMin(PIM2, PIM1)), fabs((PIM1 - PIP1) / mMin(PIM1, PIP1))),
1810 fabs((PIP1 - PIP2) / mMin(PIP1, PIP2))));
1811 return PSI;
1812}
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
Context::propertyLength
static MInt propertyLength(const MString &name, MInt solverId=m_noSolvers)
Returns the number of elements of a property.
Definition: context.cpp:538
Context::propertyExists
static MBool propertyExists(const MString &name, MInt solver=m_noSolvers)
This function checks if a property exists in general.
Definition: context.cpp:494
FvStructuredSolver2D
2D structured solver class
Definition: fvstructuredsolver2d.h:21
FvStructuredSolver2D::m_structuredBndryCnd
StructuredBndryCnd< 2 > * m_structuredBndryCnd
Definition: fvstructuredsolver2d.h:120
FvStructuredSolver2DRans::CV
std::unique_ptr< MConservativeVariables< 2 > > & CV
Definition: fvstructuredsolver2drans.h:61
FvStructuredSolver2DRans::FvStructuredSolver2DRans
FvStructuredSolver2DRans(FvStructuredSolver2D *solver)
Definition: fvstructuredsolver2drans.cpp:23
FvStructuredSolver2DRans::PV
std::unique_ptr< MPrimitiveVariables< 2 > > & PV
Definition: fvstructuredsolver2drans.h:62
FvStructuredSolver2DRans::FQ
std::unique_ptr< StructuredFQVariables > & FQ
Definition: fvstructuredsolver2drans.h:63
FvStructuredSolver2DRans::viscFluxMethod
void(FvStructuredSolver2DRans::* viscFluxMethod)()
Definition: fvstructuredsolver2drans.h:39
FvStructuredSolver2DRans::m_solver
FvStructuredSolver2D * m_solver
Definition: fvstructuredsolver2drans.h:53
FvStructuredSolver2DRans::getCellIdfromCell
MInt getCellIdfromCell(MInt origin, MInt incI, MInt incJ)
Definition: fvstructuredsolver2drans.cpp:1790
FvStructuredSolver2DRans::reconstructSurfaceData
void(FvStructuredSolver2DRans::* reconstructSurfaceData)()
Definition: fvstructuredsolver2drans.h:32
FvStructuredSolver2DRans::nDim
static constexpr const MInt nDim
Definition: fvstructuredsolver2drans.h:68
FvStructuredSolver2DRans::compTurbVisc
void(FvStructuredSolver2DRans::* compTurbVisc)()
Definition: fvstructuredsolver2drans.h:40
FvStructuredSolver< 2 >::FQ
std::unique_ptr< StructuredFQVariables > FQ
Definition: fvstructuredsolver.h:469
FvStructuredSolver::m_singularity
SingularInformation * m_singularity
Definition: fvstructuredsolver.h:835
FvStructuredSolver< 2 >::PV
std::unique_ptr< MPrimitiveVariables< nDim > > PV
Definition: fvstructuredsolver.h:333
FvStructuredSolver< 2 >::CV
std::unique_ptr< MConservativeVariables< nDim > > CV
Definition: fvstructuredsolver.h:332
Solver::m_solverId
const MInt m_solverId
a unique solver identifier
Definition: solver.h:90
StructuredBndryCnd::computeFrictionCoef
virtual void computeFrictionCoef()
Definition: fvstructuredbndrycnd.h:180
StructuredBndryCnd::distributeWallAndFPProperties
virtual void distributeWallAndFPProperties()
Definition: fvstructuredbndrycnd.h:182
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::surfaceMetricsSingularity
MFloat ** surfaceMetricsSingularity
Definition: fvstructuredcell.h:45
StructuredCell::turbTimeScale
MFloat * turbTimeScale
Definition: fvstructuredcell.h:77
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_SA
@ RANS_SA
Definition: enums.h:53
RANS_KEPSILON
@ RANS_KEPSILON
Definition: enums.h:58
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
m_log
InfoOutFile m_log
Definition: globalvariables.cpp:57
MInt
int32_t MInt
Definition: maiatypes.h:62
MUint
uint32_t MUint
Definition: maiatypes.h:63
MFloat
double MFloat
Definition: maiatypes.h:52
MBool
bool MBool
Definition: maiatypes.h:58
RANSModelConstants< RANS_KEPSILON >::C1
static constexpr MFloat C1
Definition: fvransmodelconstants.h:93
RANSModelConstants< RANS_KEPSILON >::C2
static constexpr MFloat C2
Definition: fvransmodelconstants.h:94
RANSModelConstants< RANS_KEPSILON >::C3
static constexpr MFloat C3
Definition: fvransmodelconstants.h:95
RANSModelConstants< RANS_KEPSILON >::rsigma_eps
static constexpr MFloat rsigma_eps
Definition: fvransmodelconstants.h:91
RANSModelConstants< RANS_KEPSILON >::C_mu
static constexpr MFloat C_mu
Definition: fvransmodelconstants.h:92
RANSModelConstants< RANS_KEPSILON >::C4
static constexpr MFloat C4
Definition: fvransmodelconstants.h:96
RANSModelConstants< RANS_KEPSILON >::rsigma_k
static constexpr MFloat rsigma_k
Definition: fvransmodelconstants.h:90
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