/* This file is part of the OpenLB library
*
* Copyright (C) 2008 Orestis Malaspinas, Andrea Parmigiani, 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 NAVIER_STOKES_INTO_ADVECTION_DIFFUSION_COUPLING_POST_PROCESSOR_3D_HH
#define NAVIER_STOKES_INTO_ADVECTION_DIFFUSION_COUPLING_POST_PROCESSOR_3D_HH
#include "latticeDescriptors.h"
#include "navierStokesAdvectionDiffusionCouplingPostProcessor3D.h"
#include "core/blockLattice3D.h"
#include "core/util.h"
#include "core/finiteDifference3D.h"
#include "advectionDiffusionForces.hh"
#include "advectionDiffusionForces.h"
#include "porousAdvectionDiffusionDescriptors.h"
namespace olb {
//=====================================================================================
//============== NavierStokesAdvectionDiffusionCouplingPostProcessor3D ===========
//=====================================================================================
template
NavierStokesAdvectionDiffusionCouplingPostProcessor3D::
NavierStokesAdvectionDiffusionCouplingPostProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
T gravity_, T T0_, T deltaTemp_, std::vector dir_,
std::vector partners_)
: x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_),
gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_),
dir(dir_), partners(partners_)
{
// we normalize the direction of force vector
T normDir = T();
for (unsigned iD = 0; iD < dir.size(); ++iD) {
normDir += dir[iD]*dir[iD];
}
normDir = sqrt(normDir);
for (unsigned iD = 0; iD < dir.size(); ++iD) {
dir[iD] /= normDir;
}
for (unsigned iD = 0; iD < dir.size(); ++iD) {
forcePrefactor[iD] = gravity * dir[iD];
}
tPartner = dynamic_cast> *>(partners[0]);
}
template
void NavierStokesAdvectionDiffusionCouplingPostProcessor3D::
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) {
// computation of the bousinessq force
T *force = blockLattice.get(iX,iY,iZ).template getFieldPointer();
T temperatureDifference = tPartner->get(iX,iY,iZ).computeRho() - T0;
for (unsigned iD = 0; iD < L::d; ++iD) {
force[iD] = forcePrefactor[iD] * temperatureDifference;
}
// Velocity coupling
T *u = tPartner->get(iX,iY,iZ).template getFieldPointer();
blockLattice.get(iX,iY,iZ).computeU(u);
}
}
}
}
}
template
void NavierStokesAdvectionDiffusionCouplingPostProcessor3D::
process(BlockLattice3D& blockLattice)
{
processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
}
//=====================================================================================
//============== SmagorinskyBoussinesqCouplingPostProcessor3D ===============
//=====================================================================================
template
SmagorinskyBoussinesqCouplingPostProcessor3D::
SmagorinskyBoussinesqCouplingPostProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
T gravity_, T T0_, T deltaTemp_, std::vector dir_, T PrTurb_,
std::vector partners_)
: x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_),
gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_),
dir(dir_), PrTurb(PrTurb_), partners(partners_)
{
// we normalize the direction of force vector
T normDir = T();
for (unsigned iD = 0; iD < dir.size(); ++iD) {
normDir += dir[iD]*dir[iD];
}
normDir = sqrt(normDir);
for (unsigned iD = 0; iD < dir.size(); ++iD) {
dir[iD] /= normDir;
}
for (unsigned iD = 0; iD < dir.size(); ++iD) {
forcePrefactor[iD] = gravity * dir[iD];
}
tauTurbADPrefactor = descriptors::invCs2>() / descriptors::invCs2() / PrTurb;
tPartner = dynamic_cast> *>(partners[0]);
}
template
void SmagorinskyBoussinesqCouplingPostProcessor3D::
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) {
// computation of the bousinessq force
T *force = blockLattice.get(iX,iY,iZ).template getFieldPointer();
T temperatureDifference = tPartner->get(iX,iY,iZ).computeRho() - T0;
for (unsigned iD = 0; iD < L::d; ++iD) {
force[iD] = forcePrefactor[iD] * temperatureDifference;
}
// Velocity coupling
T *u = tPartner->get(iX,iY,iZ).template getFieldPointer();
// tau coupling
T *tauNS = blockLattice.get(iX,iY,iZ).template getFieldPointer();
T *tauAD = tPartner->get(iX,iY,iZ).template getFieldPointer();
T rho, pi[util::TensorVal::n];
blockLattice.get(iX,iY,iZ).computeAllMomenta(rho, u, pi);
T PiNeqNormSqr = pi[0]*pi[0] + 2.0*pi[1]*pi[1] + pi[2]*pi[2];
if (util::TensorVal::n == 6) {
PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2*pi[4]*pi[4] +pi[5]*pi[5];
}
T PiNeqNorm = sqrt(PiNeqNormSqr);
/// Molecular realaxation time
T tau_mol_NS = 1. / blockLattice.get(iX,iY,iZ).getDynamics()->getOmega();
T tau_mol_AD = 1. / tPartner->get(iX,iY,iZ).getDynamics()->getOmega();
/// Turbulent realaxation time
T tau_turb_NS = 0.5*(sqrt(tau_mol_NS*tau_mol_NS + dynamic_cast*>(blockLattice.get(iX,iY,iZ).getDynamics())->getPreFactor()/rho*PiNeqNorm) - tau_mol_NS);
/// Effective realaxation time
tauNS[0] = tau_mol_NS+tau_turb_NS;
T tau_turb_AD = tau_turb_NS * tauTurbADPrefactor;
tauAD[0] = tau_mol_AD+tau_turb_AD;
}
}
}
}
}
template
void SmagorinskyBoussinesqCouplingPostProcessor3D::
process(BlockLattice3D& blockLattice)
{
processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
}
/// LatticeCouplingGenerator for advectionDiffusion coupling
template
SmagorinskyBoussinesqCouplingGenerator3D::
SmagorinskyBoussinesqCouplingGenerator3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
T gravity_, T T0_, T deltaTemp_, std::vector dir_, T PrTurb_)
: LatticeCouplingGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_),
gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_), dir(dir_), PrTurb(PrTurb_)
{ }
template
PostProcessor3D* SmagorinskyBoussinesqCouplingGenerator3D::generate (
std::vector partners) const
{
return new SmagorinskyBoussinesqCouplingPostProcessor3D(
this->x0,this->x1,this->y0,this->y1,this->z0,this->z1, gravity, T0, deltaTemp, dir, PrTurb, partners);
}
template
LatticeCouplingGenerator3D* SmagorinskyBoussinesqCouplingGenerator3D::clone() const
{
return new SmagorinskyBoussinesqCouplingGenerator3D(*this);
}
//=====================================================================================
//============== AdvectionDiffusionParticleCouplingPostProcessor3D ===========
//=====================================================================================
template
AdvectionDiffusionParticleCouplingPostProcessor3D::
AdvectionDiffusionParticleCouplingPostProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int iC_,
int offset_, std::vector partners_,
std::vector > > forces_)
: x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), iC(iC_), offset(offset_), forces(forces_)
{
BlockLattice3D *partnerLattice =
dynamic_cast *>(partners_[0]);
adCell = &partnerLattice->get(x0,y0,z0);
vel = &partnerLattice->get(x0,y0,z0)[offset];
vel_new = &partnerLattice->get(x0,y0,z0)[offset];
velXp = &partnerLattice->get(x0+1,y0,z0)[offset];
velXn = &partnerLattice->get(x0-1,y0,z0)[offset];
velYp = &partnerLattice->get(x0,y0+1,z0)[offset];
velYn = &partnerLattice->get(x0,y0-1,z0)[offset];
velZp = &partnerLattice->get(x0,y0,z0+1)[offset];
velZn = &partnerLattice->get(x0,y0,z0-1)[offset];
}
template
void AdvectionDiffusionParticleCouplingPostProcessor3D::
processSubDomain(BlockLattice3D& blockLattice,
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
{
int newX0, newX1, newY0, newY1, newZ0, newZ1;
int off = (par) ? 3 : 0;
int off2 = (par) ? 0 : 3;
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) {
int latticeR[4] = {iC, iX, iY, iZ};
T velGrad[3] = {0.,0.,0.};
T forceValue[3] = {0.,0.,0.};
T velF[3] = {0.,0.,0.};
nsCell = &(blockLattice.get(iX,iY,iZ));
if (forces.begin() != forces.end()) {
// calculating upwind Gradient
// vel contains velocity information on ADlattice
// velGrad contains upwind vel on ADlattice
if (vel[0+off]<0.) {
velGrad[0] = vel[0+off]*(velXp[0+off]-vel[0+off]);
velGrad[1] = vel[0+off]*(velXp[1+off]-vel[1+off]);
velGrad[2] = vel[0+off]*(velXp[2+off]-vel[2+off]);
} else {
velGrad[0] = vel[0+off]*(vel[0+off]-velXn[0+off]);
velGrad[1] = vel[0+off]*(vel[1+off]-velXn[1+off]);
velGrad[2] = vel[0+off]*(vel[2+off]-velXn[2+off]);
}
if (vel[1+off]<0.) {
velGrad[0] += vel[1+off]*(velYp[0+off]-vel[0+off]);
velGrad[1] += vel[1+off]*(velYp[1+off]-vel[1+off]);
velGrad[2] += vel[1+off]*(velYp[2+off]-vel[2+off]);
} else {
velGrad[0] += vel[1+off]*(vel[0+off]-velYn[0+off]);
velGrad[1] += vel[1+off]*(vel[1+off]-velYn[1+off]);
velGrad[2] += vel[1+off]*(vel[2+off]-velYn[2+off]);
}
if (vel[2+off]<0.) {
velGrad[0] += vel[2+off]*(velZp[0+off]-vel[0+off]);
velGrad[1] += vel[2+off]*(velZp[1+off]-vel[1+off]);
velGrad[2] += vel[2+off]*(velZp[2+off]-vel[2+off]);
} else {
velGrad[0] += vel[2+off]*(vel[0+off]-velZn[0+off]);
velGrad[1] += vel[2+off]*(vel[1+off]-velZn[1+off]);
velGrad[2] += vel[2+off]*(vel[2+off]-velZn[2+off]);
}
for (AdvectionDiffusionForce3D& f : forces) {
// writes force in forceValues, vel refers to ADlattice
f.applyForce(forceValue, nsCell, adCell, &vel[off], latticeR);
}
// compute new particle velocity
for (int i=0; i < DESCRIPTOR::d; i++) {
vel_new[i+off2] = vel[i+off] + forceValue[i] - velGrad[i];
}
} else { // set particle velocity to fluid velocity
nsCell->computeU(velF);
for (int i = 0; i < DESCRIPTOR::d; i++) {
vel_new[i + off2] = velF[i];
}
}
}
}
}
}
par = !par;
}
template
void AdvectionDiffusionParticleCouplingPostProcessor3D::
process(BlockLattice3D& blockLattice)
{
processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
}
// LatticeCouplingGenerator for advectionDiffusion coupling
template
NavierStokesAdvectionDiffusionCouplingGenerator3D::
NavierStokesAdvectionDiffusionCouplingGenerator3D(int x0_, int x1_, int y0_, int y1_,int z0_, int z1_, T gravity_, T T0_, T deltaTemp_, std::vector dir_)
: LatticeCouplingGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_),
gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_), dir(dir_)
{ }
template
PostProcessor3D* NavierStokesAdvectionDiffusionCouplingGenerator3D::generate (
std::vector partners) const
{
return new NavierStokesAdvectionDiffusionCouplingPostProcessor3D(
this->x0,this->x1,this->y0,this->y1,this->z0,this->z1, gravity, T0, deltaTemp, dir,partners);
}
template
LatticeCouplingGenerator3D* NavierStokesAdvectionDiffusionCouplingGenerator3D::clone() const
{
return new NavierStokesAdvectionDiffusionCouplingGenerator3D(*this);
}
// LatticeCouplingGenerator for one-way advectionDiffusion coupling with Stokes drag
template
AdvectionDiffusionParticleCouplingGenerator3D::
AdvectionDiffusionParticleCouplingGenerator3D(int offset_)
: LatticeCouplingGenerator3D(0, 0, 0, 0, 0, 0), offset(offset_)
{ }
template
PostProcessor3D* AdvectionDiffusionParticleCouplingGenerator3D::generate (
std::vector partners) const
{
return new AdvectionDiffusionParticleCouplingPostProcessor3D(this->x0,this->x1,this->y0,this->y1,this->z0,this->z1, this->iC, offset, partners, ADforces);
}
template
LatticeCouplingGenerator3D* AdvectionDiffusionParticleCouplingGenerator3D::clone() const
{
return new AdvectionDiffusionParticleCouplingGenerator3D(*this);
}
template
void AdvectionDiffusionParticleCouplingGenerator3D::addForce(
AdvectionDiffusionForce3D &force)
{
ADforces.push_back(force);
}
//=====================================================================================
//============== PorousNavierStokesAdvectionDiffusionCouplingPostProcessor3D ===========
//=====================================================================================
template
PorousNavierStokesAdvectionDiffusionCouplingPostProcessor3D::
PorousNavierStokesAdvectionDiffusionCouplingPostProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
T gravity_, T T0_, T deltaTemp_, std::vector dir_,
std::vector partners_)
: x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_),
gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_),
dir(dir_), partners(partners_)
{
// we normalize the direction of force vector
T normDir = T();
for (unsigned iD = 0; iD < dir.size(); ++iD) {
normDir += dir[iD]*dir[iD];
}
normDir = sqrt(normDir);
for (unsigned iD = 0; iD < dir.size(); ++iD) {
dir[iD] /= normDir;
}
}
template
void PorousNavierStokesAdvectionDiffusionCouplingPostProcessor3D::
processSubDomain(BlockLattice3D& blockLattice,
int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
{
typedef DESCRIPTOR L;
enum {x,y,z};
BlockLattice3D *tPartner =
dynamic_cast *>(partners[0]);
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) {
// Velocity coupling
T *u = tPartner->get(iX,iY,iZ).template getFieldPointer();
blockLattice.get(iX,iY,iZ).computeU(u);
//coupling between the temperature and navier stokes.
T *force = blockLattice.get(iX,iY,iZ).template getFieldPointer();
// this should return the interpolated solid-fluid temperature
T temperature = tPartner->get(iX,iY,iZ).computeRho();
T rho = blockLattice.get(iX,iY,iZ).computeRho();
for (unsigned iD = 0; iD < L::d; ++iD) {
force[iD] = gravity * rho * (temperature - T0) / deltaTemp * dir[iD];
}
}
}
}
}
}
template
void PorousNavierStokesAdvectionDiffusionCouplingPostProcessor3D::
process(BlockLattice3D& blockLattice)
{
processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
}
template
PorousNavierStokesAdvectionDiffusionCouplingGenerator3D::
PorousNavierStokesAdvectionDiffusionCouplingGenerator3D(int x0_, int x1_, int y0_, int y1_,int z0_, int z1_,
T gravity_, T T0_, T deltaTemp_, std::vector dir_)
: LatticeCouplingGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_),
gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_), dir(dir_)
{ }
template
PostProcessor3D* PorousNavierStokesAdvectionDiffusionCouplingGenerator3D::generate (
std::vector partners) const
{
return new PorousNavierStokesAdvectionDiffusionCouplingPostProcessor3D(
this->x0,this->x1,this->y0,this->y1,this->z0,this->z1, gravity, T0, deltaTemp, dir,partners);
}
template
LatticeCouplingGenerator3D* PorousNavierStokesAdvectionDiffusionCouplingGenerator3D::clone() const
{
return new PorousNavierStokesAdvectionDiffusionCouplingGenerator3D(*this);
}
} // namespace olb
#endif