dynamics.hh

open lattice boltzmann project www.openlb.org

        cell[iPop] = lbHelpers<T,Lattice>::equilibrium(iPop, rho, u, uSqr) +
                     fNeq[iPop];
    }
}

template<typename T, template<typename U> class Lattice>
void BulkMomenta<T,Lattice>::defineAllMomenta (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::d],
        const T pi[util::TensorVal<Lattice<T> >::n] )
{
    T uSqr = util::normSqr<T,Lattice<T>::d>(u);
    for (int iPop=0; iPop < Lattice<T>::q; ++iPop) {
        cell[iPop] = lbHelpers<T,Lattice>::equilibrium(iPop, rho, u, uSqr) +
                     firstOrderLbHelpers<T,Lattice>::fromPiToFneq(iPop, pi);
    }
}

////////////////////// Class ExternalVelocityMomenta //////////////////////////

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

template<typename T, template<typename U> class Lattice>
void ExternalVelocityMomenta<T,Lattice>::computeU(Cell<T,Lattice> const& cell, T u[Lattice<T>::d]) const
{
    T const* uExt = cell.getExternal(Lattice<T>::ExternalField::velocityBeginsAt);
    for (int iD=0; iD<Lattice<T>::d; ++iD) {
        u[iD] = uExt[iD];
    }
}

template<typename T, template<typename U> class Lattice>
void ExternalVelocityMomenta<T,Lattice>::computeJ(Cell<T,Lattice> const& cell, T j[Lattice<T>::d]) const
{
    T rho = computeRho(cell);
    T const* uExt = cell.getExternal(Lattice<T>::ExternalField::velocityBeginsAt);
    for (int iD=0; iD<Lattice<T>::d; ++iD) {
        j[iD] = uExt[iD]*rho;
    }
}

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

template<typename T, template<typename U> class Lattice>
void ExternalVelocityMomenta<T,Lattice>::computeRhoU (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d] ) const
{
    rho = computeRho(cell);
    computeU(cell,u);
}

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

template<typename T, template<typename U> class Lattice>
void ExternalVelocityMomenta<T,Lattice>::defineRho(Cell<T,Lattice>& cell, T rho) {
    T oldRho, u[Lattice<T>::d];
    computeRhoU(cell, oldRho, u);
    T uSqr = util::normSqr<T,Lattice<T>::d>(u);
    T fNeq[Lattice<T>::q];
    lbHelpers<T,Lattice>::computeFneq(cell, fNeq, oldRho, u);
    for (int iPop=0; iPop < Lattice<T>::q; ++iPop) {
        cell[iPop] = lbHelpers<T,Lattice>::equilibrium(iPop, rho, u, uSqr) +
                     fNeq[iPop];
    }
}

template<typename T, template<typename U> class Lattice>
void ExternalVelocityMomenta<T,Lattice>::defineU (
    Cell<T,Lattice>& cell,
    const T u[Lattice<T>::d])
{
    T* uExt = cell.getExternal(Lattice<T>::ExternalField::velocityBeginsAt);
    for (int iD=0; iD<Lattice<T>::d; ++iD) {
        uExt[iD] = u[iD];
    }
}

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

template<typename T, template<typename U> class Lattice>
void ExternalVelocityMomenta<T,Lattice>::defineAllMomenta (
        Cell<T,Lattice>& cell,
        T rho, const T u[Lattice<T>::d],
        const T pi[util::TensorVal<Lattice<T> >::n] )
{
    defineU(cell, u);
    T uSqr = util::normSqr<T,Lattice<T>::d>(u);
    for (int iPop=0; iPop < Lattice<T>::q; ++iPop) {
        cell[iPop] = lbHelpers<T,Lattice>::equilibrium(iPop, rho, u, uSqr) +
                     firstOrderLbHelpers<T,Lattice>::fromPiToFneq(iPop, pi);
    }
}


////////////////////// Class BounceBack ///////////////////////////

template<typename T, template<typename U> class Lattice>
BounceBack<T,Lattice>::BounceBack(T rho_)
    :rho(rho_)
{ }

template<typename T, template<typename U> class Lattice>
BounceBack<T,Lattice>* BounceBack<T,Lattice>::clone() const {
    return new BounceBack<T,Lattice>();
}
 
template<typename T, template<typename U> class Lattice>
T BounceBack<T,Lattice>::computeEquilibrium(int iPop, T rho, const T u[Lattice<T>::d], T uSqr) const
{
    return T();
}

template<typename T, template<typename U> class Lattice>
void BounceBack<T,Lattice>::collide (
        Cell<T,Lattice>& cell,
        LatticeStatistics<T>& statistics )
{
    for (int iPop=1; iPop <= Lattice<T>::q/2; ++iPop) {
        std::swap(cell[iPop], cell[iPop+Lattice<T>::q/2]);
    }
}

template<typename T, template<typename U> class Lattice>
void BounceBack<T,Lattice>::staticCollide (
    Cell<T,Lattice>& cell,
    const T u[Lattice<T>::d],
    LatticeStatistics<T>& statistics )
{
    this->collide(cell, statistics);
}

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

