/* 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