dynamics.h

open lattice boltzmann project www.openlb.org

    RLBdynamics(T omega_, Momenta<T,Lattice>& momenta_);
    /// Clone the object on its dynamic type.
    virtual RLBdynamics<T,Lattice>* clone() const;
    /// Compute equilibrium distribution function
    virtual T computeEquilibrium(int iPop, T rho, const T u[Lattice<T>::d], T uSqr) const;
    /// Collision step
    virtual void collide(Cell<T,Lattice>& cell,
                         LatticeStatistics<T>& statistics_);
    /// Collide with fixed velocity
    virtual void staticCollide(Cell<T,Lattice>& cell,
                               const T u[Lattice<T>::d],
                               LatticeStatistics<T>& statistics_);
    /// Get local relaxation parameter of the dynamics
    virtual T getOmega() const;
    /// Set local relaxation parameter of the dynamics
    virtual void setOmega(T omega_);
private:
    T omega;  ///< relaxation parameter
};

/// Implementation of Regularized BGK collision, followed by any Dynamics
template<typename T, template<typename U> class Lattice, typename Dynamics>
class CombinedRLBdynamics : public BasicDynamics<T,Lattice>
{
public:
    /// Constructor
    CombinedRLBdynamics(T omega_, Momenta<T,Lattice>& momenta_);
    /// Clone the object on its dynamic type.
    virtual CombinedRLBdynamics<T, Lattice, Dynamics>* clone() const;
    /// Compute equilibrium distribution function
    virtual T computeEquilibrium(int iPop, T rho, const T u[Lattice<T>::d], T uSqr) const;
    /// Collision step
    virtual void collide(Cell<T,Lattice>& cell,
                         LatticeStatistics<T>& statistics_);
    /// Collide with fixed velocity
    virtual void staticCollide(Cell<T,Lattice>& cell,
                               const T u[Lattice<T>::d],
                               LatticeStatistics<T>& statistics_);
    /// Get local relaxation parameter of the dynamics
    virtual T getOmega() const;
    /// Set local relaxation parameter of the dynamics
    virtual void setOmega(T omega_);
    /// Get local value of any parameter
    virtual T getParameter(int whichParameter) const;
    /// Set local value of any parameter
    virtual void setParameter(int whichParameter, T value);
private:
    Dynamics boundaryDynamics;
};

/// Implementation of the BGK collision step with external force
template<typename T, template<typename U> class Lattice>
class ForcedBGKdynamics : public BasicDynamics<T,Lattice> {
public:
    /// Constructor
    ForcedBGKdynamics(T omega_, Momenta<T,Lattice>& momenta_);
    /// Clone the object on its dynamic type.
    virtual ForcedBGKdynamics<T,Lattice>* clone() const;
    /// Compute equilibrium distribution function
    virtual T computeEquilibrium(int iPop, T rho, const T u[Lattice<T>::d], T uSqr) const;
    /// Collision step
    virtual void collide(Cell<T,Lattice>& cell,
                         LatticeStatistics<T>& statistics_);
    /// Collide with fixed velocity
    virtual void staticCollide(Cell<T,Lattice>& cell,
                               const T u[Lattice<T>::d],
                               LatticeStatistics<T>& statistics_);
    /// Get local relaxation parameter of the dynamics
    virtual T getOmega() const;
    /// Set local relaxation parameter of the dynamics
    virtual void setOmega(T omega_);
private:
    T omega;  ///< relaxation parameter
    static const int forceBeginsAt = Lattice<T>::ExternalField::forceBeginsAt;
    static const int sizeOfForce   = Lattice<T>::ExternalField::sizeOfForce;
};



