dynamics.hh

open lattice boltzmann project www.openlb.org

template<typename T, template<typename U> class Lattice>
T RLBdynamics<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 RLBdynamics<T,Lattice>::collide (
        Cell<T,Lattice>& cell,
        LatticeStatistics<T>& statistics )
{
    T rho, u[Lattice<T>::d], pi[util::TensorVal<Lattice<T> >::n];
    this->momenta.computeAllMomenta(cell, rho, u, pi);
    T uSqr = rlbHelpers<T,Lattice>::rlbCollision(cell, rho, u, pi, omega);
    if (cell.takesStatistics()) {
        statistics.gatherStats(rho, uSqr);
    }
}

template<typename T, template<typename U> class Lattice>
void RLBdynamics<T,Lattice>::staticCollide (
        Cell<T,Lattice>& cell,
        const T u[Lattice<T>::d],
        LatticeStatistics<T>& statistics )
{
    T rho, uDummy[Lattice<T>::d], pi[util::TensorVal<Lattice<T> >::n];
    this->momenta.computeAllMomenta(cell, rho, uDummy, pi);
    T uSqr = rlbHelpers<T,Lattice>::rlbCollision(cell, rho, u, pi, omega);
    if (cell.takesStatistics()) {
        statistics.gatherStats(rho, uSqr);
    }
}

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

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


////////////////////// Class CombinedRLBdynamics /////////////////////////

template<typename T, template<typename U> class Lattice, typename Dynamics>
CombinedRLBdynamics<T,Lattice,Dynamics>::CombinedRLBdynamics (
        T omega_, Momenta<T,Lattice>& momenta_ )
    : BasicDynamics<T,Lattice>(momenta_),
      boundaryDynamics(omega_, momenta_)
{ }

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

template<typename T, template<typename U> class Lattice, typename Dynamics>
void CombinedRLBdynamics<T,Lattice,Dynamics>::collide (
        Cell<T,Lattice>& cell,
        LatticeStatistics<T>& statistics )
{
    typedef Lattice<T> L;

    T rho, u[L::d], pi[util::TensorVal<Lattice<T> >::n];
    this->momenta.computeAllMomenta(cell,rho,u,pi);

    T uSqr = util::normSqr<T,L::d>(u);

    for (int iPop = 0; iPop < L::q; ++iPop)
    {
        cell[iPop] = computeEquilibrium(iPop,rho,u,uSqr) +
                     firstOrderLbHelpers<T,Lattice>::fromPiToFneq(iPop, pi);
    }

    boundaryDynamics.collide(cell, statistics);
}

template<typename T, template<typename U> class Lattice, typename Dynamics>
void CombinedRLBdynamics<T,Lattice,Dynamics>::staticCollide (
        Cell<T,Lattice>& cell,
        const T u[Lattice<T>::d],
        LatticeStatistics<T>& statistics )
{
    typedef Lattice<T> L;

    T rho, falseU[L::d], pi[util::TensorVal<Lattice<T> >::n];
    this->momenta.computeAllMomenta(cell, rho, falseU, pi);

    T uSqr = util::normSqr<T,L::d>(u);

    for (int iPop = 0; iPop < L::q; ++iPop)
    {
        cell[iPop] = computeEquilibrium(iPop,rho,u,uSqr) +
                     firstOrderLbHelpers<T,Lattice>::fromPiToFneq(iPop, pi);
    }

    boundaryDynamics.staticCollide(cell, u, statistics);
}

template<typename T, template<typename U> class Lattice, typename Dynamics>
T CombinedRLBdynamics<T,Lattice,Dynamics>::getOmega() const 
{
    return boundaryDynamics.getOmega();
}

template<typename T, template<typename U> class Lattice, typename Dynamics>
void CombinedRLBdynamics<T,Lattice,Dynamics>::setOmega(T omega_)
{
    boundaryDynamics.setOmega(omega_);
}

template<typename T, template<typename U> class Lattice, typename Dynamics>
T CombinedRLBdynamics<T,Lattice,Dynamics>::getParameter(int whichParameter) const {
    return boundaryDynamics.getParameter(whichParameter);
}

template<typename T, template<typename U> class Lattice, typename Dynamics>
void CombinedRLBdynamics<T,Lattice,Dynamics>::setParameter(int whichParameter, T value) {
    boundaryDynamics.setParameter(whichParameter, value);
}


////////////////////// Class ForcedBGKdynamics /////////////////////////

