OpenLB 1.8.1
Loading...
Searching...
No Matches
setTurbulentWallModel.h
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2024 Fedor Bukreev
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22*/
23
24#ifndef SET_TURBULENT_WALL_MODEL_H
25#define SET_TURBULENT_WALL_MODEL_H
26
27#include "postprocessor/turbulentWallModelPostProcessor.h"
28
29namespace olb {
30
31//======================================================================
32// ======== Dynamics for WMLES ======//
33//======================================================================
34
35template <typename T, typename DESCRIPTOR>
36using ForcedWMHRRdynamics = dynamics::Tuple<
37 T, DESCRIPTOR,
38 momenta::Tuple<
39 momenta::BulkDensity,
40 momenta::MovingPorousMomentumCombination<momenta::BulkMomentum>,
41 momenta::BulkStress,
42 momenta::DefineToNEq
43 >,
44 equilibria::ThirdOrder,
45 collision::ParameterFromCell<collision::HYBRID,collision::LocalSmagorinskyEffectiveOmega<collision::HRR>>,
46 forcing::WFGuoThirdOrder<momenta::ForcedWithStress>
47>;
48
49template <typename T, typename DESCRIPTOR>
50using WMHRRdynamics = dynamics::Tuple<
51 T, DESCRIPTOR,
52 momenta::Tuple<
53 momenta::BulkDensity,
54 momenta::MovingPorousMomentumCombination<momenta::BulkMomentum>,
55 momenta::BulkStress,
56 momenta::DefineToNEq
57 >,
58 equilibria::ThirdOrder,
59 collision::ParameterFromCell<collision::HYBRID,collision::LocalSmagorinskyEffectiveOmega<collision::HRR>>
60>;
61
62template <typename T, typename DESCRIPTOR>
63using ForcedVanDriestWMHRRdynamics = dynamics::Tuple<
64 T, DESCRIPTOR,
65 momenta::Tuple<
66 momenta::BulkDensity,
67 momenta::MovingPorousMomentumCombination<momenta::BulkMomentum>,
68 momenta::BulkStress,
69 momenta::DefineToNEq
70 >,
71 equilibria::ThirdOrder,
72 collision::ParameterFromCell<collision::HYBRID,collision::LocalVanDriestSmagorinskyEffectiveOmega<collision::HRR>>,
73 forcing::WFGuoThirdOrder<momenta::ForcedWithStress>
74>;
75
76template <typename T, typename DESCRIPTOR>
77using VanDriestWMHRRdynamics = dynamics::Tuple<
78 T, DESCRIPTOR,
79 momenta::Tuple<
80 momenta::BulkDensity,
81 momenta::MovingPorousMomentumCombination<momenta::BulkMomentum>,
82 momenta::BulkStress,
83 momenta::DefineToNEq
84 >,
85 equilibria::ThirdOrder,
86 collision::ParameterFromCell<collision::HYBRID,collision::LocalVanDriestSmagorinskyEffectiveOmega<collision::HRR>>
87>;
88
89template <typename T, typename DESCRIPTOR>
90using ForcedVanDriestExternalRhoWMHRRdynamics = dynamics::Tuple<
91 T, DESCRIPTOR,
92 momenta::Tuple<
93 momenta::BulkDensity,
94 momenta::MovingPorousMomentumCombination<momenta::BulkMomentum>,
95 momenta::BulkStress,
96 momenta::DefineToNEq
97 >,
98 equilibria::ThirdOrder,
99 collision::ParameterFromCell<collision::HYBRID,
100 collision::ParameterFromCell<collision::HYBRID_RHO,
101 collision::LocalVanDriestSmagorinskyEffectiveOmega<collision::ExternalRhoHRR>>>,
102 forcing::WFGuoThirdOrder<momenta::ForcedWithStress>
103>;
104
105template <typename T, typename DESCRIPTOR>
106using VanDriestExternalRhoWMHRRdynamics = dynamics::Tuple<
107 T, DESCRIPTOR,
108 momenta::Tuple<
109 momenta::BulkDensity,
110 momenta::MovingPorousMomentumCombination<momenta::BulkMomentum>,
111 momenta::BulkStress,
112 momenta::DefineToNEq
113 >,
114 equilibria::ThirdOrder,
115 collision::ParameterFromCell<collision::HYBRID,
116 collision::ParameterFromCell<collision::HYBRID_RHO,
117 collision::LocalVanDriestSmagorinskyEffectiveOmega<collision::ExternalRhoHRR>>>
118>;
119
120template <typename T, typename DESCRIPTOR>
121using ForcedExternalRhoWMHRRdynamics = dynamics::Tuple<
122 T, DESCRIPTOR,
123 momenta::Tuple<
124 momenta::BulkDensity,
125 momenta::MovingPorousMomentumCombination<momenta::BulkMomentum>,
126 momenta::BulkStress,
127 momenta::DefineToNEq
128 >,
129 equilibria::ThirdOrder,
130 collision::ParameterFromCell<collision::HYBRID,
131 collision::ParameterFromCell<collision::HYBRID_RHO,
132 collision::LocalSmagorinskyEffectiveOmega<collision::ExternalRhoHRR>>>,
133 forcing::WFGuoThirdOrder<momenta::ForcedWithStress>
134>;
135
136template <typename T, typename DESCRIPTOR>
137using ExternalRhoWMHRRdynamics = dynamics::Tuple<
138 T, DESCRIPTOR,
139 momenta::Tuple<
140 momenta::BulkDensity,
141 momenta::MovingPorousMomentumCombination<momenta::BulkMomentum>,
142 momenta::BulkStress,
143 momenta::DefineToNEq
144 >,
145 equilibria::ThirdOrder,
146 collision::ParameterFromCell<collision::HYBRID,
147 collision::ParameterFromCell<collision::HYBRID_RHO,
148 collision::LocalSmagorinskyEffectiveOmega<collision::ExternalRhoHRR>>>
149>;
150
151//======================================================================
152// ======== Wall model parameters ======//
153//======================================================================
154
155template <typename T>
156struct WallModelParameters{
157 /* Used method for density reconstruction
158 * 0: use local density
159 * 1: extrapolation (Guo)
160 * 2: constant (rho = 1.)
161 */
162 int rhoMethod = 0;
163
164 /* Used method for non-equilibrium population reconstruction
165 * 0: extrapolation NEQ (Guo Zhaoli)
166 * 1: first or second order finite differnce (Malaspinas)
167 * 2: equilibrium scheme (no fNeq)
168 */
169 int fNeqMethod = 1;
170
171 /* Used wall profile
172 * 0: power law profile
173 * 1: Spalding profile
174 */
175 int wallFunctionProfile = 1;
176
178 bool bodyForce = true;
179
181 bool interpolateSampleVelocity = false;
182
184 bool useVanDriest = false;
185
187 T latticeWallDistance = T(0.5);
188
190 T samplingCellDistance = T(3.);
191
192 bool movingWall = false;
193
194 bool averageVelocity = false;
195};
196
197//======================================================================
198// Set wall model dynamics depending on chosen paramertes.
199// Use before the boundary setter.
200//======================================================================
201
202template<typename T, typename DESCRIPTOR>
203void setTurbulentWallModelDynamics(SuperLattice<T, DESCRIPTOR>& sLattice,
204 FunctorPtr<SuperIndicatorF<T,DESCRIPTOR::d>>&& bulkIndicator,
205 WallModelParameters<T>& wallModelParameters)
206{
207 if( wallModelParameters.averageVelocity ) {
208 if( wallModelParameters.bodyForce ) {
209 if( wallModelParameters.useVanDriest ) {
210 if( wallModelParameters.rhoMethod != 0 ) {
211 sLattice.template defineDynamics<typename ForcedVanDriestExternalRhoWMHRRdynamics<T,DESCRIPTOR>::template wrap_collision<collision::StoreAndTrackAverageVelocity>>(std::move(bulkIndicator));
212 } else {
213 sLattice.template defineDynamics<typename ForcedVanDriestWMHRRdynamics<T,DESCRIPTOR>::template wrap_collision<collision::StoreAndTrackAverageVelocity>>(std::move(bulkIndicator));
214 }
215 }
216 else {
217 if( wallModelParameters.rhoMethod != 0 ) {
218 sLattice.template defineDynamics<typename ForcedExternalRhoWMHRRdynamics<T,DESCRIPTOR>::template wrap_collision<collision::StoreAndTrackAverageVelocity>>(std::move(bulkIndicator));
219 } else {
220 sLattice.template defineDynamics<typename ForcedWMHRRdynamics<T,DESCRIPTOR>::template wrap_collision<collision::StoreAndTrackAverageVelocity>>(std::move(bulkIndicator));
221 }
222 }
223 }
224 else {
225 if( wallModelParameters.useVanDriest ) {
226 if( wallModelParameters.rhoMethod != 0 ) {
227 sLattice.template defineDynamics<typename VanDriestExternalRhoWMHRRdynamics<T,DESCRIPTOR>::template wrap_collision<collision::StoreAndTrackAverageVelocity>>(std::move(bulkIndicator));
228 } else {
229 sLattice.template defineDynamics<typename VanDriestWMHRRdynamics<T,DESCRIPTOR>::template wrap_collision<collision::StoreAndTrackAverageVelocity>>(std::move(bulkIndicator));
230 }
231 }
232 else {
233 if( wallModelParameters.rhoMethod != 0 ) {
234 sLattice.template defineDynamics<typename ExternalRhoWMHRRdynamics<T,DESCRIPTOR>::template wrap_collision<collision::StoreAndTrackAverageVelocity>>(std::move(bulkIndicator));
235 } else {
236 sLattice.template defineDynamics<typename WMHRRdynamics<T,DESCRIPTOR>::template wrap_collision<collision::StoreAndTrackAverageVelocity>>(std::move(bulkIndicator));
237 }
238 }
239 }
240 }
241 else {
242 if( wallModelParameters.bodyForce ) {
243 if( wallModelParameters.useVanDriest ) {
244 if( wallModelParameters.rhoMethod != 0 ) {
245 sLattice.template defineDynamics<ForcedVanDriestExternalRhoWMHRRdynamics>(std::move(bulkIndicator));
246 } else {
247 sLattice.template defineDynamics<ForcedVanDriestWMHRRdynamics>(std::move(bulkIndicator));
248 }
249 }
250 else {
251 if( wallModelParameters.rhoMethod != 0 ) {
252 sLattice.template defineDynamics<ForcedExternalRhoWMHRRdynamics>(std::move(bulkIndicator));
253 } else {
254 sLattice.template defineDynamics<ForcedWMHRRdynamics>(std::move(bulkIndicator));
255 }
256 }
257 }
258 else {
259 if( wallModelParameters.useVanDriest ) {
260 if( wallModelParameters.rhoMethod != 0 ) {
261 sLattice.template defineDynamics<VanDriestExternalRhoWMHRRdynamics>(std::move(bulkIndicator));
262 } else {
263 sLattice.template defineDynamics<VanDriestWMHRRdynamics>(std::move(bulkIndicator));
264 }
265 }
266 else {
267 if( wallModelParameters.rhoMethod != 0 ) {
268 sLattice.template defineDynamics<ExternalRhoWMHRRdynamics>(std::move(bulkIndicator));
269 } else {
270 sLattice.template defineDynamics<WMHRRdynamics>(std::move(bulkIndicator));
271 }
272 }
273 }
274 }
275
276 AnalyticalConst<DESCRIPTOR::d, T, T> one(T(1.));
277 AnalyticalConst<DESCRIPTOR::d, T, T> zero(T(0.));
278 AnalyticalConst<DESCRIPTOR::d, T, T> zeroVector(Vector<T,DESCRIPTOR::d>{});
279 std::vector<T> zeroStrainVec;
280 if(DESCRIPTOR::d == 2) {
281 zeroStrainVec = {T(0), T(0), T(0)};
282 } else {
283 zeroStrainVec = {T(0), T(0), T(0), T(0), T(0), T(0)};
284 }
285 AnalyticalConst<DESCRIPTOR::d, T, T> zeroStrain(zeroStrainVec);
286 sLattice.template defineField<descriptors::TENSOR>(std::move(bulkIndicator), zeroStrain);
287 sLattice.template defineField<descriptors::OMEGA>(std::move(bulkIndicator), one);
288 sLattice.template defineField<descriptors::WMPOROSITY>(std::move(bulkIndicator), one);
289 sLattice.template defineField<descriptors::U_TAU>(std::move(bulkIndicator), zero);
290 sLattice.template defineField<collision::HYBRID>(std::move(bulkIndicator), one);
291 if( wallModelParameters.rhoMethod != 0) {
292 sLattice.template defineField<collision::HYBRID_RHO>(std::move(bulkIndicator), one);
293 sLattice.template defineField<descriptors::DENSITY>(std::move(bulkIndicator), one);
294 }
295 sLattice.template defineField<descriptors::Y1>(std::move(bulkIndicator), zeroVector);
296 sLattice.template defineField<descriptors::WMVELOCITY>(std::move(bulkIndicator), zeroVector);
297 if( wallModelParameters.useVanDriest ) {
298 sLattice.template defineField<descriptors::VISCOSITY>(std::move(bulkIndicator), zero);
299 }
300 if( wallModelParameters.movingWall ) {
301 sLattice.template defineField<descriptors::VELOCITY>(std::move(bulkIndicator), zeroVector);
302 }
303 {
304 auto& communicator = sLattice.getCommunicator(stage::PostStream());
305 communicator.template requestField<descriptors::WMVELOCITY>();
306 communicator.template requestField<descriptors::POPULATION>();
307 communicator.requestOverlap(sLattice.getOverlap());
308 communicator.exchangeRequests();
309 }
310}
311
312template<typename T, typename DESCRIPTOR>
313void setTurbulentWallModelDynamics(SuperLattice<T, DESCRIPTOR>& sLattice,
314 SuperGeometry<T,DESCRIPTOR::d>& superGeometry,
315 int materialOfBulk,
316 WallModelParameters<T>& wallModelParameters)
317{
318 // Getting the indicators by material numbers and calling the superLattice method via the indicators:
319 setTurbulentWallModelDynamics<T,DESCRIPTOR>(sLattice,
320 FunctorPtr<SuperIndicatorF<T,DESCRIPTOR::d>>(superGeometry.getMaterialIndicator(materialOfBulk)),
321 wallModelParameters);
322}
323
324//======================================================================
325// Setter of the wall model boundary.
326// Can be used with or without geometry indicator input.
327//======================================================================
328
329template<typename T, typename DESCRIPTOR>
330void setTurbulentWallModel(SuperLattice<T, DESCRIPTOR>& sLattice,
331 FunctorPtr<SuperIndicatorF<T,DESCRIPTOR::d>>&& boundaryIndicator,
332 FunctorPtr<SuperIndicatorF<T,DESCRIPTOR::d>>&& bulkIndicator,
333 WallModelParameters<T>& wallModelParameters,
334 IndicatorF<T,DESCRIPTOR::d>* indicatorAnalyticalBoundary = nullptr)
335{
336 int _overlap = 3;
337 OstreamManager clout(std::cout, "TurbulentWallModelSetter");
338
339 AnalyticalConst<DESCRIPTOR::d, T, T> one(T(1.));
340 AnalyticalConst<DESCRIPTOR::d, T, T> zero(T(0.));
341 AnalyticalConst<DESCRIPTOR::d, T, T> zeroVector(T(0.), T(0.));
342 sLattice.template defineField<descriptors::OMEGA>(std::move(boundaryIndicator), one);
343 sLattice.template defineField<descriptors::WMPOROSITY>(std::move(boundaryIndicator), zero);
344 sLattice.template defineField<descriptors::U_TAU>(std::move(boundaryIndicator), zero);
345 sLattice.template defineField<collision::HYBRID>(std::move(boundaryIndicator), one);
346 if( wallModelParameters.rhoMethod != 0) {
347 sLattice.template defineField<collision::HYBRID_RHO>(std::move(boundaryIndicator), one);
348 sLattice.template defineField<descriptors::DENSITY>(std::move(boundaryIndicator), one);
349 }
350 sLattice.template defineField<descriptors::Y1>(std::move(boundaryIndicator), zeroVector);
351 sLattice.template defineField<descriptors::WMVELOCITY>(std::move(boundaryIndicator), zeroVector);
352 if( wallModelParameters.useVanDriest ) {
353 sLattice.template defineField<descriptors::VISCOSITY>(std::move(bulkIndicator), zero);
354 }
355 if( wallModelParameters.movingWall ) {
356 sLattice.template defineField<descriptors::VELOCITY>(std::move(bulkIndicator), zeroVector);
357 }
358
359 auto& load = sLattice.getLoadBalancer();
360 for (int iC=0; iC < load.size(); ++iC) {
361 setTurbulentWallModel<T,DESCRIPTOR>(sLattice.getBlock(iC),
362 (bulkIndicator->getBlockIndicatorF(iC)).getBlockGeometry(),
363 boundaryIndicator->getBlockIndicatorF(iC),
364 bulkIndicator->getBlockIndicatorF(iC),
365 wallModelParameters,
366 indicatorAnalyticalBoundary);
367 }
368 addPoints2CommBC<T,DESCRIPTOR>(sLattice, std::forward<decltype(boundaryIndicator)>(boundaryIndicator), _overlap);
369 clout.setMultiOutput(false);
370 sLattice.template setParameter<descriptors::SAMPLING_DISTANCE>( wallModelParameters.samplingCellDistance );
371}
372
374template<typename T, typename DESCRIPTOR>
375void setTurbulentWallModel(SuperLattice<T, DESCRIPTOR>& sLattice,
376 SuperGeometry<T,DESCRIPTOR::d>& superGeometry,
377 int materialOfSolidObstacle,
378 WallModelParameters<T>& wallModelParameters,
379 IndicatorF<T,DESCRIPTOR::d>* indicatorAnalyticalBoundary = nullptr,
380 std::vector<int> bulkMaterials = std::vector<int>(1,1))
381{
382 // Getting the indicators by material numbers and calling the superLattice method via the indicators:
383 setTurbulentWallModel<T,DESCRIPTOR>(sLattice,
384 FunctorPtr<SuperIndicatorF<T,DESCRIPTOR::d>>(superGeometry.getMaterialIndicator(materialOfSolidObstacle)),
385 FunctorPtr<SuperIndicatorF<T,DESCRIPTOR::d>>(superGeometry.getMaterialIndicator(std::move(bulkMaterials))),
386 wallModelParameters,
387 indicatorAnalyticalBoundary);
388}
389
391template<typename T, typename DESCRIPTOR>
392void setTurbulentWallModel(BlockLattice<T,DESCRIPTOR>& block,
393 BlockGeometry<T,DESCRIPTOR::d>& blockGeometry,
394 BlockIndicatorF<T,DESCRIPTOR::d>& boundaryIndicator,
395 BlockIndicatorF<T,DESCRIPTOR::d>& bulkIndicator,
396 WallModelParameters<T>& wallModelParameters,
397 IndicatorF<T,DESCRIPTOR::d>* indicatorAnalyticalBoundary = nullptr)
398{
399 OstreamManager clout(std::cout, "TurbulentWallModelSetter");
400 clout.setMultiOutput(true);
401
402 // Defining boundary distance y1 and distance to exchange location y2 in wall normal direction for every boundary cell
403 const T deltaR = blockGeometry.getDeltaR();
404 // for each solid cell:
405 int kMax = 2;
406 if( wallModelParameters.fNeqMethod == 1 ) {
407 kMax = 3;
408 }
409 for(int k=1; k < kMax; k++){
410 block.forSpatialLocations([&](LatticeR<DESCRIPTOR::d> solidLatticeR) {
411 // Check if cell is solid cell
412 if (boundaryIndicator(solidLatticeR)) {
413 for (int iPop=1; iPop < DESCRIPTOR::q; ++iPop) {
414 Vector<T,DESCRIPTOR::d> boundaryLatticeR(solidLatticeR + k*descriptors::c<DESCRIPTOR>(iPop));
415 if (blockGeometry.getNeighborhoodRadius(boundaryLatticeR) >= 1) {
416 if (blockGeometry.isInside(boundaryLatticeR)) {
417 Vector<T,DESCRIPTOR::d> boundaryPhysR{ };
418 blockGeometry.getPhysR(boundaryPhysR,boundaryLatticeR);
419 // check if neighbor is fluid cell
420 if (bulkIndicator(boundaryLatticeR)) {
421 if (indicatorAnalyticalBoundary) {
422 // Calculate surface normal
423 Vector<T,DESCRIPTOR::d> normal = indicatorAnalyticalBoundary->surfaceNormal(boundaryPhysR, deltaR);
424 T y1 = 0.;
425 // Calculate boundary distance y1 in surface normal direction
426 indicatorAnalyticalBoundary->distance(y1, boundaryPhysR, normal, blockGeometry.getIcGlob());
427 if(y1 == T(0)) {
428 for ( int iD = 0; iD < DESCRIPTOR::d; iD++ ) {
429 normal[iD] *= T(-1);
430 }
431 indicatorAnalyticalBoundary->distance(y1, boundaryPhysR, normal, blockGeometry.getIcGlob());
432 }
433 y1 /= deltaR; // y1 in lattice units
434 if ( util::abs(y1) > T(k) ) {
435 y1 = util::sign(y1) * T(k);
436 }
437 block.get(boundaryLatticeR).template setField<descriptors::Y1>(-y1*normal);
438 } else {
439 Vector<T,DESCRIPTOR::d> normal;
440 for(int jPop = 0; jPop < DESCRIPTOR::q; jPop++){
441 Vector<T,DESCRIPTOR::d> boundaryLatticeR2(boundaryLatticeR + k*descriptors::c<DESCRIPTOR>(jPop));
442 if(boundaryIndicator(boundaryLatticeR2)){
443 for ( int iD = 0; iD < DESCRIPTOR::d; iD++ ) {
444 normal[iD] += k*descriptors::c<DESCRIPTOR>(jPop,iD);
445 }
446 }
447 }
448 T normalNorm = util::norm<DESCRIPTOR::d>(normal);
449 if(normalNorm != T(0)) {
450 for ( int iD = 0; iD < DESCRIPTOR::d; iD++ ) {
451 normal[iD] /= normalNorm;
452 }
453 }
454 for ( int iD = 0; iD < DESCRIPTOR::d; iD++ ) {
455 if(util::abs(normal[iD]) < T(0.4)) {
456 normal[iD] = T(0);
457 }
458 }
459 normalNorm = util::norm<DESCRIPTOR::d>(normal);
460 for ( int iD = 0; iD < DESCRIPTOR::d; iD++ ) {
461 normal[iD] /= normalNorm;
462 }
463 if(normal[0] != T(0) && normal[0] > T(0)) {
464 normal[0] /= normal[0];
465 }
466 else if(normal[0] != T(0) && normal[0] < T(0)) {
467 normal[0] /= normal[0];
468 normal[0] *= T(-1.);
469 }
470 if(normal[1] != T(0) && normal[1] > T(0)) {
471 normal[1] /= normal[1];
472 }
473 else if(normal[1] != T(0) && normal[1] < T(0)) {
474 normal[1] /= normal[1];
475 normal[1] *= T(-1.);
476 }
477 if ( DESCRIPTOR::d == 3) {
478 if(normal[2] != T(0) && normal[2] > T(0)) {
479 normal[2] /= normal[2];
480 }
481 else if(normal[2] != T(0) && normal[2] < T(0)) {
482 normal[2] /= normal[2];
483 normal[2] *= T(-1.);
484 }
485 }
486 if(normalNorm != T(0)) {
487 auto field = block.get(boundaryLatticeR).template getField<descriptors::Y1>();
488 if( util::norm<DESCRIPTOR::d>(field) == T(0)) {
489 T y1 = (k-1) + wallModelParameters.latticeWallDistance;
490 block.get(boundaryLatticeR).template setField<descriptors::Y1>(-y1*normal);
491 }
492 }
493 }
494 T pi[util::TensorVal<DESCRIPTOR>::n] {T(0)};
495 block.get(boundaryLatticeR).template setField<descriptors::TENSOR>(pi);
496 // Setting turbulent wall model post processor
497 if( wallModelParameters.movingWall) {
498 if( wallModelParameters.bodyForce ) {
499 if( wallModelParameters.interpolateSampleVelocity ) {
500 if( wallModelParameters.useVanDriest && k == 1 ) {
501 if( wallModelParameters.wallFunctionProfile == 0 ) {
502 block.addPostProcessor(typeid(stage::PostStream),
503 boundaryLatticeR,
504 meta::id<TurbulentWallModelPostProcessor<true,true,true,0,true>>{});
505 } else {
506 block.addPostProcessor(typeid(stage::PostStream),
507 boundaryLatticeR,
508 meta::id<TurbulentWallModelPostProcessor<true,true,true,1,true>>{});
509 }
510 } else {
511 if( wallModelParameters.wallFunctionProfile == 0 ) {
512 block.addPostProcessor(typeid(stage::PostStream),
513 boundaryLatticeR,
514 meta::id<TurbulentWallModelPostProcessor<true,true,false,0,true>>{});
515 } else {
516 block.addPostProcessor(typeid(stage::PostStream),
517 boundaryLatticeR,
518 meta::id<TurbulentWallModelPostProcessor<true,true,false,1,true>>{});
519 }
520 }
521 } else {
522 if( wallModelParameters.useVanDriest && k == 1 ) {
523 if( wallModelParameters.wallFunctionProfile == 0 ) {
524 block.addPostProcessor(typeid(stage::PostStream),
525 boundaryLatticeR,
526 meta::id<TurbulentWallModelPostProcessor<true,false,true,0,true>>{});
527 } else {
528 block.addPostProcessor(typeid(stage::PostStream),
529 boundaryLatticeR,
530 meta::id<TurbulentWallModelPostProcessor<true,false,true,1,true>>{});
531 }
532 } else {
533 if( wallModelParameters.wallFunctionProfile == 0 ) {
534 block.addPostProcessor(typeid(stage::PostStream),
535 boundaryLatticeR,
536 meta::id<TurbulentWallModelPostProcessor<true,false,false,0,true>>{});
537 } else {
538 block.addPostProcessor(typeid(stage::PostStream),
539 boundaryLatticeR,
540 meta::id<TurbulentWallModelPostProcessor<true,false,false,1,true>>{});
541 }
542 }
543 }
544 } else {
545 if( wallModelParameters.interpolateSampleVelocity ) {
546 if( wallModelParameters.useVanDriest && k == 1 ) {
547 if( wallModelParameters.wallFunctionProfile == 0 ) {
548 block.addPostProcessor(typeid(stage::PostStream),
549 boundaryLatticeR,
550 meta::id<TurbulentWallModelPostProcessor<false,true,true,0,true>>{});
551 } else {
552 block.addPostProcessor(typeid(stage::PostStream),
553 boundaryLatticeR,
554 meta::id<TurbulentWallModelPostProcessor<false,true,true,1,true>>{});
555 }
556 } else {
557 if( wallModelParameters.wallFunctionProfile == 0 ) {
558 block.addPostProcessor(typeid(stage::PostStream),
559 boundaryLatticeR,
560 meta::id<TurbulentWallModelPostProcessor<false,true,false,0,true>>{});
561 } else {
562 block.addPostProcessor(typeid(stage::PostStream),
563 boundaryLatticeR,
564 meta::id<TurbulentWallModelPostProcessor<false,true,false,1,true>>{});
565 }
566 }
567 } else {
568 if( wallModelParameters.useVanDriest && k == 1 ) {
569 if( wallModelParameters.wallFunctionProfile == 0 ) {
570 block.addPostProcessor(typeid(stage::PostStream),
571 boundaryLatticeR,
572 meta::id<TurbulentWallModelPostProcessor<false,false,true,0,true>>{});
573 } else {
574 block.addPostProcessor(typeid(stage::PostStream),
575 boundaryLatticeR,
576 meta::id<TurbulentWallModelPostProcessor<false,false,true,1,true>>{});
577 }
578 } else {
579 if( wallModelParameters.wallFunctionProfile == 0 ) {
580 block.addPostProcessor(typeid(stage::PostStream),
581 boundaryLatticeR,
582 meta::id<TurbulentWallModelPostProcessor<false,false,false,0,true>>{});
583 } else {
584 block.addPostProcessor(typeid(stage::PostStream),
585 boundaryLatticeR,
586 meta::id<TurbulentWallModelPostProcessor<false,false,false,1,true>>{});
587 }
588 }
589 }
590 }
591 } else {
592 if( wallModelParameters.bodyForce ) {
593 if( wallModelParameters.interpolateSampleVelocity ) {
594 if( wallModelParameters.useVanDriest && k == 1 ) {
595 if( wallModelParameters.wallFunctionProfile == 0 ) {
596 block.addPostProcessor(typeid(stage::PostStream),
597 boundaryLatticeR,
598 meta::id<TurbulentWallModelPostProcessor<true,true,true,0,false>>{});
599 } else {
600 block.addPostProcessor(typeid(stage::PostStream),
601 boundaryLatticeR,
602 meta::id<TurbulentWallModelPostProcessor<true,true,true,1,false>>{});
603 }
604 } else {
605 if( wallModelParameters.wallFunctionProfile == 0 ) {
606 block.addPostProcessor(typeid(stage::PostStream),
607 boundaryLatticeR,
608 meta::id<TurbulentWallModelPostProcessor<true,true,false,0,false>>{});
609 } else {
610 block.addPostProcessor(typeid(stage::PostStream),
611 boundaryLatticeR,
612 meta::id<TurbulentWallModelPostProcessor<true,true,false,1,false>>{});
613 }
614 }
615 } else {
616 if( wallModelParameters.useVanDriest && k == 1 ) {
617 if( wallModelParameters.wallFunctionProfile == 0 ) {
618 block.addPostProcessor(typeid(stage::PostStream),
619 boundaryLatticeR,
620 meta::id<TurbulentWallModelPostProcessor<true,false,true,0,false>>{});
621 } else {
622 block.addPostProcessor(typeid(stage::PostStream),
623 boundaryLatticeR,
624 meta::id<TurbulentWallModelPostProcessor<true,false,true,1,false>>{});
625 }
626 } else {
627 if( wallModelParameters.wallFunctionProfile == 0 ) {
628 block.addPostProcessor(typeid(stage::PostStream),
629 boundaryLatticeR,
630 meta::id<TurbulentWallModelPostProcessor<true,false,false,0,false>>{});
631 } else {
632 block.addPostProcessor(typeid(stage::PostStream),
633 boundaryLatticeR,
634 meta::id<TurbulentWallModelPostProcessor<true,false,false,1,false>>{});
635 }
636 }
637 }
638 } else {
639 if( wallModelParameters.interpolateSampleVelocity ) {
640 if( wallModelParameters.useVanDriest && k == 1 ) {
641 if( wallModelParameters.wallFunctionProfile == 0 ) {
642 block.addPostProcessor(typeid(stage::PostStream),
643 boundaryLatticeR,
644 meta::id<TurbulentWallModelPostProcessor<false,true,true,0,false>>{});
645 } else {
646 block.addPostProcessor(typeid(stage::PostStream),
647 boundaryLatticeR,
648 meta::id<TurbulentWallModelPostProcessor<false,true,true,1,false>>{});
649 }
650 } else {
651 if( wallModelParameters.wallFunctionProfile == 0 ) {
652 block.addPostProcessor(typeid(stage::PostStream),
653 boundaryLatticeR,
654 meta::id<TurbulentWallModelPostProcessor<false,true,false,0,false>>{});
655 } else {
656 block.addPostProcessor(typeid(stage::PostStream),
657 boundaryLatticeR,
658 meta::id<TurbulentWallModelPostProcessor<false,true,false,1,false>>{});
659 }
660 }
661 } else {
662 if( wallModelParameters.useVanDriest && k == 1 ) {
663 if( wallModelParameters.wallFunctionProfile == 0 ) {
664 block.addPostProcessor(typeid(stage::PostStream),
665 boundaryLatticeR,
666 meta::id<TurbulentWallModelPostProcessor<false,false,true,0,false>>{});
667 } else {
668 block.addPostProcessor(typeid(stage::PostStream),
669 boundaryLatticeR,
670 meta::id<TurbulentWallModelPostProcessor<false,false,true,1,false>>{});
671 }
672 } else {
673 if( wallModelParameters.wallFunctionProfile == 0 ) {
674 block.addPostProcessor(typeid(stage::PostStream),
675 boundaryLatticeR,
676 meta::id<TurbulentWallModelPostProcessor<false,false,false,0,false>>{});
677 } else {
678 block.addPostProcessor(typeid(stage::PostStream),
679 boundaryLatticeR,
680 meta::id<TurbulentWallModelPostProcessor<false,false,false,1,false>>{});
681 }
682 }
683 }
684 }
685 }
686 if( k == 1 ) {
687 if( wallModelParameters.fNeqMethod == 0 ) {
688 if( wallModelParameters.rhoMethod == 0 ) {
689 block.addPostProcessor(typeid(stage::PostStream),
690 boundaryLatticeR,
691 meta::id<TurbulentWallModelFneqGuoPostProcessor<0>>{});
692 }
693 else if( wallModelParameters.rhoMethod == 1 ) {
694 block.addPostProcessor(typeid(stage::PostStream),
695 boundaryLatticeR,
696 meta::id<TurbulentWallModelFneqGuoPostProcessor<1>>{});
697 }
698 else {
699 block.addPostProcessor(typeid(stage::PostStream),
700 boundaryLatticeR,
701 meta::id<TurbulentWallModelFneqGuoPostProcessor<2>>{});
702 }
703 }
704 else if( wallModelParameters.fNeqMethod == 1 ) {
705 if( wallModelParameters.rhoMethod == 0 ) {
706 block.addPostProcessor(typeid(stage::PostStream),
707 boundaryLatticeR,
708 meta::id<TurbulentWallModelFneqFDMPostProcessor<0>>{});
709 }
710 else if( wallModelParameters.rhoMethod == 1 ) {
711 block.addPostProcessor(typeid(stage::PostStream),
712 boundaryLatticeR,
713 meta::id<TurbulentWallModelFneqFDMPostProcessor<1>>{});
714 }
715 else {
716 block.addPostProcessor(typeid(stage::PostStream),
717 boundaryLatticeR,
718 meta::id<TurbulentWallModelFneqFDMPostProcessor<2>>{});
719 }
720 }
721 else {
722 if( wallModelParameters.rhoMethod == 0 ) {
723 block.addPostProcessor(typeid(stage::PostStream),
724 boundaryLatticeR,
725 meta::id<TurbulentWallModelFneqZeroPostProcessor<0>>{});
726 }
727 else if( wallModelParameters.rhoMethod == 1 ) {
728 block.addPostProcessor(typeid(stage::PostStream),
729 boundaryLatticeR,
730 meta::id<TurbulentWallModelFneqZeroPostProcessor<1>>{});
731 }
732 else {
733 block.addPostProcessor(typeid(stage::PostStream),
734 boundaryLatticeR,
735 meta::id<TurbulentWallModelFneqZeroPostProcessor<2>>{});
736 }
737 }
738 }
739 }
740 }
741 }
742 }
743 }
744 });
745 }
746}
747
748
749
750template<typename T, typename DESCRIPTOR>
751void setWallDistance(SuperLattice<T, DESCRIPTOR>& sLattice,
752 FunctorPtr<SuperIndicatorF<T,DESCRIPTOR::d>>&& boundaryIndicator,
753 FunctorPtr<SuperIndicatorF<T,DESCRIPTOR::d>>&& bulkIndicator,
754 IndicatorF<T,DESCRIPTOR::d>* indicatorAnalyticalBoundary)
755{
756 int _overlap = 3;
757 OstreamManager clout(std::cout, "WallDistanceSetter");
758
759 auto& load = sLattice.getLoadBalancer();
760 for (int iC=0; iC < load.size(); ++iC) {
761 setWallDistance<T,DESCRIPTOR>(sLattice.getBlock(iC),
762 (bulkIndicator->getBlockIndicatorF(iC)).getBlockGeometry(),
763 boundaryIndicator->getBlockIndicatorF(iC),
764 bulkIndicator->getBlockIndicatorF(iC),
765 indicatorAnalyticalBoundary);
766 }
767 addPoints2CommBC<T,DESCRIPTOR>(sLattice, std::forward<decltype(boundaryIndicator)>(boundaryIndicator), _overlap);
768 clout.setMultiOutput(false);
769}
770
772template<typename T, typename DESCRIPTOR>
773void setWallDistance(SuperLattice<T, DESCRIPTOR>& sLattice,
774 SuperGeometry<T,DESCRIPTOR::d>& superGeometry,
775 int materialOfSolidObstacle,
776 IndicatorF<T,DESCRIPTOR::d>* indicatorAnalyticalBoundary,
777 std::vector<int> bulkMaterials = std::vector<int>(1,1))
778{
779 // Getting the indicators by material numbers and calling the superLattice method via the indicators:
780 setWallDistance<T,DESCRIPTOR>(sLattice,
781 FunctorPtr<SuperIndicatorF<T,DESCRIPTOR::d>>(superGeometry.getMaterialIndicator(materialOfSolidObstacle)),
782 FunctorPtr<SuperIndicatorF<T,DESCRIPTOR::d>>(superGeometry.getMaterialIndicator(std::move(bulkMaterials))),
783 indicatorAnalyticalBoundary);
784}
785
787template<typename T, typename DESCRIPTOR>
788void setWallDistance(BlockLattice<T,DESCRIPTOR>& block,
789 BlockGeometry<T,DESCRIPTOR::d>& blockGeometry,
790 BlockIndicatorF<T,DESCRIPTOR::d>& boundaryIndicator,
791 BlockIndicatorF<T,DESCRIPTOR::d>& bulkIndicator,
792 IndicatorF<T,DESCRIPTOR::d>* indicatorAnalyticalBoundary)
793{
794 OstreamManager clout(std::cout, "WallDistanceSetter");
795 clout.setMultiOutput(true);
796
797 // Defining boundary distance y1 and distance to exchange location y2 in wall normal direction for every boundary cell
798 const T deltaR = blockGeometry.getDeltaR();
799 // for each solid cell:
800 block.forSpatialLocations([&](LatticeR<DESCRIPTOR::d> solidLatticeR) {
801 // Check if cell is solid cell
802 if (boundaryIndicator(solidLatticeR)) {
803 for (int iPop=1; iPop < DESCRIPTOR::q; ++iPop) {
804 Vector<T,DESCRIPTOR::d> boundaryLatticeR(solidLatticeR + descriptors::c<DESCRIPTOR>(iPop));
805 if (blockGeometry.getNeighborhoodRadius(boundaryLatticeR) >= 1) {
806 if (blockGeometry.isInside(boundaryLatticeR)) {
807 Vector<T,DESCRIPTOR::d> boundaryPhysR{ };
808 blockGeometry.getPhysR(boundaryPhysR,boundaryLatticeR);
809 // check if neighbor is fluid cell
810 if (bulkIndicator(boundaryLatticeR)) {
811 if (indicatorAnalyticalBoundary) {
812 // Calculate surface normal
813 Vector<T,DESCRIPTOR::d> normal = indicatorAnalyticalBoundary->surfaceNormal(boundaryPhysR, deltaR);
814 T y1 = 0.;
815 // Calculate boundary distance y1 in surface normal direction
816 indicatorAnalyticalBoundary->distance(y1, boundaryPhysR, normal, blockGeometry.getIcGlob());
817 if(y1 == T(0)) {
818 for ( int iD = 0; iD < DESCRIPTOR::d; iD++ ) {
819 normal[iD] *= T(-1);
820 }
821 indicatorAnalyticalBoundary->distance(y1, boundaryPhysR, normal, blockGeometry.getIcGlob());
822 }
823 y1 /= deltaR; // y1 in lattice units
824 if(y1 > T(1)) {
825 y1 = T(1);
826 }
827 block.get(boundaryLatticeR).template setField<descriptors::Y1>(-y1*normal);
828 }
829 }
830 }
831 }
832 }
833 }
834 });
835}
836}
837#endif
olb::AnalyticalConst
AnalyticalConst: DD -> XD, where XD is defined by value.size()
Definition analyticalF.h:101
olb::BlockGeometry
Representation of a block geometry.
Definition blockGeometryStatistics3D.h:50
olb::BlockGeometry::getPhysR
Vector< T, D > getPhysR(LatticeR< D > latticeR) const
Definition blockGeometry.h:129
olb::BlockGeometry::getIcGlob
int getIcGlob() const
Read only access to the global iC number which is given !=-1 if the block geometries are part of a su...
Definition blockGeometry.hh:148
olb::BlockGeometry::getDeltaR
T getDeltaR() const
Read only access to the voxel size given in SI units (meter)
Definition blockGeometry.hh:75
olb::BlockStructureD::getNeighborhoodRadius
CellDistance getNeighborhoodRadius(LatticeR< D > latticeR) const
Return maximum valid neighborhood sphere radius w.r.t. latticeR.
Definition blockStructure.h:223
olb::BlockStructureD::isInside
bool isInside(LatticeR< D > latticeR) const
Return whether location is valid.
Definition blockStructure.h:185
olb::FunctorPtr
Smart pointer for managing the various ways of passing functors around.
Definition functorPtr.h:60
olb::OstreamManager
class for marking output with some text
Definition ostreamManager.h:81
olb::OstreamManager::setMultiOutput
void setMultiOutput(bool b)
enable message output for all MPI processes, disabled by default
Definition ostreamManager.cpp:128
olb::SuperGeometry
Representation of a statistic for a parallel 2D geometry.
Definition superGeometryStatistics3D.h:54
olb::SuperGeometry::getMaterialIndicator
std::unique_ptr< SuperIndicatorF< T, D > > getMaterialIndicator(std::vector< int > &&materials)
Returns a material indicator using the given vector of materials.
Definition superGeometry.hh:476
olb::SuperLattice
Super class maintaining block lattices for a cuboid decomposition.
Definition superBaseF3D.h:47
olb::SuperLattice::getCommunicator
SuperCommunicator< T, SuperLattice > & getCommunicator(STAGE stage=STAGE())
Return communicator for given communication stage.
olb::SuperLattice::getBlock
BlockLattice< T, DESCRIPTOR > & getBlock(int locC)
Return BlockLattice with local index locC.
Definition superLattice.h:115
olb::SuperStructure::getOverlap
int getOverlap()
Read and write access to the overlap.
Definition superStructure.hh:57
olb::SuperStructure::getLoadBalancer
LoadBalancer< T > & getLoadBalancer()
Read and write access to the load balancer.
Definition superStructure.hh:69
olb::Vector
Plain old scalar vector.
Definition vectorHelpers.h:41
olb::descriptors::c
constexpr int c(unsigned iPop, unsigned iDim) any_platform
Definition functions.h:83
olb::util::abs
ADf< T, DIM > abs(const ADf< T, DIM > &a)
Definition aDiff.h:1019
olb::util::sign
int sign(T val) any_platform
Definition util.h:53
olb::util::norm
T norm(const std::vector< T > &a)
l2 norm of a vector of arbitrary length
Definition vectorHelpers.h:122
olb
Top level namespace for all of OpenLB.
Definition advectionDiffusionDirichlet.h:27
olb::addPoints2CommBC
void addPoints2CommBC(SuperLattice< T, DESCRIPTOR > &sLattice, FunctorPtr< SuperIndicatorF2D< T > > &&indicator, int _overlap)
Adds needed Cells to the Communicator _commBC in SuperLattice.
Definition setBoundary2D.h:58
olb::setTurbulentWallModelDynamics
void setTurbulentWallModelDynamics(SuperLattice< T, DESCRIPTOR > &sLattice, FunctorPtr< SuperIndicatorF< T, DESCRIPTOR::d > > &&bulkIndicator, WallModelParameters< T > &wallModelParameters)
Definition setTurbulentWallModel.h:203
olb::BlockIndicatorF
std::conditional_t< D==2, BlockIndicatorF2D< T >, BlockIndicatorF3D< T > > BlockIndicatorF
Definition aliases.h:207
olb::IndicatorF
std::conditional_t< D==2, IndicatorF2D< T >, IndicatorF3D< T > > IndicatorF
Definition aliases.h:247
olb::setWallDistance
void setWallDistance(SuperLattice< T, DESCRIPTOR > &sLattice, FunctorPtr< SuperIndicatorF< T, DESCRIPTOR::d > > &&boundaryIndicator, FunctorPtr< SuperIndicatorF< T, DESCRIPTOR::d > > &&bulkIndicator, IndicatorF< T, DESCRIPTOR::d > *indicatorAnalyticalBoundary)
Definition setTurbulentWallModel.h:751
olb::setTurbulentWallModel
void setTurbulentWallModel(SuperLattice< T, DESCRIPTOR > &sLattice, FunctorPtr< SuperIndicatorF< T, DESCRIPTOR::d > > &&boundaryIndicator, FunctorPtr< SuperIndicatorF< T, DESCRIPTOR::d > > &&bulkIndicator, WallModelParameters< T > &wallModelParameters, IndicatorF< T, DESCRIPTOR::d > *indicatorAnalyticalBoundary=nullptr)
Definition setTurbulentWallModel.h:330
olb::SuperIndicatorF
std::conditional_t< D==2, SuperIndicatorF2D< T >, SuperIndicatorF3D< T > > SuperIndicatorF
Definition aliases.h:197
olb::WallModelParameters::useVanDriest
bool useVanDriest
use van Driest damping function for turbulent viscosity in boundary cell
Definition setTurbulentWallModel.h:184
olb::WallModelParameters::interpolateSampleVelocity
bool interpolateSampleVelocity
interpolate sampling velocity along given normal between lattice voxels
Definition setTurbulentWallModel.h:181
olb::WallModelParameters::bodyForce
bool bodyForce
check if descriptor with body force is used
Definition setTurbulentWallModel.h:178
olb::WallModelParameters::latticeWallDistance
T latticeWallDistance
distance from cell to real wall in lattice units if no geometry indicator is given as input
Definition setTurbulentWallModel.h:187
olb::WallModelParameters::samplingCellDistance
T samplingCellDistance
distance from cell to velocity sampling point in lattice units
Definition setTurbulentWallModel.h:190
olb::collision::HYBRID_RHO
Definition collisionHRR.h:58
olb::collision::ParameterFromCell
Override COLLISION parameter PARAMETER with cell field PARAMETER.
Definition collisionModifiers.h:36
olb::dynamics::Tuple
Dynamics constructed as a tuple of momenta, equilibrium and collision.
Definition interface.h:308
olb::forcing::WFGuoThirdOrder
Dynamics combination rule implementing the forcing scheme by Guo et al. for 3rd order RLB for simulat...
Definition collisionHRR.h:431
olb::meta::id
Identity type to pass non-constructible types as value.
Definition meta.h:79
olb::momenta::BulkDensity
Standard computation for density in the bulk as zeroth moment of the population.
Definition elements.h:150
olb::momenta::BulkStress
Standard stress computation as second moment of the population.
Definition elements.h:1624
olb::momenta::DefineToNEq
When momenta are changed, the equilibrium part of the population is modified while the non-equilibriu...
Definition definitionRule.h:96
olb::stage::PostStream
Communication after propagation.
Definition stages.h:36
olb::util::TensorVal
Compute number of elements of a symmetric d-dimensional tensor.
Definition util.h:216
Generated on Sat May 17 2025 10:47:58 for OpenLB by doxygen 1.12.0