/* 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.
*/
/** \file
* A collection of dynamics classes (e.g. BGK) with which a Cell object
* can be instantiated -- generic implementation.
*/
#ifndef ADVECTION_DIFFUSION_DYNAMICS_HH
#define ADVECTION_DIFFUSION_DYNAMICS_HH
#include
#include
#include "advectionDiffusionDynamics.h"
#include "dynamics/lbHelpers.h"
namespace olb {
////////////////////// Class AdvectionDiffusionRLBdynamics //////////////////////////
/** \param omega_ relaxation parameter, related to the dynamic viscosity
* \param momenta_ a Momenta object to know how to compute velocity momenta
*/
//==================================================================//
//============= Regularized Model for Advection diffusion===========//
//==================================================================//
template
AdvectionDiffusionRLBdynamics::AdvectionDiffusionRLBdynamics (
T omega_, Momenta& momenta_ )
: BasicDynamics( momenta_ ),
omega( omega_ )
{ }
template
T AdvectionDiffusionRLBdynamics::computeEquilibrium( int iPop, T rho,
const T u[DESCRIPTOR::d], T uSqr ) const
{
return lbHelpers::equilibriumFirstOrder( iPop, rho, u );
}
template
void AdvectionDiffusionRLBdynamics::collide( Cell& cell,
LatticeStatistics& statistics )
{
T temperature = this->_momenta.computeRho( cell );
const T* u = cell.template getFieldPointer();
T uSqr = lbHelpers::
rlbCollision( cell, temperature, u, omega );
statistics.incrementStats( temperature, uSqr );
}
template
T AdvectionDiffusionRLBdynamics::getOmega() const
{
return omega;
}
template
void AdvectionDiffusionRLBdynamics::setOmega( T omega_ )
{
omega = omega_;
}
////////////////////// Class CombinedAdvectionDiffusionRLBdynamics /////////////////////////
template
CombinedAdvectionDiffusionRLBdynamics::CombinedAdvectionDiffusionRLBdynamics (
T omega, Momenta& momenta )
: BasicDynamics(momenta),
_boundaryDynamics(omega, momenta)
{ }
template
T CombinedAdvectionDiffusionRLBdynamics::
computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const
{
return lbHelpers::equilibriumFirstOrder( iPop, rho, u );
}
template
void CombinedAdvectionDiffusionRLBdynamics::collide (
Cell& cell,
LatticeStatistics& statistics )
{
typedef DESCRIPTOR L;
T temperature = this->_momenta.computeRho( cell );
const T* u = cell.template getFieldPointer();
T jNeq[DESCRIPTOR::d];
// this->_momenta.computeJ( cell, jNeq );
dynamic_cast&>(this->_momenta).computeJneq( cell, jNeq );
// cout << jNeq[0] << " " << jNeq[1] << " " << u[0] << " " << u[1] << endl;
// stripe of equilibrium part u * T
// for ( int iD = 0; iD < DESCRIPTOR::d; ++iD ) {
// jNeq[iD] -= u[iD] * temperature;
// }
// cout << jNeq[0] << " " << jNeq[1] << " " << u[0] << " " << u[1] << endl;
for (int iPop = 0; iPop < L::q; ++iPop) {
cell[iPop] = lbHelpers::equilibriumFirstOrder( iPop, temperature, u )
+ firstOrderLbHelpers::fromJneqToFneq(iPop, jNeq);
// cout << firstOrderLbHelpers::fromJneqToFneq(iPop,jNeq) << " ";
}
// cout << endl;
// lbHelpers::computeJ(cell, jNeq);
// cout << jNeq[0] << " " << jNeq[1] << " " << u[0] << " " << u[1] << endl;
_boundaryDynamics.collide(cell, statistics);
// lbHelpers::computeJ(cell, jNeq);
// cout << jNeq[0] << " " << jNeq[1] << " " << u[0] << " " << u[1] << endl;
}
template
T CombinedAdvectionDiffusionRLBdynamics::getOmega() const
{
return _boundaryDynamics.getOmega();
}
template
void CombinedAdvectionDiffusionRLBdynamics::setOmega(T omega)
{
_boundaryDynamics.setOmega(omega);
}
//==================================================================//
//============= BGK Model for Advection diffusion===========//
//==================================================================//
template
AdvectionDiffusionBGKdynamics::AdvectionDiffusionBGKdynamics (
T omega, Momenta& momenta )
: BasicDynamics( momenta ),
_omega(omega)
{ }
template
AdvectionDiffusionBGKdynamics::AdvectionDiffusionBGKdynamics (
const UnitConverter& converter, Momenta& momenta )
: BasicDynamics( momenta ),
_omega(converter.getLatticeRelaxationFrequency())
{ }
template
T AdvectionDiffusionBGKdynamics::computeEquilibrium( int iPop, T rho,
const T u[DESCRIPTOR::d], T uSqr ) const
{
return lbHelpers::equilibriumFirstOrder( iPop, rho, u );
}
template
void AdvectionDiffusionBGKdynamics::collide( Cell& cell,
LatticeStatistics& statistics )
{
T temperature = this->_momenta.computeRho( cell );
const T* u = cell.template getFieldPointer();
T uSqr = lbHelpers::
bgkCollision( cell, temperature, u, _omega );
statistics.incrementStats( temperature, uSqr );
}
template
T AdvectionDiffusionBGKdynamics::getOmega() const
{
return _omega;
}
template
void AdvectionDiffusionBGKdynamics::setOmega( T omega )
{
_omega = omega;
}
//==================================================================//
//============= TRT Model for Advection diffusion===========//
//==================================================================//
template
AdvectionDiffusionTRTdynamics::AdvectionDiffusionTRTdynamics (
T omega, Momenta& momenta, T magicParameter )
: AdvectionDiffusionBGKdynamics( omega, momenta ),
_omega2(1/(magicParameter/(1/omega-0.5)+0.5)), _magicParameter(magicParameter)
{ }
template
void AdvectionDiffusionTRTdynamics::collide( Cell& cell,
LatticeStatistics& statistics )
{
T temperature = this->_momenta.computeRho( cell );
const T* u = cell.template getFieldPointer();
T fPlus[DESCRIPTOR::q], fMinus[DESCRIPTOR::q];
T fEq[DESCRIPTOR::q], fEqPlus[DESCRIPTOR::q], fEqMinus[DESCRIPTOR::q];
for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
fPlus[iPop] = 0.5 * ( cell[iPop] + cell[descriptors::opposite(iPop)] );
fMinus[iPop] = 0.5 * ( cell[iPop] - cell[descriptors::opposite(iPop)] );
fEq[iPop] = lbHelpers::equilibriumFirstOrder(iPop, temperature, u);
}
for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
fEqPlus[iPop] = 0.5 * ( fEq[iPop] + fEq[descriptors::opposite(iPop)] );
fEqMinus[iPop] = 0.5 * ( fEq[iPop] - fEq[descriptors::opposite(iPop)] );
}
for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
cell[iPop] -= _omega2 * (fPlus[iPop] - fEqPlus[iPop]) + this->_omega * (fMinus[iPop] - fEqMinus[iPop]);
}
statistics.incrementStats( temperature, 0. );
}
//==================================================================//
//============= BGK Model for Advection diffusion with source ===========//
//==================================================================//
template
SourcedAdvectionDiffusionBGKdynamics::SourcedAdvectionDiffusionBGKdynamics (
T omega, Momenta& momenta )
: AdvectionDiffusionBGKdynamics( omega, momenta ), _omegaMod(1. - 0.5 * omega)
{
OLB_PRECONDITION( DESCRIPTOR::template provides() );
}
template
void SourcedAdvectionDiffusionBGKdynamics::collide( Cell& cell,
LatticeStatistics& statistics )
{
const T* u = cell.template getFieldPointer();
const T* source = cell.template getFieldPointer();
const T temperature = this->_momenta.computeRho( cell ) + 0.5 * source[0];
T uSqr = lbHelpers::
bgkCollision( cell, temperature, u, this->_omega );
// Q / t_i = ( 1 - 0.5*_omega) * q
const T sourceMod = source[0] * _omegaMod;
for ( int iPop = 0; iPop < DESCRIPTOR::q; iPop++ ) {
cell[iPop] += sourceMod * descriptors::t(iPop);
}
statistics.incrementStats( temperature, uSqr );
}
template
T SourcedAdvectionDiffusionBGKdynamics::computeRho(Cell const& cell) const
{
const T* source = cell.template getFieldPointer();
return this->_momenta.computeRho( cell ) + 0.5 * source[0];
}
template
void SourcedAdvectionDiffusionBGKdynamics::computeRhoU (
Cell const& cell,
T& rho, T u[DESCRIPTOR::d]) const
{
this->_momenta.computeRhoU( cell, rho, u );
const T* source = cell.template getFieldPointer();
rho += 0.5 * source[0] * _omegaMod;
}
//==================================================================================//
//=========== BGK Model for Advection diffusion with Stokes drag and Smagorinsky====//
//==================================================================================//
template
SmagorinskyParticleAdvectionDiffusionBGKdynamics::SmagorinskyParticleAdvectionDiffusionBGKdynamics (
T omega_, Momenta& momenta_, T smagoConst_, T dx_, T dt_)
: AdvectionDiffusionBGKdynamics(omega_,momenta_), smagoConst(smagoConst_), preFactor(computePreFactor(omega_,smagoConst_, dx_, dt_) )
{ }
template
void SmagorinskyParticleAdvectionDiffusionBGKdynamics::collide(Cell& cell, LatticeStatistics& statistics )
{
T temperature, uad[DESCRIPTOR::d], pi[util::TensorVal::n];
this->_momenta.computeAllMomenta(cell, temperature, uad, pi);
const T* u = (statistics.getTime() % 2 == 0) ? cell.template getFieldPointer() : cell.template getFieldPointer();
T newOmega = computeOmega(this->getOmega(), preFactor, temperature, pi);
T uSqr = lbHelpers::bgkCollision(cell, temperature, u, newOmega);
statistics.incrementStats(temperature, uSqr);
}
template
T SmagorinskyParticleAdvectionDiffusionBGKdynamics::getSmagorinskyOmega(Cell& cell)
{
T temperature, uTemp[DESCRIPTOR::d], pi[util::TensorVal::n];
this->_momenta.computeAllMomenta(cell, temperature, uTemp, pi);
T newOmega = computeOmega(this->getOmega(), preFactor, temperature, pi);
return newOmega;
}
template
void SmagorinskyParticleAdvectionDiffusionBGKdynamics::setOmega(T omega_)
{
preFactor = computePreFactor(omega_, smagoConst, dx, dt);
}
template
T SmagorinskyParticleAdvectionDiffusionBGKdynamics::computePreFactor(T omega_, T smagoConst_, T dx_, T dt_)
{
return (T)(smagoConst_*smagoConst_*dx_*dx_)*descriptors::invCs2()/dt_*4*sqrt(2);
}
template
T SmagorinskyParticleAdvectionDiffusionBGKdynamics::computeOmega(T omega0, T preFactor_, T rho, T pi[util::TensorVal::n] )
{
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 = 1. /omega0;
/// Turbulent realaxation time
T tau_turb = 0.5*(sqrt(tau_mol*tau_mol+(preFactor_*tau_eff*PiNeqNorm))-tau_mol);
/// Effective realaxation time
tau_eff = tau_mol+tau_turb;
T omega_new= 1./tau_eff;
return omega_new;
}
//==================================================================//
//=========== BGK Model for Advection diffusion with Stokes Drag ====//
//==================================================================//
template
ParticleAdvectionDiffusionBGKdynamics::ParticleAdvectionDiffusionBGKdynamics (
T omega_, Momenta& momenta_ )
: AdvectionDiffusionBGKdynamics(omega_,momenta_), omega( omega_ )
{ }
template
void ParticleAdvectionDiffusionBGKdynamics::collide( Cell& cell,
LatticeStatistics& statistics )
{
T temperature = this->_momenta.computeRho( cell );
const T* u = (statistics.getTime() % 2 == 0) ? cell.template getFieldPointer() : cell.template getFieldPointer();
T uSqr = lbHelpers::
bgkCollision( cell, temperature, u, omega );
statistics.incrementStats( temperature, uSqr );
}
//==================================================================//
//================= MRT Model for Advection diffusion ==============//
//==================================================================//
template
AdvectionDiffusionMRTdynamics::AdvectionDiffusionMRTdynamics(
T omega, Momenta& momenta) :
BasicDynamics(momenta), _omega(omega)
{
T rt[DESCRIPTOR::q]; // relaxation times vector.
for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
rt[iPop] = descriptors::s(iPop);
}
for (int iPop = 0; iPop < descriptors::shearIndexes(); ++iPop) {
rt[descriptors::shearViscIndexes(iPop)] = omega;
}
for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
for (int jPop = 0; jPop < DESCRIPTOR::q; ++jPop) {
invM_S[iPop][jPop] = T();
for (int kPop = 0; kPop < DESCRIPTOR::q; ++kPop) {
if (kPop == jPop) {
invM_S[iPop][jPop] += descriptors::invM(iPop,kPop) * rt[kPop];
}
}
}
}
}
template
T AdvectionDiffusionMRTdynamics::computeEquilibrium(int iPop, T rho,
const T u[DESCRIPTOR::d], T uSqr) const
{
return lbHelpers::equilibrium(iPop, rho, u, uSqr);
}
template
void AdvectionDiffusionMRTdynamics::collide(Cell& cell,
LatticeStatistics& statistics)
{
T temperature = this->_momenta.computeRho(cell);
const T* u = cell.template getFieldPointer();
T uSqr = lbHelpers::mrtCollision(cell, temperature, u, invM_S);
statistics.incrementStats(temperature, uSqr);
}
template
T AdvectionDiffusionMRTdynamics::getOmega() const
{
return _omega;
}
template
void AdvectionDiffusionMRTdynamics::setOmega(T omega)
{
_omega = omega;
}
} // namespace olb
#endif