/** \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>
ForcedBGKdynamics<T,Lattice>::ForcedBGKdynamics (
        T omega_, Momenta<T,Lattice>& momenta_ )
    : BasicDynamics<T,Lattice>(momenta_),
      omega(omega_)
{
    // This ensures both that the constant sizeOfForce is defined in
    // ExternalField and that it has the proper size
    OLB_PRECONDITION( Lattice<T>::d == Lattice<T>::ExternalField::sizeOfForce );
}

template<typename T, template<typename U> class Lattice>
ForcedBGKdynamics<T,Lattice>* ForcedBGKdynamics<T,Lattice>::clone() const {
    return new ForcedBGKdynamics<T,Lattice>(*this);
}
 
template<typename T, template<typename U> class Lattice>
T ForcedBGKdynamics<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 ForcedBGKdynamics<T,Lattice>::collide (
        Cell<T,Lattice>& cell,
        LatticeStatistics<T>& statistics )
{
    T rho, u[Lattice<T>::d];
    this->momenta.computeRhoU(cell, rho, u);
    T* force = cell.getExternal(forceBeginsAt);
    for (int iVel=0; iVel<Lattice<T>::d; ++iVel) {
        u[iVel] += force[iVel] / (T)2.;
    }
    T uSqr = lbHelpers<T,Lattice>::bgkCollision(cell, rho, u, omega);
    lbHelpers<T,Lattice>::addExternalForce(cell, u, omega);
    if (cell.takesStatistics()) {
        statistics.gatherStats(rho, uSqr);
    }
}

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

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

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



////////////////////// Class D3Q13dynamics /////////////////////////

/** \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>
D3Q13dynamics<T,Lattice>::D3Q13dynamics (
        T omega_, Momenta<T,Lattice>& momenta_ )
    : BasicDynamics<T,Lattice>(momenta_)
{
    setOmega(omega_);
}

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

template<typename T, template<typename U> class Lattice>
T D3Q13dynamics<T,Lattice>::computeEquilibrium(int iPop, T rho, const T u[Lattice<T>::d], T uSqr) const
{
    // To get at the equilibrium, execute collision with relaxation parameters 1
    Cell<T,Lattice> tmp[Lattice<T>::q];
    for (int pop=0; pop<Lattice<T>::q; ++pop) {
        tmp[pop] = Lattice<T>::t[pop];
    }
    d3q13Helpers<T>::collision(tmp, rho, u, (T)1, (T)1);
    return tmp[iPop];
}


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

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

template<typename T, template<typename U> class Lattice>
T D3Q13dynamics<T,Lattice>::getOmega() const {
    return (T)4 / ( (T)3/lambda_nu + (T)1/(T)2 );
}

template<typename T, template<typename U> class Lattice>
void D3Q13dynamics<T,Lattice>::setOmega(T omega_) {
    lambda_nu = (T)3 / ( (T)4/omega_ - (T)1/(T)2 );
    lambda_nu_prime = (T)3 / ( (T)2/omega_ + (T)1/(T)2 );
}

////////////////////// Class Momenta //////////////////////////////

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

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

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

////////////////////// Class BulkMomenta //////////////////////////

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

template<typename T, template<typename U> class Lattice>
void BulkMomenta<T,Lattice>::computeU(Cell<T,Lattice> const& cell, T u[Lattice<T>::d]) const
{
    T dummyRho;
    lbHelpers<T,Lattice>::computeRhoU(cell, dummyRho, u);
}

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

template<typename T, template<typename U> class Lattice>
void BulkMomenta<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 BulkMomenta<T,Lattice>::computeRhoU (
        Cell<T,Lattice> const& cell,
        T& rho, T u[Lattice<T>::d] ) const
{
    lbHelpers<T,Lattice>::computeRhoU(cell, rho,u);
}

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

template<typename T, template<typename U> class Lattice>
void BulkMomenta<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 BulkMomenta<T,Lattice>::defineU (
    Cell<T,Lattice>& cell,
    const T u[Lattice<T>::d])
{
    T rho, oldU[Lattice<T>::d];
    computeRhoU(cell, rho, oldU);
    T uSqr = util::normSqr<T,Lattice<T>::d>(u);
    T fNeq[Lattice<T>::q];
    lbHelpers<T,Lattice>::computeFneq(cell, fNeq, rho, oldU);
    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 BulkMomenta<T,Lattice>::defineRhoU (
    Cell<T,Lattice>& cell,
    T rho, const T u[Lattice<T>::d])
{
    T oldRho, oldU[Lattice<T>::d];
    computeRhoU(cell, oldRho, oldU);
    T uSqr = util::normSqr<T,Lattice<T>::d>(u);
    T fNeq[Lattice<T>::q];
    lbHelpers<T,Lattice>::computeFneq(cell, fNeq, oldRho, oldU);
    for (int iPop=0; iPop < Lattice<T>::q; ++iPop) {