/// 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, template<typename U> class Lattice>
class D3Q13dynamics : public BasicDynamics<T,Lattice> {
public:
    /// Constructor
    D3Q13dynamics(T omega_, Momenta<T,Lattice>& momenta_);
    /// Clone the object on its dynamic type.
    virtual D3Q13dynamics<T,Lattice>* clone() const;
    /// Compute equilibrium distribution function
    virtual T computeEquilibrium(int iPop, T rho, const T u[Lattice<T>::d], T uSqr) const;
    /// Collision step
    virtual void collide(Cell<T,Lattice>& cell,
                         LatticeStatistics<T>& statistics_);
    /// Collide with fixed velocity
    virtual void staticCollide(Cell<T,Lattice>& cell,
                               const T u[Lattice<T>::d],
                               LatticeStatistics<T>& statistics_);
    /// Get local relaxation parameter of the dynamics
    virtual T getOmega() const;
    /// Set local relaxation parameter of the dynamics
    virtual void setOmega(T omega_);
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, template<typename U> class Lattice>
struct BulkMomenta : public Momenta<T,Lattice> {
    /// Compute particle density on the cell.
    virtual T computeRho(Cell<T,Lattice> const& cell) const;
    /// Compute fluid velocity on the cell.
    virtual void computeU (
        Cell<T,Lattice> const& cell,
        T u[Lattice<T>::d] ) const;
    /// Compute fluid momentum on the cell.
    virtual void computeJ (
        Cell<T,Lattice> const& cell,
        T j[Lattice<T>::d] ) const;
    /// Compute components of the stress tensor on the cell.
    virtual void computeStress (
        Cell<T,Lattice> const& cell, 
        T rho, const T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const;
    /// Compute fluid velocity and particle density on the cell.
    virtual void computeRhoU (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d]) const;
    /// Compute all momenta on the cell, up to second order.
    virtual void computeAllMomenta (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const;
    /// Set particle density on the cell.
    virtual void defineRho(Cell<T,Lattice>& cell, T rho);
    /// Set fluid velocity on the cell.
    virtual void defineU(Cell<T,Lattice>& cell,
                         const T u[Lattice<T>::d]);
    /// Define fluid velocity and particle density on the cell.
    virtual void defineRhoU (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::d]);
    /// Define all momenta on the cell, up to second order.
    virtual void defineAllMomenta (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::d],
        const T pi[util::TensorVal<Lattice<T> >::n] );
};

/// Velocity is stored in external scalar (and computed e.g. in a PostProcessor)
template<typename T, template<typename U> class Lattice>
struct ExternalVelocityMomenta : public Momenta<T,Lattice> {
    /// Compute particle density on the cell.
    virtual T computeRho(Cell<T,Lattice> const& cell) const;
    /// Compute fluid velocity on the cell.
    virtual void computeU (
        Cell<T,Lattice> const& cell,
        T u[Lattice<T>::d] ) const;
    /// Compute fluid momentum on the cell.
    virtual void computeJ (
        Cell<T,Lattice> const& cell,
        T j[Lattice<T>::d] ) const;
    /// Compute components of the stress tensor on the cell.
    virtual void computeStress (
        Cell<T,Lattice> const& cell, 
        T rho, const T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const;
    /// Compute fluid velocity and particle density on the cell.
    virtual void computeRhoU (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d]) const;
    /// Compute all momenta on the cell, up to second order.
    virtual void computeAllMomenta (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const;
    /// Set particle density on the cell.
    virtual void defineRho(Cell<T,Lattice>& cell, T rho);
    /// Set fluid velocity on the cell.
    virtual void defineU(Cell<T,Lattice>& cell,
                         const T u[Lattice<T>::d]);
    /// Define fluid velocity and particle density on the cell.
    virtual void defineRhoU (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::d]);
    /// Define all momenta on the cell, up to second order.
    virtual void defineAllMomenta (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::d],
        const T pi[util::TensorVal<Lattice<T> >::n] );
};

