From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Mon, 24 Jun 2019 14:43:36 +0200 Subject: Initialize at openlb-1-3 --- src/dynamics/dynamics.h | 995 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 995 insertions(+) create mode 100644 src/dynamics/dynamics.h (limited to 'src/dynamics/dynamics.h') diff --git a/src/dynamics/dynamics.h b/src/dynamics/dynamics.h new file mode 100644 index 0000000..40970b5 --- /dev/null +++ b/src/dynamics/dynamics.h @@ -0,0 +1,995 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2006, 2007 Jonas Latt, Mathias J. Krause + * 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 -- header file. + */ +#ifndef LB_DYNAMICS_H +#define LB_DYNAMICS_H + +#include "latticeDescriptors.h" +#include "core/util.h" +#include "core/postProcessing.h" +#include "core/latticeStatistics.h" + +namespace olb { + + +template class Cell; + +/// Interface for the dynamics classes +template +struct Dynamics { + /// Destructor: virtual to enable inheritance + virtual ~Dynamics() { } + /// Implementation of the collision step + virtual void collide(Cell& cell, + LatticeStatistics& statistics_) =0; + /// Compute equilibrium distribution function + virtual T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const; + /// Initialize cell at equilibrium distribution + void iniEquilibrium(Cell& cell, T rho, const T u[DESCRIPTOR::d]); + /// Compute particle density on the cell. + /** \return particle density + */ + virtual T computeRho(Cell const& cell) const =0; + /// Compute fluid velocity on the cell. + /** \param u fluid velocity + */ + virtual void computeU( Cell const& cell, + T u[DESCRIPTOR::d] ) const =0; + /// Compute fluid momentum (j=rho*u) on the cell. + /** \param j fluid momentum + */ + virtual void computeJ( Cell const& cell, + T j[DESCRIPTOR::d] ) const =0; + /// Compute components of the stress tensor on the cell. + /** \param pi stress tensor */ + virtual void computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const =0; + /// Compute fluid velocity and particle density on the cell. + /** \param rho particle density + * \param u fluid velocity + */ + virtual void computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const =0; + /// Compute all momenta on the cell, up to second order. + /** \param rho particle density + * \param u fluid velocity + * \param pi stress tensor + */ + virtual void computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const =0; + /// Set particle density on the cell. + /** \param rho particle density + */ + virtual void defineRho(Cell& cell, T rho) =0; + virtual void defineRho(int iPop, T rho); + /// Set fluid velocity on the cell. + /** \param u fluid velocity + */ + virtual void defineU(Cell& cell, + const T u[DESCRIPTOR::d]) =0; + /// Functions for offLattice Velocity boundary conditions + virtual void setBoundaryIntersection(int iPop, T distance); + virtual bool getBoundaryIntersection(int iPop, T point[DESCRIPTOR::d]); + virtual void defineU(const T u[DESCRIPTOR::d]); + virtual void defineU(int iPop, const T u[DESCRIPTOR::d]); + virtual T getVelocityCoefficient(int iPop); + + /// Define fluid velocity and particle density on the cell. + /** \param rho particle density + * \param u fluid velocity + */ + virtual void defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) =0; + /// Define all momenta on the cell, up to second order. + /** \param rho particle density + * \param u fluid velocity + * \param pi stress tensor + */ + virtual void defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) =0; + /// Get local relaxation parameter of the dynamics + virtual T getOmega() const =0; + /// Set local relaxation parameter of the dynamics + virtual void setOmega(T omega) =0; +}; + +/// Interface for classes that compute velocity momenta +/** This class is useful for example to distinguish between bulk and + * boundary nodes, given that on the boundaries, a particular strategy + * must be applied to compute the velocity momenta. + */ +template +struct Momenta { + /// Destructor: virtual to enable inheritance + virtual ~Momenta() { } + /// Compute particle density on the cell. + virtual T computeRho(Cell const& cell) const =0; + /// Compute fluid velocity on the cell. + virtual void computeU ( + Cell const& cell, + T u[DESCRIPTOR::d] ) const =0; + /// Compute fluid momentum on the cell. + virtual void computeJ ( + Cell const& cell, + T j[DESCRIPTOR::d] ) const =0; + /// Compute components of the stress tensor on the cell. + virtual void computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const =0; + /// Compute fluid velocity and particle density on the cell. + virtual void computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const; + /// Compute all momenta on the cell, up to second order. + virtual void computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const; + /// Set particle density on the cell. + virtual void defineRho(Cell& cell, T rho) =0; + /// Set fluid velocity on the cell. + virtual void defineU(Cell& cell, + const T u[DESCRIPTOR::d]) =0; + /// Define fluid velocity and particle density on the cell. + virtual void defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]); + /// Define all momenta on the cell, up to second order. + virtual void defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) =0; +}; + +/// Abstract base for dynamics classes +/** In this version of the Dynamics classes, the computation of the + * velocity momenta is taken care of by an object of type Momenta. + */ +template +class BasicDynamics : public Dynamics { +public: + /// Must be contructed with an object of type Momenta + BasicDynamics(Momenta& momenta); + /// Implemented via the Momenta object + T computeRho(Cell const& cell) const override; + /// Implemented via the Momenta object + void computeU ( + Cell const& cell, + T u[DESCRIPTOR::d] ) const override; + /// Implemented via the Momenta object + void computeJ ( + Cell const& cell, + T j[DESCRIPTOR::d] ) const override; + /// Implemented via the Momenta object + void computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const override; + /// Implemented via the Momenta object + void computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const override; + /// Implemented via the Momenta object + void computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const override; + /// Implemented via the Momenta object + void defineRho(Cell& cell, T rho) override; + /// Implemented via the Momenta object + void defineU(Cell& cell, + const T u[DESCRIPTOR::d]) override; + /// Implemented via the Momenta object + void defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) override; + /// Implemented via the Momenta object + void defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) override; +protected: + Momenta& _momenta; ///< computation of velocity momenta +}; + +/// Implementation of a generic dynamics to realize a pressure drop at a periodic boundary +template +class PeriodicPressureDynamics : public BaseDynamics { + +public: + /// Constructor + PeriodicPressureDynamics(BaseDynamics& baseDynamics, T densityOffset, int nx, int ny, int nz=0) : BaseDynamics(baseDynamics), _densityOffset(densityOffset), _nx(nx), _ny(ny), _nz(nz) {}; + + /// Implementation of the collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override{ + BaseDynamics::collide(cell,statistics_); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + if ( ((_nx==1 || _nx==-1) && descriptors::c(iPop,0)==_nx) || ((_ny==1 || _ny==-1) && descriptors::c(iPop,1)==_ny) || (DESCRIPTOR::d==3 && !_nz && descriptors::c(iPop,2)==_nz) ) { + cell[iPop] += (cell[iPop] + descriptors::t(iPop))*_densityOffset; + } + } +}; + +private: + T _densityOffset; + int _nx, _ny, _nz; +}; + +/// Implementation of the BGK collision step +template +class BGKdynamics : public BasicDynamics { +public: + /// Constructor + BGKdynamics(T omega, Momenta& momenta); + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Get local relaxation parameter of the dynamics + T getOmega() const override; + /// Set local relaxation parameter of the dynamics + void setOmega(T omega) override; +private: + T _omega; ///< relaxation parameter +}; + +/// Implementation of the TRT collision step +template +class TRTdynamics : public BasicDynamics { +public: + /// Constructor + TRTdynamics(T omega, Momenta& momenta, T magicParameter); + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Get local relaxation parameter of the dynamics + T getOmega() const override; + /// Set local relaxation parameter of the dynamics + void setOmega(T omega) override; +private: + T _omega; ///< relaxation parameter + T _omega2; /// relaxation parameter for odd moments + T _magicParameter; +}; + +/// Implementation of the pressure-corrected BGK collision step +template +class ConstRhoBGKdynamics : public BasicDynamics { +public: + /// Constructor + ConstRhoBGKdynamics(T omega, Momenta& momenta); + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Get local relaxation parameter of the dynamics + T getOmega() const override; + /// Set local relaxation parameter of the dynamics + void setOmega(T omega) override; +private: + T _omega; ///< relaxation parameter +}; + +/// Implementation of the so-called incompressible collision step +template +class IncBGKdynamics : public BasicDynamics { +public: + /// Constructor + IncBGKdynamics(T omega, Momenta& momenta); + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Get local relaxation parameter of the dynamics + T getOmega() const override; + /// Set local relaxation parameter of the dynamics + void setOmega(T omega) override; +private: + T _omega; ///< relaxation parameter +}; + + + +/// Implementation of the Regularized BGK collision step +/** This model is substantially more stable than plain BGK, and has roughly + * the same efficiency. However, it cuts out the modes at higher Knudsen + * numbers and can not be used in the regime of rarefied gases. + */ +template +class RLBdynamics : public BasicDynamics { +public: + /// Constructor + RLBdynamics(T omega, Momenta& momenta); + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Get local relaxation parameter of the dynamics + T getOmega() const override; + /// Set local relaxation parameter of the dynamics + void setOmega(T omega) override; +private: + T _omega; ///< relaxation parameter +}; + +/// Implementation of Regularized BGK collision, followed by any Dynamics +template +class CombinedRLBdynamics : public BasicDynamics { +public: + /// Constructor + CombinedRLBdynamics(T omega, Momenta& momenta); + /// Compute equilibrium distribution function + T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const override; + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Get local relaxation parameter of the dynamics + T getOmega() const override; + /// Set local relaxation parameter of the dynamics + void setOmega(T omega) override; +private: + Dynamics _boundaryDynamics; +}; + +/// Implementation of the BGK collision step with external force +template +class ForcedBGKdynamics : public BasicDynamics { +public: + /// Constructor + ForcedBGKdynamics(T omega, Momenta& momenta); + /// Compute fluid velocity on the cell. + void computeU ( + Cell const& cell, + T u[DESCRIPTOR::d] ) const override; + /// Compute fluid velocity and particle density on the cell. + void computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const override; + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Get local relaxation parameter of the dynamics + T getOmega() const override; + /// Set local relaxation parameter of the dynamics + void setOmega(T omega) override; +protected: + T _omega; ///< relaxation parameter +}; + +/// Implementation of the BGK collision step with external force +template +class ForcedKupershtokhBGKdynamics : public BasicDynamics { +public: + /// Constructor + ForcedKupershtokhBGKdynamics(T omega, Momenta& momenta); + /// Compute fluid velocity on the cell. + void computeU ( + Cell const& cell, + T u[DESCRIPTOR::d] ) const override; + /// Compute fluid velocity and particle density on the cell. + void computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const override; + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Get local relaxation parameter of the dynamics + T getOmega() const override; + /// Set local relaxation parameter of the dynamics + void setOmega(T omega) override; +protected: + T _omega; ///< relaxation parameter +}; + +// no momenta call, thus implement all function from BasicDynamics +template +class PoissonDynamics : public BasicDynamics { +public: + /// Constructor + PoissonDynamics(T omega, Momenta& momenta, T sink); + // TODO AM zerothOrderEquilibrium in lbHelpers + T computeEquilibrium( int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr ) const override; + T computeRho(Cell const& cell) const override; + void computeU(Cell const& cell, T u[DESCRIPTOR::d]) const override; + void computeJ(Cell const& cell, T j[DESCRIPTOR::d]) const override; + void computeStress(Cell const& cell, T rho, + const T u[DESCRIPTOR::d], T pi[util::TensorVal::n] ) const override; + void computeRhoU( Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const override; + void computeAllMomenta( Cell const& cell, T &rho, + T u[DESCRIPTOR::q], T pi[util::TensorVal::n] ) const override; + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + T getOmega() const override; + void setOmega(T omega) override; +protected: + T _omega; + T _sink; +}; + + +template +class P1Dynamics : public BasicDynamics { +public: + /// Constructor + P1Dynamics(T omega, Momenta& momenta, T absorption, T scattering); + // TODO AM zerothOrderEquilibrium in lbHelpers + T computeEquilibrium( int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr ) const override; + T computeRho(Cell const& cell) const override; + void computeU(Cell const& cell, T u[DESCRIPTOR::d]) const override; + void computeJ(Cell const& cell, T j[DESCRIPTOR::d]) const override; + void computeStress(Cell const& cell, T rho, + const T u[DESCRIPTOR::d], T pi[util::TensorVal::n] ) const override; + void computeRhoU( Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const override; + void computeAllMomenta( Cell const& cell, T &rho, + T u[DESCRIPTOR::q], T pi[util::TensorVal::n] ) const override; + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + T getOmega() const override; + void setOmega(T omega) override; +protected: + T _omega; + T _absorption; + T _scattering; +}; +/// Implementation of the BGK collision step with external force +template +class ResettingForcedBGKdynamics : public ForcedBGKdynamics { +public: + ResettingForcedBGKdynamics(T omega, Momenta& momenta); + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + inline void setForce(T force[3]) + { +// _frc[0] = force[0]; +// _frc[1] = force[1]; +// _frc[2] = force[2]; + _frc[0] = 0.0; + _frc[1] = 0.0; + _frc[2] = 0.0; + } +private: + T _frc[3]; +}; + +/// Other Implementation of the BGK collision step with external force +template +class ForcedShanChenBGKdynamics : public ForcedBGKdynamics { +public: + /// Constructor + ForcedShanChenBGKdynamics(T omega, Momenta& momenta); + /// Compute fluid velocity on the cell. + void computeU ( + Cell const& cell, + T u[DESCRIPTOR::d] ) const override; + /// Compute fluid velocity and particle density on the cell. + void computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const override; + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; +}; + +/// Implementation of the TRT collision step with external force +template +class ForcedTRTdynamics : public BasicDynamics { +public: + /// Constructor + ForcedTRTdynamics(T omega, Momenta& momenta, T magicParameter); + /// Compute fluid velocity on the cell. + void computeU ( + Cell const& cell, + T u[DESCRIPTOR::d] ) const override; + /// Compute fluid velocity and particle density on the cell. + void computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const override; + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Get local relaxation parameter of the dynamics + T getOmega() const override; + /// Set local relaxation parameter of the dynamics + void setOmega(T omega) override; +protected: + T _omega; ///< relaxation parameter + T _omega2; /// relaxation parameter for odd moments + T _magicParameter; +}; + +/// Implementation of the 3D D3Q13 dynamics +/** This is (so far) the minimal existing 3D model, with only 13 + * directions. Three different relaxation times are used to achieve + * asymptotic hydrodynamics, isotropy and galilean invariance. + */ +template +class D3Q13dynamics : public BasicDynamics { +public: + /// Constructor + D3Q13dynamics(T omega, Momenta& momenta); + /// Compute equilibrium distribution function + T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const override; + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Get local relaxation parameter of the dynamics + T getOmega() const override; + /// Set local relaxation parameter of the dynamics + void setOmega(T omega) override; +private: + T lambda_nu; ///< first relaxation parameter + T lambda_nu_prime; ///< second relaxation parameter +}; + +/// Standard computation of velocity momenta in the bulk +template +struct BulkMomenta : public Momenta { + /// Compute particle density on the cell. + T computeRho(Cell const& cell) const override; + /// Compute fluid velocity on the cell. + void computeU ( + Cell const& cell, + T u[DESCRIPTOR::d] ) const override; + /// Compute fluid momentum on the cell. + void computeJ ( + Cell const& cell, + T j[DESCRIPTOR::d] ) const override; + /// Compute components of the stress tensor on the cell. + void computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const override; + /// Compute fluid velocity and particle density on the cell. + void computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const override; + /// Compute all momenta on the cell, up to second order. + void computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const override; + /// Set particle density on the cell. + void defineRho(Cell& cell, T rho) override; + /// Set fluid velocity on the cell. + void defineU(Cell& cell, + const T u[DESCRIPTOR::d]) override; + /// Define fluid velocity and particle density on the cell. + void defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) override; + /// Define all momenta on the cell, up to second order. + void defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) override; +}; + +/// Velocity is stored in external scalar (and computed e.g. in a PostProcessor) +template +struct ExternalVelocityMomenta : public Momenta { + /// Compute particle density on the cell. + T computeRho(Cell const& cell) const override; + /// Compute fluid velocity on the cell. + void computeU ( + Cell const& cell, + T u[DESCRIPTOR::d] ) const override; + /// Compute fluid momentum on the cell. + void computeJ ( + Cell const& cell, + T j[DESCRIPTOR::d] ) const override; + /// Compute components of the stress tensor on the cell. + void computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const override; + /// Compute fluid velocity and particle density on the cell. + void computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const override; + /// Compute all momenta on the cell, up to second order. + void computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const override; + /// Set particle density on the cell. + void defineRho(Cell& cell, T rho) override; + /// Set fluid velocity on the cell. + void defineU(Cell& cell, + const T u[DESCRIPTOR::d]) override; + /// Define fluid velocity and particle density on the cell. + void defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) override; + /// Define all momenta on the cell, up to second order. + void defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) override; +}; + +/// Implementation of "bounce-back" dynamics +/** This is a very popular way to implement no-slip boundary conditions, + * because the dynamics are independent of the orientation of the boundary. + * It is a special case, because it implements no usual LB dynamics. + * For that reason, it derives directly from the class Dynamics. + * + * The code works for both 2D and 3D lattices. + */ +template +class BounceBack : public Dynamics { +public: + /// A fictitious density value on bounce-back in not fixed on nodes via this constructor. + BounceBack(); + /// You may fix a fictitious density value on bounce-back nodes via this constructor. + BounceBack(T rho); + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Yields 1; + T computeRho(Cell const& cell) const override; + /// Yields 0; + void computeU ( + Cell const& cell, + T u[DESCRIPTOR::d] ) const override; + /// Yields 0; + void computeJ ( + Cell const& cell, + T j[DESCRIPTOR::d] ) const override; + /// Yields NaN + void computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const override; + void computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const override; + void computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const override; + /// Does nothing + void defineRho(Cell& cell, T rho) override; + /// Does nothing + void defineU(Cell& cell, + const T u[DESCRIPTOR::d]) override; + /// Does nothing + void defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) override; + /// Does nothing + void defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) override; + /// Yields NaN + T getOmega() const override; + /// Does nothing + void setOmega(T omega) override; +private: + T _rho; + bool _rhoFixed; +}; + +/// Implementation of "bounce-back velocity" dynamics +/** This is a very popular way to implement no-slip boundary conditions, + * because the dynamics are independent of the orientation of the boundary. + * It is a special case, because it implements no usual LB dynamics. + * For that reason, it derives directly from the class Dynamics. It + * fixes the velociy to a given velocity _u. + * + * The code works for both 2D and 3D lattices. + */ +template +class BounceBackVelocity : public Dynamics { +public: + /// A fictitious density value on bounce-back in not fixed on nodes via this constructor. + BounceBackVelocity(const T u[DESCRIPTOR::d]); + /// You may fix a fictitious density value on bounce-back nodes via this constructor. + BounceBackVelocity(const T rho, const T u[DESCRIPTOR::d]); + /// Collision step, bounce back with a fixed velocity _u + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Retuns rho (if defined else zero) + T computeRho(Cell const& cell) const override; + /// Retuns _u + void computeU ( + Cell const& cell, + T u[DESCRIPTOR::d] ) const override; + /// Retuns rho (if defined else zero) times _u + void computeJ ( + Cell const& cell, + T j[DESCRIPTOR::d] ) const override; + /// Yields NaN + void computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const override; + /// Retuns rho (if defined else zero) and _u + void computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const override; + void computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const override; + /// Devines the velocity rho + void defineRho(Cell& cell, T rho) override; + /// Devines the velocity _u + void defineU(Cell& cell, + const T u[DESCRIPTOR::d]) override; + /// Devines rho and _u + void defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) override; + /// Does nothing + void defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) override; + /// Yields NaN + T getOmega() const override; + /// Does nothing + void setOmega(T omega) override; +private: + T _rho; + bool _rhoFixed; + T _u[DESCRIPTOR::d]; +}; + +/// Implementation of "bounce-back anti" dynamics +/** This is a way to implement a Dirichlet rho/pressure boundary conditions, + * because the dynamics are independent of the orientation of the boundary. + * It is a special case, because it implements no usual LB dynamics. + * For that reason, it derives directly from the class Dynamics. It + * fixes the rho to a given _rho. + * + * The code works for both 2D and 3D lattices. + */ +template +class BounceBackAnti : public Dynamics { +public: + /// A fictitious density value on bounce-back in not fixed on nodes via this constructor. + BounceBackAnti(); + /// You may fix a fictitious density value on bounce-back nodes via this constructor. + BounceBackAnti(T rho); + /// Collision step, bounce back with a fixed velocity _u + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Retuns rho (if defined else zero) + T computeRho(Cell const& cell) const override; + /// Retuns _u + void computeU ( + Cell const& cell, + T u[DESCRIPTOR::d] ) const override; + /// Retuns rho (if defined else zero) times _u + void computeJ ( + Cell const& cell, + T j[DESCRIPTOR::d] ) const override; + /// Yields NaN + void computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const override; + /// Retuns rho (if defined else zero) and _u + void computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const override; + void computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const override; + /// Devines the velocity rho + void defineRho(Cell& cell, T rho) override; + /// Devines the velocity _u + void defineU(Cell& cell, + const T u[DESCRIPTOR::d]) override; + /// Devines rho and _u + void defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) override; + /// Does nothing + void defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) override; + /// Yields NaN + T getOmega() const override; + /// Does nothing + void setOmega(T omega) override; +private: + T _rho; + bool _rhoFixed; + T _u[DESCRIPTOR::d]; +}; + + +/** Corresponds to macro Robin boundary, micro Fresnel surface + * Motivated by Hiorth et al. 2008; doi 10.1002/fld.1822 + */ +template +class PartialBounceBack final: public BounceBack { +public: + PartialBounceBack(T rf); + T computeEquilibrium( int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr ) const override; + /// Collision step + void collide(Cell& cell, LatticeStatistics& statistics_) override; +private: + T _rf; +}; + + +/// Implementation of a "dead cell" that does nothing +template +class NoDynamics : public Dynamics { +public: + /// You may fix a fictitious density value on no dynamics node via this constructor. + NoDynamics(T rho = T(1) ); + /// Yields 0; + T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const override; + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Yields 1; + T computeRho(Cell const& cell) const override; + /// Yields 0; + void computeU ( + Cell const& cell, + T u[DESCRIPTOR::d] ) const override; + /// Yields 0; + void computeJ ( + Cell const& cell, + T j[DESCRIPTOR::d] ) const override; + /// Yields NaN + void computeStress ( + Cell const& cell, + T rho, const T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const override; + void computeRhoU ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d]) const override; + void computeAllMomenta ( + Cell const& cell, + T& rho, T u[DESCRIPTOR::d], + T pi[util::TensorVal::n] ) const override; + /// Does nothing + void defineRho(Cell& cell, T rho) override; + /// Does nothing + void defineU(Cell& cell, + const T u[DESCRIPTOR::d]) override; + /// Does nothing + void defineRhoU ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d]) override; + /// Does nothing + void defineAllMomenta ( + Cell& cell, + T rho, const T u[DESCRIPTOR::d], + const T pi[util::TensorVal::n] ) override; + /// Yields NaN + T getOmega() const override; + /// Does nothing + void setOmega(T omega) override; + +private: + /// Default rho=1 + T _rho; +}; + +/// Dynamics for offLattice boundary conditions +/// OffDynamics are basically NoDynamics with the additional functionality +/// to store given velocities exactly at boundary links. +template +class OffDynamics : public NoDynamics { +public: + /// Constructor + OffDynamics(const T _location[DESCRIPTOR::d]); + /// Constructor + OffDynamics(const T _location[DESCRIPTOR::d], T _distances[DESCRIPTOR::q]); + /// Returns local stored rho which is updated if the bc is used as velocity!=0 condition + T computeRho(Cell const& cell) const override; + /// Returns an average of the locally stored u + void computeU(Cell const& cell, T u[DESCRIPTOR::d] ) const override; + /// Set Intersection of the link and the boundary + void setBoundaryIntersection(int iPop, T distance) override; + /// Get Intersection of the link and the boundary + bool getBoundaryIntersection(int iPop, T intersection[DESCRIPTOR::d]) override; + /// Set particle density on the cell. + void defineRho(Cell& cell, T rho) override; + /// Set single velocity + void defineRho(int iPop, T rho) override; + /// Set fluid velocity on the cell. + void defineU(Cell& cell, const T u[DESCRIPTOR::d]) override; + /// Set constant velocity + void defineU(const T u[DESCRIPTOR::d]) override; + /// Set single velocity + void defineU(int iPop, const T u[DESCRIPTOR::d]) override; + /// Get VelocitySummand for Bouzidi-Boundary Condition + T getVelocityCoefficient(int iPop) override; + +private: + T _rho; + T _u[DESCRIPTOR::q][DESCRIPTOR::d]; + T location[DESCRIPTOR::d]; + T distances[DESCRIPTOR::q]; + T boundaryIntersection[DESCRIPTOR::q][DESCRIPTOR::d]; + T velocityCoefficient[DESCRIPTOR::q]; +}; + +/// Implementation of density sink by setting a zero distribution on the cell +template +class ZeroDistributionDynamics : public NoDynamics { +public: + /// Constructor. + ZeroDistributionDynamics(); + /// Collision step + void collide(Cell& cell, + LatticeStatistics& statistics_) override; + /// Yields 1 + T computeRho(Cell const& cell) const override; +}; + + +namespace instances { + +template +BulkMomenta& getBulkMomenta(); + +template +ExternalVelocityMomenta& getExternalVelocityMomenta(); + +template +BounceBack& getBounceBack(); + +template +PartialBounceBack& getPartialBounceBack(const double rf); + +template +BounceBackVelocity& getBounceBackVelocity(const double rho, const double u[DESCRIPTOR::d]); + +template +BounceBackAnti& getBounceBackAnti(const double rho); + +template +NoDynamics& getNoDynamics(T rho = T(1) ); + +template +ZeroDistributionDynamics& getZeroDistributionDynamics(); + +} + +} + +#endif -- cgit v1.2.3