24#ifndef WALLFUNCTION_BOUNDARY_POST_PROCESSORS_3D_HH
25#define WALLFUNCTION_BOUNDARY_POST_PROCESSORS_3D_HH
36template <
typename T,
typename S>
42template <
typename T,
typename S>
48 T b = 0.119760479041916168;
49 T c = 0.488023952095808383;
51 T e = 3.50727901936264842;
58 output[0] = y_plus *
util::sqrt(tau_w[0]/_rho);
64template <
typename T,
typename S>
65PowerLawProfile<T,S>::PowerLawProfile(T nu, T u2, T y2, T y1, T rho) :
AnalyticalF<1,T,S>(2), _nu(nu), _u2(u2), _y2(y2), _y1(y1), _rho(rho)
67 this->
getName() =
"PowerLawProfile";
70template <
typename T,
typename S>
75 T y_plus = _y1 * u_tau / _nu;
79 output[1] = u_tau * 8.3 *
util::pow(y_plus, 1./7.);
81 else if (y_plus < 30.0 && y_plus > 5.) {
83 output[1] = u_tau * (-2.81742 + 4.85723 *
util::log(y_plus));
86 output[0] = 2.*_u2 * _rho * _nu / _y2;
91 y_plus = _y1 * u_tau / _nu;
92 output[1] = u_tau * y_plus;
99template<
typename T,
typename DESCRIPTOR>
102 std::vector<int> discreteNormal, std::vector<int> missingIndices,
105 : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_),
106 _blockGeometryStructure(blockGeometryStructure),
107 _discreteNormal(discreteNormal), _missingIndices(missingIndices),
110 this->
getName() =
"WallFunctionBoundaryProcessor3D";
112 discreteNormalX = _discreteNormal[0];
113 discreteNormalY = _discreteNormal[1];
114 discreteNormalZ = _discreteNormal[2];
116 int Type_BC = discreteNormalX*discreteNormalX + discreteNormalY*discreteNormalY + discreteNormalZ*discreteNormalZ;
119 if (discreteNormalX != 0) {
120 orientation = discreteNormalX;
123 else if (discreteNormalY != 0) {
124 orientation = discreteNormalY;
127 else if (discreteNormalZ != 0) {
128 orientation = discreteNormalZ;
132 else if (Type_BC == 2) {
136 else if (Type_BC == 3) {
141 if (geoIndicator == NULL) {
144 if (_wallFunctionParam.latticeWalldistance > 0.) {
145 y_1 = _wallFunctionParam.latticeWalldistance;
150 unit_normal[0] = discreteNormalX / normal_norm;
151 unit_normal[1] = discreteNormalY / normal_norm;
152 unit_normal[2] = discreteNormalZ / normal_norm;
155 calculateWallDistances(geoIndicator);
156 if (y_1 > 2. || std::isnan(unit_normal[0]) || std::isnan(unit_normal[1]) || std::isnan(unit_normal[2])) {
158 if (_wallFunctionParam.latticeWalldistance > 0.) {
159 y_1 = _wallFunctionParam.latticeWalldistance;
164 unit_normal[0] = discreteNormalX / normal_norm;
165 unit_normal[1] = discreteNormalY / normal_norm;
166 unit_normal[2] = discreteNormalZ / normal_norm;
169 y_2 = y_1 + normal_norm;
172 getIndices(direction, 0, onWallIndices);
173 getIndices(direction, orientation, normalIndices);
178 getIndices(direction, -orientation, normalInwardsIndices);
181template<
typename T,
typename DESCRIPTOR>
184 int newX0, newX1, newY0, newY1, newZ0, newZ1;
186 x0, x1, y0, y1, z0, z1,
187 x0_, x1_, y0_, y1_, z0_, z1_,
188 newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) {
191#ifdef PARALLEL_MODE_OMP
192 #pragma omp parallel for
194 for (iX=newX0; iX<=newX1; ++iX) {
195 for (
int iY=newY0; iY<=newY1; ++iY) {
196 for (
int iZ=newZ0; iZ<=newZ1; ++iZ) {
197 ComputeWallFunction(blockLattice, iX,iY,iZ);
204template<
typename T,
typename DESCRIPTOR>
207 for (
int iVel=0; iVel<DESCRIPTOR::q; ++iVel) {
208 if (descriptors::c<DESCRIPTOR>(iVel,index)==value) {
209 indices.push_back(iVel);
214template<
typename T,
typename DESCRIPTOR>
215void WallFunctionBoundaryProcessor3D<T,DESCRIPTOR>::calculateWallDistances(IndicatorF3D<T>* indicator)
217 using S = BaseType<T>;
223 T scaling = _converter.getConversionFactorLength() * 0.1;
224 _blockGeometryStructure.getPhysR(physR,{iX, iY, iZ});
225 Vector<T,3> origin(physR[0],physR[1],physR[2]);
226 Vector<T,3> normal(0.,0.,0.);
228 Vector<T,3> direction(0.,0.,0.);
229 int smallestDistance_i = 0;
230 T smallestDistance = 0.;
231 bool firstHit =
true;
232 origin[0] = physR[0];
233 origin[1] = physR[1];
234 origin[2] = physR[2];
235 int discreteDirection[6][3];
236 for (
int i=0; i < 6; i++) {
237 for (
int j=0; j < 3; j++) {
238 discreteDirection[i][j] = 0;
241 discreteDirection[0][0] = 1;
242 discreteDirection[1][0] = -1;
243 discreteDirection[2][1] = 1;
244 discreteDirection[3][1] = -1;
245 discreteDirection[4][2] = 1;
246 discreteDirection[5][2] = -1;
247 for (
int i=0; i < 6; i++) {
248 direction[0] = discreteDirection[i][0] * scaling;
249 direction[1] = discreteDirection[i][1] * scaling;
250 direction[2] = discreteDirection[i][2] * scaling;
251 if (indicator->distance(distance, origin, direction)) {
254 smallestDistance_i = i;
258 if (distance < smallestDistance) {
260 smallestDistance_i = i;
266 direction[0] = discreteDirection[smallestDistance_i][0] * scaling;
267 direction[1] = discreteDirection[smallestDistance_i][1] * scaling;
268 direction[2] = discreteDirection[smallestDistance_i][2] * scaling;
270 Vector<T,3> direction2(direction);
271 Vector<T,3> direction3(direction);
272 if (smallestDistance_i == 0 || smallestDistance_i == 1 ) {
273 direction2[1] = direction2[0] * scaling;
274 direction3[2] = direction3[0] * scaling;
276 else if (smallestDistance_i == 2 || smallestDistance_i == 3 ) {
277 direction2[0] = direction2[1] * scaling;
278 direction3[2] = direction3[1] * scaling;
281 direction2[0] = direction2[2] * scaling;
282 direction3[1] = direction3[2] * scaling;
285 Vector<S,3> directionN = direction*(1/
norm(direction));
286 Vector<S,3> POS(origin + smallestDistance*directionN);
288 indicator->distance(distance, origin, direction2);
289 Vector<S,3> direction2N = direction2*(1/
norm(direction2));
290 Vector<S,3> POS2(origin + distance*direction2N);
292 indicator->distance(distance, origin, direction3);
294 Vector<S,3> direction3N = direction3*(1/
norm(direction3));
295 Vector<S,3> POS3(origin + distance*direction3N);
297 Vector<S,3> vec1 (POS - POS2);
298 Vector<S,3> vec2 (POS - POS3);
300 normal[0] = -(vec1[1]*vec2[2] - vec1[2]*vec2[1]);
301 normal[1] = -(vec1[2]*vec2[0] - vec1[0]*vec2[2]);
302 normal[2] = -(vec1[0]*vec2[1] - vec1[1]*vec2[0]);
304 T normalMagnitude =
util::sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
305 normal[0] /= normalMagnitude;
306 normal[1] /= normalMagnitude;
307 normal[2] /= normalMagnitude;
309 unit_normal[0] = normal[0];
310 unit_normal[1] = normal[0];
311 unit_normal[2] = normal[0];
313 direction[0] = normal[0] * scaling;
314 direction[1] = normal[1] * scaling;
315 direction[2] = normal[2] * scaling;
317 indicator->distance(distance, origin, direction);
318 y_1 =
distance / _converter.getConversionFactorLength();
321template<
typename T,
typename DESCRIPTOR>
322void WallFunctionBoundaryProcessor3D<T,DESCRIPTOR>::VelGradFromSecondOrderFD(
bool NormalGradient, T Vel_BC[DESCRIPTOR::d], T Vel_1[DESCRIPTOR::d], T Vel_2[DESCRIPTOR::d], T VelGrad[DESCRIPTOR::d])
324 if (NormalGradient) {
325 if (orientation == 1) {
326 for (
int Dim=0; Dim<DESCRIPTOR::d; Dim++) {
327 VelGrad[Dim] =
fd::BSGradient(Vel_BC[Dim], Vel_1[Dim], Vel_2[Dim]);
331 for (
int Dim=0; Dim<DESCRIPTOR::d; Dim++) {
332 VelGrad[Dim] =
fd::FSGradient(Vel_BC[Dim], Vel_1[Dim], Vel_2[Dim]);
337 if (orientation == 1) {
338 for (
int Dim=0; Dim<DESCRIPTOR::d; Dim++) {
343 for (
int Dim=0; Dim<DESCRIPTOR::d; Dim++) {
350template<
typename T,
typename DESCRIPTOR>
351void WallFunctionBoundaryProcessor3D<T,DESCRIPTOR>::computeVelocityGradientTensor(T u_bc[DESCRIPTOR::d], T u_x1[DESCRIPTOR::d], T u_x2[DESCRIPTOR::d], T u_y1[DESCRIPTOR::d], T u_y2[DESCRIPTOR::d], T u_z1[DESCRIPTOR::d], T u_z2[DESCRIPTOR::d], T VelGrad[DESCRIPTOR::d][DESCRIPTOR::d])
353 T dx_u[DESCRIPTOR::d], dy_u[DESCRIPTOR::d], dz_u[DESCRIPTOR::d];
354 VelGradFromSecondOrderFD(direction == 0, u_bc, u_x1, u_x2, dx_u);
355 VelGradFromSecondOrderFD(direction == 1, u_bc, u_y1, u_y2, dy_u);
356 VelGradFromSecondOrderFD(direction == 2, u_bc, u_z1, u_z2, dz_u);
361 for (
int Dim = 0; Dim < DESCRIPTOR::d; Dim++) {
362 VelGrad[Dim][0] = dx_u[Dim];
363 VelGrad[Dim][1] = dy_u[Dim];
364 VelGrad[Dim][2] = dz_u[Dim];
368template<
typename T,
typename DESCRIPTOR>
369void WallFunctionBoundaryProcessor3D<T,DESCRIPTOR>::computeNeighborsU(BlockLattice<T,DESCRIPTOR>& blockLattice,
int x,
int y,
int z,
370 T u_x1[DESCRIPTOR::d], T u_x2[DESCRIPTOR::d], T u_y1[DESCRIPTOR::d], T u_y2[DESCRIPTOR::d], T u_z1[DESCRIPTOR::d], T u_z2[DESCRIPTOR::d])
374 if (direction == 0) {
375 blockLattice.get(x - discreteNormalX, y, z).computeU(u_x1);
376 blockLattice.get(x - 2*discreteNormalX, y, z).computeU(u_x2);
379 ComputeUWall(blockLattice, x + discreteNormalY + discreteNormalZ, y, z, u_x1);
380 ComputeUWall(blockLattice, x - discreteNormalY - discreteNormalZ, y, z, u_x2);
382 if (direction == 1) {
383 blockLattice.get(x, y - discreteNormalY, z).computeU(u_y1);
384 blockLattice.get(x, y - 2*discreteNormalY, z).computeU(u_y2);
387 ComputeUWall(blockLattice, x, y + discreteNormalX + discreteNormalZ, z, u_y1);
388 ComputeUWall(blockLattice, x, y - discreteNormalX - discreteNormalZ, z, u_y2);
390 if (direction == 2) {
391 blockLattice.get(x, y, z - discreteNormalZ).computeU(u_z1);
392 blockLattice.get(x, y, z - 2*discreteNormalZ).computeU(u_z2);
395 ComputeUWall(blockLattice, x, y, z + discreteNormalX + discreteNormalY, u_z1);
396 ComputeUWall(blockLattice, x, y, z - discreteNormalX - discreteNormalY, u_z2);
400template<
typename T,
typename DESCRIPTOR>
404 T u_x1[3], u_x2[3], u_y1[3], u_y2[3], u_z1[3], u_z2[3];
405 computeNeighborsU(blockLattice, x, y, z, u_x1, u_x2, u_y1, u_y2, u_z1, u_z2);
407 computeVelocityGradientTensor(u_bc, u_x1, u_x2, u_y1, u_y2, u_z1, u_z2, VelGrad);
410template<
typename T,
typename DESCRIPTOR>
412 T u_x1[DESCRIPTOR::d], T u_x2[DESCRIPTOR::d], T u_y1[DESCRIPTOR::d], T u_y2[DESCRIPTOR::d], T u_z1[DESCRIPTOR::d], T u_z2[DESCRIPTOR::d],
413 T& rho_x1, T& rho_x2, T& rho_y1, T& rho_y2, T& rho_z1, T& rho_z2)
417 if (direction == 0) {
418 rho_x1 = blockLattice.
get(x - discreteNormalX, y, z).computeRho();
419 rho_x2 = blockLattice.
get(x - 2*discreteNormalX, y, z).computeRho();
422 ComputeRhoWall(blockLattice, blockLattice.
get(x + discreteNormalY + discreteNormalZ, y, z), x + discreteNormalY + discreteNormalZ, y, z, u_x1, rho_x1);
423 ComputeRhoWall(blockLattice, blockLattice.
get(x - discreteNormalY - discreteNormalZ, y, z), x - discreteNormalY - discreteNormalZ, y, z, u_x2, rho_x2);
425 if (direction == 1) {
426 rho_y1 = blockLattice.
get(x, y - discreteNormalY, z).computeRho();
427 rho_y2 = blockLattice.
get(x, y - 2*discreteNormalY, z).computeRho();
430 ComputeRhoWall(blockLattice, blockLattice.
get(x, y - discreteNormalX - discreteNormalZ, z), x, y - discreteNormalX - discreteNormalZ, z, u_y1, rho_y1);
431 ComputeRhoWall(blockLattice, blockLattice.
get(x, y - discreteNormalX - discreteNormalZ, z), x, y - discreteNormalX - discreteNormalZ, z, u_y2, rho_y2);
433 if (direction == 2) {
434 rho_z1 = blockLattice.
get(x, y, z - discreteNormalZ).computeRho();
435 rho_z2 = blockLattice.
get(x, y, z - 2*discreteNormalZ).computeRho();
438 ComputeRhoWall(blockLattice, blockLattice.
get(x, y, z + discreteNormalX + discreteNormalY), x, y, z + discreteNormalX + discreteNormalY, u_z1, rho_z1);
439 ComputeRhoWall(blockLattice, blockLattice.
get(x, y, z - discreteNormalX - discreteNormalY), x, y, z - discreteNormalX - discreteNormalY, u_z2, rho_z2);
443template<
typename T,
typename DESCRIPTOR>
445 T u_x1[DESCRIPTOR::d], T u_x2[DESCRIPTOR::d], T u_y1[DESCRIPTOR::d], T u_y2[DESCRIPTOR::d], T u_z1[DESCRIPTOR::d], T u_z2[DESCRIPTOR::d],
446 T& rho_x1, T& rho_x2, T& rho_y1, T& rho_y2, T& rho_z1, T& rho_z2)
449 computeNeighborsU(blockLattice, x, y, z, u_x1, u_x2, u_y1, u_y2, u_z1, u_z2);
451 computeNeighborsRho(blockLattice, x, y, z, u_x1, u_x2, u_y1, u_y2, u_z1, u_z2, rho_x1, rho_x2, rho_y1, rho_y2, rho_z1, rho_z2);
456template<
typename T,
typename DESCRIPTOR>
464 T u_2[DESCRIPTOR::d];
466 blockLattice.
get(x - discreteNormalX, y - discreteNormalY, z - discreteNormalZ).computeRhoU(rho_2,u_2);
469 T u_2_dot_unit_normal = u_2[0] * unit_normal[0] +
470 u_2[1] * unit_normal[1] +
471 u_2[2] * unit_normal[2];
473 u_2_parallel[0] = u_2[0] - (u_2_dot_unit_normal * unit_normal[0]);
474 u_2_parallel[1] = u_2[1] - (u_2_dot_unit_normal * unit_normal[1]);
475 u_2_parallel[2] = u_2[2] - (u_2_dot_unit_normal * unit_normal[2]);
477 T u_2_parallel_norm =
util::sqrt(u_2_parallel[0] * u_2_parallel[0] +
478 u_2_parallel[1] * u_2_parallel[1] +
479 u_2_parallel[2] * u_2_parallel[2]);
481 T e_x_loc[3] = {0., 0., 0.};
482 if (u_2_parallel_norm != 0) {
483 T inv_u_2_parallel_norm = (1. /u_2_parallel_norm);
484 e_x_loc[0] = u_2_parallel[0] * inv_u_2_parallel_norm;
485 e_x_loc[1] = u_2_parallel[1] * inv_u_2_parallel_norm;
486 e_x_loc[2] = u_2_parallel[2] * inv_u_2_parallel_norm;
488 T u2 = e_x_loc[0] * u_2[0] +
489 e_x_loc[1] * u_2[1] +
498 T rho_phy = _converter.getPhysDensity();
499 T rho_lat = _converter.getLatticeDensity(rho_phy);
500 T Mol_Lat_Nu = _converter.getLatticeViscosity();
502 if (_wallFunctionParam.wallProfile == 0) {
504 Musker<T,T> musker_u2_Lat_MolVis(Mol_Lat_Nu, y_2, rho_lat);
506 T tau_w_guess = cell.template getField<descriptors::TAU_W>();
507 tau_w[0] = newton.solve(tau_w_guess,
false);
509 if (std::isnan(tau_w[0])) {
513 tau_w[0] = newton.solve(0.,
false);
516 if (std::isnan(tau_w[0])) {
518 clout <<
"Position = [" << x <<
", " << y <<
"," << z <<
"]" << std::endl;
519 clout <<
"Tau guess = " << tau_w_guess << std::endl;
524 Musker<T,T> musker_u1_Lat_MolVis(Mol_Lat_Nu, y_1, rho_lat);
525 musker_u1_Lat_MolVis(u1,tau_w);
534 tau_w[0] = output[0];
542 cell.template setField<descriptors::TAU_W>(&(tau_w[0]));
544 u[0] = e_x_loc[0] * u1[0];
545 u[1] = e_x_loc[1] * u1[0];
546 u[2] = e_x_loc[2] * u1[0];
548 if (_wallFunctionParam.bodyForce) {
549 auto force = cell.template getFieldPointer<descriptors::FORCE>();
550 for (
int iDim=0; iDim<DESCRIPTOR::d; ++iDim) {
551 u[iDim] -= force[iDim] / 2.;
557template<
typename T,
typename DESCRIPTOR>
561 T rho_phy = _converter.getPhysDensity();
562 T rho = _converter.getLatticeDensity(rho_phy);
563 T nu_mol = _converter.getLatticeViscosity();
568 T vanDriest = 1 -
util::exp(-y_plus / 26.);
570 T nu_turb =
util::pow(_wallFunctionParam.vonKarman * y_bc * vanDriest, 2.) * uxdz_abs;
572 T nu_eff = nu_turb + nu_mol;
574 tau_eff = (3. * nu_eff) + 0.5;
578template<
typename T,
typename DESCRIPTOR>
584 blockLattice.
get(x - discreteNormalX, y - discreteNormalY, z - discreteNormalZ).computeU(u_z);
585 blockLattice.
get(x - 2*discreteNormalX, y - 2*discreteNormalY, z - 2*discreteNormalZ).computeU(u_z2);
588 if (_wallFunctionParam.latticeWalldistance > 0.) {
589 y_bc = _wallFunctionParam.latticeWalldistance;
591 T normal_norm =
util::sqrt(discreteNormalX * discreteNormalX +
592 discreteNormalY * discreteNormalY +
593 discreteNormalZ * discreteNormalZ);
595 T tau_w = cell.template getField<descriptors::TAU_W>();
599 if (_wallFunctionParam.curved ==
true || orientation == 0) {
600 T inv_normal_norm = 1. / normal_norm;
602 normal[0] = discreteNormalX * (inv_normal_norm);
603 normal[1] = discreteNormalY * (inv_normal_norm);
604 normal[2] = discreteNormalZ * (inv_normal_norm);
607 T sP_normal_u2 = normal[0] * u_z[0] +
611 u_2_parallel[0] = u_z[0] - sP_normal_u2 * normal[0];
612 u_2_parallel[1] = u_z[1] - sP_normal_u2 * normal[1];
613 u_2_parallel[2] = u_z[2] - sP_normal_u2 * normal[2];
615 T u_2_parallel_norm =
util::sqrt(u_2_parallel[0] * u_2_parallel[0] +
616 u_2_parallel[1] * u_2_parallel[1] +
617 u_2_parallel[2] * u_2_parallel[2]);
620 T inv_u_2_parallel_norm = (1. /u_2_parallel_norm );
622 if (u_2_parallel_norm != 0) {
623 e_x_loc[0] = u_2_parallel[0] * inv_u_2_parallel_norm;
624 e_x_loc[1] = u_2_parallel[1] * inv_u_2_parallel_norm;
625 e_x_loc[2] = u_2_parallel[2] * inv_u_2_parallel_norm;
628 e_x_loc[0] = T(-discreteNormalX);
629 e_x_loc[1] = T(-discreteNormalY);
630 e_x_loc[2] = T(-discreteNormalZ);
636 u1 = e_x_loc[0] * u_bc[0] +
637 e_x_loc[1] * u_bc[1] +
638 e_x_loc[2] * u_bc[2];
639 u2 = e_x_loc[0] * u_z[0] +
640 e_x_loc[1] * u_z[1] +
642 u3 = e_x_loc[0] * u_z2[0] +
643 e_x_loc[1] * u_z2[1] +
644 e_x_loc[2] * u_z2[2];
646 computeVanDriestTauEff(y_bc, tau_w, u1, u2, u3, tau_eff);
649 computeVanDriestTauEff(y_bc, tau_w, u_bc[direction], u_z[direction], u_z2[direction], tau_eff);
651 cell.template setField<descriptors::TAU_EFF>(&(tau_eff));
655template<
typename T,
typename DESCRIPTOR>
660 if (_wallFunctionParam.rhoMethod == 0) {
663 T u_bc_Perpendicular = 0.;
664 if (_wallFunctionParam.bodyForce) {
665 T u_bc_tmp[DESCRIPTOR::d];
666 auto force = cell.template getFieldPointer<descriptors::FORCE>();
667 for (
int iDim = 0; iDim<DESCRIPTOR::d; ++iDim) {
668 u_bc_tmp[iDim] = u_bc[iDim] - force[iDim] / 2.;
670 for (
int iDim = 0; iDim<DESCRIPTOR::d; ++iDim) {
671 u_bc_Perpendicular += u_bc_tmp[iDim]*unit_normal[iDim];
675 for (
int iDim = 0; iDim<DESCRIPTOR::d; ++iDim) {
676 u_bc_Perpendicular += u_bc[iDim]*unit_normal[iDim];
681 for (
unsigned fIndex=0; fIndex<onWallIndices.size(); ++fIndex) {
682 rhoOnWall += cell[onWallIndices[fIndex]];
686 for (
unsigned fIndex=0; fIndex<normalIndices.size(); ++fIndex) {
687 rhoNormal += cell[normalIndices[fIndex]];
690 rho_bc = ((T)2*rhoNormal+rhoOnWall+(T)1) /
691 ((T)1+u_bc_Perpendicular);
694 else if (_wallFunctionParam.rhoMethod == 1) {
696 rho_bc = blockLattice.
get(x - discreteNormalX, y - discreteNormalY, z - discreteNormalZ).computeRho();
699 else if (_wallFunctionParam.rhoMethod == 2) {
705template<
typename T,
typename DESCRIPTOR>
708 T sum_cell_fneq = 0.;
709 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
710 sum_cell_fneq += fneq_bc[iPop];
714 for (
int iAlpha=0; iAlpha < DESCRIPTOR::d; ++iAlpha) {
715 for (
int iBeta=iAlpha; iBeta < DESCRIPTOR::d; ++iBeta) {
717 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
718 pi_bc[iPi] += descriptors::c<DESCRIPTOR>(iPop)[iAlpha]*descriptors::c<DESCRIPTOR>(iPop)[iBeta]*fneq_bc[iPop];
721 pi_bc[iPi] -= (1./descriptors::invCs2<T,DESCRIPTOR>())*sum_cell_fneq;
727 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
732template<
typename T,
typename DESCRIPTOR>
736 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
737 fneq_bc[iPop] = cell[iPop] - dynamics -> computeEquilibrium(iPop,rho_bc,u_bc);
739 for (
unsigned fIndex=0; fIndex<normalInwardsIndices.size(); ++fIndex) {
740 fneq_bc[normalInwardsIndices[fIndex]] = fneq_bc[descriptors::opposite<DESCRIPTOR>(normalInwardsIndices[fIndex])];
743 computeRFneqfromFneq(fneq_bc);
746template<
typename T,
typename DESCRIPTOR>
747void WallFunctionBoundaryProcessor3D<T,DESCRIPTOR>::computeFneqENeq(
BlockLattice<T,DESCRIPTOR>& blockLattice,
Cell<T,DESCRIPTOR>& cell,
int x,
int y,
int z, T u_bc[DESCRIPTOR::d], T rho_bc, T fneq_bc[DESCRIPTOR::q])
749 Cell<T,DESCRIPTOR> cell_fluid = blockLattice.
get(x - discreteNormalX, y - discreteNormalY, z - discreteNormalZ);
751 T rho_fluid, u_fluid[DESCRIPTOR::d];
753 for (
int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
754 fneq_bc[iPop] = cell_fluid[iPop] - dynamics_fluid -> computeEquilibrium(iPop,rho_fluid,u_fluid);
758template<
typename T,
typename DESCRIPTOR>
759void WallFunctionBoundaryProcessor3D<T,DESCRIPTOR>::computeFneqRSOFD(
BlockLattice<T,DESCRIPTOR>& blockLattice,
Cell<T,DESCRIPTOR>& cell,
int x,
int y,
int z, T u_bc[DESCRIPTOR::d], T rho_bc, T fneq_bc[DESCRIPTOR::q])
762 T Velocity_Grad[DESCRIPTOR::d][DESCRIPTOR::d];
763 computeVelocityGradient(blockLattice, x, y, z, u_bc, Velocity_Grad);
765 T tau_eff = cell.template getField<descriptors::TAU_EFF>();
766 T Factor = -1.*tau_eff*rho_bc/descriptors::invCs2<T,DESCRIPTOR>();
768 for (
int Alpha=0; Alpha<DESCRIPTOR::d; ++Alpha) {
769 for (
int Beta=Alpha; Beta<DESCRIPTOR::d; ++Beta) {
770 pi_bc[iPi] = Factor*(Velocity_Grad[Alpha][Beta] + Velocity_Grad[Beta][Alpha]);
775 if (_wallFunctionParam.bodyForce) {
778 auto F = cell.template getField<descriptors::FORCE>();
781 for (
int Alpha=0; Alpha<DESCRIPTOR::d; ++Alpha) {
782 for (
int Beta=Alpha; Beta<DESCRIPTOR::d; ++Beta) {
783 pi_bc[iPi] += -1.*(rho_bc/2.)*(F[Alpha]*u_bc[Beta] + u_bc[Alpha]*F[Beta]);
789 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
797template<
typename T,
typename DESCRIPTOR>
798void WallFunctionBoundaryProcessor3D<T,DESCRIPTOR>::ComputeFneqWall(
BlockLattice<T,DESCRIPTOR>& blockLattice,
Cell<T,DESCRIPTOR>& cell,
int x,
int y,
int z, T u_bc[DESCRIPTOR::d], T rho_bc, T fneq_bc[DESCRIPTOR::q])
802 if (_wallFunctionParam.fneqMethod == 0) {
803 computeFneqRNEBB(cell, u_bc, rho_bc, fneq_bc);
807 else if (_wallFunctionParam.fneqMethod == 1) {
808 computeFneqENeq(blockLattice, cell, x, y, z, u_bc, rho_bc, fneq_bc);
812 else if (_wallFunctionParam.fneqMethod == 2) {
813 computeFneqRSOFD(blockLattice, cell, x, y, z, u_bc, rho_bc, fneq_bc);
817 else if (_wallFunctionParam.fneqMethod == 3) {
818 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
825template<
typename T,
typename DESCRIPTOR>
830 T u_bc[DESCRIPTOR::d];
831 T fneq_bc[DESCRIPTOR::q];
834 ComputeUWall(blockLattice, x, y, z, u_bc);
836 if (_wallFunctionParam.useVanDriest) {
838 ComputeTauEff(blockLattice, cell_bc, x, y, z, u_bc);
842 ComputeRhoWall(blockLattice, cell_bc, x, y, z, u_bc, rho_bc);
845 ComputeFneqWall(blockLattice, cell_bc, x, y, z, u_bc, rho_bc, fneq_bc);
850 for (
int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
851 cell_bc[iPop] = dynamics_bc -> computeEquilibrium(iPop,rho_bc,u_bc) + fneq_bc[iPop];
852 if (std::isnan(cell_bc[iPop])) {
854 clout <<
"Musker Profile Computation" << std::endl;
855 clout <<
"Position = [" << x <<
", " << y <<
"," << z <<
"]" << std::endl;
856 clout <<
"Normal outwards = [" << discreteNormalX <<
", " << discreteNormalY <<
"," << discreteNormalZ <<
"]" << std::endl;
857 clout <<
"Velocity at boundary u_bc = [" << u_bc[0] <<
", " << u_bc[1] <<
"," << u_bc[2] <<
"]" << std::endl;
858 clout <<
"Density at boundary rho_bc = " << rho_bc << std::endl;
865template<
typename T,
typename DESCRIPTOR>
868 processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
873template<
typename T,
typename DESCRIPTOR>
875 std::vector<int> discreteNormal, std::vector<int> missingIndices,
879 _blockGeometryStructure(blockGeometryStructure),
880 _discreteNormal(discreteNormal), _missingIndices(missingIndices),
881 _converter(converter), _wallFunctionParam(
wallFunctionParam), _geoIndicator(geoIndicator)
884template<
typename T,
typename DESCRIPTOR>
889 _blockGeometryStructure, _discreteNormal, _missingIndices,
890 _converter, _wallFunctionParam, _geoIndicator);
893template<
typename T,
typename DESCRIPTOR>
898 _blockGeometryStructure, _discreteNormal, _missingIndices,
899 _converter, _wallFunctionParam, _geoIndicator);
AnalyticalF are applications from DD to XD, where X is set by the constructor.
Representation of a block geometry.
Platform-abstracted block lattice for external access and inter-block interaction.
Cell< T, DESCRIPTOR > get(CellID iCell)
Get Cell interface for index iCell.
Highest-level interface to Cell data.
Dynamics< T, DESCRIPTOR > * getDynamics()
Get a pointer to the dynamics.
void computeRhoU(T &rho, T u[descriptors::d< DESCRIPTOR >()]) const
Compute fluid velocity and particle density on the cell.
std::string & getName()
read and write access to name
IndicatorF3D is an application from .
bool operator()(T output[], const S tau_w[])
has to be implemented for 'every' derived class
class for marking output with some text
std::string & getName()
read and write access to name
PowerLawProfile(T nu, T u2, T y2, T y1, T rho)
bool operator()(T output[], const S tau_w[])
has to be implemented for 'every' derived class
Conversion between physical and lattice units, as well as discretization.
WallFunctionBoundaryProcessor3D(int x0, int x1, int y0, int y1, int z0, int z1, BlockGeometry< T, 3 > &blockGeometryStructure, std::vector< int > discreteNormal, std::vector< int > missingIndices, UnitConverter< T, DESCRIPTOR > const &converter, wallFunctionParam< T > const &wallFunctionParam, IndicatorF3D< T > *geoIndicator)
virtual void processSubDomain(BlockLattice< T, DESCRIPTOR > &blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
Execute post-processing step on a sublattice.
virtual void ComputeWallFunction(BlockLattice< T, DESCRIPTOR > &blockLattice, int x, int y, int z)
virtual void process(BlockLattice< T, DESCRIPTOR > &blockLattice)
Execute post-processing step.
PostProcessorGenerator3D< T, DESCRIPTOR > * clone() const override
WallFunctionBoundaryProcessorGenerator3D(int x0, int x1, int y0, int y1, int z0, int z1, BlockGeometry< T, 3 > &blockGeometryStructure, std::vector< int > discreteNormal, std::vector< int > missingIndices, UnitConverter< T, DESCRIPTOR > const &converter, wallFunctionParam< T > const &wallFunctionParam, IndicatorF3D< T > *geoIndicator)
PostProcessor3D< T, DESCRIPTOR > * generate() const override
constexpr T BSGradient(T u_0, T u_1, T u_2)
Backward second-order first derivative.
constexpr T centralGradient(T u_p1, T u_m1) any_platform
Second-order central gradient (u_p1 = u(x+1))
constexpr T FSGradient(T u_0, T u_1, T u_2)
Forward second-order first derivative.
constexpr T boundaryGradient(T u_0, T u_1, T u_2) any_platform
Second-order asymmetric gradient (u_1 = u(x+1))
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
ADf< T, DIM > abs(const ADf< T, DIM > &a)
bool distance(S &distance, const Vector< S, D > &origin, const Vector< S, D > &direction, S precision, S pitch, F1 isInside, F2 isInsideBoundingBox)
ADf< T, DIM > log(const ADf< T, DIM > &a)
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
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)
ADf< T, DIM > atan(const ADf< T, DIM > &a)
ADf< T, DIM > exp(const ADf< T, DIM > &a)
Top level namespace for all of OpenLB.
constexpr T norm(const ScalarVector< T, D, IMPL > &a)
Euclidean vector norm.
Interface for per-cell dynamics.
Compute number of elements of a symmetric d-dimensional tensor.
Set of functions commonly used in LB computations – header file.