dynamics.hh

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 -- generic implementation.
 */
#ifndef LB_DYNAMICS_HH
#define LB_DYNAMICS_HH

#include <algorithm>
#include <limits>
#include "dynamics.h"
#include "cell.h"
#include "lbHelpers.h"
#include "firstOrderLbHelpers.h"
#include "d3q13Helpers.h"

namespace olb {

////////////////////// Class Dynamics ////////////////////////

template<typename T, template<typename U> class Lattice>
void Dynamics<T,Lattice>::computePopulations(Cell<T,Lattice> const& cell, T* f) const {
    for (int iPop=0; iPop<Lattice<T>::q; ++iPop) {
        f[iPop] = cell[iPop];
    }
}

template<typename T, template<typename U> class Lattice>
void Dynamics<T,Lattice>::iniEquilibrium(Cell<T,Lattice>& cell, T rho, const T u[Lattice<T>::d]) {
    T uSqr = util::normSqr<T,Lattice<T>::d>(u);
    for (int iPop=0; iPop<Lattice<T>::q; ++iPop) {
        cell[iPop] = computeEquilibrium(iPop, rho, u, uSqr);
    }
}

template<typename T, template<typename U> class Lattice>
void Dynamics<T,Lattice>::computeExternalField (
        Cell<T,Lattice> const& cell, int pos, int size, T* ext) const {
    OLB_PRECONDITION(pos+size <= Lattice<T>::ExternalField::numScalars);
    T const* externalData = cell.getExternal(pos);
    for (int iExt=0; iExt<size; ++iExt) {
        ext[iExt] = externalData[iExt];
    }
}

template<typename T, template<typename U> class Lattice>
void Dynamics<T,Lattice>::definePopulations(Cell<T,Lattice>& cell, const T* f)
{
    for (int iPop=0; iPop<Lattice<T>::q; ++iPop) {
        cell[iPop] = f[iPop];
    }
}

template<typename T, template<typename U> class Lattice>
void Dynamics<T,Lattice>::defineExternalField (
        Cell<T,Lattice>& cell, int pos, int size, const T* ext)
{
    OLB_PRECONDITION(pos+size <= Lattice<T>::ExternalField::numScalars);
    T* externalData = cell.getExternal(pos);
    for (int iExt=0; iExt<size; ++iExt) {
        externalData[iExt] = ext[iExt];
    }
}

template<typename T, template<typename U> class Lattice>
T Dynamics<T,Lattice>::getParameter(int whichParameter) const {
    if (whichParameter == dynamicParams::omega_shear) {
        return getOmega();
    }
    return 0.;
}

template<typename T, template<typename U> class Lattice>
void Dynamics<T,Lattice>::setParameter(int whichParameter, T value) {
    if (whichParameter == dynamicParams::omega_shear) {
        setOmega(value);
    }
}

////////////////////// Class BasicDynamics ////////////////////////

template<typename T, template<typename U> class Lattice>
BasicDynamics<T,Lattice>::BasicDynamics(Momenta<T,Lattice>& momenta_)
    : momenta(momenta_)
{ }

template<typename T, template<typename U> class Lattice>
T BasicDynamics<T,Lattice>::computeRho(Cell<T,Lattice> const& cell) const
{
    return momenta.computeRho(cell);
}

template<typename T, template<typename U> class Lattice>
void BasicDynamics<T,Lattice>::computeU (
        Cell<T,Lattice> const& cell,
        T u[Lattice<T>::d]) const
{
    momenta.computeU(cell, u);
}

template<typename T, template<typename U> class Lattice>
void BasicDynamics<T,Lattice>::computeJ (
        Cell<T,Lattice> const& cell,
        T j[Lattice<T>::d]) const
{
    momenta.computeJ(cell, j);
}

template<typename T, template<typename U> class Lattice>
void BasicDynamics<T,Lattice>::computeStress (
        Cell<T,Lattice> const& cell,
        T rho, const T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const
{
    momenta.computeStress(cell, rho, u, pi);
}

template<typename T, template<typename U> class Lattice>
void BasicDynamics<T,Lattice>::computeRhoU (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d]) const
{
    momenta.computeRhoU(cell, rho, u);
}

template<typename T, template<typename U> class Lattice>
void BasicDynamics<T,Lattice>::computeAllMomenta (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const
{
    momenta.computeAllMomenta(cell, rho, u, pi);
}

template<typename T, template<typename U> class Lattice>
void BasicDynamics<T,Lattice>::defineRho(Cell<T,Lattice>& cell, T rho) {
    momenta.defineRho(cell, rho);
}

template<typename T, template<typename U> class Lattice>
void BasicDynamics<T,Lattice>::defineU (
        Cell<T,Lattice>& cell,
        const T u[Lattice<T>::d])
{
    momenta.defineU(cell, u);
}

template<typename T, template<typename U> class Lattice>
void BasicDynamics<T,Lattice>::defineRhoU (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::d])
{
    momenta.defineRhoU(cell, rho, u);
}

template<typename T, template<typename U> class Lattice>
void BasicDynamics<T,Lattice>::defineAllMomenta (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::d],
        const T pi[util::TensorVal<Lattice<T> >::n] )
{
    momenta.defineAllMomenta(cell, rho, u, pi);
}


