dynamics.h

open lattice boltzmann project www.openlb.org

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2006, 2007 Jonas Latt
 *  Address: Rue General Dufour 24,  1211 Geneva 4, Switzerland 
 *  E-mail: [e-mail]
 *
 *  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 "util.h"
#include "postProcessing.h"

namespace olb {

namespace dynamicParams {
    // Use 0-99 for relaxation parameters
    const int omega_shear = 0;
    const int omega_bulk  = 1;

    // Use 100-199 for material constants
    const int sqrSpeedOfSound = 100; // Speed of sound squared

    // Use 1000 and higher for custom user-defined constants
}

template<typename T, template<typename U> class Lattice> class Cell;

/// Interface for the dynamics classes
template<typename T, template<typename U> class Lattice>
struct Dynamics {
    /// Destructor: virtual to enable inheritance
    virtual ~Dynamics() { }
    /// Clone the object on its dynamic type.
    virtual Dynamics<T,Lattice>* clone() const =0;
    /// Implementation of the collision step
    virtual void collide(Cell<T,Lattice>& cell,
                         LatticeStatistics<T>& statistics_) =0;
    /// Collide with fixed velocity
    virtual void staticCollide(Cell<T,Lattice>& cell,
                               const T u[Lattice<T>::d],
                               LatticeStatistics<T>& statistics_) =0;
    /// Compute equilibrium distribution function
    virtual T computeEquilibrium(int iPop, T rho, const T u[Lattice<T>::d], T uSqr) const =0;
    /// Initialize cell at equilibrium distribution
    virtual void iniEquilibrium(Cell<T,Lattice>& cell, T rho, const T u[Lattice<T>::d]);
    /// Compute particle density on the cell.
    /** \return particle density
     */
    virtual T computeRho(Cell<T,Lattice> const& cell) const =0;
    /// Compute fluid velocity on the cell.
    /** \param u fluid velocity
     */
    virtual void computeU( Cell<T,Lattice> const& cell,
                           T u[Lattice<T>::d] ) const =0;
    /// Compute fluid momentum (j=rho*u) on the cell.
    /** \param j fluid momentum
     */
    virtual void computeJ( Cell<T,Lattice> const& cell,
                           T j[Lattice<T>::d] ) const =0;
    /// Compute components of the stress tensor on the cell.
    /** \param pi stress tensor */
    virtual void computeStress (
        Cell<T,Lattice> const& cell, 
        T rho, const T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::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,Lattice> const& cell,
        T& rho, T u[Lattice<T>::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,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const =0;
    /// Access particle populations through the dynamics object.
    /** Default implementation: access cell directly.
     */
    virtual void computePopulations(Cell<T,Lattice> const& cell, T* f) const;
    /// Access external fields through the dynamics object.
    /** Default implementation: access cell directly.
     */
    virtual void computeExternalField (
            Cell<T,Lattice> const& cell, int pos, int size, T* ext ) const;
    /// Set particle density on the cell.
    /** \param rho particle density
     */
    virtual void defineRho(Cell<T,Lattice>& cell, T rho) =0;
    /// Set fluid velocity on the cell.
    /** \param u fluid velocity
     */
    virtual void defineU(Cell<T,Lattice>& cell,
                         const T u[Lattice<T>::d]) =0;
    /// Define fluid velocity and particle density on the cell.
    /** \param rho particle density
     *  \param u fluid velocity
     */
    virtual void defineRhoU (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::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,Lattice>& cell,
        T rho, const T u[Lattice<T>::d],
        const T pi[util::TensorVal<Lattice<T> >::n] ) =0;
    /// Define particle populations through the dynamics object.
    /** Default implementation: access cell directly.
     */
    virtual void definePopulations(Cell<T,Lattice>& cell, const T* f);
    /// Define external fields through the dynamics object.
    /** Default implementation: access cell directly.
     */
    virtual void defineExternalField (
            Cell<T,Lattice>& cell, int pos, int size, const T* ext);
    /// 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;
    /// 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);
};

/// 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, template<typename U> class Lattice>
struct Momenta {
    /// Destructor: virtual to enable inheritance
    virtual ~Momenta() { }
    /// Compute particle density on the cell.
    virtual T computeRho(Cell<T,Lattice> const& cell) const =0;
    /// Compute fluid velocity on the cell.
    virtual void computeU (
        Cell<T,Lattice> const& cell,
        T u[Lattice<T>::d] ) const =0;
    /// Compute fluid momentum on the cell.
    virtual void computeJ (
        Cell<T,Lattice> const& cell,
        T j[Lattice<T>::d] ) const =0;
    /// 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 =0;
    /// 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) =0;
    /// Set fluid velocity on the cell.
    virtual void defineU(Cell<T,Lattice>& cell,
                         const T u[Lattice<T>::d]) =0;
    /// 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] ) =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, template<typename U> class Lattice>
class BasicDynamics : public Dynamics<T,Lattice> {
public:
    /// Must be contructed with an object of type Momenta
    BasicDynamics(Momenta<T,Lattice>& momenta_);
    /// Implemented via the Momenta object
    virtual T computeRho(Cell<T,Lattice> const& cell) const;
    /// Implemented via the Momenta object
    virtual void computeU (
        Cell<T,Lattice> const& cell,
        T u[Lattice<T>::d] ) const;
    /// Implemented via the Momenta object
    virtual void computeJ (
        Cell<T,Lattice> const& cell,
        T j[Lattice<T>::d] ) const;
    /// Implemented via the Momenta object
    virtual void computeStress (
        Cell<T,Lattice> const& cell, 
        T rho, const T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const;
    /// Implemented via the Momenta object
    virtual void computeRhoU (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d]) const;
    /// Implemented via the Momenta object
    virtual void computeAllMomenta (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const;
    /// Implemented via the Momenta object
    virtual void defineRho(Cell<T,Lattice>& cell, T rho);
    /// Implemented via the Momenta object
    virtual void defineU(Cell<T,Lattice>& cell,
                         const T u[Lattice<T>::d]);
    /// Implemented via the Momenta object
    virtual void defineRhoU (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::d]);
    /// Implemented via the Momenta object
    virtual void defineAllMomenta (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::d],
        const T pi[util::TensorVal<Lattice<T> >::n] );
protected:
    Momenta<T,Lattice>& momenta;  ///< computation of velocity momenta
};

/// Implementation of the BGK collision step
template<typename T, template<typename U> class Lattice>
class BGKdynamics : public BasicDynamics<T,Lattice> {
public:
    /// Constructor
    BGKdynamics(T omega_, Momenta<T,Lattice>& momenta_);
    /// Clone the object on its dynamic type.
    virtual BGKdynamics<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 the pressure-corrected BGK collision step
template<typename T, template<typename U> class Lattice>
class ConstRhoBGKdynamics : public BasicDynamics<T,Lattice> {
public:
    /// Constructor
    ConstRhoBGKdynamics(T omega_, Momenta<T,Lattice>& momenta_);
    /// Clone the object on its dynamic type.
    virtual ConstRhoBGKdynamics<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 the so-called incompressible collision step
template<typename T, template<typename U> class Lattice>
class IncBGKdynamics : public BasicDynamics<T,Lattice> {
public:
    /// Constructor
    IncBGKdynamics(T omega_, Momenta<T,Lattice>& momenta_);
    /// Clone the object on its dynamic type.
    virtual IncBGKdynamics<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 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, template<typename U> class Lattice>
class RLBdynamics : public BasicDynamics<T,Lattice> {
public:
    /// Constructor