/* This file is part of the OpenLB library
*
* Copyright (C) 2008 Orestis Malaspinas, Andrea Parmigiani
* 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 ADVECTION_DIFFUSION_BOUNDARIES_HH
#define ADVECTION_DIFFUSION_BOUNDARIES_HH
#include "advectionDiffusionBoundaries.h"
#include "dynamics/latticeDescriptors.h"
#include "core/util.h"
#include "dynamics/lbHelpers.h"
namespace olb {
//==================================================================================================
//==================== For regularized Advection Diffusion Boundary Condition ======================
//============================================================================================
// For flat Walls
template
AdvectionDiffusionBoundariesDynamics::
AdvectionDiffusionBoundariesDynamics( T omega_, Momenta& momenta_)
: BasicDynamics(momenta_), boundaryDynamics(omega_, momenta_)
{
}
template
T AdvectionDiffusionBoundariesDynamics::
computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const
{
return lbHelpers::equilibriumFirstOrder( iPop, rho, u );
}
template
void AdvectionDiffusionBoundariesDynamics::
collide(Cell& cell,LatticeStatistics& statistics)
{
typedef DESCRIPTOR L;
typedef lbHelpers lbH;
T dirichletTemperature = this->_momenta.computeRho(cell);
T* u = cell.template getFieldPointer();
std::vector unknownIndexes = util::subIndexOutgoing();
std::vector knownIndexes = util::remainingIndexes(unknownIndexes);
int missingNormal = 0;
if ((L::d == 3 && L::q == 7)||(L::d == 2 && L::q == 5)) {
T sum = T();
for (unsigned i = 0; i < knownIndexes.size(); ++i) {
sum += cell[knownIndexes[i]];
}
T difference = dirichletTemperature - (T) 1 - sum; // on cell there are non-shiftet values -> temperature has to be changed
// here I know all missing and non missing f_i
for (unsigned i = 0; i < unknownIndexes.size(); ++i) {
int numOfNonNullComp = 0;
for (int iDim = 0; iDim < L::d; ++iDim) {
numOfNonNullComp += abs(descriptors::c(unknownIndexes[i],iDim));
}
if (numOfNonNullComp == 1) {
missingNormal = unknownIndexes[i];
// here missing diagonal directions are erased
// just the normal direction stays (D3Q7)
unknownIndexes.erase(unknownIndexes.begin() + i);
break;
}
}
cell[missingNormal] = difference; // on cell there are non-shiftet values -> temperature has to be changed
boundaryDynamics.collide(cell, statistics); // only for D3Q7
} else {
// part for q=19 copied from AdvectionDiffusionEdgesDynamics.collide()
// but here just all missing directions, even at border of inlet area
// has to be checked!
for (unsigned iteratePop = 0; iteratePop < unknownIndexes.size();
++iteratePop) {
cell[unknownIndexes[iteratePop]] =
lbH::equilibriumFirstOrder(unknownIndexes[iteratePop], dirichletTemperature, u)
- (cell[util::opposite(unknownIndexes[iteratePop])]
- lbH::equilibriumFirstOrder(
util::opposite(unknownIndexes[iteratePop]),
dirichletTemperature, u));
}
}
}
template
T AdvectionDiffusionBoundariesDynamics::
getOmega() const
{
return boundaryDynamics.getOmega();
}
template
void AdvectionDiffusionBoundariesDynamics::
setOmega(T omega_)
{
boundaryDynamics.setOmega(omega_);
}
//=================================================================
// For 2D Corners with regularized Dynamic ==============================================
//=================================================================
template
AdvectionDiffusionCornerDynamics2D::AdvectionDiffusionCornerDynamics2D(
T omega_, Momenta& momenta_)
: BasicDynamics(momenta_),
boundaryDynamics(omega_, momenta_)
{
}
template
T AdvectionDiffusionCornerDynamics2D::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const
{
return lbHelpers::equilibriumFirstOrder( iPop, rho, u );
}
template
void AdvectionDiffusionCornerDynamics2D::collide(Cell& cell,LatticeStatistics& statistics)
{
typedef DESCRIPTOR L;
typedef lbHelpers lbH;
T temperature = this->_momenta.computeRho(cell);
T* u = cell.template getFieldPointer();
// I need to get Missing information on the corners !!!!
std::vector unknownIndexes = util::subIndexOutgoing2DonCorners();
// here I know all missing and non missing f_i
// The collision procedure for D2Q5 and D3Q7 lattice is the same ...
// Given the rule f_i_neq = -f_opposite(i)_neq
// I have the right number of equations for the number of unknowns using these lattices
for (unsigned iPop = 0; iPop < unknownIndexes.size(); ++iPop) {
cell[unknownIndexes[iPop]] = lbH::equilibriumFirstOrder(unknownIndexes[iPop], temperature, u)
-(cell[util::opposite(unknownIndexes[iPop])]
- lbH::equilibriumFirstOrder(util::opposite(unknownIndexes[iPop]), temperature, u) ) ;
}
// Once all the f_i are known, I can call the collision for the Regularized Model.
boundaryDynamics.collide(cell, statistics);
}
template
T AdvectionDiffusionCornerDynamics2D::getOmega() const
{
return boundaryDynamics.getOmega();
}
template
void AdvectionDiffusionCornerDynamics2D::setOmega(T omega_)
{
boundaryDynamics.setOmega(omega_);
}
//=================================================================
// For 3D Corners with regularized Dynamic ==============================================
//=================================================================
template
AdvectionDiffusionCornerDynamics3D::AdvectionDiffusionCornerDynamics3D(
T omega_, Momenta& momenta_)
: BasicDynamics(momenta_),
boundaryDynamics(omega_, momenta_)
{
}
template
T AdvectionDiffusionCornerDynamics3D::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const
{
return lbHelpers::equilibriumFirstOrder( iPop, rho, u );
}
template
void AdvectionDiffusionCornerDynamics3D::collide(Cell& cell,LatticeStatistics& statistics)
{
typedef DESCRIPTOR L;
typedef lbHelpers lbH;
T temperature = this->_momenta.computeRho(cell);
T* u = cell.template getFieldPointer();
// I need to get Missing information on the corners !!!!
std::vector unknownIndexes = util::subIndexOutgoing3DonCorners();
// here I know all missing and non missing f_i
// The collision procedure for D2Q5 and D3Q7 lattice is the same ...
// Given the rule f_i_neq = -f_opposite(i)_neq
// I have the right number of equations for the number of unknowns using these lattices
for (unsigned iPop = 0; iPop < unknownIndexes.size(); ++iPop) {
cell[unknownIndexes[iPop]] = lbH::equilibriumFirstOrder(unknownIndexes[iPop], temperature, u)
-(cell[util::opposite(unknownIndexes[iPop])]
- lbH::equilibriumFirstOrder(util::opposite(unknownIndexes[iPop]), temperature, u) ) ;
}
// Once all the f_i are known, I can call the collision for the Regularized Model.
boundaryDynamics.collide(cell, statistics);
}
template
T AdvectionDiffusionCornerDynamics3D::getOmega() const
{
return boundaryDynamics.getOmega();
}
template
void AdvectionDiffusionCornerDynamics3D::setOmega(T omega_)
{
boundaryDynamics.setOmega(omega_);
}
//=================================================================
// For 3D Edges with regularized Dynamic ==============================================
//=================================================================
template
AdvectionDiffusionEdgesDynamics::AdvectionDiffusionEdgesDynamics(
T omega_, Momenta& momenta_)
: BasicDynamics(momenta_),
boundaryDynamics(omega_, momenta_)
{
}
template
T AdvectionDiffusionEdgesDynamics::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const
{
return lbHelpers::equilibriumFirstOrder( iPop, rho, u );
}
template
void AdvectionDiffusionEdgesDynamics::collide(Cell& cell,LatticeStatistics& statistics)
{
typedef DESCRIPTOR L;
typedef lbHelpers lbH;
T temperature = this->_momenta.computeRho(cell);
T* u = cell.template getFieldPointer();
// I need to get Missing information on the corners !!!!
std::vector unknownIndexes = util::subIndexOutgoing3DonEdges();
// here I know all missing and non missing f_i
// The collision procedure for D2Q5 and D3Q7 lattice is the same ...
// Given the rule f_i_neq = -f_opposite(i)_neq
// I have the right number of equations for the number of unknowns using these lattices
for (unsigned iPop = 0; iPop < unknownIndexes.size(); ++iPop) {
cell[unknownIndexes[iPop]] = lbH::equilibriumFirstOrder(unknownIndexes[iPop], temperature, u)
-(cell[util::opposite(unknownIndexes[iPop])]
- lbH::equilibriumFirstOrder(util::opposite(unknownIndexes[iPop]), temperature, u) ) ;
}
// Once all the f_i are known, I can call the collision for the Regularized Model.
boundaryDynamics.collide(cell, statistics);
}
template
T AdvectionDiffusionEdgesDynamics::getOmega() const
{
return boundaryDynamics.getOmega();
}
template
void AdvectionDiffusionEdgesDynamics::setOmega(T omega_)
{
boundaryDynamics.setOmega(omega_);
}
} // namespace olb
#endif