/* This file is part of the OpenLB library
*
* Copyright (C) 2018 Sam Avis, Robin Trunk
* 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 FREE_ENERGY_POST_PROCESSOR_3D_HH
#define FREE_ENERGY_POST_PROCESSOR_3D_HH
#include "freeEnergyPostProcessor3D.h"
#include "core/blockLattice3D.h"
namespace olb {
//////// FreeEnergyChemicalPotentialCoupling3D ///////////////////////////////////
template
FreeEnergyChemicalPotentialCoupling3D ::FreeEnergyChemicalPotentialCoupling3D (
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, T alpha_, T kappa1_, T kappa2_, T kappa3_,
std::vector partners_)
: x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), alpha(alpha_),
kappa1(kappa1_), kappa2(kappa2_), kappa3(kappa3_), partners(partners_)
{ }
template
FreeEnergyChemicalPotentialCoupling3D ::FreeEnergyChemicalPotentialCoupling3D (
T alpha_, T kappa1_, T kappa2_, T kappa3_,
std::vector partners_)
: x0(0), x1(0), y0(0), y1(0), z0(0), z1(0), alpha(alpha_),
kappa1(kappa1_), kappa2(kappa2_), kappa3(kappa3_), partners(partners_)
{ }
template
void FreeEnergyChemicalPotentialCoupling3D::processSubDomain (
BlockLattice3D& blockLattice,
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_ )
{
// If partners.size() == 1: two fluid components
// If partners.size() == 2: three fluid components
BlockLattice3D *partnerLattice1 = dynamic_cast *>(partners[0]);
BlockLattice3D *partnerLattice2 = 0;
if (partners.size() > 1) {
partnerLattice2 = dynamic_cast *>(partners[1]);
}
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 nx = newX1-newX0+3; // include a one-cell boundary
int ny = newY1-newY0+3; // include a one-cell boundary
int nz = newZ1-newZ0+3; // include a one-cell boundary
int offsetX = newX0-1;
int offsetY = newY0-1;
int offsetZ = newZ0-1;
// compute the density fields for each lattice
BlockData3D rhoField1(nx, ny, nz);
BlockData3D rhoField2(nx, ny, nz);
BlockData3D rhoField3(nx, ny, nz);
for (int iX=newX0-1; iX<=newX1+1; ++iX)
for (int iY=newY0-1; iY<=newY1+1; ++iY)
for (int iZ=newZ0-1; iZ<=newZ1+1; ++iZ) {
rhoField1.get(iX-offsetX, iY-offsetY, iZ-offsetZ) = blockLattice.get(iX,iY,iZ).computeRho();
}
for (int iX=newX0-1; iX<=newX1+1; ++iX)
for (int iY=newY0-1; iY<=newY1+1; ++iY)
for (int iZ=newZ0-1; iZ<=newZ1+1; ++iZ) {
rhoField2.get(iX-offsetX, iY-offsetY, iZ-offsetZ) = partnerLattice1->get(iX,iY,iZ).computeRho();
}
if (partners.size() > 1) {
for (int iX=newX0-1; iX<=newX1+1; ++iX)
for (int iY=newY0-1; iY<=newY1+1; ++iY)
for (int iZ=newZ0-1; iZ<=newZ1+1; ++iZ) {
rhoField3.get(iX-offsetX, iY-offsetY, iZ-offsetZ) = partnerLattice2->get(iX,iY,iZ).computeRho();
}
}
// calculate chemical potential
for (int iX=newX0; iX<=newX1; ++iX) {
for (int iY=newY0; iY<=newY1; ++iY) {
for (int iZ=newZ0; iZ<=newZ1; ++iZ) {
T densitySum = rhoField1.get(iX-offsetX, iY-offsetY, iZ-offsetZ)
+ rhoField2.get(iX-offsetX, iY-offsetY, iZ-offsetZ);
T densityDifference = rhoField1.get(iX-offsetX, iY-offsetY, iZ-offsetZ)
- rhoField2.get(iX-offsetX, iY-offsetY, iZ-offsetZ);
if (partners.size() > 1) {
densitySum -= rhoField3.get(iX-offsetX, iY-offsetY, iZ-offsetZ);
densityDifference -= rhoField3.get(iX-offsetX, iY-offsetY, iZ-offsetZ);
}
T term1 = 0.125 * kappa1 * (densitySum)
* (densitySum-1.) * (densitySum-2.);
T term2 = 0.125 * kappa2 * (densityDifference)
* (densityDifference-1.) * (densityDifference-2.);
T term3 = 0.;
if (partners.size() > 1) {
T rho3 = rhoField3.get(iX-offsetX, iY-offsetY, iZ-offsetZ);
term3 = kappa3 * rho3 * (rho3 - 1.) * (2.*rho3 - 1.);
}
T laplaceRho1 = 1.0 / 6.0 * (
rhoField1.get(iX-offsetX, iY-offsetY-1, iZ-offsetZ-1)
+ rhoField1.get(iX-offsetX-1, iY-offsetY, iZ-offsetZ-1)
+ 2. * rhoField1.get(iX-offsetX, iY-offsetY, iZ-offsetZ-1)
+ rhoField1.get(iX-offsetX+1, iY-offsetY, iZ-offsetZ-1)
+ rhoField1.get(iX-offsetX, iY-offsetY+1, iZ-offsetZ-1)
+ rhoField1.get(iX-offsetX-1, iY-offsetY-1, iZ-offsetZ)
+ 2. * rhoField1.get(iX-offsetX, iY-offsetY-1, iZ-offsetZ)
+ rhoField1.get(iX-offsetX+1, iY-offsetY-1, iZ-offsetZ)
+ 2. * rhoField1.get(iX-offsetX-1, iY-offsetY, iZ-offsetZ)
-24. * rhoField1.get(iX-offsetX, iY-offsetY, iZ-offsetZ)
+ 2. * rhoField1.get(iX-offsetX+1, iY-offsetY, iZ-offsetZ)
+ rhoField1.get(iX-offsetX-1, iY-offsetY+1, iZ-offsetZ)
+ 2. * rhoField1.get(iX-offsetX, iY-offsetY+1, iZ-offsetZ)
+ rhoField1.get(iX-offsetX+1, iY-offsetY+1, iZ-offsetZ)
+ rhoField1.get(iX-offsetX, iY-offsetY-1, iZ-offsetZ+1)
+ rhoField1.get(iX-offsetX-1, iY-offsetY, iZ-offsetZ+1)
+ 2. * rhoField1.get(iX-offsetX, iY-offsetY, iZ-offsetZ+1)
+ rhoField1.get(iX-offsetX+1, iY-offsetY, iZ-offsetZ+1)
+ rhoField1.get(iX-offsetX, iY-offsetY+1, iZ-offsetZ+1)
);
T laplaceRho2 = 1.0 / 6.0 * (
rhoField2.get(iX-offsetX, iY-offsetY-1, iZ-offsetZ-1)
+ rhoField2.get(iX-offsetX-1, iY-offsetY, iZ-offsetZ-1)
+ 2. * rhoField2.get(iX-offsetX, iY-offsetY, iZ-offsetZ-1)
+ rhoField2.get(iX-offsetX+1, iY-offsetY, iZ-offsetZ-1)
+ rhoField2.get(iX-offsetX, iY-offsetY+1, iZ-offsetZ-1)
+ rhoField2.get(iX-offsetX-1, iY-offsetY-1, iZ-offsetZ)
+ 2. * rhoField2.get(iX-offsetX, iY-offsetY-1, iZ-offsetZ)
+ rhoField2.get(iX-offsetX+1, iY-offsetY-1, iZ-offsetZ)
+ 2. * rhoField2.get(iX-offsetX-1, iY-offsetY, iZ-offsetZ)
-24. * rhoField2.get(iX-offsetX, iY-offsetY, iZ-offsetZ)
+ 2. * rhoField2.get(iX-offsetX+1, iY-offsetY, iZ-offsetZ)
+ rhoField2.get(iX-offsetX-1, iY-offsetY+1, iZ-offsetZ)
+ 2. * rhoField2.get(iX-offsetX, iY-offsetY+1, iZ-offsetZ)
+ rhoField2.get(iX-offsetX+1, iY-offsetY+1, iZ-offsetZ)
+ rhoField2.get(iX-offsetX, iY-offsetY-1, iZ-offsetZ+1)
+ rhoField2.get(iX-offsetX-1, iY-offsetY, iZ-offsetZ+1)
+ 2. * rhoField2.get(iX-offsetX, iY-offsetY, iZ-offsetZ+1)
+ rhoField2.get(iX-offsetX+1, iY-offsetY, iZ-offsetZ+1)
+ rhoField2.get(iX-offsetX, iY-offsetY+1, iZ-offsetZ+1)
);
T laplaceRho3 = 0.;
if (partners.size() > 1) {
laplaceRho3 = 1.0 / 6.0 * (
rhoField3.get(iX-offsetX, iY-offsetY-1, iZ-offsetZ-1)
+ rhoField3.get(iX-offsetX-1, iY-offsetY, iZ-offsetZ-1)
+ 2. * rhoField3.get(iX-offsetX, iY-offsetY, iZ-offsetZ-1)
+ rhoField3.get(iX-offsetX+1, iY-offsetY, iZ-offsetZ-1)
+ rhoField3.get(iX-offsetX, iY-offsetY+1, iZ-offsetZ-1)
+ rhoField3.get(iX-offsetX-1, iY-offsetY-1, iZ-offsetZ)
+ 2. * rhoField3.get(iX-offsetX, iY-offsetY-1, iZ-offsetZ)
+ rhoField3.get(iX-offsetX+1, iY-offsetY-1, iZ-offsetZ)
+ 2. * rhoField3.get(iX-offsetX-1, iY-offsetY, iZ-offsetZ)
-24. * rhoField3.get(iX-offsetX, iY-offsetY, iZ-offsetZ)
+ 2. * rhoField3.get(iX-offsetX+1, iY-offsetY, iZ-offsetZ)
+ rhoField3.get(iX-offsetX-1, iY-offsetY+1, iZ-offsetZ)
+ 2. * rhoField3.get(iX-offsetX, iY-offsetY+1, iZ-offsetZ)
+ rhoField3.get(iX-offsetX+1, iY-offsetY+1, iZ-offsetZ)
+ rhoField3.get(iX-offsetX, iY-offsetY-1, iZ-offsetZ+1)
+ rhoField3.get(iX-offsetX-1, iY-offsetY, iZ-offsetZ+1)
+ 2. * rhoField3.get(iX-offsetX, iY-offsetY, iZ-offsetZ+1)
+ rhoField3.get(iX-offsetX+1, iY-offsetY, iZ-offsetZ+1)
+ rhoField3.get(iX-offsetX, iY-offsetY+1, iZ-offsetZ+1)
);
}
// setting chemical potential to the respective lattices
blockLattice.get(iX,iY,iZ).template setField(
term1 + term2
+ 0.25*alpha*alpha*( (kappa2 - kappa1) * laplaceRho2
+(kappa2 + kappa1) * (laplaceRho3 - laplaceRho1) )
);
partnerLattice1->get(iX, iY, iZ).template setField(
term1 - term2
+ 0.25*alpha*alpha*( (kappa2 - kappa1) * (laplaceRho1 - laplaceRho3)
-(kappa2 + kappa1) * laplaceRho2 )
);
if (partners.size() > 1) {
partnerLattice2->get(iX, iY, iZ).template setField(
- term1 - term2 + term3
+ 0.25*alpha*alpha*( (kappa2 + kappa1) * laplaceRho1
-(kappa2 - kappa1) * laplaceRho2
-(kappa2 + kappa1 + 4.*kappa3) * laplaceRho3 )
);
}
}
}
}
}
}
template
void FreeEnergyChemicalPotentialCoupling3D::process (
BlockLattice3D& blockLattice)
{
processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
}
//////// FreeEnergyForceCoupling3D ///////////////////////////////////
template
FreeEnergyForceCoupling3D ::FreeEnergyForceCoupling3D (
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
std::vector partners_)
: x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), partners(partners_)
{ }
template
FreeEnergyForceCoupling3D ::FreeEnergyForceCoupling3D (
std::vector partners_)
: x0(0), x1(0), y0(0), y1(0), z0(0), z1(0), partners(partners_)
{ }
template
void FreeEnergyForceCoupling3D::processSubDomain (
BlockLattice3D& blockLattice,
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_ )
{
// If partners.size() == 1: two fluid components
// If partners.size() == 2: three fluid components
BlockLattice3D *partnerLattice1 = dynamic_cast *>(partners[0]);
BlockLattice3D *partnerLattice2 = 0;
if (partners.size() > 1) {
partnerLattice2 = dynamic_cast *>(partners[1]);
}
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 phi = blockLattice.get(iX,iY,iZ).computeRho();
T rho = partnerLattice1->get(iX,iY,iZ).computeRho();
T gradMuPhiX = 1./12. * ( -blockLattice.get(iX-1,iY,iZ-1).template getField()
- blockLattice.get(iX-1,iY,iZ+1).template getField()
- blockLattice.get(iX-1,iY-1,iZ ).template getField()
- 2.* blockLattice.get(iX-1,iY,iZ ).template getField()
- blockLattice.get(iX-1,iY+1,iZ ).template getField()
+ blockLattice.get(iX+1,iY-1,iZ ).template getField()
+ 2.* blockLattice.get(iX+1,iY,iZ ).template getField()
+ blockLattice.get(iX+1,iY+1,iZ ).template getField()
+ blockLattice.get(iX+1,iY,iZ-1).template getField()
+ blockLattice.get(iX+1,iY,iZ+1).template getField() );
T gradMuPhiY = 1./12. * ( -blockLattice.get(iX-1,iY-1,iZ ).template getField()
- blockLattice.get(iX+1,iY-1,iZ ).template getField()
- blockLattice.get(iX,iY-1,iZ-1).template getField()
- 2.* blockLattice.get(iX,iY-1,iZ ).template getField()
- blockLattice.get(iX,iY-1,iZ+1).template getField()
+ blockLattice.get(iX,iY+1,iZ-1).template getField()
+ 2.* blockLattice.get(iX,iY+1,iZ ).template getField()
+ blockLattice.get(iX,iY+1,iZ+1).template getField()
+ blockLattice.get(iX-1,iY+1,iZ ).template getField()
+ blockLattice.get(iX+1,iY+1,iZ ).template getField() );
T gradMuPhiZ = 1./12. * ( -blockLattice.get(iX,iY-1,iZ-1).template getField()
- blockLattice.get(iX,iY+1,iZ-1).template getField()
- blockLattice.get(iX-1,iY,iZ-1).template getField()
- 2.* blockLattice.get(iX,iY,iZ-1).template getField()
- blockLattice.get(iX+1,iY,iZ-1).template getField()
+ blockLattice.get(iX-1,iY,iZ+1).template getField()
+ 2.* blockLattice.get(iX,iY,iZ+1).template getField()
+ blockLattice.get(iX+1,iY,iZ+1).template getField()
+ blockLattice.get(iX,iY-1,iZ+1).template getField()
+ blockLattice.get(iX,iY+1,iZ+1).template getField() );
T gradMuRhoX = 1./12. * ( -partnerLattice1->get(iX-1,iY,iZ-1).template getField()
- partnerLattice1->get(iX-1,iY,iZ+1).template getField()
- partnerLattice1->get(iX-1,iY-1,iZ ).template getField()
- 2.* partnerLattice1->get(iX-1,iY,iZ ).template getField()
- partnerLattice1->get(iX-1,iY+1,iZ ).template getField()
+ partnerLattice1->get(iX+1,iY-1,iZ ).template getField()
+ 2.* partnerLattice1->get(iX+1,iY,iZ ).template getField()
+ partnerLattice1->get(iX+1,iY+1,iZ ).template getField()
+ partnerLattice1->get(iX+1,iY,iZ-1).template getField()
+ partnerLattice1->get(iX+1,iY,iZ+1).template getField() );
T gradMuRhoY = 1./12. * ( -partnerLattice1->get(iX-1,iY-1,iZ ).template getField()
- partnerLattice1->get(iX+1,iY-1,iZ ).template getField()
- partnerLattice1->get(iX,iY-1,iZ-1).template getField()
- 2.* partnerLattice1->get(iX,iY-1,iZ ).template getField()
- partnerLattice1->get(iX,iY-1,iZ+1).template getField()
+ partnerLattice1->get(iX,iY+1,iZ-1).template getField()
+ 2.* partnerLattice1->get(iX,iY+1,iZ ).template getField()
+ partnerLattice1->get(iX,iY+1,iZ+1).template getField()
+ partnerLattice1->get(iX-1,iY+1,iZ ).template getField()
+ partnerLattice1->get(iX+1,iY+1,iZ ).template getField() );
T gradMuRhoZ = 1./12. * ( -partnerLattice1->get(iX,iY-1,iZ-1).template getField()
- partnerLattice1->get(iX,iY+1,iZ-1).template getField()
- partnerLattice1->get(iX-1,iY,iZ-1).template getField()
- 2.* partnerLattice1->get(iX,iY,iZ-1).template getField()
- partnerLattice1->get(iX+1,iY,iZ-1).template getField()
+ partnerLattice1->get(iX-1,iY,iZ+1).template getField()
+ 2.* partnerLattice1->get(iX,iY,iZ+1).template getField()
+ partnerLattice1->get(iX+1,iY,iZ+1).template getField()
+ partnerLattice1->get(iX,iY-1,iZ+1).template getField()
+ partnerLattice1->get(iX,iY+1,iZ+1).template getField() );
T psi = 0.;
T gradMuPsiX = 0.;
T gradMuPsiY = 0.;
T gradMuPsiZ = 0.;
if (partners.size() > 1) {
psi = partnerLattice2->get(iX,iY,iZ).computeRho();
gradMuPsiX = 1./12. * ( -partnerLattice2->get(iX-1,iY,iZ-1).template getField()
- partnerLattice2->get(iX-1,iY,iZ+1).template getField()
- partnerLattice2->get(iX-1,iY-1,iZ ).template getField()
- 2.* partnerLattice2->get(iX-1,iY,iZ ).template getField()
- partnerLattice2->get(iX-1,iY+1,iZ ).template getField()
+ partnerLattice2->get(iX+1,iY-1,iZ ).template getField()
+ 2.* partnerLattice2->get(iX+1,iY,iZ ).template getField()
+ partnerLattice2->get(iX+1,iY+1,iZ ).template getField()
+ partnerLattice2->get(iX+1,iY,iZ-1).template getField()
+ partnerLattice2->get(iX+1,iY,iZ+1).template getField() );
gradMuPsiY = 1./12. * ( -partnerLattice2->get(iX-1,iY-1,iZ ).template getField()
- partnerLattice2->get(iX+1,iY-1,iZ ).template getField()
- partnerLattice2->get(iX,iY-1,iZ-1).template getField()
- 2.* partnerLattice2->get(iX,iY-1,iZ ).template getField()
- partnerLattice2->get(iX,iY-1,iZ+1).template getField()
+ partnerLattice2->get(iX,iY+1,iZ-1).template getField()
+ 2.* partnerLattice2->get(iX,iY+1,iZ ).template getField()
+ partnerLattice2->get(iX,iY+1,iZ+1).template getField()
+ partnerLattice2->get(iX-1,iY+1,iZ ).template getField()
+ partnerLattice2->get(iX+1,iY+1,iZ ).template getField() );
gradMuPsiZ = 1./12. * ( -partnerLattice2->get(iX,iY-1,iZ-1).template getField()
- partnerLattice2->get(iX,iY+1,iZ-1).template getField()
- partnerLattice2->get(iX-1,iY,iZ-1).template getField()
- 2.* partnerLattice2->get(iX,iY,iZ-1).template getField()
- partnerLattice2->get(iX+1,iY,iZ-1).template getField()
+ partnerLattice2->get(iX-1,iY,iZ+1).template getField()
+ 2.* partnerLattice2->get(iX,iY,iZ+1).template getField()
+ partnerLattice2->get(iX+1,iY,iZ+1).template getField()
+ partnerLattice2->get(iX,iY-1,iZ+1).template getField()
+ partnerLattice2->get(iX,iY+1,iZ+1).template getField() );
}
T forceX = -rho*gradMuRhoX - phi*gradMuPhiX - psi*gradMuPsiX;
T forceY = -rho*gradMuRhoY - phi*gradMuPhiY - psi*gradMuPsiY;
T forceZ = -rho*gradMuRhoZ - phi*gradMuPhiZ - psi*gradMuPsiZ;
partnerLattice1->get(iX,iY,iZ).template setField({forceX, forceY, forceZ});
T u[3];
partnerLattice1->get(iX,iY,iZ).computeU(u);
blockLattice.get(iX,iY,iZ).template setField(u);
if (partners.size() > 1) {
partnerLattice2->get(iX,iY,iZ).template setField(u);
}
}
}
}
}
}
template
void FreeEnergyForceCoupling3D::process (
BlockLattice3D& blockLattice)
{
processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
}
//////// FreeEnergyInletOutletCoupling3D ///////////////////////////////////
template
FreeEnergyInletOutletCoupling3D ::FreeEnergyInletOutletCoupling3D (
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
std::vector partners_)
: x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), partners(partners_)
{ }
template
FreeEnergyInletOutletCoupling3D ::FreeEnergyInletOutletCoupling3D (
std::vector partners_)
: x0(0), x1(0), y0(0), y1(0), z0(0), z1(0), partners(partners_)
{ }
template
void FreeEnergyInletOutletCoupling3D::processSubDomain (
BlockLattice3D& blockLattice,
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_ )
{
// If partners.size() == 1: two fluid components
// If partners.size() == 2: three fluid components
BlockLattice3D *partnerLattice1 = dynamic_cast *>(partners[0]);
BlockLattice3D *partnerLattice2 = 0;
if (partners.size() > 1) {
partnerLattice2 = dynamic_cast *>(partners[1]);
}
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 u[DESCRIPTOR::d];
partnerLattice1->get(iX,iY,iZ).computeU(u);
blockLattice.get(iX,iY,iZ).defineU(u);
if (partners.size() > 1) {
partnerLattice2->get(iX,iY,iZ).defineU(u);
}
}
}
}
}
}
template
void FreeEnergyInletOutletCoupling3D::process(
BlockLattice3D& blockLattice)
{
processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
}
//////// FreeEnergyDensityOutletCoupling3D ///////////////////////////////////
template
FreeEnergyDensityOutletCoupling3D ::FreeEnergyDensityOutletCoupling3D (
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, T rho_,
std::vector partners_)
: x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_),
rho(rho_), partners(partners_)
{ }
template
FreeEnergyDensityOutletCoupling3D ::FreeEnergyDensityOutletCoupling3D (
T rho_, std::vector partners_)
: x0(0), x1(0), y0(0), y1(0), z0(0), z1(0), rho(rho_), partners(partners_)
{ }
template
void FreeEnergyDensityOutletCoupling3D::processSubDomain (
BlockLattice3D& blockLattice,
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_ )
{
// If partners.size() == 1: two fluid components
// If partners.size() == 2: three fluid components
BlockLattice3D *partnerLattice1 = dynamic_cast *>(partners[0]);
BlockLattice3D *partnerLattice2 = 0;
if (partners.size() > 1) {
partnerLattice2 = dynamic_cast *>(partners[1]);
}
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 rho0, phi, psi;
rho0 = blockLattice.get(iX,iY,iZ).computeRho();
phi = partnerLattice1->get(iX,iY,iZ).computeRho();
blockLattice.get(iX,iY,iZ).defineRho(rho);
partnerLattice1->get(iX,iY,iZ).defineRho(phi * rho / rho0);
if (partners.size() > 1) {
psi = partnerLattice2->get(iX,iY,iZ).computeRho();
partnerLattice2->get(iX,iY,iZ).defineRho(psi * rho / rho0);
}
}
}
}
}
}
template
void FreeEnergyDensityOutletCoupling3D::process(
BlockLattice3D& blockLattice)
{
processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
}
//////// FreeEnergyChemicalPotentialGenerator3D ///////////////////////////////////
template
FreeEnergyChemicalPotentialGenerator3D::FreeEnergyChemicalPotentialGenerator3D (
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
T alpha_, T kappa1_, T kappa2_ )
: LatticeCouplingGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_), alpha(alpha_),
kappa1(kappa1_), kappa2(kappa2_), kappa3(0)
{ }
template
FreeEnergyChemicalPotentialGenerator3D::FreeEnergyChemicalPotentialGenerator3D (
T alpha_, T kappa1_, T kappa2_ )
: LatticeCouplingGenerator3D(0, 0, 0, 0, 0, 0), alpha(alpha_),
kappa1(kappa1_), kappa2(kappa2_), kappa3(0)
{ }
template
FreeEnergyChemicalPotentialGenerator3D::FreeEnergyChemicalPotentialGenerator3D (
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
T alpha_, T kappa1_, T kappa2_, T kappa3_ )
: LatticeCouplingGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_), alpha(alpha_),
kappa1(kappa1_), kappa2(kappa2_), kappa3(kappa3_)
{ }
template
FreeEnergyChemicalPotentialGenerator3D::FreeEnergyChemicalPotentialGenerator3D (
T alpha_, T kappa1_, T kappa2_, T kappa3_ )
: LatticeCouplingGenerator3D(0, 0, 0, 0, 0, 0), alpha(alpha_),
kappa1(kappa1_), kappa2(kappa2_), kappa3(kappa3_)
{ }
template
PostProcessor3D* FreeEnergyChemicalPotentialGenerator3D::generate (
std::vector partners) const
{
return new FreeEnergyChemicalPotentialCoupling3D(
this->x0, this->x1, this->y0, this->y1, this->z0, this->z1,
alpha, kappa1, kappa2, kappa3, partners);
}
template
LatticeCouplingGenerator3D*
FreeEnergyChemicalPotentialGenerator3D::clone() const
{
return new FreeEnergyChemicalPotentialGenerator3D(*this);
}
//////// FreeEnergyForceGenerator3D ///////////////////////////////////
template
FreeEnergyForceGenerator3D::FreeEnergyForceGenerator3D (
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
: LatticeCouplingGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_)
{ }
template
FreeEnergyForceGenerator3D::FreeEnergyForceGenerator3D ( )
: LatticeCouplingGenerator3D(0, 0, 0, 0, 0, 0)
{ }
template
PostProcessor3D* FreeEnergyForceGenerator3D::generate (
std::vector partners) const
{
return new FreeEnergyForceCoupling3D(
this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, partners);
}
template
LatticeCouplingGenerator3D* FreeEnergyForceGenerator3D::clone() const
{
return new FreeEnergyForceGenerator3D(*this);
}
//////// FreeEnergyInletOutletGenerator3D ///////////////////////////////////
template
FreeEnergyInletOutletGenerator3D::FreeEnergyInletOutletGenerator3D (
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
: LatticeCouplingGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_)
{ }
template
FreeEnergyInletOutletGenerator3D::FreeEnergyInletOutletGenerator3D ( )
: LatticeCouplingGenerator3D(0, 0, 0, 0, 0, 0)
{ }
template
PostProcessor3D* FreeEnergyInletOutletGenerator3D::generate (
std::vector partners) const
{
return new FreeEnergyInletOutletCoupling3D(
this->x0, this->x1, this->y0, this->y1, this->z0, this->z1, partners);
}
template
LatticeCouplingGenerator3D* FreeEnergyInletOutletGenerator3D::clone() const
{
return new FreeEnergyInletOutletGenerator3D(*this);
}
//////// FreeEnergyDensityOutletGenerator3D ///////////////////////////////////
template
FreeEnergyDensityOutletGenerator3D::FreeEnergyDensityOutletGenerator3D (
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, T rho_)
: LatticeCouplingGenerator3D(x0_, x1_, y0_, y1_), rho(rho_)
{ }
template
FreeEnergyDensityOutletGenerator3D::FreeEnergyDensityOutletGenerator3D (
T rho_)
: LatticeCouplingGenerator3D(0, 0, 0, 0, 0, 0), rho(rho_)
{ }
template
PostProcessor3D* FreeEnergyDensityOutletGenerator3D::generate (
std::vector partners) const
{
return new FreeEnergyDensityOutletCoupling3D(
this->x0,this->x1,this->y0,this->y1,this->z0,this->z1, rho, partners);
}
template
LatticeCouplingGenerator3D* FreeEnergyDensityOutletGenerator3D::clone() const
{
return new FreeEnergyDensityOutletGenerator3D(*this);
}
} // namespace olb
#endif