diff options
Diffstat (limited to 'src/dynamics/dynamics.hh')
| -rw-r--r-- | src/dynamics/dynamics.hh | 2026 |
1 files changed, 2026 insertions, 0 deletions
diff --git a/src/dynamics/dynamics.hh b/src/dynamics/dynamics.hh new file mode 100644 index 0000000..12b77bf --- /dev/null +++ b/src/dynamics/dynamics.hh @@ -0,0 +1,2026 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2006-2015 Jonas Latt, Mathias J. Krause + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * <http://www.openlb.net/> + * + * 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 LB_DYNAMICS_HH +#define LB_DYNAMICS_HH + +#include <type_traits> +#include <numeric> +#include "dynamics.h" +#include "core/cell.h" +#include "lbHelpers.h" +#include "firstOrderLbHelpers.h" +#include "d3q13Helpers.h" + +namespace olb { + +////////////////////// Class Dynamics //////////////////////// + +template<typename T, typename DESCRIPTOR> +T Dynamics<T,DESCRIPTOR>::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const +{ + return lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr); +} + +template<typename T, typename DESCRIPTOR> +void Dynamics<T,DESCRIPTOR>::iniEquilibrium(Cell<T,DESCRIPTOR>& cell, T rho, const T u[DESCRIPTOR::d]) +{ + T uSqr = util::normSqr<T,DESCRIPTOR::d>(u); + for (int iPop=0; iPop<DESCRIPTOR::q; ++iPop) { + cell[iPop] = computeEquilibrium(iPop, rho, u, uSqr); + } +} + +template<typename T, typename DESCRIPTOR> +void Dynamics<T,DESCRIPTOR>::setBoundaryIntersection(int iPop, T distance) +{ } + +template<typename T, typename DESCRIPTOR> +bool Dynamics<T,DESCRIPTOR>::getBoundaryIntersection(int iPop, T point[DESCRIPTOR::d]) +{ + return 0; +} + +template<typename T, typename DESCRIPTOR> +void Dynamics<T,DESCRIPTOR>::defineRho(int iPop, T rho) +{ } + +template<typename T, typename DESCRIPTOR> +void Dynamics<T,DESCRIPTOR>::defineU(const T u[DESCRIPTOR::d]) +{ } + +template<typename T, typename DESCRIPTOR> +void Dynamics<T,DESCRIPTOR>::defineU(int iPop, const T u[DESCRIPTOR::d]) +{ } + +template<typename T, typename DESCRIPTOR> +T Dynamics<T,DESCRIPTOR>::getVelocityCoefficient(int iPop) +{ + return 0; +} + +////////////////////// Class BasicDynamics //////////////////////// + +template<typename T, typename DESCRIPTOR> +BasicDynamics<T,DESCRIPTOR>::BasicDynamics(Momenta<T,DESCRIPTOR>& momenta) + : _momenta(momenta) +{ } + +template<typename T, typename DESCRIPTOR> +T BasicDynamics<T,DESCRIPTOR>::computeRho(Cell<T,DESCRIPTOR> const& cell) const +{ + return _momenta.computeRho(cell); +} + +template<typename T, typename DESCRIPTOR> +void BasicDynamics<T,DESCRIPTOR>::computeU ( + Cell<T,DESCRIPTOR> const& cell, + T u[DESCRIPTOR::d]) const +{ + _momenta.computeU(cell, u); +} + +template<typename T, typename DESCRIPTOR> +void BasicDynamics<T,DESCRIPTOR>::computeJ ( + Cell<T,DESCRIPTOR> const& cell, + T j[DESCRIPTOR::d]) const +{ + _momenta.computeJ(cell, j); +} + +template<typename T, typename DESCRIPTOR> +void BasicDynamics<T,DESCRIPTOR>::computeStress ( + Cell<T,DESCRIPTOR> const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal<DESCRIPTOR >::n] ) const +{ + _momenta.computeStress(cell, rho, u, pi); +} + +template<typename T, typename DESCRIPTOR> +void BasicDynamics<T,DESCRIPTOR>::computeRhoU ( + Cell<T,DESCRIPTOR> const& cell, + T& rho, T u[DESCRIPTOR::d]) const +{ + _momenta.computeRhoU(cell, rho, u); +} + +template<typename T, typename DESCRIPTOR> +void BasicDynamics<T,DESCRIPTOR>::computeAllMomenta ( + Cell<T,DESCRIPTOR> const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal<DESCRIPTOR >::n] ) const +{ + this->computeRhoU(cell, rho, u); + this->computeStress(cell, rho, u, pi); +} + +template<typename T, typename DESCRIPTOR> +void BasicDynamics<T,DESCRIPTOR>::defineRho(Cell<T,DESCRIPTOR>& cell, T rho) +{ + _momenta.defineRho(cell, rho); +} + +template<typename T, typename DESCRIPTOR> +void BasicDynamics<T,DESCRIPTOR>::defineU ( + Cell<T,DESCRIPTOR>& cell, + const T u[DESCRIPTOR::d]) +{ + _momenta.defineU(cell, u); +} + +template<typename T, typename DESCRIPTOR> +void BasicDynamics<T,DESCRIPTOR>::defineRhoU ( + Cell<T,DESCRIPTOR>& cell, + T rho, const T u[DESCRIPTOR::d]) +{ + _momenta.defineRhoU(cell, rho, u); +} + +template<typename T, typename DESCRIPTOR> +void BasicDynamics<T,DESCRIPTOR>::defineAllMomenta ( + Cell<T,DESCRIPTOR>& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal<DESCRIPTOR >::n] ) +{ + _momenta.defineAllMomenta(cell, rho, u, pi); +} + + +////////////////////// Class BGKdynamics ////////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template<typename T, typename DESCRIPTOR> +BGKdynamics<T,DESCRIPTOR>::BGKdynamics ( + T omega, Momenta<T,DESCRIPTOR>& momenta ) + : BasicDynamics<T,DESCRIPTOR>(momenta), + _omega(omega) +{ } + +template<typename T, typename DESCRIPTOR> +void BGKdynamics<T,DESCRIPTOR>::collide ( + Cell<T,DESCRIPTOR>& cell, + LatticeStatistics<T>& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + T uSqr = lbHelpers<T,DESCRIPTOR>::bgkCollision(cell, rho, u, _omega); + statistics.incrementStats(rho, uSqr); +} + +template<typename T, typename DESCRIPTOR> +T BGKdynamics<T,DESCRIPTOR>::getOmega() const +{ + return _omega; +} + +template<typename T, typename DESCRIPTOR> +void BGKdynamics<T,DESCRIPTOR>::setOmega(T omega) +{ + _omega = omega; +} + + +////////////////////// Class TRTdynamics ////////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template<typename T, typename DESCRIPTOR> +TRTdynamics<T,DESCRIPTOR>::TRTdynamics ( + T omega, Momenta<T,DESCRIPTOR>& momenta, T magicParameter ) + : BasicDynamics<T,DESCRIPTOR>(momenta), + _omega(omega), _omega2(1/(magicParameter/(1/omega-0.5)+0.5)), _magicParameter(magicParameter) +{ } + +template<typename T, typename DESCRIPTOR> +void TRTdynamics<T,DESCRIPTOR>::collide ( + Cell<T,DESCRIPTOR>& cell, + LatticeStatistics<T>& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + // T uSqr = lbHelpers<T,DESCRIPTOR>::bgkCollision(cell, rho, u, _omega); + const T uSqr = util::normSqr<T,DESCRIPTOR::d>(u); + 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<DESCRIPTOR>(iPop)] ); + fMinus[iPop] = 0.5 * ( cell[iPop] - cell[descriptors::opposite<DESCRIPTOR>(iPop)] ); + fEq[iPop] = lbHelpers<T, DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr); + } + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fEqPlus[iPop] = 0.5 * ( fEq[iPop] + fEq[descriptors::opposite<DESCRIPTOR>(iPop)] ); + fEqMinus[iPop] = 0.5 * ( fEq[iPop] - fEq[descriptors::opposite<DESCRIPTOR>(iPop)] ); + } + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] -= _omega * (fPlus[iPop] - fEqPlus[iPop]) + _omega2 * (fMinus[iPop] - fEqMinus[iPop]); + } + statistics.incrementStats(rho, uSqr); +} + +template<typename T, typename DESCRIPTOR> +T TRTdynamics<T,DESCRIPTOR>::getOmega() const +{ + return _omega; +} + +template<typename T, typename DESCRIPTOR> +void TRTdynamics<T,DESCRIPTOR>::setOmega(T omega) +{ + _omega = omega; +} + +////////////////////// Class ConstRhoBGKdynamics ////////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template<typename T, typename DESCRIPTOR> +ConstRhoBGKdynamics<T,DESCRIPTOR>::ConstRhoBGKdynamics ( + T omega, Momenta<T,DESCRIPTOR>& momenta ) + : BasicDynamics<T,DESCRIPTOR>(momenta), + _omega(omega) +{ } + +template<typename T, typename DESCRIPTOR> +void ConstRhoBGKdynamics<T,DESCRIPTOR>::collide ( + Cell<T,DESCRIPTOR>& cell, + LatticeStatistics<T>& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + + T deltaRho = (T)1 - (statistics).getAverageRho(); + T ratioRho = (T)1 + deltaRho/rho; + + T uSqr = lbHelpers<T,DESCRIPTOR>::constRhoBgkCollision ( + cell, rho, u, ratioRho, _omega ); + statistics.incrementStats(rho+deltaRho, uSqr); +} + +template<typename T, typename DESCRIPTOR> +T ConstRhoBGKdynamics<T,DESCRIPTOR>::getOmega() const +{ + return _omega; +} + +template<typename T, typename DESCRIPTOR> +void ConstRhoBGKdynamics<T,DESCRIPTOR>::setOmega(T omega) +{ + _omega = omega; +} + +////////////////////// Class IncBGKdynamics ////////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template<typename T, typename DESCRIPTOR> +IncBGKdynamics<T,DESCRIPTOR>::IncBGKdynamics ( + T omega, Momenta<T,DESCRIPTOR>& momenta ) + : BasicDynamics<T,DESCRIPTOR>(momenta), _omega(omega) +{ } + +template<typename T, typename DESCRIPTOR> +void IncBGKdynamics<T,DESCRIPTOR>::collide ( + Cell<T,DESCRIPTOR>& cell, + LatticeStatistics<T>& statistics ) +{ + T rho = this->_momenta.computeRho(cell); + T p = rho / descriptors::invCs2<T,DESCRIPTOR>(); + T j[DESCRIPTOR::d]; + this->_momenta.computeJ(cell, j); + T uSqr = lbHelpers<T,DESCRIPTOR>::incBgkCollision(cell, p, j, _omega); + statistics.incrementStats(rho, uSqr); +} + +template<typename T, typename DESCRIPTOR> +T IncBGKdynamics<T,DESCRIPTOR>::getOmega() const +{ + return _omega; +} + +template<typename T, typename DESCRIPTOR> +void IncBGKdynamics<T,DESCRIPTOR>::setOmega(T omega) +{ + _omega = omega; +} + + + +////////////////////// Class RLBdynamics ///////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template<typename T, typename DESCRIPTOR> +RLBdynamics<T,DESCRIPTOR>::RLBdynamics ( + T omega, Momenta<T,DESCRIPTOR>& momenta ) + : BasicDynamics<T,DESCRIPTOR>(momenta), + _omega(omega) +{ } + +template<typename T, typename DESCRIPTOR> +void RLBdynamics<T,DESCRIPTOR>::collide ( + Cell<T,DESCRIPTOR>& cell, + LatticeStatistics<T>& statistics ) +{ + T rho, u[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR >::n]; + this->_momenta.computeAllMomenta(cell, rho, u, pi); + T uSqr = rlbHelpers<T,DESCRIPTOR>::rlbCollision(cell, rho, u, pi, _omega); + statistics.incrementStats(rho, uSqr); +} + +template<typename T, typename DESCRIPTOR> +T RLBdynamics<T,DESCRIPTOR>::getOmega() const +{ + return _omega; +} + +template<typename T, typename DESCRIPTOR> +void RLBdynamics<T,DESCRIPTOR>::setOmega(T omega) +{ + _omega = omega; +} + +////////////////////// Class CombinedRLBdynamics ///////////////////////// + +template<typename T, typename DESCRIPTOR, typename Dynamics> +CombinedRLBdynamics<T,DESCRIPTOR,Dynamics>::CombinedRLBdynamics ( + T omega, Momenta<T,DESCRIPTOR>& momenta ) + : BasicDynamics<T,DESCRIPTOR>(momenta), + _boundaryDynamics(omega, momenta) +{ } + +template<typename T, typename DESCRIPTOR, typename Dynamics> +T CombinedRLBdynamics<T,DESCRIPTOR,Dynamics>:: +computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const +{ + return _boundaryDynamics.computeEquilibrium(iPop, rho, u, uSqr); +} + +template<typename T, typename DESCRIPTOR, typename Dynamics> +void CombinedRLBdynamics<T,DESCRIPTOR,Dynamics>::collide ( + Cell<T,DESCRIPTOR>& cell, + LatticeStatistics<T>& statistics ) +{ + typedef DESCRIPTOR L; + + T rho, u[L::d], pi[util::TensorVal<DESCRIPTOR >::n]; + this->_momenta.computeAllMomenta(cell,rho,u,pi); + + T uSqr = util::normSqr<T,L::d>(u); + + for (int iPop = 0; iPop < L::q; ++iPop) { + cell[iPop] = computeEquilibrium(iPop,rho,u,uSqr) + + firstOrderLbHelpers<T,DESCRIPTOR>::fromPiToFneq(iPop, pi); + } + + _boundaryDynamics.collide(cell, statistics); +} + +template<typename T, typename DESCRIPTOR, typename Dynamics> +T CombinedRLBdynamics<T,DESCRIPTOR,Dynamics>::getOmega() const +{ + return _boundaryDynamics.getOmega(); +} + +template<typename T, typename DESCRIPTOR, typename Dynamics> +void CombinedRLBdynamics<T,DESCRIPTOR,Dynamics>::setOmega(T omega) +{ + _boundaryDynamics.setOmega(omega); +} + + +////////////////////// Class ForcedBGKdynamics ///////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta Momenta object to know how to compute velocity momenta + */ +template<typename T, typename DESCRIPTOR> +ForcedBGKdynamics<T,DESCRIPTOR>::ForcedBGKdynamics ( + T omega, Momenta<T,DESCRIPTOR>& momenta ) + : BasicDynamics<T,DESCRIPTOR>(momenta), _omega(omega) +{ + OLB_PRECONDITION( DESCRIPTOR::template provides<descriptors::FORCE>() ); +} + +template<typename T, typename DESCRIPTOR> +void ForcedBGKdynamics<T,DESCRIPTOR>::computeU (Cell<T,DESCRIPTOR> const& cell, T u[DESCRIPTOR::d] ) const +{ + T rho; + this->_momenta.computeRhoU(cell, rho, u); + const T* force = cell.template getFieldPointer<descriptors::FORCE>(); + for (int iVel=0; iVel < DESCRIPTOR::d; ++iVel) { + u[iVel] += force[iVel] / (T)2.; + } +} + +template<typename T, typename DESCRIPTOR> +void ForcedBGKdynamics<T,DESCRIPTOR>::computeRhoU (Cell<T,DESCRIPTOR> const& cell, T& rho, T u[DESCRIPTOR::d] ) const +{ + this->_momenta.computeRhoU(cell, rho, u); + const T* force = cell.template getFieldPointer<descriptors::FORCE>(); + for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) { + u[iVel] += force[iVel] / (T)2.; + } +} + + +template<typename T, typename DESCRIPTOR> +void ForcedBGKdynamics<T,DESCRIPTOR>::collide ( + Cell<T,DESCRIPTOR>& cell, + LatticeStatistics<T>& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + const T* force = cell.template getFieldPointer<descriptors::FORCE>(); + for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) { + u[iVel] += force[iVel] / (T)2.; + } + T uSqr = lbHelpers<T,DESCRIPTOR>::bgkCollision(cell, rho, u, _omega); + lbHelpers<T,DESCRIPTOR>::addExternalForce(cell, u, _omega, rho); + statistics.incrementStats(rho, uSqr); +} + +template<typename T, typename DESCRIPTOR> +T ForcedBGKdynamics<T,DESCRIPTOR>::getOmega() const +{ + return _omega; +} + +template<typename T, typename DESCRIPTOR> +void ForcedBGKdynamics<T,DESCRIPTOR>::setOmega(T omega) +{ + _omega = omega; +} + +////////////////////// Class ForcedKupershtokhBGKdynamics ///////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta Momenta object to know how to compute velocity momenta + */ +template<typename T, typename DESCRIPTOR> +ForcedKupershtokhBGKdynamics<T,DESCRIPTOR>::ForcedKupershtokhBGKdynamics ( + T omega, Momenta<T,DESCRIPTOR>& momenta ) + : BasicDynamics<T,DESCRIPTOR>(momenta), _omega(omega) +{ + OLB_PRECONDITION( DESCRIPTOR::template provides<descriptors::FORCE>() ); +} + +template<typename T, typename DESCRIPTOR> +void ForcedKupershtokhBGKdynamics<T,DESCRIPTOR>::computeU (Cell<T,DESCRIPTOR> const& cell, T u[DESCRIPTOR::d] ) const +{ + T rho; + this->_momenta.computeRhoU(cell, rho, u); + for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) { + u[iVel] += cell.template getFieldPointer<descriptors::FORCE>()[iVel] / (T)2.; + } +} + +template<typename T, typename DESCRIPTOR> +void ForcedKupershtokhBGKdynamics<T,DESCRIPTOR>::computeRhoU (Cell<T,DESCRIPTOR> const& cell, T& rho, T u[DESCRIPTOR::d] ) const +{ + this->_momenta.computeRhoU(cell, rho, u); + for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) { + u[iVel] += cell.template getFieldPointer<descriptors::FORCE>()[iVel] / (T)2.; + } +} + + +template<typename T, typename DESCRIPTOR> +void ForcedKupershtokhBGKdynamics<T,DESCRIPTOR>::collide ( + Cell<T,DESCRIPTOR>& cell, + LatticeStatistics<T>& statistics ) +{ + T rho, u[DESCRIPTOR::d], uPlusDeltaU[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + T* force = cell.template getFieldPointer<descriptors::FORCE>(); + T uSqr = lbHelpers<T,DESCRIPTOR>::bgkCollision(cell, rho, u, _omega); + + for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) { + uPlusDeltaU[iVel] = u[iVel] + force[iVel]; + } + const T uPlusDeltaUSqr = util::normSqr<T,DESCRIPTOR::d>(uPlusDeltaU); + + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] += lbHelpers<T, DESCRIPTOR>::equilibrium(iPop, rho, uPlusDeltaU, uPlusDeltaUSqr) + - lbHelpers<T, DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr); + } + + statistics.incrementStats(rho, uSqr); +} + +template<typename T, typename DESCRIPTOR> +T ForcedKupershtokhBGKdynamics<T,DESCRIPTOR>::getOmega() const +{ + return _omega; +} + +template<typename T, typename DESCRIPTOR> +void ForcedKupershtokhBGKdynamics<T,DESCRIPTOR>::setOmega(T omega) +{ + _omega = omega; +} + + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta Momenta object to know how to compute velocity momenta + * \param sink counterpart of a source term + */ +template<typename T, typename DESCRIPTOR> +PoissonDynamics<T,DESCRIPTOR>::PoissonDynamics ( + T omega, Momenta<T,DESCRIPTOR>& momenta , T sink) + : BasicDynamics<T,DESCRIPTOR>(momenta), _omega(omega), _sink(sink) +{ +} + +template<typename T, typename DESCRIPTOR> +T PoissonDynamics<T,DESCRIPTOR>::computeEquilibrium( int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr ) const +{ + return descriptors::t<T,DESCRIPTOR>(iPop) * rho - descriptors::t<T,DESCRIPTOR>(iPop); +} + + +template<typename T, typename DESCRIPTOR> +T PoissonDynamics<T,DESCRIPTOR>::computeRho(Cell<T,DESCRIPTOR> const& cell) const +{ + return lbHelpers<T,DESCRIPTOR>::computeRho(cell); +} + +template<typename T, typename DESCRIPTOR> +void PoissonDynamics<T,DESCRIPTOR>::computeU(Cell<T,DESCRIPTOR> const& cell, T u[DESCRIPTOR::d]) const +{ + for ( int iDim = 0; iDim < DESCRIPTOR::d; iDim++ ) { + u[iDim] = T(); + } +} + +template<typename T, typename DESCRIPTOR> +void PoissonDynamics<T,DESCRIPTOR>::computeJ(Cell<T,DESCRIPTOR> const& cell, T j[DESCRIPTOR::d]) const +{ + lbHelpers<T,DESCRIPTOR>::computeJ(cell, j); +} + +template<typename T, typename DESCRIPTOR> +void PoissonDynamics<T,DESCRIPTOR>::computeStress(Cell<T,DESCRIPTOR> const& cell, T rho, + const T u[DESCRIPTOR::d], T pi[util::TensorVal<DESCRIPTOR>::n] ) const +{ + for ( int iDim = 0; iDim < DESCRIPTOR::d; iDim++ ) { + pi[iDim] = T(); + } +} + +template<typename T, typename DESCRIPTOR> +void PoissonDynamics<T,DESCRIPTOR>::computeRhoU( Cell<T,DESCRIPTOR> const& cell, T& rho, T u[DESCRIPTOR::d]) const +{ + rho = computeRho(cell); + computeU(cell, u); +} + +template<typename T, typename DESCRIPTOR> +void PoissonDynamics<T,DESCRIPTOR>::computeAllMomenta( Cell<T,DESCRIPTOR> const& cell, T &rho, + T u[DESCRIPTOR::q], T pi[util::TensorVal<DESCRIPTOR>::n] ) const +{ + rho = computeRho(cell); + computeU(cell, u); + computeStress(cell, rho, u, pi); +} + + +template<typename T, typename DESCRIPTOR> +void PoissonDynamics<T,DESCRIPTOR>::collide( Cell<T,DESCRIPTOR>& cell, LatticeStatistics<T>& statistics ) +{ + T rho = computeRho(cell); + + for ( int iPop = 0; iPop < DESCRIPTOR::q; ++iPop ) { + cell[iPop] = (1 - _omega) * (cell[iPop] + descriptors::t<T,DESCRIPTOR>(iPop)) + + _omega * descriptors::t<T,DESCRIPTOR>(iPop) * rho + - _sink*(cell[iPop] + descriptors::t<T,DESCRIPTOR>(iPop)) - descriptors::t<T,DESCRIPTOR>(iPop); + } + +// // add spherical harmonic definition f_i = 3*pi/4 *density + pi/4 *flux \cdot e_i +// T density_post_collision = lbDynamicsHelpers<T,DESCRIPTOR>::computeRho(cell); +// T flux_post_collision[DESCRIPTOR::q]; +// lbDynamicsHelpers<T,DESCRIPTOR>::computeJ(cell,flux_post_collision); +// T s = flux_post_collision[0]*descriptors::t<T,DESCRIPTOR>(0)+ +// flux_post_collision[1]*descriptors::t<T,DESCRIPTOR>(1)+ +// flux_post_collision[2]*descriptors::t<T,DESCRIPTOR>(2); +// std::cout << "s " << s << std::endl; +// for ( int iPop = 0; iPop < DESCRIPTOR::q; ++iPop ) { +// cell[iPop] = 3*M_PI/4. *density_post_collision +// + M_PI/4. *s; +// } + statistics.incrementStats(rho, 0); +} + +template<typename T, typename DESCRIPTOR> +T PoissonDynamics<T,DESCRIPTOR>::getOmega() const +{ + return _omega; +} + +template<typename T, typename DESCRIPTOR> +void PoissonDynamics<T,DESCRIPTOR>::setOmega(T omega) +{ + _omega = omega; +} + + +template<typename T, typename DESCRIPTOR> +P1Dynamics<T,DESCRIPTOR>::P1Dynamics ( + T omega, Momenta<T,DESCRIPTOR>& momenta , T absorption, T scattering) + : BasicDynamics<T,DESCRIPTOR>(momenta), _omega(omega), _absorption(absorption), _scattering(scattering) +{ +} + +template<typename T, typename DESCRIPTOR> +T P1Dynamics<T,DESCRIPTOR>::computeEquilibrium( int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr ) const +{ + return descriptors::t<T,DESCRIPTOR>(iPop) * rho - descriptors::t<T,DESCRIPTOR>(iPop); +} + + +template<typename T, typename DESCRIPTOR> +T P1Dynamics<T,DESCRIPTOR>::computeRho(Cell<T,DESCRIPTOR> const& cell) const +{ + return lbHelpers<T,DESCRIPTOR>::computeRho(cell); +} + +template<typename T, typename DESCRIPTOR> +void P1Dynamics<T,DESCRIPTOR>::computeU(Cell<T,DESCRIPTOR> const& cell, T u[DESCRIPTOR::d]) const +{ + for ( int iDim = 0; iDim < DESCRIPTOR::d; iDim++ ) { + u[iDim] = T(); + } +} + +template<typename T, typename DESCRIPTOR> +void P1Dynamics<T,DESCRIPTOR>::computeJ(Cell<T,DESCRIPTOR> const& cell, T j[DESCRIPTOR::d]) const +{ + lbHelpers<T,DESCRIPTOR>::computeJ(cell, j); +} + +template<typename T, typename DESCRIPTOR> +void P1Dynamics<T,DESCRIPTOR>::computeStress(Cell<T,DESCRIPTOR> const& cell, T rho, + const T u[DESCRIPTOR::d], T pi[util::TensorVal<DESCRIPTOR>::n] ) const +{ + for ( int iDim = 0; iDim < DESCRIPTOR::d; iDim++ ) { + pi[iDim] = T(); + } +} + +template<typename T, typename DESCRIPTOR> +void P1Dynamics<T,DESCRIPTOR>::computeRhoU( Cell<T,DESCRIPTOR> const& cell, T& rho, T u[DESCRIPTOR::d]) const +{ + rho = computeRho(cell); + computeU(cell, u); +} + +template<typename T, typename DESCRIPTOR> +void P1Dynamics<T,DESCRIPTOR>::computeAllMomenta( Cell<T,DESCRIPTOR> const& cell, T &rho, + T u[DESCRIPTOR::q], T pi[util::TensorVal<DESCRIPTOR>::n] ) const +{ + rho = computeRho(cell); + computeU(cell, u); + computeStress(cell, rho, u, pi); +} + + +template<typename T, typename DESCRIPTOR> +void P1Dynamics<T,DESCRIPTOR>::collide( Cell<T,DESCRIPTOR>& cell, LatticeStatistics<T>& statistics ) +{ + T rho = computeRho(cell); + + // compute equilibrium f^eq_i eq. (5.32) from DOI:10.1002/9780470177013 + + // ----move to computeEq + // give me a std::array, please + std::array<T,DESCRIPTOR::q> cellShifted; + for (int iPop = 0; iPop <DESCRIPTOR::q; ++iPop) { + cellShifted[iPop] = cell[iPop] + descriptors::t<T,DESCRIPTOR>(iPop); + } + double moment0 = std::accumulate(cellShifted.begin(),cellShifted.end(), T(0)); + std::array<T,DESCRIPTOR::d> moment1; + moment1.fill( T(0) ); + // sum_j v_j f_j + for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) { + for (int iDim = 0; iDim < DESCRIPTOR::d; ++iDim) { + moment1[iDim] += descriptors::c<DESCRIPTOR>(iPop,iDim)*cellShifted[iPop]; + } + } + // do scalar product in eq. sum_j v_j f_j * v_i + std::array<T,DESCRIPTOR::q> moment1_v; + moment1_v.fill( T(0) ); + for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) { + // scalar product + for (int iDim = 0; iDim < DESCRIPTOR::d; ++iDim) { + moment1_v[iPop] += moment1[iDim]*descriptors::c<DESCRIPTOR>(iPop,iDim); + } + } + std::array<T,DESCRIPTOR::q> fEq; + for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) { + fEq[iPop] = descriptors::t<T,DESCRIPTOR>(iPop) *(moment0 +moment1_v[iPop]); + //fEq[iPop] = 1./(4*M_PI) *moment0 +3./(4*M_PI) * moment1_v[iPop] -descriptors::t<T,DESCRIPTOR>(iPop); + } + + + std::array<T,DESCRIPTOR::q> relaxCell; + for (int iPop = 0; iPop <DESCRIPTOR::q; ++iPop) { + relaxCell[iPop] = + cellShifted[iPop] + - _scattering*(cellShifted[iPop] - fEq[iPop]) + - (_absorption+_scattering)*descriptors::t<T,DESCRIPTOR>(iPop)*cellShifted[iPop]; + } + + for (int iPop = 0; iPop <DESCRIPTOR::q; ++iPop) { + cell[iPop] = relaxCell[iPop] - descriptors::t<T,DESCRIPTOR>(iPop); + } + + statistics.incrementStats(rho, 0); +} + +template<typename T, typename DESCRIPTOR> +T P1Dynamics<T,DESCRIPTOR>::getOmega() const +{ + return _omega; +} + +template<typename T, typename DESCRIPTOR> +void P1Dynamics<T,DESCRIPTOR>::setOmega(T omega) +{ + _omega = omega; +} + +////////////////////// Class ResettingForcedBGKdynamics ///////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template<typename T, typename DESCRIPTOR> +ResettingForcedBGKdynamics<T,DESCRIPTOR>::ResettingForcedBGKdynamics ( + T omega, Momenta<T,DESCRIPTOR>& momenta ) + : ForcedBGKdynamics<T,DESCRIPTOR>(omega, momenta) +{ + OLB_PRECONDITION( DESCRIPTOR::template provides<descriptors::FORCE>() ); +} + +template<typename T, typename DESCRIPTOR> +void ResettingForcedBGKdynamics<T,DESCRIPTOR>::collide ( + Cell<T,DESCRIPTOR>& cell, + LatticeStatistics<T>& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + T* force = cell.template getFieldPointer<descriptors::FORCE>(); + if ( !util::nearZero(force[0]) || !util::nearZero(force[1]) || !util::nearZero(force[2]) ) // TODO: unnecessary?? + for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) { + u[iVel] += force[iVel] / (T)2.; + } +// if (force[2] != 0) +// std::cout << force[0] << " " << force[1] << " " << force[2] << std::endl; + T uSqr = lbHelpers<T,DESCRIPTOR>::bgkCollision(cell, rho, u, this->_omega); + lbHelpers<T,DESCRIPTOR>::addExternalForce(cell, u, this->_omega, rho); + statistics.incrementStats(rho, uSqr); + + force[0] = _frc[0]; + force[1] = _frc[1]; + force[2] = _frc[2]; +// force[0] = 0.; +// force[1] = 0.; +// force[2] = 0.; +} + +////////////////////// Class ForcedShanChenBGKdynamics ///////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template<typename T, typename DESCRIPTOR> +ForcedShanChenBGKdynamics<T,DESCRIPTOR>::ForcedShanChenBGKdynamics ( + T omega, Momenta<T,DESCRIPTOR>& momenta ) + : ForcedBGKdynamics<T,DESCRIPTOR>(omega, momenta ) +{ + OLB_PRECONDITION( DESCRIPTOR::template provides<descriptors::FORCE>() ); +} + +template<typename T, typename DESCRIPTOR> +void ForcedShanChenBGKdynamics<T,DESCRIPTOR>::computeU (Cell<T,DESCRIPTOR> const& cell, T u[DESCRIPTOR::d] ) const +{ + T rho; + this->_momenta.computeRhoU(cell, rho, u); + for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) { + u[iVel] += cell.template getFieldPointer<descriptors::FORCE>()[iVel] / (T)2.; + } +} + +template<typename T, typename DESCRIPTOR> +void ForcedShanChenBGKdynamics<T,DESCRIPTOR>::computeRhoU (Cell<T,DESCRIPTOR> const& cell, T& rho, T u[DESCRIPTOR::d] ) const +{ + this->_momenta.computeRhoU(cell, rho, u); + for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) { + u[iVel] += cell.template getFieldPointer<descriptors::FORCE>()[iVel] / (T)2.; + } +} + +template<typename T, typename DESCRIPTOR> +void ForcedShanChenBGKdynamics<T,DESCRIPTOR>::collide ( + Cell<T,DESCRIPTOR>& cell, + LatticeStatistics<T>& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + T* force = cell.template getFieldPointer<descriptors::FORCE>(); + for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) { + u[iVel] += force[iVel] / this->getOmega(); + } + T uSqr = lbHelpers<T,DESCRIPTOR>::bgkCollision(cell, rho, u, this->getOmega() ); + uSqr=0.; + for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) { + u[iVel] += force[iVel] / (T)2.; + u[iVel] -= force[iVel] / this->getOmega(); + uSqr += u[iVel]*u[iVel]; + } + statistics.incrementStats(rho, uSqr); +} + +////////////////////// Class ForcedTRTdynamics ///////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta Momenta object to know how to compute velocity momenta + */ +template<typename T, typename DESCRIPTOR> +ForcedTRTdynamics<T,DESCRIPTOR>::ForcedTRTdynamics ( + T omega, Momenta<T,DESCRIPTOR>& momenta, T magicParameter ) + : BasicDynamics<T,DESCRIPTOR>(momenta), + _omega(omega), _omega2(1/(magicParameter/(1/omega-0.5)+0.5)), _magicParameter(magicParameter) +{ + OLB_PRECONDITION( DESCRIPTOR::template provides<descriptors::FORCE>() ); +} + +template<typename T, typename DESCRIPTOR> +void ForcedTRTdynamics<T,DESCRIPTOR>::computeU (Cell<T,DESCRIPTOR> const& cell, T u[DESCRIPTOR::d] ) const +{ + T rho; + this->_momenta.computeRhoU(cell, rho, u); + for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) { + u[iVel] += cell.template getFieldPointer<descriptors::FORCE>()[iVel] / (T)2.; + } +} + +template<typename T, typename DESCRIPTOR> +void ForcedTRTdynamics<T,DESCRIPTOR>::computeRhoU (Cell<T,DESCRIPTOR> const& cell, T& rho, T u[DESCRIPTOR::d] ) const +{ + this->_momenta.computeRhoU(cell, rho, u); + for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) { + u[iVel] += cell.template getFieldPointer<descriptors::FORCE>()[iVel] / (T)2.; + } +} + + +template<typename T, typename DESCRIPTOR> +void ForcedTRTdynamics<T,DESCRIPTOR>::collide ( + Cell<T,DESCRIPTOR>& cell, + LatticeStatistics<T>& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + T* force = cell.template getFieldPointer<descriptors::FORCE>(); + for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) { + u[iVel] += force[iVel] / (T)2.; + } + const T uSqr = util::normSqr<T,DESCRIPTOR::d>(u); + 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<DESCRIPTOR>(iPop)] ); + fMinus[iPop] = 0.5 * ( cell[iPop] - cell[descriptors::opposite<DESCRIPTOR>(iPop)] ); + fEq[iPop] = lbHelpers<T, DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr); + } + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fEqPlus[iPop] = 0.5 * ( fEq[iPop] + fEq[descriptors::opposite<DESCRIPTOR>(iPop)] ); + fEqMinus[iPop] = 0.5 * ( fEq[iPop] - fEq[descriptors::opposite<DESCRIPTOR>(iPop)] ); + } + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] -= _omega * (fPlus[iPop] - fEqPlus[iPop]) + _omega2 * (fMinus[iPop] - fEqMinus[iPop]); + } + lbHelpers<T,DESCRIPTOR>::addExternalForce(cell, u, _omega, rho); + statistics.incrementStats(rho, uSqr); +} + +template<typename T, typename DESCRIPTOR> +T ForcedTRTdynamics<T,DESCRIPTOR>::getOmega() const +{ + return _omega; +} + +template<typename T, typename DESCRIPTOR> +void ForcedTRTdynamics<T,DESCRIPTOR>::setOmega(T omega) +{ + _omega = omega; +} +////////////////////// Class D3Q13dynamics ///////////////////////// + +/** \param omega relaxation parameter, related to the dynamic viscosity + * \param momenta a Momenta object to know how to compute velocity momenta + */ +template<typename T, typename DESCRIPTOR> +D3Q13dynamics<T,DESCRIPTOR>::D3Q13dynamics ( + T omega, Momenta<T,DESCRIPTOR>& momenta ) + : BasicDynamics<T,DESCRIPTOR>(momenta) +{ + setOmega(omega); +} + +template<typename T, typename DESCRIPTOR> +T D3Q13dynamics<T,DESCRIPTOR>::computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const +{ + // To get at the equilibrium, execute collision with relaxation parameters 1 + Cell<T,DESCRIPTOR> tmp; + for (int pop=0; pop<DESCRIPTOR::q; ++pop) { + tmp[pop] = descriptors::t<T,DESCRIPTOR>(pop); + } + d3q13Helpers<T>::collision(tmp, rho, u, (T)1, (T)1); + return tmp[iPop]; +} + +template<typename T, typename DESCRIPTOR> +void D3Q13dynamics<T,DESCRIPTOR>::collide ( + Cell<T,DESCRIPTOR>& cell, + LatticeStatistics<T>& statistics ) +{ + T rho, u[DESCRIPTOR::d]; + this->_momenta.computeRhoU(cell, rho, u); + T uSqr = d3q13Helpers<T>::collision ( + cell, rho, u, lambda_nu, lambda_nu_prime ); + statistics.incrementStats(rho, uSqr); +} + +template<typename T, typename DESCRIPTOR> +T D3Q13dynamics<T,DESCRIPTOR>::getOmega() const +{ + return (T)4 / ( (T)3/lambda_nu + (T)1/(T)2 ); +} + +template<typename T, typename DESCRIPTOR> +void D3Q13dynamics<T,DESCRIPTOR>::setOmega(T omega) +{ + lambda_nu = (T)3 / ( (T)4/omega - (T)1/(T)2 ); + lambda_nu_prime = (T)3 / ( (T)2/omega + (T)1/(T)2 ); +} + +////////////////////// Class Momenta ////////////////////////////// + +template<typename T, typename DESCRIPTOR> +void Momenta<T,DESCRIPTOR>::computeRhoU ( + Cell<T,DESCRIPTOR> const& cell, + T& rho, T u[DESCRIPTOR::d]) const +{ + rho = this->computeRho(cell); + this->computeU(cell, u); + +} + +template<typename T, typename DESCRIPTOR> +void Momenta<T,DESCRIPTOR>::computeAllMomenta ( + Cell<T,DESCRIPTOR> const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal<DESCRIPTOR >::n] ) const +{ + this->computeRhoU(cell, rho, u); + this->computeStress(cell, rho, u, pi); +} + +template<typename T, typename DESCRIPTOR> +void Momenta<T,DESCRIPTOR>::defineRhoU ( + Cell<T,DESCRIPTOR>& cell, + T rho, const T u[DESCRIPTOR::d]) +{ + this->defineRho(cell, rho); + this->defineU(cell, u); + +} + +////////////////////// Class BulkMomenta ////////////////////////// + +template<typename T, typename DESCRIPTOR> +T BulkMomenta<T,DESCRIPTOR>::computeRho(Cell<T,DESCRIPTOR> const |