////////////////////// Class BGKdynamics //////////////////////////

/** \param omega_ relaxation parameter, related to the dynamic viscosity
 *  \param momenta_ a Momenta object to know how to compute velocity momenta
 */
template<typename T, template<typename U> class Lattice>
BGKdynamics<T,Lattice>::BGKdynamics (
        T omega_, Momenta<T,Lattice>& momenta_ )
    : BasicDynamics<T,Lattice>(momenta_),
      omega(omega_)
{ }

template<typename T, template<typename U> class Lattice>
BGKdynamics<T,Lattice>* BGKdynamics<T,Lattice>::clone() const {
    return new BGKdynamics<T,Lattice>(*this);
}

template<typename T, template<typename U> class Lattice>
T BGKdynamics<T,Lattice>::computeEquilibrium(int iPop, T rho, const T u[Lattice<T>::d], T uSqr) const
{
    return lbHelpers<T,Lattice>::equilibrium(iPop, rho, u, uSqr);
}

template<typename T, template<typename U> class Lattice>
void BGKdynamics<T,Lattice>::collide (
        Cell<T,Lattice>& cell,
        LatticeStatistics<T>& statistics )
{
    T rho, u[Lattice<T>::d];
    this->momenta.computeRhoU(cell, rho, u);
    T uSqr = lbHelpers<T,Lattice>::bgkCollision(cell, rho, u, omega);
    if (cell.takesStatistics()) {
        statistics.gatherStats(rho, uSqr);
    }
}

template<typename T, template<typename U> class Lattice>
void BGKdynamics<T,Lattice>::staticCollide (
        Cell<T,Lattice>& cell,
        const T u[Lattice<T>::d],
        LatticeStatistics<T>& statistics )
{
    T rho;
    rho = this->momenta.computeRho(cell);
    T uSqr = lbHelpers<T,Lattice>::bgkCollision(cell, rho, u, omega);
    if (cell.takesStatistics()) {
        statistics.gatherStats(rho, uSqr);
    }
}

template<typename T, template<typename U> class Lattice>
T BGKdynamics<T,Lattice>::getOmega() const {
    return omega;
}

template<typename T, template<typename U> class Lattice>
void BGKdynamics<T,Lattice>::setOmega(T omega_) {
    omega = omega_;
}


////////////////////// Class ConstRhoBGKdynamics //////////////////////////

/** \param omega_ relaxation parameter, related to the dynamic viscosity
 *  \param momenta_ a Momenta object to know how to compute velocity momenta
 */
template<typename T, template<typename U> class Lattice>
ConstRhoBGKdynamics<T,Lattice>::ConstRhoBGKdynamics (
        T omega_, Momenta<T,Lattice>& momenta_ )
    : BasicDynamics<T,Lattice>(momenta_),
      omega(omega_)
{ }

template<typename T, template<typename U> class Lattice>
ConstRhoBGKdynamics<T,Lattice>* ConstRhoBGKdynamics<T,Lattice>::clone()
    const
{
    return new ConstRhoBGKdynamics<T,Lattice>(*this);
}

template<typename T, template<typename U> class Lattice>
T ConstRhoBGKdynamics<T,Lattice>::computeEquilibrium(int iPop, T rho, const T u[Lattice<T>::d], T uSqr) const
{
    return lbHelpers<T,Lattice>::equilibrium(iPop, rho, u, uSqr);
}

