From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Mon, 24 Jun 2019 14:43:36 +0200 Subject: Initialize at openlb-1-3 --- src/boundary/boundaryPostProcessors3D.hh | 994 +++++++++++++++++++++++++++++++ 1 file changed, 994 insertions(+) create mode 100644 src/boundary/boundaryPostProcessors3D.hh (limited to 'src/boundary/boundaryPostProcessors3D.hh') diff --git a/src/boundary/boundaryPostProcessors3D.hh b/src/boundary/boundaryPostProcessors3D.hh new file mode 100644 index 0000000..2eae9cd --- /dev/null +++ b/src/boundary/boundaryPostProcessors3D.hh @@ -0,0 +1,994 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2006, 2007 Jonas Latt + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program; if not, write to the Free + * Software Foundation, Inc., 51 Franklin Street, Fifth Floor, + * Boston, MA 02110-1301, USA. +*/ + +#ifndef BOUNDARY_POST_PROCESSORS_3D_HH +#define BOUNDARY_POST_PROCESSORS_3D_HH + +#include "boundaryPostProcessors3D.h" +#include "core/finiteDifference3D.h" +#include "core/blockLattice3D.h" +#include "dynamics/firstOrderLbHelpers.h" +#include "core/util.h" +#include "utilities/vectorHelpers.h" + +namespace olb { + +//////// PlaneFdBoundaryProcessor3D /////////////////////////////////// + +template +PlaneFdBoundaryProcessor3D:: +PlaneFdBoundaryProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_) +{ + OLB_PRECONDITION(x0==x1 || y0==y1 || z0==z1); +} + +template +void PlaneFdBoundaryProcessor3D:: +processSubDomain(BlockLattice3D& blockLattice, + int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) +{ + using namespace olb::util::tensorIndices3D; + + int newX0, newX1, newY0, newY1, newZ0, newZ1; + if ( util::intersect ( + x0, x1, y0, y1, z0, z1, + x0_, x1_, y0_, y1_, z0_, z1_, + newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) { + int iX; + +#ifdef PARALLEL_MODE_OMP + #pragma omp parallel for +#endif + for (iX=newX0; iX<=newX1; ++iX) { + T dx_u[DESCRIPTOR::d], dy_u[DESCRIPTOR::d], dz_u[DESCRIPTOR::d]; + for (int iY=newY0; iY<=newY1; ++iY) { + for (int iZ=newZ0; iZ<=newZ1; ++iZ) { + Cell& cell = blockLattice.get(iX,iY,iZ); + Dynamics* dynamics = blockLattice.getDynamics(iX, iY, iZ); + T rho, u[DESCRIPTOR::d]; + cell.computeRhoU(rho,u); + + interpolateGradients<0> ( blockLattice, dx_u, iX, iY, iZ ); + interpolateGradients<1> ( blockLattice, dy_u, iX, iY, iZ ); + interpolateGradients<2> ( blockLattice, dz_u, iX, iY, iZ ); + T dx_ux = dx_u[0]; + T dy_ux = dy_u[0]; + T dz_ux = dz_u[0]; + T dx_uy = dx_u[1]; + T dy_uy = dy_u[1]; + T dz_uy = dz_u[1]; + T dx_uz = dx_u[2]; + T dy_uz = dy_u[2]; + T dz_uz = dz_u[2]; + T omega = dynamics->getOmega(); + T sToPi = - rho / descriptors::invCs2() / omega; + T pi[util::TensorVal::n]; + pi[xx] = (T)2 * dx_ux * sToPi; + pi[yy] = (T)2 * dy_uy * sToPi; + pi[zz] = (T)2 * dz_uz * sToPi; + pi[xy] = (dx_uy + dy_ux) * sToPi; + pi[xz] = (dx_uz + dz_ux) * sToPi; + pi[yz] = (dy_uz + dz_uy) * sToPi; + + // Computation of the particle distribution functions + // according to the regularized formula + T uSqr = util::normSqr(u); + for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) + cell[iPop] = dynamics -> computeEquilibrium(iPop,rho,u,uSqr) + + firstOrderLbHelpers::fromPiToFneq(iPop, pi); + } + } + } + } +} + +template +void PlaneFdBoundaryProcessor3D:: +process(BlockLattice3D& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1); +} + + +template +template +void PlaneFdBoundaryProcessor3D:: +interpolateGradients(BlockLattice3D const& blockLattice, T velDeriv[DESCRIPTOR::d], + int iX, int iY, int iZ) const +{ + fd::DirectedGradients3D:: + interpolateVector(velDeriv, blockLattice, iX, iY, iZ); +} + + +//////// PlaneFdBoundaryProcessorGenerator3D /////////////////////////////// + +template +PlaneFdBoundaryProcessorGenerator3D:: +PlaneFdBoundaryProcessorGenerator3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) + : PostProcessorGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_) +{ } + +template +PostProcessor3D* PlaneFdBoundaryProcessorGenerator3D:: +generate() const +{ + return new PlaneFdBoundaryProcessor3D + ( this->x0, this->x1, this->y0, this->y1, this->z0, this->z1 ); +} + +template +PostProcessorGenerator3D* +PlaneFdBoundaryProcessorGenerator3D::clone() const +{ + return new PlaneFdBoundaryProcessorGenerator3D + (this->x0, this->x1, this->y0, this->y1, this->z0, this->z1); +} + + +//////// StraightConvectionBoundaryProcessorGenerator3D //////////////////////////////// + +template +StraightConvectionBoundaryProcessor3D:: +StraightConvectionBoundaryProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, T* uAv_) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), uAv(uAv_) +{ + OLB_PRECONDITION(x0==x1 || y0==y1 || z0==z1); + + saveCell = new T*** [(size_t)(x1_-x0_+1)]; + for (int iX=0; iX<=x1_-x0_; ++iX) { + saveCell[iX] = new T** [(size_t)(y1_-y0_+1)]; + for (int iY=0; iY<=y1_-y0_; ++iY) { + saveCell[iX][iY] = new T* [(size_t)(z1_-z0_+1)]; + for (int iZ=0; iZ<=z1_-z0_; ++iZ) { + saveCell[iX][iY][iZ] = new T [(size_t)(DESCRIPTOR::q)]; + for (int iPop=0; iPop +StraightConvectionBoundaryProcessor3D:: +~StraightConvectionBoundaryProcessor3D() +{ + for (int iX=0; iX<=x1-x0; ++iX) { + for (int iY=0; iY<=y1-y0; ++iY) { + for (int iZ=0; iZ<=z1-z0; ++iZ) { + delete [] saveCell[iX][iY][iZ]; + } + delete [] saveCell[iX][iY]; + } + delete [] saveCell[iX]; + } + delete [] saveCell; +} + +template +void StraightConvectionBoundaryProcessor3D:: +processSubDomain(BlockLattice3D& blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) +{ + int newX0, newX1, newY0, newY1, newZ0, newZ1; + if ( util::intersect ( + x0, x1, y0, y1, z0, z1, + x0_, x1_, y0_, y1_, z0_, z1_, + newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) { + + int iX; +#ifdef PARALLEL_MODE_OMP + #pragma omp parallel for +#endif + for (iX=newX0; iX<=newX1; ++iX) { + for (int iY=newY0; iY<=newY1; ++iY) { + for (int iZ=newZ0; iZ<=newZ1; ++iZ) { + Cell& cell = blockLattice.get(iX,iY,iZ); + for (int iPop = 0; iPop < DESCRIPTOR::q ; ++iPop) { + if (descriptors::c(iPop,direction)==-orientation) { + // using default -1 avoids wrong first call + if (!util::nearZero(1 + saveCell[iX-newX0][iY-newY0][iZ-newZ0][iPop]) ) { + cell[iPop] = saveCell[iX-newX0][iY-newY0][iZ-newZ0][iPop]; + } + } + } + + T rho0, u0[3]; + T rho1, u1[3]; + T rho2, u2[3]; + if (direction==0) { + blockLattice.get(iX,iY,iZ).computeRhoU(rho0,u0); + blockLattice.get(iX-orientation,iY,iZ).computeRhoU(rho1,u1); + blockLattice.get(iX-orientation*2,iY,iZ).computeRhoU(rho2,u2); + } + else if (direction==1) { + blockLattice.get(iX,iY,iZ).computeRhoU(rho0,u0); + blockLattice.get(iX,iY-orientation,iZ).computeRhoU(rho1,u1); + blockLattice.get(iX,iY-orientation*2,iZ).computeRhoU(rho2,u2); + } + else { + blockLattice.get(iX,iY,iZ).computeRhoU(rho0,u0); + blockLattice.get(iX,iY,iZ-orientation).computeRhoU(rho1,u1); + blockLattice.get(iX,iY,iZ-orientation*2).computeRhoU(rho2,u2); + } + + // rho0 = T(1); rho1 = T(1); rho2 = T(1); + + T uDelta[3]; + T uAverage = rho0*u0[direction]; + if (uAv!=nullptr) { + uAverage = *uAv * rho0; + } + uDelta[0]=-uAverage*0.5*(3*rho0*u0[0]-4*rho1*u1[0]+rho2*u2[0]); + uDelta[1]=-uAverage*0.5*(3*rho0*u0[1]-4*rho1*u1[1]+rho2*u2[1]); + uDelta[2]=-uAverage*0.5*(3*rho0*u0[2]-4*rho1*u1[2]+rho2*u2[2]); + + for (int iPop = 0; iPop < DESCRIPTOR::q ; ++iPop) { + if (descriptors::c(iPop,direction) == -orientation) { + saveCell[iX-newX0][iY-newY0][iZ-newZ0][iPop] = cell[iPop] + descriptors::invCs2()*descriptors::t(iPop)*(uDelta[0]*descriptors::c(iPop,0)+uDelta[1]*descriptors::c(iPop,1)+uDelta[2]*descriptors::c(iPop,2)); + } + } + } + } + } + } +} + +template +void StraightConvectionBoundaryProcessor3D:: +process(BlockLattice3D& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1); +} + + +//////// StraightConvectionBoundaryProcessorGenerator2D //////////////////////////////// + +template +StraightConvectionBoundaryProcessorGenerator3D:: +StraightConvectionBoundaryProcessorGenerator3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, T* uAv_) + : PostProcessorGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_), uAv(uAv_) +{ } + +template +PostProcessor3D* +StraightConvectionBoundaryProcessorGenerator3D::generate() const +{ + return new StraightConvectionBoundaryProcessor3D + ( this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, uAv); +} + +template +PostProcessorGenerator3D* +StraightConvectionBoundaryProcessorGenerator3D::clone() const +{ + return new StraightConvectionBoundaryProcessorGenerator3D + (this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, uAv); +} + + +//////// OuterVelocityEdgeProcessor3D /////////////////////////////////// + +template +OuterVelocityEdgeProcessor3D:: +OuterVelocityEdgeProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_) +{ + OLB_PRECONDITION ( + (plane==2 && x0==x1 && y0==y1) || + (plane==1 && x0==x1 && z0==z1) || + (plane==0 && y0==y1 && z0==z1) ); + +} + +template +void OuterVelocityEdgeProcessor3D:: +processSubDomain(BlockLattice3D& blockLattice, + int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) +{ + using namespace olb::util::tensorIndices3D; + + int newX0, newX1, newY0, newY1, newZ0, newZ1; + if ( util::intersect ( x0, x1, y0, y1, z0, z1, + x0_, x1_, y0_, y1_, z0_, z1_, + newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) { + int iX; + +#ifdef PARALLEL_MODE_OMP + #pragma omp parallel for +#endif + for (iX=newX0; iX<=newX1; ++iX) { + for (int iY=newY0; iY<=newY1; ++iY) { + for (int iZ=newZ0; iZ<=newZ1; ++iZ) { + Cell& cell = blockLattice.get(iX,iY,iZ); + Dynamics* dynamics = blockLattice.getDynamics(iX, iY, iZ); + + T rho10 = getNeighborRho(iX,iY,iZ,1,0, blockLattice); + T rho01 = getNeighborRho(iX,iY,iZ,0,1, blockLattice); + T rho20 = getNeighborRho(iX,iY,iZ,2,0, blockLattice); + T rho02 = getNeighborRho(iX,iY,iZ,0,2, blockLattice); + T rho = (T)2/(T)3*(rho01+rho10)-(T)1/(T)6*(rho02+rho20); + + T dA_uB_[3][3]; + interpolateGradients ( blockLattice, dA_uB_[0], iX, iY, iZ ); + interpolateGradients ( blockLattice, dA_uB_[1], iX, iY, iZ ); + interpolateGradients ( blockLattice, dA_uB_[2], iX, iY, iZ ); + T dA_uB[3][3]; + for (int iBeta=0; iBeta<3; ++iBeta) { + dA_uB[plane][iBeta] = dA_uB_[0][iBeta]; + dA_uB[direction1][iBeta] = dA_uB_[1][iBeta]; + dA_uB[direction2][iBeta] = dA_uB_[2][iBeta]; + } + T omega = dynamics -> getOmega(); + T sToPi = - rho / descriptors::invCs2() / omega; + T pi[util::TensorVal::n]; + pi[xx] = (T)2 * dA_uB[0][0] * sToPi; + pi[yy] = (T)2 * dA_uB[1][1] * sToPi; + pi[zz] = (T)2 * dA_uB[2][2] * sToPi; + pi[xy] = (dA_uB[0][1]+dA_uB[1][0]) * sToPi; + pi[xz] = (dA_uB[0][2]+dA_uB[2][0]) * sToPi; + pi[yz] = (dA_uB[1][2]+dA_uB[2][1]) * sToPi; + + // Computation of the particle distribution functions + // according to the regularized formula + T u[DESCRIPTOR::d]; + cell.computeU(u); + T uSqr = util::normSqr(u); + + for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] = dynamics -> computeEquilibrium(iPop,rho,u,uSqr) + + firstOrderLbHelpers::fromPiToFneq(iPop, pi); + } + } + } + } + } +} + +template +void OuterVelocityEdgeProcessor3D:: +process(BlockLattice3D& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1); +} + +template +T OuterVelocityEdgeProcessor3D:: +getNeighborRho(int x, int y, int z, int step1, int step2, BlockLattice3D const& blockLattice) +{ + int coords[3] = {x, y, z}; + coords[direction1] += -normal1*step1; + coords[direction2] += -normal2*step2; + return blockLattice.get(coords[0], coords[1], coords[2]).computeRho(); +} + +template +template +void OuterVelocityEdgeProcessor3D:: +interpolateGradients(BlockLattice3D const& blockLattice, + T velDeriv[DESCRIPTOR::d], + int iX, int iY, int iZ) const +{ + fd::DirectedGradients3D:: + interpolateVector(velDeriv, blockLattice, iX, iY, iZ); +} + +//////// OuterVelocityEdgeProcessorGenerator3D /////////////////////////////// + +template +OuterVelocityEdgeProcessorGenerator3D:: +OuterVelocityEdgeProcessorGenerator3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) + : PostProcessorGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_) +{ } + +template +PostProcessor3D* +OuterVelocityEdgeProcessorGenerator3D:: +generate() const +{ + return new OuterVelocityEdgeProcessor3D < T,DESCRIPTOR, plane,normal1,normal2 > + ( this->x0, this->x1, this->y0, this->y1, this->z0, this->z1); +} + +template +PostProcessorGenerator3D* +OuterVelocityEdgeProcessorGenerator3D::clone() const +{ + return new OuterVelocityEdgeProcessorGenerator3D + (this->x0, this->x1, this->y0, this->y1, this->z0, this->z1); +} + +/////////// OuterVelocityCornerProcessor3D ///////////////////////////////////// + +template +OuterVelocityCornerProcessor3D:: +OuterVelocityCornerProcessor3D ( int x_, int y_, int z_ ) + : x(x_), y(y_), z(z_) +{ } + +template +void OuterVelocityCornerProcessor3D:: +process(BlockLattice3D& blockLattice) +{ + using namespace olb::util::tensorIndices3D; + Cell& cell = blockLattice.get(x,y,z); + Dynamics* dynamics = blockLattice.getDynamics(x, y, z); + + T rho100 = blockLattice.get(x - 1*xNormal, y - 0*yNormal, z - 0*zNormal).computeRho(); + T rho010 = blockLattice.get(x - 0*xNormal, y - 1*yNormal, z - 0*zNormal).computeRho(); + T rho001 = blockLattice.get(x - 0*xNormal, y - 0*yNormal, z - 1*zNormal).computeRho(); + T rho200 = blockLattice.get(x - 2*xNormal, y - 0*yNormal, z - 0*zNormal).computeRho(); + T rho020 = blockLattice.get(x - 0*xNormal, y - 2*yNormal, z - 0*zNormal).computeRho(); + T rho002 = blockLattice.get(x - 0*xNormal, y - 0*yNormal, z - 2*zNormal).computeRho(); + T rho = (T)4/(T)9 * (rho001 + rho010 + rho100) - (T)1/(T)9 * (rho002 + rho020 + rho200); + + T dx_u[DESCRIPTOR::d], dy_u[DESCRIPTOR::d], dz_u[DESCRIPTOR::d]; + fd::DirectedGradients3D::interpolateVector(dx_u, blockLattice, x,y,z); + fd::DirectedGradients3D::interpolateVector(dy_u, blockLattice, x,y,z); + fd::DirectedGradients3D::interpolateVector(dz_u, blockLattice, x,y,z); + + T dx_ux = dx_u[0]; + T dy_ux = dy_u[0]; + T dz_ux = dz_u[0]; + T dx_uy = dx_u[1]; + T dy_uy = dy_u[1]; + T dz_uy = dz_u[1]; + T dx_uz = dx_u[2]; + T dy_uz = dy_u[2]; + T dz_uz = dz_u[2]; + T omega = dynamics -> getOmega(); + T sToPi = - rho / descriptors::invCs2() / omega; + T pi[util::TensorVal::n]; + pi[xx] = (T)2 * dx_ux * sToPi; + pi[yy] = (T)2 * dy_uy * sToPi; + pi[zz] = (T)2 * dz_uz * sToPi; + pi[xy] = (dx_uy + dy_ux) * sToPi; + pi[xz] = (dx_uz + dz_ux) * sToPi; + pi[yz] = (dy_uz + dz_uy) * sToPi; + + // Computation of the particle distribution functions + // according to the regularized formula + T u[DESCRIPTOR::d]; + cell.computeU(u); + T uSqr = util::normSqr(u); + + for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] = dynamics -> computeEquilibrium(iPop,rho,u,uSqr) + + firstOrderLbHelpers::fromPiToFneq(iPop, pi); + } + +} + +template +void OuterVelocityCornerProcessor3D:: +processSubDomain(BlockLattice3D& blockLattice, + int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) +{ + if (util::contained(x, y, z, x0_, x1_, y0_, y1_, z0_, z1_)) { + process(blockLattice); + } +} + +//////// OuterVelocityCornerProcessorGenerator3D /////////////////////////////// + +template +OuterVelocityCornerProcessorGenerator3D:: +OuterVelocityCornerProcessorGenerator3D(int x_, int y_, int z_) + : PostProcessorGenerator3D(x_,x_, y_,y_, z_,z_) +{ } + +template +PostProcessor3D* +OuterVelocityCornerProcessorGenerator3D:: +generate() const +{ + return new OuterVelocityCornerProcessor3D + ( this->x0, this->y0, this->z0 ); +} + +template +PostProcessorGenerator3D* +OuterVelocityCornerProcessorGenerator3D::clone() const +{ + return new OuterVelocityCornerProcessorGenerator3D + (this->x0, this->y0, this->z0); +} + + +//////// SlipBoundaryProcessor3D //////////////////////////////// + +template +SlipBoundaryProcessor3D:: +SlipBoundaryProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX, int discreteNormalY, int discreteNormalZ) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_) +{ + OLB_PRECONDITION(x0==x1 || y0==y1 || z0==z1); + int mirrorDirection0; + int mirrorDirection1; + int mirrorDirection2; + int mult = 2 / (discreteNormalX*discreteNormalX + discreteNormalY*discreteNormalY + discreteNormalZ*discreteNormalZ); + reflectionPop[0] = 0; + for (int iPop = 1; iPop < DESCRIPTOR::q; iPop++) { + reflectionPop[iPop] = 0; + // iPop are the directions which pointing into the fluid, discreteNormal is pointing outwarts + int scalarProduct = descriptors::c(iPop,0)*discreteNormalX + descriptors::c(iPop,1)*discreteNormalY + descriptors::c(iPop,2)*discreteNormalZ; + if ( scalarProduct < 0) { + // bounce back for the case discreteNormalX = discreteNormalY = discreteNormalZ = 1, that is mult=0 + if (mult == 0) { + mirrorDirection0 = -descriptors::c(iPop,0); + mirrorDirection1 = -descriptors::c(iPop,1); + mirrorDirection2 = -descriptors::c(iPop,2); + } + else { + mirrorDirection0 = descriptors::c(iPop,0) - mult*scalarProduct*discreteNormalX; + mirrorDirection1 = descriptors::c(iPop,1) - mult*scalarProduct*discreteNormalY; + mirrorDirection2 = descriptors::c(iPop,2) - mult*scalarProduct*discreteNormalZ; + } + + // run through all lattice directions and look for match of direction + for (int i = 1; i < DESCRIPTOR::q; i++) { + if (descriptors::c(i,0)==mirrorDirection0 + && descriptors::c(i,1)==mirrorDirection1 + && descriptors::c(i,2)==mirrorDirection2) { + reflectionPop[iPop] = i; + break; + } + } + } + } +} + +template +void SlipBoundaryProcessor3D:: +processSubDomain(BlockLattice3D& blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) +{ + int newX0, newX1, newY0, newY1, newZ0, newZ1; + if ( util::intersect ( + x0, x1, y0, y1, z0, z1, + x0_, x1_, y0_, y1_, z0_, z1_, + newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) { + + int iX; +#ifdef PARALLEL_MODE_OMP + #pragma omp parallel for +#endif + for (iX=newX0; iX<=newX1; ++iX) { + for (int iY=newY0; iY<=newY1; ++iY) { + for (int iZ=newZ0; iZ<=newZ1; ++iZ) { + for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) { + if (reflectionPop[iPop]!=0) { + //do reflection + blockLattice.get(iX,iY,iZ)[iPop] = blockLattice.get(iX,iY,iZ)[reflectionPop[iPop]]; + } + } + } + } + } + } +} + +template +void SlipBoundaryProcessor3D:: +process(BlockLattice3D& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1); +} + +//////// SlipBoundaryProcessorGenerator3D //////////////////////////////// + +template +SlipBoundaryProcessorGenerator3D:: +SlipBoundaryProcessorGenerator3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX_, int discreteNormalY_, int discreteNormalZ_) + : PostProcessorGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_), discreteNormalX(discreteNormalX_), discreteNormalY(discreteNormalY_), discreteNormalZ(discreteNormalZ_) +{ } + +template +PostProcessor3D* +SlipBoundaryProcessorGenerator3D::generate() const +{ + return new SlipBoundaryProcessor3D + ( this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, discreteNormalX, discreteNormalY, discreteNormalZ); +} + +template +PostProcessorGenerator3D* +SlipBoundaryProcessorGenerator3D::clone() const +{ + return new SlipBoundaryProcessorGenerator3D + (this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, discreteNormalX, discreteNormalY, discreteNormalZ); +} + +//////// PartialSlipBoundaryProcessor3D //////////////////////////////// + +template +PartialSlipBoundaryProcessor3D:: +PartialSlipBoundaryProcessor3D(T tuner_, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX, int discreteNormalY, int discreteNormalZ) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), tuner(tuner_) +{ + OLB_PRECONDITION(x0==x1 || y0==y1 || z0==z1); + reflectionPop[0] = 0; + for (int iPop = 1; iPop < DESCRIPTOR::q; iPop++) { + reflectionPop[iPop] = 0; + // iPop are the directions which pointing into the fluid, discreteNormal is pointing outwarts + if (descriptors::c(iPop,0)*discreteNormalX + descriptors::c(iPop,1)*discreteNormalY + descriptors::c(iPop,2)*discreteNormalZ < 0) { + //std::cout << "-----" <(iPop,0) - mult*(descriptors::c(iPop,0)*discreteNormalX + descriptors::c(iPop,1)*discreteNormalY + descriptors::c(iPop,2)*discreteNormalZ)*discreteNormalX); + mirrorDirection1 = (descriptors::c(iPop,1) - mult*(descriptors::c(iPop,0)*discreteNormalX + descriptors::c(iPop,1)*discreteNormalY + descriptors::c(iPop,2)*discreteNormalZ)*discreteNormalY); + mirrorDirection2 = (descriptors::c(iPop,2) - mult*(descriptors::c(iPop,0)*discreteNormalX + descriptors::c(iPop,1)*discreteNormalY + descriptors::c(iPop,2)*discreteNormalZ)*discreteNormalZ); + + // bounce back for the case discreteNormalX = discreteNormalY = discreteNormalZ = 1, that is mult=0 + if (mult == 0) { + mirrorDirection0 = -descriptors::c(iPop,0); + mirrorDirection1 = -descriptors::c(iPop,1); + mirrorDirection2 = -descriptors::c(iPop,2); + } + + // computes mirror jPop + for (reflectionPop[iPop] = 1; reflectionPop[iPop] < DESCRIPTOR::q ; reflectionPop[iPop]++) { + if (descriptors::c(reflectionPop[iPop],0)==mirrorDirection0 && descriptors::c(reflectionPop[iPop],1)==mirrorDirection1 && descriptors::c(reflectionPop[iPop],2)==mirrorDirection2) { + break; + } + } + //std::cout < +void PartialSlipBoundaryProcessor3D:: +processSubDomain(BlockLattice3D& blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) +{ + int newX0, newX1, newY0, newY1, newZ0, newZ1; + if ( util::intersect ( + x0, x1, y0, y1, z0, z1, + x0_, x1_, y0_, y1_, z0_, z1_, + newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) { + + int iX; +#ifdef PARALLEL_MODE_OMP + #pragma omp parallel for +#endif + for (iX=newX0; iX<=newX1; ++iX) { + for (int iY=newY0; iY<=newY1; ++iY) { + for (int iZ=newZ0; iZ<=newZ1; ++iZ) { + for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) { + if (reflectionPop[iPop]!=0) { + //do reflection + blockLattice.get(iX,iY,iZ)[iPop] = tuner*blockLattice.get(iX,iY,iZ)[reflectionPop[iPop]]; + } + } + for (int iPop = 1; iPop < DESCRIPTOR::q/2 ; ++iPop) { + T provv = blockLattice.get(iX,iY,iZ)[descriptors::opposite(iPop)]; + blockLattice.get(iX,iY,iZ)[descriptors::opposite(iPop)] += + (1.-tuner)*blockLattice.get(iX,iY,iZ)[iPop]; + blockLattice.get(iX,iY,iZ)[iPop] += (1.-tuner)*provv; + } + } + } + } + } +} + +template +void PartialSlipBoundaryProcessor3D:: +process(BlockLattice3D& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1); +} + +//////// PartialSlipBoundaryProcessorGenerator3D //////////////////////////////// + +template +PartialSlipBoundaryProcessorGenerator3D:: +PartialSlipBoundaryProcessorGenerator3D(T tuner_, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX_, int discreteNormalY_, int discreteNormalZ_) + : PostProcessorGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_), discreteNormalX(discreteNormalX_), discreteNormalY(discreteNormalY_), discreteNormalZ(discreteNormalZ_), tuner(tuner_) +{ } + +template +PostProcessor3D* +PartialSlipBoundaryProcessorGenerator3D::generate() const +{ + return new PartialSlipBoundaryProcessor3D + (tuner, this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, discreteNormalX, discreteNormalY, discreteNormalZ); +} + +template +PostProcessorGenerator3D* +PartialSlipBoundaryProcessorGenerator3D::clone() const +{ + return new PartialSlipBoundaryProcessorGenerator3D + (tuner, this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, discreteNormalX, discreteNormalY, discreteNormalZ); +} + +//////// FreeEnergyWallProcessor3D //////////////////////////////// + +template +FreeEnergyWallProcessor3D::FreeEnergyWallProcessor3D( + int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, + int discreteNormalX_, int discreteNormalY_, int discreteNormalZ_, T addend_) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), discreteNormalX(discreteNormalX_), + discreteNormalY(discreteNormalY_), discreteNormalZ(discreteNormalZ_), addend(addend_) +{ } + +template +void FreeEnergyWallProcessor3D:: +processSubDomain(BlockLattice3D& blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) +{ + int newX0, newX1, newY0, newY1, newZ0, newZ1; + if ( util::intersect ( + x0, x1, y0, y1, z0, z1, + x0_, x1_, y0_, y1_, z0_, z1_, + newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) { + + for (int iX=newX0; iX<=newX1; ++iX) { + for (int iY=newY0; iY<=newY1; ++iY) { + for (int iZ=newZ0; iZ<=newZ1; ++iZ) { + T rhoBulk = blockLattice.get(iX-discreteNormalX, iY-discreteNormalY, iZ-discreteNormalZ).computeRho(); + T rhoTmp = 0; + for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) { + rhoTmp += blockLattice.get(iX,iY,iZ)[iPop]; + } + T rho = rhoBulk + addend; + blockLattice.get(iX,iY,iZ)[0] = rho - rhoTmp - 1; + } + } + } + } +} + +template +void FreeEnergyWallProcessor3D:: +process(BlockLattice3D& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1); +} + +//////// FreeEnergyWallProcessorGenerator3D //////////////////////////////// + +template +FreeEnergyWallProcessorGenerator3D::FreeEnergyWallProcessorGenerator3D( + int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, + int discreteNormalX_, int discreteNormalY_, int discreteNormalZ_, T addend_) + : PostProcessorGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_), discreteNormalX(discreteNormalX_), + discreteNormalY(discreteNormalY_), discreteNormalZ(discreteNormalZ_), addend(addend_) +{ } + +template +PostProcessor3D* +FreeEnergyWallProcessorGenerator3D::generate() const +{ + return new FreeEnergyWallProcessor3D + (this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, + discreteNormalX, discreteNormalY, discreteNormalZ, addend); +} + +template +PostProcessorGenerator3D* +FreeEnergyWallProcessorGenerator3D::clone() const +{ + return new FreeEnergyWallProcessorGenerator3D + (this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, + discreteNormalX, discreteNormalY, discreteNormalZ, addend); +} + + +//////// FreeEnergyChemPotBoundaryProcessor3D //////////////////////////////// + +template +FreeEnergyChemPotBoundaryProcessor3D:: +FreeEnergyChemPotBoundaryProcessor3D( + int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX_, + int discreteNormalY_, int discreteNormalZ_, int latticeNumber_) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), discreteNormalX(discreteNormalX_), + discreteNormalY(discreteNormalY_), discreteNormalZ(discreteNormalZ_), latticeNumber(latticeNumber_) +{ } + +template +void FreeEnergyChemPotBoundaryProcessor3D:: +processSubDomain( + BlockLattice3D& blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) +{ + int newX0, newX1, newY0, newY1, newZ0, newZ1; + if ( util::intersect ( + x0, x1, y0, y1, z0, z1, + x0_, x1_, y0_, y1_, z0_, z1_, + newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) { + + for (int iX=newX0; iX<=newX1; ++iX) { + for (int iY=newY0; iY<=newY1; ++iY) { + for (int iZ=newZ0; iZ<=newZ1; ++iZ) { + blockLattice.get(iX,iY,iZ).template setField( + blockLattice.get(iX-discreteNormalX, iY-discreteNormalY, iZ-discreteNormalZ).template getField() + ); + if (latticeNumber == 1) { + T rho0 = blockLattice.get(iX,iY,iZ).computeRho(); + T rho1 = blockLattice.get(iX-discreteNormalX, iY-discreteNormalY, iZ-discreteNormalZ).computeRho(); + *(blockLattice.get(iX,iY,iZ).template getFieldPointer()) += (rho1 / rho0 - 1) / descriptors::invCs2(); + } + } + } + } + } +} + +template +void FreeEnergyChemPotBoundaryProcessor3D:: +process(BlockLattice3D& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1); +} + +//////// FreeEnergyChemPotBoundaryProcessorGenerator3D //////////////////////////////// + +template +FreeEnergyChemPotBoundaryProcessorGenerator3D:: +FreeEnergyChemPotBoundaryProcessorGenerator3D( + int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, + int discreteNormalX_, int discreteNormalY_, int discreteNormalZ_, int latticeNumber_) + : PostProcessorGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_), discreteNormalX(discreteNormalX_), + discreteNormalY(discreteNormalY_), discreteNormalZ(discreteNormalZ_), latticeNumber(latticeNumber_) +{ } + +template +PostProcessor3D* +FreeEnergyChemPotBoundaryProcessorGenerator3D::generate() const +{ + return new FreeEnergyChemPotBoundaryProcessor3D + (this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, + discreteNormalX, discreteNormalY, discreteNormalZ, latticeNumber); +} + +template +PostProcessorGenerator3D* +FreeEnergyChemPotBoundaryProcessorGenerator3D::clone() const +{ + return new FreeEnergyChemPotBoundaryProcessorGenerator3D + (this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, + discreteNormalX, discreteNormalY, discreteNormalZ, latticeNumber); +} + + +//////// FreeEnergyConvectiveProcessor3D //////////////////////////////// + +template +FreeEnergyConvectiveProcessor3D:: +FreeEnergyConvectiveProcessor3D( + int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, + int discreteNormalX_, int discreteNormalY_, int discreteNormalZ_) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), discreteNormalX(discreteNormalX_), + discreteNormalY(discreteNormalY_), discreteNormalZ(discreteNormalZ_) +{ } + +template +void FreeEnergyConvectiveProcessor3D:: +processSubDomain( + BlockLattice3D& blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) +{ + int newX0, newX1, newY0, newY1, newZ0, newZ1; + if ( util::intersect ( + x0, x1, y0, y1, z0, z1, + x0_, x1_, y0_, y1_, z0_, z1_, + newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) { + + for (int iX=newX0; iX<=newX1; ++iX) { + for (int iY=newY0; iY<=newY1; ++iY) { + for (int iZ=newZ0; iZ<=newZ1; ++iZ) { + T rho, uPerp, rho0, rho1, u[3]; + rho0 = blockLattice.get(iX,iY,iZ).computeRho(); + blockLattice.get(iX-discreteNormalX,iY-discreteNormalY,iZ-discreteNormalZ).computeRhoU(rho1, u); + + if (discreteNormalZ == 0) { + if (discreteNormalY == 0) { + if (discreteNormalX < 0) { + uPerp = u[0]; + } + else { + uPerp = -u[0]; + } + } + else if (discreteNormalX == 0) { + if (discreteNormalY < 0) { + uPerp = u[1]; + } + else { + uPerp = -u[1]; + } + } + else { + uPerp = sqrt(u[0] * u[0] + u[1] * u[1]); + } + } + else if (discreteNormalY == 0) { + if (discreteNormalX == 0) { + if (discreteNormalZ < 0) { + uPerp = u[2]; + } + else { + uPerp = -u[2]; + } + } + else { + uPerp = sqrt(u[0] * u[0] + u[2] * u[2]); + } + } + else if (discreteNormalX == 0) { + uPerp = sqrt(u[1] * u[1] + u[2] * u[2]); + } + else { + uPerp = sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]); + } + + rho = (rho0 + uPerp * rho1) / (1. + uPerp); + blockLattice.get(iX,iY,iZ).defineRho(rho); + } + } + } + } +} + +template +void FreeEnergyConvectiveProcessor3D:: +process(BlockLattice3D& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1); +} + +//////// FreeEnergyConvectiveProcessorGenerator3D //////////////////////////////// + +template +FreeEnergyConvectiveProcessorGenerator3D:: +FreeEnergyConvectiveProcessorGenerator3D( + int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, + int discreteNormalX_, int discreteNormalY_, int discreteNormalZ_) + : PostProcessorGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_), discreteNormalX(discreteNormalX_), + discreteNormalY(discreteNormalY_), discreteNormalZ(discreteNormalZ_) +{ } + +template +PostProcessor3D* +FreeEnergyConvectiveProcessorGenerator3D::generate() const +{ + return new FreeEnergyConvectiveProcessor3D + (this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, + discreteNormalX, discreteNormalY, discreteNormalZ); +} + +template +PostProcessorGenerator3D* +FreeEnergyConvectiveProcessorGenerator3D::clone() const +{ + return new FreeEnergyConvectiveProcessorGenerator3D + (this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, + discreteNormalX, discreteNormalY, discreteNormalZ); +} + +} // namespace olb + +#endif -- cgit v1.2.3