/// 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, template<typename U> class Lattice>
class BounceBack : public Dynamics<T,Lattice> {
public:
    /// You may fix a fictitious density value on bounce-back nodes via the constructor.
    BounceBack(T rho_=T());
    /// Clone the object on its dynamic type.
    virtual BounceBack<T,Lattice>* clone() const;
    /// Yields 0;
    virtual T computeEquilibrium(int iPop, T rho, const T u[Lattice<T>::d], T uSqr) const;
    /// Collision step
    virtual void collide(Cell<T,Lattice>& cell,
                         LatticeStatistics<T>& statistics_);
    /// Collide with fixed velocity
    virtual void staticCollide(Cell<T,Lattice>& cell,
                               const T u[Lattice<T>::d],
                               LatticeStatistics<T>& statistics_);
    /// Yields 1;
    virtual T computeRho(Cell<T,Lattice> const& cell) const;
    /// Yields 0;
    virtual void computeU (
        Cell<T,Lattice> const& cell,
        T u[Lattice<T>::d] ) const;
    /// Yields 0;
    virtual void computeJ (
        Cell<T,Lattice> const& cell,
        T j[Lattice<T>::d] ) const;
    /// Yields NaN
    virtual void computeStress (
        Cell<T,Lattice> const& cell, 
        T rho, const T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const;
    virtual void computeRhoU (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d]) const;
    virtual void computeAllMomenta (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const;
    /// Does nothing
    virtual void defineRho(Cell<T,Lattice>& cell, T rho);
    /// Does nothing
    virtual void defineU(Cell<T,Lattice>& cell,
                         const T u[Lattice<T>::d]);
    /// Does nothing
    virtual void defineRhoU (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::d]);
    /// Does nothing
    virtual void defineAllMomenta (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::d],
        const T pi[util::TensorVal<Lattice<T> >::n] );
    /// Yields NaN
    virtual T getOmega() const;
    /// Does nothing
    virtual void setOmega(T omega_);
private:
    T rho;
};

/// Implementation of a "dead cell" that does nothing
template<typename T, template<typename U> class Lattice>
class NoDynamics : public Dynamics<T,Lattice> {
public:
    /// Clone the object on its dynamic type.
    virtual NoDynamics<T,Lattice>* clone() const;
    /// Yields 0;
    virtual T computeEquilibrium(int iPop, T rho, const T u[Lattice<T>::d], T uSqr) const;
    /// Collision step
    virtual void collide(Cell<T,Lattice>& cell,
                         LatticeStatistics<T>& statistics_);
    /// Collide with fixed velocity
    virtual void staticCollide(Cell<T,Lattice>& cell,
                               const T u[Lattice<T>::d],
                               LatticeStatistics<T>& statistics_);
    /// Yields 1;
    virtual T computeRho(Cell<T,Lattice> const& cell) const;
    /// Yields 0;
    virtual void computeU (
        Cell<T,Lattice> const& cell,
        T u[Lattice<T>::d] ) const;
    /// Yields 0;
    virtual void computeJ (
        Cell<T,Lattice> const& cell,
        T j[Lattice<T>::d] ) const;
    /// Yields NaN
    virtual void computeStress (
        Cell<T,Lattice> const& cell, 
        T rho, const T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const;
    virtual void computeRhoU (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d]) const;
    virtual void computeAllMomenta (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const;
    /// Does nothing
    virtual void defineRho(Cell<T,Lattice>& cell, T rho);
    /// Does nothing
    virtual void defineU(Cell<T,Lattice>& cell,
                         const T u[Lattice<T>::d]);
    /// Does nothing
    virtual void defineRhoU (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::d]);
    /// Does nothing
    virtual void defineAllMomenta (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::d],
        const T pi[util::TensorVal<Lattice<T> >::n] );
    /// Yields NaN
    virtual T getOmega() const;
    /// Does nothing
    virtual void setOmega(T omega_);
};

namespace instances {

    template<typename T, template<typename U> class Lattice>
    BulkMomenta<T,Lattice>& getBulkMomenta();

    template<typename T, template<typename U> class Lattice>
    ExternalVelocityMomenta<T,Lattice>& getExternalVelocityMomenta();

    template<typename T, template<typename U> class Lattice>
    BounceBack<T,Lattice>& getBounceBack();

    template<typename T, template<typename U> class Lattice>
    NoDynamics<T,Lattice>& getNoDynamics();

}

}

#endif