template<typename T, template<typename U> class Lattice>
void BounceBack<T,Lattice>::computeU (
        Cell<T,Lattice> const& cell,
        T u[Lattice<T>::d]) const
{
    for (int iD=0; iD<Lattice<T>::d; ++iD) {
        u[iD] = T();
    }
}

template<typename T, template<typename U> class Lattice>
void BounceBack<T,Lattice>::computeJ (
        Cell<T,Lattice> const& cell,
        T j[Lattice<T>::d]) const
{
    for (int iD=0; iD<Lattice<T>::d; ++iD) {
        j[iD] = T();
    }
}

template<typename T, template<typename U> class Lattice>
void BounceBack<T,Lattice>::computeStress (
        Cell<T,Lattice> const& cell,
        T rho, const T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const
{
    for (int iPi=0; iPi<util::TensorVal<Lattice<T> >::n; ++iPi) {
        pi[iPi] = std::numeric_limits<T>::signaling_NaN();
    }
}

template<typename T, template<typename U> class Lattice>
void BounceBack<T,Lattice>::computeRhoU (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d]) const
{
    rho = computeRho(cell);
    computeU(cell, u);
}

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

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

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

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

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

template<typename T, template<typename U> class Lattice>
T BounceBack<T,Lattice>::getOmega() const {
    return std::numeric_limits<T>::signaling_NaN();
}

template<typename T, template<typename U> class Lattice>
void BounceBack<T,Lattice>::setOmega(T omega_)
{ }

////////////////////// Class NoDynamics ///////////////////////////

template<typename T, template<typename U> class Lattice>
NoDynamics<T,Lattice>* NoDynamics<T,Lattice>::clone() const {
    return new NoDynamics<T,Lattice>();
}
 
template<typename T, template<typename U> class Lattice>
T NoDynamics<T,Lattice>::computeEquilibrium(int iPop, T rho, const T u[Lattice<T>::d], T uSqr) const
{
    return T();
}

template<typename T, template<typename U> class Lattice>
void NoDynamics<T,Lattice>::collide (
        Cell<T,Lattice>& cell,
        LatticeStatistics<T>& statistics )
{ }

template<typename T, template<typename U> class Lattice>
void NoDynamics<T,Lattice>::staticCollide (
    Cell<T,Lattice>& cell,
    const T u[Lattice<T>::d],
    LatticeStatistics<T>& statistics )
{ }

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

template<typename T, template<typename U> class Lattice>
void NoDynamics<T,Lattice>::computeU (
        Cell<T,Lattice> const& cell,
        T u[Lattice<T>::d]) const
{
    for (int iD=0; iD<Lattice<T>::d; ++iD) {
        u[iD] = T();
    }
}

template<typename T, template<typename U> class Lattice>
void NoDynamics<T,Lattice>::computeJ (
        Cell<T,Lattice> const& cell,
        T j[Lattice<T>::d]) const
{
    for (int iD=0; iD<Lattice<T>::d; ++iD) {
        j[iD] = T();
    }
}

template<typename T, template<typename U> class Lattice>
void NoDynamics<T,Lattice>::computeStress (
        Cell<T,Lattice> const& cell,
        T rho, const T u[Lattice<T>::d],
        T pi[util::TensorVal<Lattice<T> >::n] ) const
{
    for (int iPi=0; iPi<util::TensorVal<Lattice<T> >::n; ++iPi) {
        pi[iPi] = std::numeric_limits<T>::signaling_NaN();
    }
}

template<typename T, template<typename U> class Lattice>
void NoDynamics<T,Lattice>::computeRhoU (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d]) const
{
    rho = computeRho(cell);
    computeU(cell, u);
}

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

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

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

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

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

template<typename T, template<typename U> class Lattice>
T NoDynamics<T,Lattice>::getOmega() const {
    return std::numeric_limits<T>::signaling_NaN();
}

template<typename T, template<typename U> class Lattice>
void NoDynamics<T,Lattice>::setOmega(T omega_)
{ }


/////////////// Singletons //////////////////////////////////

namespace instances {

    template<typename T, template<typename U> class Lattice>
    BulkMomenta<T,Lattice>& getBulkMomenta() {
        static BulkMomenta<T,Lattice> bulkMomentaSingleton;
        return bulkMomentaSingleton;
    }

    template<typename T, template<typename U> class Lattice>
    ExternalVelocityMomenta<T,Lattice>& getExternalVelocityMomenta() {
        static ExternalVelocityMomenta<T,Lattice> externalVelocityMomentaSingleton;
        return externalVelocityMomentaSingleton;
    }

    template<typename T, template<typename U> class Lattice>
    BounceBack<T,Lattice>& getBounceBack() {
        static BounceBack<T,Lattice> bounceBackSingleton;
        return bounceBackSingleton;
    }

    template<typename T, template<typename U> class Lattice>
    NoDynamics<T,Lattice>& getNoDynamics() {
        static NoDynamics<T,Lattice> noDynamicsSingleton;
        return noDynamicsSingleton;
    }

}

}

#endif