Initialize at openlb-1-3

1 files changed, 995 insertions, 0 deletions

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
+ *  <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 -- 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<typename T, typename DESCRIPTOR> class Cell;
+
+/// Interface for the dynamics classes
+template<typename T, typename DESCRIPTOR>
+struct Dynamics {
+  /// Destructor: virtual to enable inheritance
+  virtual ~Dynamics() { }
+  /// Implementation of the collision step
+  virtual void collide(Cell<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& 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<T,DESCRIPTOR>& cell, T rho, const T u[DESCRIPTOR::d]);
+  /// Compute particle density on the cell.
+  /** \return particle density
+   */
+  virtual T computeRho(Cell<T,DESCRIPTOR> const& cell) const =0;
+  /// Compute fluid velocity on the cell.
+  /** \param u fluid velocity
+   */
+  virtual void computeU( Cell<T,DESCRIPTOR> 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<T,DESCRIPTOR> 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<T,DESCRIPTOR> const& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const =0;
+  /// Compute fluid velocity and particle density on the cell.
+  /** \param rho particle density
+   *  \param u fluid velocity
+   */
+  virtual void computeRhoU (
+    Cell<T,DESCRIPTOR> 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<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const =0;
+  /// Set particle density on the cell.
+  /** \param rho particle density
+   */
+  virtual void defineRho(Cell<T,DESCRIPTOR>& 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<T,DESCRIPTOR>& 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<T,DESCRIPTOR>& 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<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    const T pi[util::TensorVal<DESCRIPTOR >::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<typename T, typename DESCRIPTOR>
+struct Momenta {
+  /// Destructor: virtual to enable inheritance
+  virtual ~Momenta() { }
+  /// Compute particle density on the cell.
+  virtual T computeRho(Cell<T,DESCRIPTOR> const& cell) const =0;
+  /// Compute fluid velocity on the cell.
+  virtual void computeU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T u[DESCRIPTOR::d] ) const =0;
+  /// Compute fluid momentum on the cell.
+  virtual void computeJ (
+    Cell<T,DESCRIPTOR> const& cell,
+    T j[DESCRIPTOR::d] ) const =0;
+  /// Compute components of the stress tensor on the cell.
+  virtual void computeStress (
+    Cell<T,DESCRIPTOR> const& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const =0;
+  /// Compute fluid velocity and particle density on the cell.
+  virtual void computeRhoU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d]) const;
+  /// Compute all momenta on the cell, up to second order.
+  virtual void computeAllMomenta (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const;
+  /// Set particle density on the cell.
+  virtual void defineRho(Cell<T,DESCRIPTOR>& cell, T rho) =0;
+  /// Set fluid velocity on the cell.
+  virtual void defineU(Cell<T,DESCRIPTOR>& cell,
+                       const T u[DESCRIPTOR::d]) =0;
+  /// Define fluid velocity and particle density on the cell.
+  virtual void defineRhoU (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d]);
+  /// Define all momenta on the cell, up to second order.
+  virtual void defineAllMomenta (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    const T pi[util::TensorVal<DESCRIPTOR >::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<typename T, typename DESCRIPTOR>
+class BasicDynamics : public Dynamics<T,DESCRIPTOR> {
+public:
+  /// Must be contructed with an object of type Momenta
+  BasicDynamics(Momenta<T,DESCRIPTOR>& momenta);
+  /// Implemented via the Momenta object
+  T computeRho(Cell<T,DESCRIPTOR> const& cell) const override;
+  /// Implemented via the Momenta object
+  void computeU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T u[DESCRIPTOR::d] ) const override;
+  /// Implemented via the Momenta object
+  void computeJ (
+    Cell<T,DESCRIPTOR> const& cell,
+    T j[DESCRIPTOR::d] ) const override;
+  /// Implemented via the Momenta object
+  void computeStress (
+    Cell<T,DESCRIPTOR> const& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
+  /// Implemented via the Momenta object
+  void computeRhoU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d]) const override;
+  /// Implemented via the Momenta object
+  void computeAllMomenta (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
+  /// Implemented via the Momenta object
+  void defineRho(Cell<T,DESCRIPTOR>& cell, T rho) override;
+  /// Implemented via the Momenta object
+  void defineU(Cell<T,DESCRIPTOR>& cell,
+                       const T u[DESCRIPTOR::d]) override;
+  /// Implemented via the Momenta object
+  void defineRhoU (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d]) override;
+  /// Implemented via the Momenta object
+  void defineAllMomenta (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    const T pi[util::TensorVal<DESCRIPTOR >::n] ) override;
+protected:
+  Momenta<T,DESCRIPTOR>& _momenta;  ///< computation of velocity momenta
+};
+
+/// Implementation of a generic dynamics to realize a pressure drop at a periodic boundary
+template<typename T, typename DESCRIPTOR, typename BaseDynamics>
+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<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& statistics_) override{
+    BaseDynamics::collide(cell,statistics_);
+    for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+      if ( ((_nx==1 || _nx==-1) && descriptors::c<DESCRIPTOR>(iPop,0)==_nx) || ((_ny==1 || _ny==-1) && descriptors::c<DESCRIPTOR>(iPop,1)==_ny) || (DESCRIPTOR::d==3 && !_nz && descriptors::c<DESCRIPTOR>(iPop,2)==_nz) ) {
+        cell[iPop] += (cell[iPop] + descriptors::t<T,DESCRIPTOR>(iPop))*_densityOffset;
+      }
+    }
+};
+
+private:
+  T _densityOffset;
+  int _nx, _ny, _nz;
+};
+
+/// Implementation of the BGK collision step
+template<typename T, typename DESCRIPTOR>
+class BGKdynamics : public BasicDynamics<T,DESCRIPTOR> {
+public:
+  /// Constructor
+  BGKdynamics(T omega, Momenta<T,DESCRIPTOR>& momenta);
+  /// Collision step
+  void collide(Cell<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& 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<typename T, typename DESCRIPTOR>
+class TRTdynamics : public BasicDynamics<T,DESCRIPTOR> {
+public:
+  /// Constructor
+  TRTdynamics(T omega, Momenta<T,DESCRIPTOR>& momenta, T magicParameter);
+  /// Collision step
+  void collide(Cell<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& 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<typename T, typename DESCRIPTOR>
+class ConstRhoBGKdynamics : public BasicDynamics<T,DESCRIPTOR> {
+public:
+  /// Constructor
+  ConstRhoBGKdynamics(T omega, Momenta<T,DESCRIPTOR>& momenta);
+  /// Collision step
+  void collide(Cell<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& 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<typename T, typename DESCRIPTOR>
+class IncBGKdynamics : public BasicDynamics<T,DESCRIPTOR> {
+public:
+  /// Constructor
+  IncBGKdynamics(T omega, Momenta<T,DESCRIPTOR>& momenta);
+  /// Collision step
+  void collide(Cell<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& 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<typename T, typename DESCRIPTOR>
+class RLBdynamics : public BasicDynamics<T,DESCRIPTOR> {
+public:
+  /// Constructor
+  RLBdynamics(T omega, Momenta<T,DESCRIPTOR>& momenta);
+  /// Collision step
+  void collide(Cell<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& 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<typename T, typename DESCRIPTOR, typename Dynamics>
+class CombinedRLBdynamics : public BasicDynamics<T,DESCRIPTOR> {
+public:
+  /// Constructor
+  CombinedRLBdynamics(T omega, Momenta<T,DESCRIPTOR>& 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<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& 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<typename T, typename DESCRIPTOR>
+class ForcedBGKdynamics : public BasicDynamics<T,DESCRIPTOR> {
+public:
+  /// Constructor
+  ForcedBGKdynamics(T omega, Momenta<T,DESCRIPTOR>& momenta);
+  ///  Compute fluid velocity on the cell.
+  void computeU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T u[DESCRIPTOR::d] ) const override;
+  /// Compute fluid velocity and particle density on the cell.
+  void computeRhoU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d]) const override;
+  /// Collision step
+  void collide(Cell<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& 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<typename T, typename DESCRIPTOR>
+class ForcedKupershtokhBGKdynamics : public BasicDynamics<T,DESCRIPTOR> {
+public:
+  /// Constructor
+  ForcedKupershtokhBGKdynamics(T omega, Momenta<T,DESCRIPTOR>& momenta);
+  ///  Compute fluid velocity on the cell.
+  void computeU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T u[DESCRIPTOR::d] ) const override;
+  /// Compute fluid velocity and particle density on the cell.
+  void computeRhoU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d]) const override;
+  /// Collision step
+  void collide(Cell<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& 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<typename T, typename DESCRIPTOR>
+class PoissonDynamics : public BasicDynamics<T,DESCRIPTOR> {
+public:
+  /// Constructor
+  PoissonDynamics(T omega, Momenta<T,DESCRIPTOR>& 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<T,DESCRIPTOR> const& cell) const override;
+  void computeU(Cell<T,DESCRIPTOR> const& cell, T u[DESCRIPTOR::d]) const override;
+  void computeJ(Cell<T,DESCRIPTOR> const& cell, T j[DESCRIPTOR::d]) const override;
+  void computeStress(Cell<T,DESCRIPTOR> const& cell, T rho,
+                     const T u[DESCRIPTOR::d], T pi[util::TensorVal<DESCRIPTOR>::n] ) const override;
+  void computeRhoU( Cell<T,DESCRIPTOR> const& cell,
+                    T& rho, T u[DESCRIPTOR::d]) const override;
+  void computeAllMomenta( Cell<T,DESCRIPTOR> const& cell, T &rho,
+                          T u[DESCRIPTOR::q], T pi[util::TensorVal<DESCRIPTOR>::n] ) const override;
+  void collide(Cell<T,DESCRIPTOR>& cell,
+               LatticeStatistics<T>& statistics_) override;
+  T getOmega() const override;
+  void setOmega(T omega) override;
+protected:
+  T _omega;
+  T _sink;
+};
+
+
+template<typename T, typename DESCRIPTOR>
+class P1Dynamics : public BasicDynamics<T,DESCRIPTOR> {
+public:
+  /// Constructor
+  P1Dynamics(T omega, Momenta<T,DESCRIPTOR>& 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<T,DESCRIPTOR> const& cell) const override;
+  void computeU(Cell<T,DESCRIPTOR> const& cell, T u[DESCRIPTOR::d]) const override;
+  void computeJ(Cell<T,DESCRIPTOR> const& cell, T j[DESCRIPTOR::d]) const override;
+  void computeStress(Cell<T,DESCRIPTOR> const& cell, T rho,
+                     const T u[DESCRIPTOR::d], T pi[util::TensorVal<DESCRIPTOR>::n] ) const override;
+  void computeRhoU( Cell<T,DESCRIPTOR> const& cell,
+                    T& rho, T u[DESCRIPTOR::d]) const override;
+  void computeAllMomenta( Cell<T,DESCRIPTOR> const& cell, T &rho,
+                          T u[DESCRIPTOR::q], T pi[util::TensorVal<DESCRIPTOR>::n] ) const override;
+  void collide(Cell<T,DESCRIPTOR>& cell,
+               LatticeStatistics<T>& 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<typename T, typename DESCRIPTOR>
+class ResettingForcedBGKdynamics : public ForcedBGKdynamics<T,DESCRIPTOR> {
+public:
+  ResettingForcedBGKdynamics(T omega, Momenta<T,DESCRIPTOR>& momenta);
+  /// Collision step
+  void collide(Cell<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& 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<typename T, typename DESCRIPTOR>
+class ForcedShanChenBGKdynamics : public ForcedBGKdynamics<T,DESCRIPTOR> {
+public:
+  /// Constructor
+  ForcedShanChenBGKdynamics(T omega, Momenta<T,DESCRIPTOR>& momenta);
+  ///  Compute fluid velocity on the cell.
+  void computeU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T u[DESCRIPTOR::d] ) const override;
+  /// Compute fluid velocity and particle density on the cell.
+  void computeRhoU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d]) const override;
+  /// Collision step
+  void collide(Cell<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& statistics_) override;
+};
+
+/// Implementation of the TRT collision step with external force
+template<typename T, typename DESCRIPTOR>
+class ForcedTRTdynamics : public BasicDynamics<T,DESCRIPTOR> {
+public:
+  /// Constructor
+  ForcedTRTdynamics(T omega, Momenta<T,DESCRIPTOR>& momenta, T magicParameter);
+  ///  Compute fluid velocity on the cell.
+  void computeU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T u[DESCRIPTOR::d] ) const override;
+  /// Compute fluid velocity and particle density on the cell.
+  void computeRhoU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d]) const override;
+  /// Collision step
+  void collide(Cell<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& 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<typename T, typename DESCRIPTOR>
+class D3Q13dynamics : public BasicDynamics<T,DESCRIPTOR> {
+public:
+  /// Constructor
+  D3Q13dynamics(T omega, Momenta<T,DESCRIPTOR>& 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<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& 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<typename T, typename DESCRIPTOR>
+struct BulkMomenta : public Momenta<T,DESCRIPTOR> {
+  /// Compute particle density on the cell.
+  T computeRho(Cell<T,DESCRIPTOR> const& cell) const override;
+  /// Compute fluid velocity on the cell.
+  void computeU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T u[DESCRIPTOR::d] ) const override;
+  /// Compute fluid momentum on the cell.
+  void computeJ (
+    Cell<T,DESCRIPTOR> const& cell,
+    T j[DESCRIPTOR::d] ) const override;
+  /// Compute components of the stress tensor on the cell.
+  void computeStress (
+    Cell<T,DESCRIPTOR> const& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
+  /// Compute fluid velocity and particle density on the cell.
+  void computeRhoU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d]) const override;
+  /// Compute all momenta on the cell, up to second order.
+  void computeAllMomenta (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
+  /// Set particle density on the cell.
+  void defineRho(Cell<T,DESCRIPTOR>& cell, T rho) override;
+  /// Set fluid velocity on the cell.
+  void defineU(Cell<T,DESCRIPTOR>& cell,
+                       const T u[DESCRIPTOR::d]) override;
+  /// Define fluid velocity and particle density on the cell.
+  void defineRhoU (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d]) override;
+  /// Define all momenta on the cell, up to second order.
+  void defineAllMomenta (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    const T pi[util::TensorVal<DESCRIPTOR >::n] ) override;
+};
+
+/// Velocity is stored in external scalar (and computed e.g. in a PostProcessor)
+template<typename T, typename DESCRIPTOR>
+struct ExternalVelocityMomenta : public Momenta<T,DESCRIPTOR> {
+  /// Compute particle density on the cell.
+  T computeRho(Cell<T,DESCRIPTOR> const& cell) const override;
+  /// Compute fluid velocity on the cell.
+  void computeU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T u[DESCRIPTOR::d] ) const override;
+  /// Compute fluid momentum on the cell.
+  void computeJ (
+    Cell<T,DESCRIPTOR> const& cell,
+    T j[DESCRIPTOR::d] ) const override;
+  /// Compute components of the stress tensor on the cell.
+  void computeStress (
+    Cell<T,DESCRIPTOR> const& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
+  /// Compute fluid velocity and particle density on the cell.
+  void computeRhoU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d]) const override;
+  /// Compute all momenta on the cell, up to second order.
+  void computeAllMomenta (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
+  /// Set particle density on the cell.
+  void defineRho(Cell<T,DESCRIPTOR>& cell, T rho) override;
+  /// Set fluid velocity on the cell.
+  void defineU(Cell<T,DESCRIPTOR>& cell,
+                       const T u[DESCRIPTOR::d]) override;
+  /// Define fluid velocity and particle density on the cell.
+  void defineRhoU (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d]) override;
+  /// Define all momenta on the cell, up to second order.
+  void defineAllMomenta (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    const T pi[util::TensorVal<DESCRIPTOR >::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<typename T, typename DESCRIPTOR>
+class BounceBack : public Dynamics<T,DESCRIPTOR> {
+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<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& statistics_) override;
+  /// Yields 1;
+  T computeRho(Cell<T,DESCRIPTOR> const& cell) const override;
+  /// Yields 0;
+  void computeU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T u[DESCRIPTOR::d] ) const override;
+  /// Yields 0;
+  void computeJ (
+    Cell<T,DESCRIPTOR> const& cell,
+    T j[DESCRIPTOR::d] ) const override;
+  /// Yields NaN
+  void computeStress (
+    Cell<T,DESCRIPTOR> const& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
+  void computeRhoU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d]) const override;
+  void computeAllMomenta (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
+  /// Does nothing
+  void defineRho(Cell<T,DESCRIPTOR>& cell, T rho) override;
+  /// Does nothing
+  void defineU(Cell<T,DESCRIPTOR>& cell,
+                       const T u[DESCRIPTOR::d]) override;
+  /// Does nothing
+  void defineRhoU (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d]) override;
+  /// Does nothing
+  void defineAllMomenta (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    const T pi[util::TensorVal<DESCRIPTOR >::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<typename T, typename DESCRIPTOR>
+class BounceBackVelocity : public Dynamics<T,DESCRIPTOR> {
+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<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& statistics_) override;
+  /// Retuns rho (if defined else zero)
+  T computeRho(Cell<T,DESCRIPTOR> const& cell) const override;
+  /// Retuns _u
+  void computeU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T u[DESCRIPTOR::d] ) const override;
+  /// Retuns rho (if defined else zero) times _u
+  void computeJ (
+    Cell<T,DESCRIPTOR> const& cell,
+    T j[DESCRIPTOR::d] ) const override;
+  /// Yields NaN
+  void computeStress (
+    Cell<T,DESCRIPTOR> const& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
+  /// Retuns rho (if defined else zero) and _u
+  void computeRhoU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d]) const override;
+  void computeAllMomenta (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
+  /// Devines the velocity rho
+  void defineRho(Cell<T,DESCRIPTOR>& cell, T rho) override;
+  /// Devines the velocity _u
+  void defineU(Cell<T,DESCRIPTOR>& cell,
+                       const T u[DESCRIPTOR::d]) override;
+  /// Devines rho and _u
+  void defineRhoU (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d]) override;
+  /// Does nothing
+  void defineAllMomenta (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    const T pi[util::TensorVal<DESCRIPTOR >::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<typename T, typename DESCRIPTOR>
+class BounceBackAnti : public Dynamics<T,DESCRIPTOR> {
+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<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& statistics_) override;
+  /// Retuns rho (if defined else zero)
+  T computeRho(Cell<T,DESCRIPTOR> const& cell) const override;
+  /// Retuns _u
+  void computeU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T u[DESCRIPTOR::d] ) const override;
+  /// Retuns rho (if defined else zero) times _u
+  void computeJ (
+    Cell<T,DESCRIPTOR> const& cell,
+    T j[DESCRIPTOR::d] ) const override;
+  /// Yields NaN
+  void computeStress (
+    Cell<T,DESCRIPTOR> const& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
+  /// Retuns rho (if defined else zero) and _u
+  void computeRhoU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d]) const override;
+  void computeAllMomenta (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
+  /// Devines the velocity rho
+  void defineRho(Cell<T,DESCRIPTOR>& cell, T rho) override;
+  /// Devines the velocity _u
+  void defineU(Cell<T,DESCRIPTOR>& cell,
+                       const T u[DESCRIPTOR::d]) override;
+  /// Devines rho and _u
+  void defineRhoU (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d]) override;
+  /// Does nothing
+  void defineAllMomenta (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    const T pi[util::TensorVal<DESCRIPTOR >::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<typename T, typename DESCRIPTOR>
+class PartialBounceBack final: public BounceBack<T,DESCRIPTOR> {
+public:
+  PartialBounceBack(T rf);
+  T computeEquilibrium( int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr ) const override;
+  /// Collision step
+  void collide(Cell<T,DESCRIPTOR>& cell, LatticeStatistics<T>& statistics_) override;
+private:
+  T _rf;
+};
+
+
+/// Implementation of a "dead cell" that does nothing
+template<typename T, typename DESCRIPTOR>
+class NoDynamics : public Dynamics<T,DESCRIPTOR> {
+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<T,DESCRIPTOR>& cell,
+                       LatticeStatistics<T>& statistics_) override;
+  /// Yields 1;
+  T computeRho(Cell<T,DESCRIPTOR> const& cell) const override;
+  /// Yields 0;
+  void computeU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T u[DESCRIPTOR::d] ) const override;
+  /// Yields 0;
+  void computeJ (
+    Cell<T,DESCRIPTOR> const& cell,
+    T j[DESCRIPTOR::d] ) const override;
+  /// Yields NaN
+  void computeStress (
+    Cell<T,DESCRIPTOR> const& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
+  void computeRhoU (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d]) const override;
+  void computeAllMomenta (
+    Cell<T,DESCRIPTOR> const& cell,
+    T& rho, T u[DESCRIPTOR::d],
+    T pi[util::TensorVal<DESCRIPTOR >::n] ) const override;
+  /// Does nothing
+  void defineRho(Cell<T,DESCRIPTOR>& cell, T rho) override;
+  /// Does nothing
+  void defineU(Cell<T,DESCRIPTOR>& cell,
+                       const T u[DESCRIPTOR::d]) override;
+  /// Does nothing
+  void defineRhoU (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d]) override;
+  /// Does nothing
+  void defineAllMomenta (
+    Cell<T,DESCRIPTOR>& cell,
+    T rho, const T u[DESCRIPTOR::d],
+    const T pi[util::TensorVal<DESCRIPTOR >::n] ) override;
+  /// Yields NaN
+  T getOmega() const override;
+  /// Does nothing
+  void setOmega(T omega) override;