/* This file is part of the OpenLB library
*
* Copyright (C) 2007 Orestis Malaspinas, 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 EXTENDED_FINITE_DIFFERENCE_BOUNDARY_3D_HH
#define EXTENDED_FINITE_DIFFERENCE_BOUNDARY_3D_HH
#include "extendedFiniteDifferenceBoundary3D.h"
#include "core/finiteDifference3D.h"
#include "core/blockLattice3D.h"
#include "core/util.h"
#include "dynamics/lbHelpers.h"
#include "dynamics/firstOrderLbHelpers.h"
#include "boundaryInstantiator3D.h"
#include "momentaOnBoundaries3D.h"
namespace olb {
//////// ExtendedFdPlaneBoundaryPostProcessor3D ///////////////////////////////////
template
ExtendedFdPlaneBoundaryPostProcessor3D ::
ExtendedFdPlaneBoundaryPostProcessor3D(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 ExtendedFdPlaneBoundaryPostProcessor3D::
processSubDomain(BlockLattice3D& blockLattice,
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
{
typedef DESCRIPTOR L;
using namespace olb::util::tensorIndices3D;
typedef lbHelpers lbH;
enum {x,y,z};
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) {
Cell& cell = blockLattice.get(iX,iY,iZ);
T rho, u[L::d];
cell.computeRhoU(rho,u);
T dx_U[DESCRIPTOR::d], dy_U[DESCRIPTOR::d], dz_U[DESCRIPTOR::d];
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 rhoGradU[L::d][L::d];
rhoGradU[x][x] = rho *dx_U[x];
rhoGradU[x][y] = rho *dx_U[y];
rhoGradU[x][z] = rho *dx_U[z];
rhoGradU[y][x] = rho *dy_U[x];
rhoGradU[y][y] = rho *dy_U[y];
rhoGradU[y][z] = rho *dy_U[z];
rhoGradU[z][x] = rho *dz_U[x];
rhoGradU[z][y] = rho *dz_U[y];
rhoGradU[z][z] = rho *dz_U[z];
T omega = blockLattice.getDynamics(iX, iY, iZ) -> getOmega();
T sToPi = - (T)1 / descriptors::invCs2() / omega;
T pi[util::TensorVal::n];
pi[xx] = (T)2 * rhoGradU[x][x] * sToPi;
pi[yy] = (T)2 * rhoGradU[y][y] * sToPi;
pi[zz] = (T)2 * rhoGradU[z][z] * sToPi;
pi[xy] = (rhoGradU[x][y] + rhoGradU[y][x]) * sToPi;
pi[xz] = (rhoGradU[x][z] + rhoGradU[z][x]) * sToPi;
pi[yz] = (rhoGradU[y][z] + rhoGradU[z][y]) * sToPi;
// here ends the "regular" fdBoudaryCondition
// implemented in OpenLB
T uSqr = util::normSqr(u);
// first we compute the term
// (c_{i\alpha} \nabla_\beta)(rho*u_\alpha*u_\beta)
T dx_rho, dy_rho, dz_rho;
interpolateGradients<0>(blockLattice, dx_rho, iX, iY, iZ);
interpolateGradients<1>(blockLattice, dy_rho, iX, iY, iZ);
interpolateGradients<2>(blockLattice, dz_rho, iX, iY, iZ);
for (int iPop = 0; iPop < L::q; ++iPop) {
T cGradRhoUU = T();
for (int iAlpha=0; iAlpha < L::d; ++iAlpha) {
cGradRhoUU += descriptors::c(iPop,iAlpha) * (
dx_rho*u[iAlpha]*u[x] +
dx_U[iAlpha]*rho*u[x] +
dx_U[x]*rho*u[iAlpha] + //end of dx derivative
dy_rho*u[iAlpha]*u[y] +
dy_U[iAlpha]*rho*u[y] +
dy_U[y]*rho*u[iAlpha] +//end of dy derivative
dz_rho*u[iAlpha]*u[z] +
dz_U[iAlpha]*rho*u[z] +
dz_U[z]*rho*u[iAlpha]);
}
// then we compute the term
// c_{i\gamma}\nabla_{\gamma}(\rho*u_\alpha * u_\beta)
T cDivRhoUU[L::d][L::d]; //first step towards QcdivRhoUU
for (int iAlpha = 0; iAlpha < L::d; ++iAlpha) {
for (int iBeta = 0; iBeta < L::d; ++iBeta) {
cDivRhoUU[iAlpha][iBeta] = descriptors::c(iPop,x)*
(dx_rho*u[iAlpha]*u[iBeta] +
dx_U[iAlpha]*rho*u[iBeta] +
dx_U[iBeta]*rho*u[iAlpha])
+descriptors::c(iPop,y)*
(dy_rho*u[iAlpha]*u[iBeta] +
dy_U[iAlpha]*rho*u[iBeta] +
dy_U[iBeta]*rho*u[iAlpha])
+descriptors::c(iPop,z)*
(dz_rho*u[iAlpha]*u[iBeta] +
dz_U[iAlpha]*rho*u[iBeta] +
dz_U[iBeta]*rho*u[iAlpha]);
}
}
//Finally we can compute
// Q_{i\alpha\beta}c_{i\gamma}\nabla_{\gamma}(\rho*u_\alpha * u_\beta)
// and Q_{i\alpha\beta}\rho\nabla_{\alpha}u_\beta
T qCdivRhoUU = T();
T qRhoGradU = T();
for (int iAlpha = 0; iAlpha < L::d; ++iAlpha) {
for (int iBeta = 0; iBeta < L::d; ++iBeta) {
int ci_ci = descriptors::c(iPop,iAlpha)*descriptors::c(iPop,iBeta);
qCdivRhoUU += ci_ci * cDivRhoUU[iAlpha][iBeta];
qRhoGradU += ci_ci * rhoGradU[iAlpha][iBeta];
if (iAlpha == iBeta) {
qCdivRhoUU -= cDivRhoUU[iAlpha][iBeta]/descriptors::invCs2();
qRhoGradU -= rhoGradU[iAlpha][iBeta]/descriptors::invCs2();
}
}
}
// we then can reconstruct the value of the populations
// according to the complete C-E expansion term
cell[iPop] = lbH::equilibrium(iPop,rho,u,uSqr)
- descriptors::t(iPop) * descriptors::invCs2() / omega
* (qRhoGradU - cGradRhoUU + 0.5*descriptors::invCs2()*qCdivRhoUU);
}
}
}
}
}
}
template
void ExtendedFdPlaneBoundaryPostProcessor3D::
process(BlockLattice3D& blockLattice)
{
processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
}
template
template
void ExtendedFdPlaneBoundaryPostProcessor3D::
interpolateGradients(BlockLattice3D const& blockLattice,
T velDeriv[DESCRIPTOR::d],
int iX, int iY, int iZ) const
{
fd::DirectedGradients3D::
interpolateVector(velDeriv, blockLattice, iX, iY, iZ);
}
template
template
void ExtendedFdPlaneBoundaryPostProcessor3D::
interpolateGradients(BlockLattice3D const& blockLattice,
T& rhoDeriv, int iX, int iY, int iZ) const
{
fd::DirectedGradients3D::
interpolateScalar(rhoDeriv, blockLattice, iX, iY, iZ);
}
//////// ExtendedFdPlaneBoundaryProcessorGenertor3D ///////////////////////////////
template
ExtendedFdPlaneBoundaryProcessorGenerator3D::
ExtendedFdPlaneBoundaryProcessorGenerator3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
: PostProcessorGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_)
{ }
template
PostProcessor3D*
ExtendedFdPlaneBoundaryProcessorGenerator3D::generate() const
{
return new ExtendedFdPlaneBoundaryPostProcessor3D
(this->x0, this->x1, this->y0, this->y1, this->z0, this->z1);
}
template
PostProcessorGenerator3D*
ExtendedFdPlaneBoundaryProcessorGenerator3D::clone() const
{
return new ExtendedFdPlaneBoundaryProcessorGenerator3D
(this->x0, this->x1, this->y0, this->y1, this->z0, this->z1);
}
////////// Class ExtendedFdBoundaryManager3D /////////////////////////////////////////
template
class ExtendedFdBoundaryManager3D {
public:
template static Momenta*
getVelocityBoundaryMomenta();
template static Dynamics*
getVelocityBoundaryDynamics(T omega, Momenta& momenta);
template static PostProcessorGenerator3D*
getVelocityBoundaryProcessor(int x0, int x1, int y0, int y1, int z0, int z1);
template static Momenta*
getPressureBoundaryMomenta();
template static Dynamics*
getPressureBoundaryDynamics(T omega, Momenta& momenta);
template static PostProcessorGenerator3D*
getPressureBoundaryProcessor(int x0, int x1, int y0, int y1, int z0, int z1);
template static PostProcessorGenerator3D*
getConvectionBoundaryProcessor(int x0, int x1, int y0, int y1, int z0, int z1, T* uAv=NULL);
template static Momenta*
getExternalVelocityEdgeMomenta();
template static Dynamics*
getExternalVelocityEdgeDynamics(T omega, Momenta& momenta);
template static PostProcessorGenerator3D*
getExternalVelocityEdgeProcessor(int x0, int x1, int y0, int y1, int z0, int z1);
template static Momenta*
getInternalVelocityEdgeMomenta();
template static Dynamics*
getInternalVelocityEdgeDynamics(T omega, Momenta& momenta);
template static PostProcessorGenerator3D*
getInternalVelocityEdgeProcessor(int x0, int x1, int y0, int y1, int z0, int z1);
template static Momenta*
getExternalVelocityCornerMomenta();
template static Dynamics*
getExternalVelocityCornerDynamics(T omega, Momenta& momenta);
template static PostProcessorGenerator3D*
getExternalVelocityCornerProcessor(int x, int y, int z);
template static Momenta*
getInternalVelocityCornerMomenta();
template static Dynamics*
getInternalVelocityCornerDynamics(T omega, Momenta& momenta);
template static PostProcessorGenerator3D*
getInternalVelocityCornerProcessor(int x, int y, int z);
};
template
template
Momenta*
ExtendedFdBoundaryManager3D::getVelocityBoundaryMomenta()
{
return new BasicDirichletBM;
}
template
template
Dynamics* ExtendedFdBoundaryManager3D::
getVelocityBoundaryDynamics(T omega, Momenta& momenta)
{
return new MixinDynamics(omega, momenta);
}
template
template
PostProcessorGenerator3D*
ExtendedFdBoundaryManager3D::
getVelocityBoundaryProcessor(int x0, int x1, int y0, int y1, int z0, int z1)
{
return new ExtendedFdPlaneBoundaryProcessorGenerator3D
(x0,x1, y0,y1, z0,z1);
}
template
template
Momenta*
ExtendedFdBoundaryManager3D::getPressureBoundaryMomenta()
{
return new BasicDirichletBM;
}
template
template
Dynamics* ExtendedFdBoundaryManager3D::
getPressureBoundaryDynamics(T omega, Momenta& momenta)
{
return new MixinDynamics(omega, momenta);
}
template
template
PostProcessorGenerator3D*
ExtendedFdBoundaryManager3D::
getPressureBoundaryProcessor(int x0, int x1, int y0, int y1, int z0, int z1)
{
return new ExtendedFdPlaneBoundaryProcessorGenerator3D
(x0,x1, y0,y1, z0,z1);
}
template
template
PostProcessorGenerator3D*
ExtendedFdBoundaryManager3D::
getConvectionBoundaryProcessor(int x0, int x1, int y0, int y1, int z0, int z1, T* uAv)
{
return nullptr;
}
template
template
Momenta*
ExtendedFdBoundaryManager3D::getExternalVelocityEdgeMomenta()
{
return new FixedVelocityBM;
}
template
template
Dynamics*
ExtendedFdBoundaryManager3D::
getExternalVelocityEdgeDynamics(T omega, Momenta& momenta)
{
return new MixinDynamics(omega, momenta);
}
template
template
PostProcessorGenerator3D*
ExtendedFdBoundaryManager3D::
getExternalVelocityEdgeProcessor(int x0, int x1, int y0, int y1, int z0, int z1)
{
return new OuterVelocityEdgeProcessorGenerator3D(x0,x1, y0,y1, z0,z1);
}
template
template
Momenta*
ExtendedFdBoundaryManager3D::getInternalVelocityEdgeMomenta()
{
return new InnerEdgeVelBM3D;
}
template
template
Dynamics*
ExtendedFdBoundaryManager3D::
getInternalVelocityEdgeDynamics(T omega, Momenta& momenta)
{
return new CombinedRLBdynamics(omega, momenta);
}
template
template
PostProcessorGenerator3D*
ExtendedFdBoundaryManager3D::
getInternalVelocityEdgeProcessor(int x0, int x1, int y0, int y1, int z0, int z1)
{
return nullptr;
}
template
template
Momenta*
ExtendedFdBoundaryManager3D::getExternalVelocityCornerMomenta()
{
return new FixedVelocityBM;
}
template
template
Dynamics*
ExtendedFdBoundaryManager3D::
getExternalVelocityCornerDynamics(T omega, Momenta& momenta)
{
return new MixinDynamics(omega, momenta);
}
template
template
PostProcessorGenerator3D*
ExtendedFdBoundaryManager3D::
getExternalVelocityCornerProcessor(int x, int y, int z)
{
return new OuterVelocityCornerProcessorGenerator3D (x,y,z);
}
template
template
Momenta*
ExtendedFdBoundaryManager3D::getInternalVelocityCornerMomenta()
{
return new InnerCornerVelBM3D;
}
template
template
Dynamics*
ExtendedFdBoundaryManager3D::
getInternalVelocityCornerDynamics(T omega, Momenta& momenta)
{
return new CombinedRLBdynamics(omega, momenta);
}
template
template
PostProcessorGenerator3D*
ExtendedFdBoundaryManager3D::
getInternalVelocityCornerProcessor(int x, int y, int z)
{
return nullptr;
}
////////// Factory function //////////////////////////////////////////////////
template
OnLatticeBoundaryCondition3D*
createExtendedFdBoundaryCondition3D(BlockLatticeStructure3D& block)
{
return new BoundaryConditionInstantiator3D <
T, DESCRIPTOR,
ExtendedFdBoundaryManager3D > (block);
}
} // namespace olb
#endif