OpenLB 1.8.1
Loading...
Searching...
No Matches
boundaryPostProcessors3D.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2006, 2007 Jonas Latt
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 BOUNDARY_POST_PROCESSORS_3D_HH
25#define BOUNDARY_POST_PROCESSORS_3D_HH
26
27#include "boundaryPostProcessors3D.h"
28
29#include "utilities/finiteDifference3D.h"
30#include "core/util.h"
31
32#include "dynamics/dynamics.h"
33#include "dynamics/lbm.h"
34
35namespace olb {
36
38
39template <typename T, typename DESCRIPTOR, int direction, int orientation>
40template <concepts::DynamicCell CELL>
41void PlaneFdBoundaryProcessor3D<T,DESCRIPTOR,direction,orientation>::apply(CELL& cell)
42{
43 using V = typename CELL::value_t;
44 using namespace olb::util::tensorIndices3D;
45
46 V dx_u[DESCRIPTOR::d], dy_u[DESCRIPTOR::d], dz_u[DESCRIPTOR::d];
47 V rho, u[DESCRIPTOR::d];
48
49 auto& dynamics = cell.getDynamics();
50
51 cell.computeRhoU(rho,u);
52
53 interpolateGradients<0>(cell, dx_u);
54 interpolateGradients<1>(cell, dy_u);
55 interpolateGradients<2>(cell, dz_u);
56
57 V dx_ux = dx_u[0];
58 V dy_ux = dy_u[0];
59 V dz_ux = dz_u[0];
60 V dx_uy = dx_u[1];
61 V dy_uy = dy_u[1];
62 V dz_uy = dz_u[1];
63 V dx_uz = dx_u[2];
64 V dy_uz = dy_u[2];
65 V dz_uz = dz_u[2];
66 V omega = dynamics.getOmegaOrFallback(std::numeric_limits<V>::signaling_NaN());
67 V sToPi = - rho / descriptors::invCs2<V,DESCRIPTOR>() / omega;
68 V pi[util::TensorVal<DESCRIPTOR >::n];
69 pi[xx] = (V)2 * dx_ux * sToPi;
70 pi[yy] = (V)2 * dy_uy * sToPi;
71 pi[zz] = (V)2 * dz_uz * sToPi;
72 pi[xy] = (dx_uy + dy_ux) * sToPi;
73 pi[xz] = (dx_uz + dz_ux) * sToPi;
74 pi[yz] = (dy_uz + dz_uy) * sToPi;
75
76 // Computation of the particle distribution functions
77 // according to the regularized formula
78 V fEq[DESCRIPTOR::q] { };
79 dynamics.computeEquilibrium(cell, rho, u, fEq);
80 for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
81 cell[iPop] = fEq[iPop] + equilibrium<DESCRIPTOR>::template fromPiToFneq<V>(iPop, pi);
82 }
83}
84
85template <typename T, typename DESCRIPTOR, int direction, int orientation>
86template <int deriveDirection, typename CELL, typename V>
87void PlaneFdBoundaryProcessor3D<T,DESCRIPTOR,direction,orientation>::interpolateGradients(
88 CELL& cell, V velDeriv[DESCRIPTOR::d]) const
89{
90 fd::DirectedGradients3D<V,DESCRIPTOR,direction,orientation,deriveDirection,direction==deriveDirection>
91 ::interpolateVector(velDeriv, cell);
92}
93
94template <typename T, typename DESCRIPTOR, int direction, int orientation>
95template <concepts::DynamicCell CELL>
96void StraightConvectionBoundaryProcessor3D<T,DESCRIPTOR,direction,orientation>::initialize(CELL& cell)
97{
98 constexpr auto missing = util::populationsContributingToVelocity<DESCRIPTOR,direction,-orientation>();
99 auto prevCell = cell.template getFieldPointer<PREV_CELL>();
100 for (unsigned i=0; i < missing.size(); ++i) {
101 prevCell[i] = cell[missing[i]];
102 }
103}
104
105template <typename T, typename DESCRIPTOR, int direction, int orientation>
106template <concepts::DynamicCell CELL>
107void StraightConvectionBoundaryProcessor3D<T,DESCRIPTOR,direction,orientation>::apply(CELL& cell)
108{
109 using V = typename CELL::value_t;
110 constexpr auto missing = util::populationsContributingToVelocity<DESCRIPTOR,direction,-orientation>();
111
112 auto prevCell = cell.template getField<PREV_CELL>();
113
114 for (unsigned i=0; i < missing.size(); ++i) {
115 cell[missing[i]] = prevCell[i];
116 }
117
118 V rho0, u0[3];
119 V rho1, u1[3];
120 V rho2, u2[3];
121
122 cell.computeRhoU(rho0, u0);
123
124 static_assert(direction == 0 || direction == 1 || direction ==2,
125 "Direction must be one of 0, 1 or 2");
126 if constexpr (direction == 0) {
127 cell.neighbor({-orientation ,0,0}).computeRhoU(rho1, u1);
128 cell.neighbor({-orientation*2,0,0}).computeRhoU(rho2, u2);
129 }
130 else if constexpr (direction == 1) {
131 cell.neighbor({0,-orientation ,0}).computeRhoU(rho1, u1);
132 cell.neighbor({0,-orientation*2,0}).computeRhoU(rho2, u2);
133 }
134 else if constexpr (direction == 2) {
135 cell.neighbor({0,0,-orientation }).computeRhoU(rho1, u1);
136 cell.neighbor({0,0,-orientation*2}).computeRhoU(rho2, u2);
137 }
138
139 V uDelta[3];
140 V uAverage = rho0*u0[direction];
141 //if (uAv!=nullptr) {
142 // rho0 = V{1};
143 // rho1 = V{1};
144 // rho2 = V{1};
145 // uAverage = *uAv * rho0;
146 //}
147
148 uDelta[0] = -uAverage*0.5*(3*rho0*u0[0]-4*rho1*u1[0]+rho2*u2[0]);
149 uDelta[1] = -uAverage*0.5*(3*rho0*u0[1]-4*rho1*u1[1]+rho2*u2[1]);
150 uDelta[2] = -uAverage*0.5*(3*rho0*u0[2]-4*rho1*u1[2]+rho2*u2[2]);
151
152 for (unsigned i=0; i < missing.size(); ++i) {
153 auto iPop = missing[i];
154 prevCell[i] = cell[iPop]
155 + descriptors::invCs2<V,DESCRIPTOR>()*descriptors::t<V,DESCRIPTOR>(iPop)
156 * ( uDelta[0]*descriptors::c<DESCRIPTOR>(iPop,0)
157 + uDelta[1]*descriptors::c<DESCRIPTOR>(iPop,1)
158 + uDelta[2]*descriptors::c<DESCRIPTOR>(iPop,2));
159 }
160
161 cell.template setField<PREV_CELL>(prevCell);
162}
163
165
166template <typename T, typename DESCRIPTOR, int plane, int normal1, int normal2>
167template <concepts::DynamicCell CELL>
168void OuterVelocityEdgeProcessor3D<T,DESCRIPTOR,plane,normal1,normal2>::apply(CELL& cell)
169{
170 using V = typename CELL::value_t;
171 using namespace olb::util::tensorIndices3D;
172
173 constexpr auto direction1 = (plane+1)%3;
174 constexpr auto direction2 = (plane+2)%3;
175
176 auto& dynamics = cell.getDynamics();
177
178 V rho10 = getNeighborRho(cell, 1,0);
179 V rho01 = getNeighborRho(cell, 0,1);
180 V rho20 = getNeighborRho(cell, 2,0);
181 V rho02 = getNeighborRho(cell, 0,2);
182 V rho = (V)2/(V)3*(rho01+rho10)-(V)1/(V)6*(rho02+rho20);
183
184 V dA_uB_[3][3];
185
186 interpolateGradients<plane,0> (cell, dA_uB_[0]);
187 interpolateGradients<direction1,normal1>(cell, dA_uB_[1]);
188 interpolateGradients<direction2,normal2>(cell, dA_uB_[2]);
189
190 V dA_uB[3][3];
191 for (int iBeta=0; iBeta<3; ++iBeta) {
192 dA_uB[plane][iBeta] = dA_uB_[0][iBeta];
193 dA_uB[direction1][iBeta] = dA_uB_[1][iBeta];
194 dA_uB[direction2][iBeta] = dA_uB_[2][iBeta];
195 }
196 V omega = dynamics.getOmegaOrFallback(std::numeric_limits<V>::signaling_NaN());
197 V sToPi = - rho / descriptors::invCs2<V,DESCRIPTOR>() / omega;
198 V pi[util::TensorVal<DESCRIPTOR>::n];
199 pi[xx] = (V)2 * dA_uB[0][0] * sToPi;
200 pi[yy] = (V)2 * dA_uB[1][1] * sToPi;
201 pi[zz] = (V)2 * dA_uB[2][2] * sToPi;
202 pi[xy] = (dA_uB[0][1]+dA_uB[1][0]) * sToPi;
203 pi[xz] = (dA_uB[0][2]+dA_uB[2][0]) * sToPi;
204 pi[yz] = (dA_uB[1][2]+dA_uB[2][1]) * sToPi;
205
206 // Computation of the particle distribution functions
207 // according to the regularized formula
208 V u[DESCRIPTOR::d];
209 cell.computeU(u);
210
211 V fEq[DESCRIPTOR::q] { };
212 dynamics.computeEquilibrium(cell, rho, u, fEq);
213 for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
214 cell[iPop] = fEq[iPop] + equilibrium<DESCRIPTOR>::template fromPiToFneq<V>(iPop, pi);
215 }
216}
217
218template <typename T, typename DESCRIPTOR, int plane, int normal1, int normal2>
219template <typename CELL, typename V>
220V OuterVelocityEdgeProcessor3D<T,DESCRIPTOR,plane,normal1,normal2>
221 ::getNeighborRho(CELL& cell, int step1, int step2)
222{
223 constexpr auto direction1 = (plane+1)%3;
224 constexpr auto direction2 = (plane+2)%3;
225 int coords[3] { };
226 coords[direction1] = -normal1*step1;
227 coords[direction2] = -normal2*step2;
228 return cell.neighbor(coords).computeRho();
229}
230
231template <typename T, typename DESCRIPTOR, int plane, int normal1, int normal2>
232template <int deriveDirection, int orientation, typename CELL, typename V>
233void OuterVelocityEdgeProcessor3D<T,DESCRIPTOR,plane,normal1,normal2>
234 ::interpolateGradients(CELL& cell, V velDeriv[DESCRIPTOR::d]) const
235{
236 fd::DirectedGradients3D<V,DESCRIPTOR,deriveDirection,orientation,deriveDirection,deriveDirection!=plane>
237 ::interpolateVector(velDeriv, cell);
238}
239
241
242template <typename T, typename DESCRIPTOR, int xNormal, int yNormal, int zNormal>
243template <concepts::DynamicCell CELL>
244void OuterVelocityCornerProcessor3D<T,DESCRIPTOR,xNormal,yNormal,zNormal>::apply(CELL& cell)
245{
246 using V = typename CELL::value_t;
247 using namespace olb::util::tensorIndices3D;
248
249 auto& dynamics = cell.getDynamics();
250
251 V rho100 = cell.neighbor({-1*xNormal, -0*yNormal, -0*zNormal}).computeRho();
252 V rho010 = cell.neighbor({-0*xNormal, -1*yNormal, -0*zNormal}).computeRho();
253 V rho001 = cell.neighbor({-0*xNormal, -0*yNormal, -1*zNormal}).computeRho();
254 V rho200 = cell.neighbor({-2*xNormal, -0*yNormal, -0*zNormal}).computeRho();
255 V rho020 = cell.neighbor({-0*xNormal, -2*yNormal, -0*zNormal}).computeRho();
256 V rho002 = cell.neighbor({-0*xNormal, -0*yNormal, -2*zNormal}).computeRho();
257
258 V rho = (V)4/(V)9 * (rho001 + rho010 + rho100) - (V)1/(V)9 * (rho002 + rho020 + rho200);
259
260 V dx_u[DESCRIPTOR::d], dy_u[DESCRIPTOR::d], dz_u[DESCRIPTOR::d];
261 fd::DirectedGradients3D<V, DESCRIPTOR, 0, xNormal, 0, true>::interpolateVector(dx_u, cell);
262 fd::DirectedGradients3D<V, DESCRIPTOR, 1, yNormal, 0, true>::interpolateVector(dy_u, cell);
263 fd::DirectedGradients3D<V, DESCRIPTOR, 2, zNormal, 0, true>::interpolateVector(dz_u, cell);
264
265 V dx_ux = dx_u[0];
266 V dy_ux = dy_u[0];
267 V dz_ux = dz_u[0];
268 V dx_uy = dx_u[1];
269 V dy_uy = dy_u[1];
270 V dz_uy = dz_u[1];
271 V dx_uz = dx_u[2];
272 V dy_uz = dy_u[2];
273 V dz_uz = dz_u[2];
274 V omega = dynamics.getOmegaOrFallback(std::numeric_limits<V>::signaling_NaN());
275 V sToPi = - rho / descriptors::invCs2<V,DESCRIPTOR>() / omega;
276 V pi[util::TensorVal<DESCRIPTOR >::n];
277 pi[xx] = (V)2 * dx_ux * sToPi;
278 pi[yy] = (V)2 * dy_uy * sToPi;
279 pi[zz] = (V)2 * dz_uz * sToPi;
280 pi[xy] = (dx_uy + dy_ux) * sToPi;
281 pi[xz] = (dx_uz + dz_ux) * sToPi;
282 pi[yz] = (dy_uz + dz_uy) * sToPi;
283
284 // Computation of the particle distribution functions
285 // according to the regularized formula
286 V u[DESCRIPTOR::d];
287 cell.computeU(u);
288
289 V fEq[DESCRIPTOR::q] { };
290 dynamics.computeEquilibrium(cell, rho, u, fEq);
291 for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
292 cell[iPop] = fEq[iPop] + equilibrium<DESCRIPTOR>::template fromPiToFneq<V>(iPop, pi);
293 }
294}
295
297
298template<typename T, typename DESCRIPTOR>
299SlipBoundaryProcessor3D<T,DESCRIPTOR>::
300SlipBoundaryProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX, int discreteNormalY, int discreteNormalZ)
301 : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_)
302{
303 this->getName() = "SlipBoundaryProcessor3D";
304 OLB_PRECONDITION(x0==x1 || y0==y1 || z0==z1);
305 int mirrorDirection0;
306 int mirrorDirection1;
307 int mirrorDirection2;
308 int mult = 2 / (discreteNormalX*discreteNormalX + discreteNormalY*discreteNormalY + discreteNormalZ*discreteNormalZ);
309 reflectionPop[0] = 0;
310 for (int iPop = 1; iPop < DESCRIPTOR::q; iPop++) {
311 reflectionPop[iPop] = 0;
312 // iPop are the directions which pointing into the fluid, discreteNormal is pointing outwarts
313 int scalarProduct = descriptors::c<DESCRIPTOR>(iPop,0)*discreteNormalX + descriptors::c<DESCRIPTOR>(iPop,1)*discreteNormalY + descriptors::c<DESCRIPTOR>(iPop,2)*discreteNormalZ;
314 if ( scalarProduct < 0) {
315 // bounce back for the case discreteNormalX = discreteNormalY = discreteNormalZ = 1, that is mult=0
316 if (mult == 0) {
317 mirrorDirection0 = -descriptors::c<DESCRIPTOR>(iPop,0);
318 mirrorDirection1 = -descriptors::c<DESCRIPTOR>(iPop,1);
319 mirrorDirection2 = -descriptors::c<DESCRIPTOR>(iPop,2);
320 }
321 else {
322 mirrorDirection0 = descriptors::c<DESCRIPTOR>(iPop,0) - mult*scalarProduct*discreteNormalX;
323 mirrorDirection1 = descriptors::c<DESCRIPTOR>(iPop,1) - mult*scalarProduct*discreteNormalY;
324 mirrorDirection2 = descriptors::c<DESCRIPTOR>(iPop,2) - mult*scalarProduct*discreteNormalZ;
325 }
326
327 // run through all lattice directions and look for match of direction
328 for (int i = 1; i < DESCRIPTOR::q; i++) {
329 if (descriptors::c<DESCRIPTOR>(i,0)==mirrorDirection0
330 && descriptors::c<DESCRIPTOR>(i,1)==mirrorDirection1
331 && descriptors::c<DESCRIPTOR>(i,2)==mirrorDirection2) {
332 reflectionPop[iPop] = i;
333 break;
334 }
335 }
336 }
337 }
338}
339
340template<typename T, typename DESCRIPTOR>
341void SlipBoundaryProcessor3D<T,DESCRIPTOR>::
342processSubDomain(BlockLattice<T,DESCRIPTOR>& blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
343{
344 int newX0, newX1, newY0, newY1, newZ0, newZ1;
345 if ( util::intersect (
346 x0, x1, y0, y1, z0, z1,
347 x0_, x1_, y0_, y1_, z0_, z1_,
348 newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) {
349
350 int iX;
351#ifdef PARALLEL_MODE_OMP
352 #pragma omp parallel for
353#endif
354 for (iX=x0; iX<=x1; ++iX) {
355 for (int iY=y0; iY<=y1; ++iY) {
356 for (int iZ=z0; iZ<=z1; ++iZ) {
357 for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) {
358 if (reflectionPop[iPop]!=0) {
359 //do reflection
360 blockLattice.get(iX,iY,iZ)[iPop] = blockLattice.get(iX,iY,iZ)[reflectionPop[iPop]];
361 }
362 }
363 }
364 }
365 }
366 }
367}
368
369template<typename T, typename DESCRIPTOR>
370void SlipBoundaryProcessor3D<T,DESCRIPTOR>::
371process(BlockLattice<T,DESCRIPTOR>& blockLattice)
372{
373 processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
374}
375
377
378template<typename T, typename DESCRIPTOR>
379SlipBoundaryProcessorGenerator3D<T,DESCRIPTOR>::
380SlipBoundaryProcessorGenerator3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX_, int discreteNormalY_, int discreteNormalZ_)
381 : PostProcessorGenerator3D<T,DESCRIPTOR>(x0_, x1_, y0_, y1_, z0_, z1_), discreteNormalX(discreteNormalX_), discreteNormalY(discreteNormalY_), discreteNormalZ(discreteNormalZ_)
382{ }
383
384template<typename T, typename DESCRIPTOR>
385PostProcessor3D<T,DESCRIPTOR>*
386SlipBoundaryProcessorGenerator3D<T,DESCRIPTOR>::generate() const
387{
388 return new SlipBoundaryProcessor3D<T,DESCRIPTOR>
389 ( this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, discreteNormalX, discreteNormalY, discreteNormalZ);
390}
391
392template<typename T, typename DESCRIPTOR>
393PostProcessorGenerator3D<T,DESCRIPTOR>*
394SlipBoundaryProcessorGenerator3D<T,DESCRIPTOR>::clone() const
395{
396 return new SlipBoundaryProcessorGenerator3D<T,DESCRIPTOR>
397 (this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, discreteNormalX, discreteNormalY, discreteNormalZ);
398}
399
401
402template<typename T, typename DESCRIPTOR>
403PartialSlipBoundaryProcessor3D<T,DESCRIPTOR>::
404PartialSlipBoundaryProcessor3D(T tuner_, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX, int discreteNormalY, int discreteNormalZ)
405 : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), tuner(tuner_)
406{
407 this->getName() = "PartialSlipBoundaryProcessor3D";
408 OLB_PRECONDITION(x0==x1 || y0==y1 || z0==z1);
409 reflectionPop[0] = 0;
410 for (int iPop = 1; iPop < DESCRIPTOR::q; iPop++) {
411 reflectionPop[iPop] = 0;
412 // iPop are the directions which pointing into the fluid, discreteNormal is pointing outwarts
413 if (descriptors::c<DESCRIPTOR>(iPop,0)*discreteNormalX + descriptors::c<DESCRIPTOR>(iPop,1)*discreteNormalY + descriptors::c<DESCRIPTOR>(iPop,2)*discreteNormalZ < 0) {
414 //std::cout << "-----" <<std::endl;
415 int mirrorDirection0;
416 int mirrorDirection1;
417 int mirrorDirection2;
418 int mult = 2 / (discreteNormalX*discreteNormalX + discreteNormalY*discreteNormalY + discreteNormalZ*discreteNormalZ);
419
420 mirrorDirection0 = (descriptors::c<DESCRIPTOR>(iPop,0) - mult*(descriptors::c<DESCRIPTOR>(iPop,0)*discreteNormalX + descriptors::c<DESCRIPTOR>(iPop,1)*discreteNormalY + descriptors::c<DESCRIPTOR>(iPop,2)*discreteNormalZ)*discreteNormalX);
421
422 mirrorDirection1 = (descriptors::c<DESCRIPTOR>(iPop,1) - mult*(descriptors::c<DESCRIPTOR>(iPop,0)*discreteNormalX + descriptors::c<DESCRIPTOR>(iPop,1)*discreteNormalY + descriptors::c<DESCRIPTOR>(iPop,2)*discreteNormalZ)*discreteNormalY);
423 mirrorDirection2 = (descriptors::c<DESCRIPTOR>(iPop,2) - mult*(descriptors::c<DESCRIPTOR>(iPop,0)*discreteNormalX + descriptors::c<DESCRIPTOR>(iPop,1)*discreteNormalY + descriptors::c<DESCRIPTOR>(iPop,2)*discreteNormalZ)*discreteNormalZ);
424
425 // bounce back for the case discreteNormalX = discreteNormalY = discreteNormalZ = 1, that is mult=0
426 if (mult == 0) {
427 mirrorDirection0 = -descriptors::c<DESCRIPTOR>(iPop,0);
428 mirrorDirection1 = -descriptors::c<DESCRIPTOR>(iPop,1);
429 mirrorDirection2 = -descriptors::c<DESCRIPTOR>(iPop,2);
430 }
431
432 // computes mirror jPop
433 for (reflectionPop[iPop] = 1; reflectionPop[iPop] < DESCRIPTOR::q ; reflectionPop[iPop]++) {
434 if (descriptors::c<DESCRIPTOR>(reflectionPop[iPop],0)==mirrorDirection0 && descriptors::c<DESCRIPTOR>(reflectionPop[iPop],1)==mirrorDirection1 && descriptors::c<DESCRIPTOR>(reflectionPop[iPop],2)==mirrorDirection2) {
435 break;
436 }
437 }
438 //std::cout <<iPop << " to "<< jPop <<" for discreteNormal= "<< discreteNormalX << "/"<<discreteNormalY <<std::endl;
439 }
440 }
441}
442
443template<typename T, typename DESCRIPTOR>
444void PartialSlipBoundaryProcessor3D<T,DESCRIPTOR>::
445processSubDomain(BlockLattice<T,DESCRIPTOR>& blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
446{
447 int newX0, newX1, newY0, newY1, newZ0, newZ1;
448
449 if ( util::intersect (
450 x0, x1, y0, y1, z0, z1,
451 x0_, x1_, y0_, y1_, z0_, z1_,
452 newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) {
453
454 int iX;
455#ifdef PARALLEL_MODE_OMP
456 #pragma omp parallel for
457#endif
458 for (iX=newX0; iX<=newX1; ++iX) {
459 for (int iY=newY0; iY<=newY1; ++iY) {
460 for (int iZ=newZ0; iZ<=newZ1; ++iZ) {
461 for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) {
462 if (reflectionPop[iPop]!=0) {
463 //do reflection
464 blockLattice.get(iX,iY,iZ)[iPop] = tuner*blockLattice.get(iX,iY,iZ)[reflectionPop[iPop]];
465 }
466 }
467 for (int iPop = 1; iPop < DESCRIPTOR::q/2 ; ++iPop) {
468 T provv = blockLattice.get(iX,iY,iZ)[descriptors::opposite<DESCRIPTOR>(iPop)];
469 blockLattice.get(iX,iY,iZ)[descriptors::opposite<DESCRIPTOR>(iPop)] +=
470 (1.-tuner)*blockLattice.get(iX,iY,iZ)[iPop];
471 blockLattice.get(iX,iY,iZ)[iPop] += (1.-tuner)*provv;
472 }
473 }
474 }
475 }
476 }
477}
478
479template<typename T, typename DESCRIPTOR>
480void PartialSlipBoundaryProcessor3D<T,DESCRIPTOR>::
481process(BlockLattice<T,DESCRIPTOR>& blockLattice)
482{
483 processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
484}
485
487
488template<typename T, typename DESCRIPTOR>
489PartialSlipBoundaryProcessorGenerator3D<T,DESCRIPTOR>::
490PartialSlipBoundaryProcessorGenerator3D(T tuner_, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX_, int discreteNormalY_, int discreteNormalZ_)
491 : PostProcessorGenerator3D<T,DESCRIPTOR>(x0_, x1_, y0_, y1_, z0_, z1_), discreteNormalX(discreteNormalX_), discreteNormalY(discreteNormalY_), discreteNormalZ(discreteNormalZ_), tuner(tuner_)
492{ }
493
494template<typename T, typename DESCRIPTOR>
495PostProcessor3D<T,DESCRIPTOR>*
496PartialSlipBoundaryProcessorGenerator3D<T,DESCRIPTOR>::generate() const
497{
498 return new PartialSlipBoundaryProcessor3D<T,DESCRIPTOR>
499 (tuner, this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, discreteNormalX, discreteNormalY, discreteNormalZ);
500}
501
502template<typename T, typename DESCRIPTOR>
503PostProcessorGenerator3D<T,DESCRIPTOR>*
504PartialSlipBoundaryProcessorGenerator3D<T,DESCRIPTOR>::clone() const
505{
506 return new PartialSlipBoundaryProcessorGenerator3D<T,DESCRIPTOR>
507 (tuner, this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, discreteNormalX, discreteNormalY, discreteNormalZ);
508}
509
511template<typename T, typename DESCRIPTOR, int NORMAL_X, int NORMAL_Y, int NORMAL_Z>
512template <concepts::DynamicCell CELL, typename PARAMETERS>
513void FreeEnergyWallProcessor3D<T,DESCRIPTOR,NORMAL_X,NORMAL_Y,NORMAL_Z>::apply(CELL& cell, PARAMETERS& vars) {
514
515 T addend = cell.template getField<descriptors::ADDEND>();
516
517 T rhoBulk = cell.neighbor({-NORMAL_X, -NORMAL_Y, -NORMAL_Z}).computeRho();
518 T rhoTmp = 0.;
519
520 for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) {
521 rhoTmp += cell[iPop];
522 }
523
524 T rho = rhoBulk + addend;
525 rho -= rhoTmp;
526
527 cell[0] = rho - 1.;
528
529}
530
531
533template<typename T, typename DESCRIPTOR, int NORMAL_X, int NORMAL_Y, int NORMAL_Z>
534template <concepts::DynamicCell CELL>
535void FreeEnergyChemPotBoundaryProcessor3DA<T,DESCRIPTOR,NORMAL_X,NORMAL_Y,NORMAL_Z>::apply(CELL& cell) {
536
537 cell.template setField<descriptors::CHEM_POTENTIAL>(
538 cell.neighbor({-NORMAL_X,-NORMAL_Y,-NORMAL_Z}).template getField<descriptors::CHEM_POTENTIAL>());
539
540 T rho0 = cell.computeRho();
541 T rho1 = cell.neighbor({-NORMAL_X,-NORMAL_Y,-NORMAL_Z}).computeRho();
542
543 cell.template setField<descriptors::CHEM_POTENTIAL>(
544 cell.template getField<descriptors::CHEM_POTENTIAL>() + (rho1 / rho0 - 1) / descriptors::invCs2<T,DESCRIPTOR>());
545}
546
547template<typename T, typename DESCRIPTOR, int NORMAL_X, int NORMAL_Y, int NORMAL_Z>
548template <concepts::DynamicCell CELL>
549void FreeEnergyChemPotBoundaryProcessor3DB<T,DESCRIPTOR,NORMAL_X,NORMAL_Y,NORMAL_Z>::apply(CELL& cell) {
550
551 cell.template setField<descriptors::CHEM_POTENTIAL>(
552 cell.neighbor({-NORMAL_X,-NORMAL_Y,-NORMAL_Z}).template getField<descriptors::CHEM_POTENTIAL>());
553}
554
555
557template<typename T, typename DESCRIPTOR, int NORMAL_X, int NORMAL_Y, int NORMAL_Z>
558template <concepts::DynamicCell CELL>
559void FreeEnergyConvectiveProcessor3D<T,DESCRIPTOR,NORMAL_X,NORMAL_Y,NORMAL_Z>::apply(CELL& cell) {
560
561 T rho, rho0, rho1, u[3];
562
563 rho0 = cell.computeRho();
564
565 cell.neighbor({-NORMAL_X,-NORMAL_Y,-NORMAL_Z}).computeRhoU(rho1, u);
566
567 T uPerp = 0;
568
569 Vector<T, 3> normalVec({NORMAL_X,NORMAL_Y,NORMAL_Z});
570
571 if (normalVec[2] == 0) {
572 if (normalVec[1] == 0) {
573 if (normalVec[0] < 0) {
574 uPerp = -u[0];
575 } else {
576 uPerp = u[0];
577 }
578 } else if (normalVec[0] == 0) {
579 if (normalVec[1] < 0) {
580 uPerp = -u[1];
581 } else {
582 uPerp = u[1];
583 }
584 } else {
585 uPerp = util::sqrt(u[0] * u[0] + u[1] * u[1]);
586 }
587 } else if (normalVec[1] == 0) {
588 if (normalVec[0] == 0) {
589 if (normalVec[2] < 0) {
590 uPerp = -u[2];
591 } else {
592 uPerp = u[2];
593 }
594 } else {
595 uPerp = util::sqrt(u[0] * u[0] + u[2] * u[2]);
596 }
597 } else if (normalVec[0] == 0) {
598 uPerp = util::sqrt(u[1] * u[1] + u[2] * u[2]);
599 } else {
600 uPerp = util::sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]);
601 }
602
603 rho = (rho0 + uPerp * rho1) / (1. + uPerp);
604 cell.defineRho(rho);
605}
606
607
608} // namespace olb
609
610#endif
olb::FreeEnergyWallProcessor3D::apply
void apply(CELL &cell, PARAMETERS &parameters) any_platform
Definition boundaryPostProcessors3D.hh:513
olb::OuterVelocityEdgeProcessor3D
This class computes the skordos BC on a convex edge wall in 3D but with a limited number of terms add...
Definition boundaryPostProcessors3D.h:88
olb::PartialSlipBoundaryProcessor3D
This class computes a partial slip BC in 3D.
Definition boundaryPostProcessors3D.h:164
olb::PartialSlipBoundaryProcessor3D::PartialSlipBoundaryProcessor3D
PartialSlipBoundaryProcessor3D(T tuner_, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX_, int discreteNormalY_, int discreteNormalZ_)
Definition boundaryPostProcessors3D.hh:404
olb::PartialSlipBoundaryProcessor3D::processSubDomain
void processSubDomain(BlockLattice< T, DESCRIPTOR > &blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) override
Execute post-processing step on a sublattice.
Definition boundaryPostProcessors3D.hh:445
olb::PartialSlipBoundaryProcessor3D::process
void process(BlockLattice< T, DESCRIPTOR > &blockLattice) override
Execute post-processing step.
Definition boundaryPostProcessors3D.hh:481
olb::PartialSlipBoundaryProcessorGenerator3D::PartialSlipBoundaryProcessorGenerator3D
PartialSlipBoundaryProcessorGenerator3D(T tuner_, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX_, int discreteNormalY_, int discreteNormalZ_)
Definition boundaryPostProcessors3D.hh:490
olb::PartialSlipBoundaryProcessorGenerator3D::generate
PostProcessor3D< T, DESCRIPTOR > * generate() const override
Definition boundaryPostProcessors3D.hh:496
olb::PartialSlipBoundaryProcessorGenerator3D::clone
PostProcessorGenerator3D< T, DESCRIPTOR > * clone() const override
Definition boundaryPostProcessors3D.hh:504
olb::PlaneFdBoundaryProcessor3D
This class computes the skordos BC on a plane wall in 3D but with a limited number of terms added to ...
Definition boundaryPostProcessors3D.h:42
olb::PlaneFdBoundaryProcessor3D::apply
void apply(CELL &cell) any_platform
Definition boundaryPostProcessors3D.hh:41
olb::PostProcessor3D::getName
std::string & getName()
read and write access to name
Definition postProcessing.hh:192
olb::SlipBoundaryProcessor3D
This class computes a slip BC in 3D.
Definition boundaryPostProcessors3D.h:127
olb::SlipBoundaryProcessor3D::SlipBoundaryProcessor3D
SlipBoundaryProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX_, int discreteNormalY_, int discreteNormalZ_)
Definition boundaryPostProcessors3D.hh:300
olb::SlipBoundaryProcessor3D::processSubDomain
void processSubDomain(BlockLattice< T, DESCRIPTOR > &blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) override
Execute post-processing step on a sublattice.
Definition boundaryPostProcessors3D.hh:342
olb::SlipBoundaryProcessor3D::process
void process(BlockLattice< T, DESCRIPTOR > &blockLattice) override
Execute post-processing step.
Definition boundaryPostProcessors3D.hh:371
olb::SlipBoundaryProcessorGenerator3D::clone
PostProcessorGenerator3D< T, DESCRIPTOR > * clone() const override
Definition boundaryPostProcessors3D.hh:394
olb::SlipBoundaryProcessorGenerator3D::SlipBoundaryProcessorGenerator3D
SlipBoundaryProcessorGenerator3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX_, int discreteNormalY_, int discreteNormalZ_)
Definition boundaryPostProcessors3D.hh:380
olb::SlipBoundaryProcessorGenerator3D::generate
PostProcessor3D< T, DESCRIPTOR > * generate() const override
Definition boundaryPostProcessors3D.hh:386
olb::Vector
Plain old scalar vector.
Definition vectorHelpers.h:41
olb::descriptors::invCs2
constexpr T invCs2() any_platform
Definition functions.h:107
olb::descriptors::t
constexpr T t(unsigned iPop, tag::CUM) any_platform
Definition cum.h:108
olb::descriptors::c
constexpr int c(unsigned iPop, unsigned iDim) any_platform
Definition functions.h:83
olb::descriptors::opposite
constexpr int opposite(unsigned iPop) any_platform
Definition functions.h:95
olb::util::sqrt
Expr sqrt(Expr x)
Definition expr.cpp:225
olb::util::populationsContributingToVelocity
constexpr auto populationsContributingToVelocity() any_platform
Return array of population indices where c[iVel] == value.
Definition util.h:222
olb::util::intersect
bool intersect(int x0, int x1, int y0, int y1, int x0_, int x1_, int y0_, int y1_, int &newX0, int &newX1, int &newY0, int &newY1)
Definition util.h:85
olb
Top level namespace for all of OpenLB.
Definition advectionDiffusionDirichlet.h:27
OLB_PRECONDITION
#define OLB_PRECONDITION(COND)
Definition olbDebug.h:46
olb::fd::DirectedGradients3D::interpolateVector
static void interpolateVector(T velDeriv[DESCRIPTOR::d], BlockLattice< T, DESCRIPTOR > const &blockLattice, int iX, int iY, int iZ)
olb::util::TensorVal
Compute number of elements of a symmetric d-dimensional tensor.
Definition util.h:216
util.h
Set of functions commonly used in LB computations – header file.
Generated on Sat May 17 2025 10:47:58 for OpenLB by doxygen 1.12.0