template<typename T, template<typename U> class Lattice>
void ConstRhoBGKdynamics<T,Lattice>::collide (
        Cell<T,Lattice>& cell,
        LatticeStatistics<T>& statistics )
{
    T rho, u[Lattice<T>::d];
    this->momenta.computeRhoU(cell, rho, u);

    T deltaRho = (T)1 - (statistics).getAverageRho();
    T ratioRho = (T)1 + deltaRho/rho;

    T uSqr = lbHelpers<T,Lattice>::constRhoBgkCollision (
                cell, rho, u, ratioRho, omega );
    if (cell.takesStatistics()) {
        statistics.gatherStats(rho+deltaRho, uSqr);
    }
}

template<typename T, template<typename U> class Lattice>
void ConstRhoBGKdynamics<T,Lattice>::staticCollide (
        Cell<T,Lattice>& cell,
        const T u[Lattice<T>::d],
        LatticeStatistics<T>& statistics )
{
    T rho = this->momenta.computeRho(cell);
    T uSqr = lbHelpers<T,Lattice>::bgkCollision(cell, rho, u, omega);
    if (cell.takesStatistics()) {
        statistics.gatherStats(rho, uSqr);
    }
}

template<typename T, template<typename U> class Lattice>
T ConstRhoBGKdynamics<T,Lattice>::getOmega() const {
    return omega;
}

template<typename T, template<typename U> class Lattice>
void ConstRhoBGKdynamics<T,Lattice>::setOmega(T omega_) {
    omega = omega_;
}

////////////////////// Class IncBGKdynamics //////////////////////////

/** \param omega_ relaxation parameter, related to the dynamic viscosity
 *  \param momenta_ a Momenta object to know how to compute velocity momenta
 */
template<typename T, template<typename U> class Lattice>
IncBGKdynamics<T,Lattice>::IncBGKdynamics (
        T omega_, Momenta<T,Lattice>& momenta_ )
    : BasicDynamics<T,Lattice>(momenta_),
      omega(omega_)
{ }

template<typename T, template<typename U> class Lattice>
IncBGKdynamics<T,Lattice>* IncBGKdynamics<T,Lattice>::clone() const {
    return new IncBGKdynamics<T,Lattice>(*this);
}

template<typename T, template<typename U> class Lattice>
T IncBGKdynamics<T,Lattice>::computeEquilibrium(int iPop, T rho, const T u[Lattice<T>::d], T uSqr) const
{
    return lbHelpers<T,Lattice>::equilibrium(iPop, rho, u, uSqr);
}

template<typename T, template<typename U> class Lattice>
void IncBGKdynamics<T,Lattice>::collide (
        Cell<T,Lattice>& cell,
        LatticeStatistics<T>& statistics )
{
    T rho = this->momenta.computeRho(cell);
    T p = rho / Lattice<T>::invCs2;
    T j[Lattice<T>::d];
    this->momenta.computeJ(cell, j);
    T uSqr = lbHelpers<T,Lattice>::incBgkCollision(cell, p, j, omega);
    if (cell.takesStatistics()) {
        statistics.gatherStats(rho, uSqr);
    }
}

template<typename T, template<typename U> class Lattice>
void IncBGKdynamics<T,Lattice>::staticCollide (
         Cell<T,Lattice>& cell,
         const T j[Lattice<T>::d],
         LatticeStatistics<T>& statistics )
{
    T rho = this->momenta.computeRho(cell);
    T p = rho / Lattice<T>::invCs2;
    T uSqr = lbHelpers<T,Lattice>::incBgkCollision(cell, p, j, omega);
    if (cell.takesStatistics()) {
        statistics.gatherStats(rho, uSqr);
    }
}

template<typename T, template<typename U> class Lattice>
T IncBGKdynamics<T,Lattice>::getOmega() const {
    return omega;
}

template<typename T, template<typename U> class Lattice>
void IncBGKdynamics<T,Lattice>::setOmega(T omega_) {
    omega = omega_;
}



////////////////////// Class RLBdynamics /////////////////////////

/** \param omega_ relaxation parameter, related to the dynamic viscosity
 *  \param momenta_ a Momenta object to know how to compute velocity momenta
 */
template<typename T, template<typename U> class Lattice>
RLBdynamics<T,Lattice>::RLBdynamics (
        T omega_, Momenta<T,Lattice>& momenta_ )
    : BasicDynamics<T,Lattice>(momenta_),
      omega(omega_)
{ }

template<typename T, template<typename U> class Lattice>
RLBdynamics<T,Lattice>* RLBdynamics<T,Lattice>::clone() const {
    return new RLBdynamics<T,Lattice>(*this);
}