dynamicsTemplates.h

Go to the documentation of this file.

/* This file is part of the Palabos library.
 *
 * Copyright (C) 2011 FlowKit Sarl
 * Avenue de Chailly 23
 * 1012 Lausanne, Switzerland
 * E-mail contact: contact@flowkit.com
 *
 * The most recent release of Palabos can be downloaded at 
 * <http://www.palabos.org/>
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * The library 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 Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef DYNAMICS_TEMPLATES_H
#define DYNAMICS_TEMPLATES_H

#include "core/globalDefs.h"
#include "core/cell.h"
#include "core/util.h"
#include "latticeBoltzmann/offEquilibriumTemplates.h"

namespace plb {

template<typename T, class Descriptor> struct dynamicsTemplatesImpl;

template<typename T, template<typename U> class Descriptor>
struct dynamicsTemplates {

static T bgk_ma2_equilibrium(plint iPop, T rhoBar, T invRho, Array<T,Descriptor<T>::d> const& j, T jSqr) {
    return dynamicsTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
        ::bgk_ma2_equilibrium(iPop, rhoBar, invRho, j, jSqr);
}

static void bgk_ma2_equilibria( T rhoBar, T invRho, Array<T,Descriptor<T>::d> const& j,
                                T jSqr, Array<T,Descriptor<T>::q>& eqPop )
{
    dynamicsTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
        ::bgk_ma2_equilibria(rhoBar, invRho, j, jSqr, eqPop);
}

static T bgk_ma2_collision(Cell<T,Descriptor>& cell, T rhoBar, Array<T,Descriptor<T>::d> const& j, T omega)
{
    return dynamicsTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
        ::bgk_ma2_collision(cell.getRawPopulations(), rhoBar, j, omega);
}

static T bgk_inc_collision(Cell<T,Descriptor>& cell, T rhoBar, Array<T,Descriptor<T>::d> const& j, T omega)
{
    return dynamicsTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
        ::bgk_inc_collision(cell.getRawPopulations(), rhoBar, j, omega);
}

static T rlb_collision(Cell<T,Descriptor>& cell, T rhoBar, T invRho, Array<T,Descriptor<T>::d> const& j,
                       Array<T,SymmetricTensor<T,Descriptor>::n> const& PiNeq, T omega )
{
    return dynamicsTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
        ::rlb_collision(cell.getRawPopulations(), rhoBar, invRho, j, PiNeq, omega);
}

static T bgk_ma2_constRho_collision(Cell<T,Descriptor>& cell,
                                T rhoBar, Array<T,Descriptor<T>::d> const& j, T ratioRho, T omega)
{
    return dynamicsTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
               ::bgk_ma2_constRho_collision(cell.getRawPopulations(), rhoBar, j, ratioRho, omega);
}

static T precond_bgk_ma2_equilibrium(plint iPop, T rhoBar, T invRho, Array<T,Descriptor<T>::d> const& j, T jSqr, T invGamma)
{
    return dynamicsTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
        ::precond_bgk_ma2_equilibrium(iPop, rhoBar, invRho, j, jSqr, invGamma);
}

static T precond_bgk_ma2_collision(Cell<T,Descriptor>& cell, T rhoBar, Array<T,Descriptor<T>::d> const& j, T omega, T invGamma)
{
    return dynamicsTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
        ::precond_bgk_ma2_collision(cell.getRawPopulations(), rhoBar, j, omega, invGamma);
}

static T bgk_ma4_equilibrium(plint iPop, T rhoBar, T invRho, Array<T,Descriptor<T>::d> const& j, T jSqr, T rhoThetaBar) {
    return dynamicsTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
        ::bgk_ma4_equilibrium(iPop, rhoBar, invRho, j, jSqr, rhoThetaBar);
}

static T bgk_ma4_collision(Cell<T,Descriptor>& cell, T rhoBar, Array<T,Descriptor<T>::d> const& j, T rhoThetaBar, T omega)
{
    return dynamicsTemplatesImpl<T,typename Descriptor<T>::BaseDescriptor>
        ::bgk_ma4_collision(cell.getRawPopulations(), rhoBar, j, rhoThetaBar, omega);
}

};  // struct dynamicsTemplates


template<typename T, class Descriptor>
struct dynamicsTemplatesImpl {

static T bgk_ma2_equilibrium(plint iPop, T rhoBar, T invRho, Array<T,Descriptor::d> const& j, T jSqr) {
    T c_j = Descriptor::c[iPop][0]*j[0];
    for (int iD=1; iD < Descriptor::d; ++iD) {
       c_j += Descriptor::c[iPop][iD]*j[iD];
    }
    return Descriptor::t[iPop] * (
           rhoBar + Descriptor::invCs2 * c_j +
           Descriptor::invCs2/(T)2 * invRho * (
               Descriptor::invCs2 * c_j*c_j - jSqr )
       );
}

static void bgk_ma2_equilibria( T rhoBar, T invRho, Array<T,Descriptor::d>
        const& j,
                                T jSqr, Array<T,Descriptor::q>& eqPop )
{
    for (int iPop=0; iPop<Descriptor::d; ++iPop) {
        eqPop[iPop] = bgk_ma2_equilibrium(iPop, rhoBar, invRho, j, jSqr);
    }
}

static T bgk_ma2_collision(Array<T,Descriptor::q>& f, T rhoBar, Array<T,Descriptor::d> const& j, T omega) {
    T invRho = Descriptor::invRho(rhoBar);
    const T jSqr = VectorTemplateImpl<T,Descriptor::d>::normSqr(j);
    for (plint iPop=0; iPop < Descriptor::q; ++iPop) {
        f[iPop] *= (T)1-omega;
        f[iPop] += omega * dynamicsTemplatesImpl<T,Descriptor>::bgk_ma2_equilibrium (
                                iPop, rhoBar, invRho, j, jSqr );
    }
    return jSqr*invRho*invRho;
}

static T bgk_inc_collision(Array<T,Descriptor::q>& f, T rhoBar, Array<T,Descriptor::d> const& j, T omega) {
    // In incompressible BGK, the Ma^2 term is preceeded by 1 instead of 1/rho.
    T invRho = (T)1;
    const T jSqr = VectorTemplateImpl<T,Descriptor::d>::normSqr(j);
    for (plint iPop=0; iPop < Descriptor::q; ++iPop) {
        f[iPop] *= (T)1-omega;
        f[iPop] += omega * dynamicsTemplatesImpl<T,Descriptor>::bgk_ma2_equilibrium (
                                iPop, rhoBar, invRho, j, jSqr );
    }
    return jSqr;
}

static T rlb_collision (
    Array<T,Descriptor::q>& f, T rhoBar, T invRho, Array<T,Descriptor::d> const& j,
    Array<T,SymmetricTensorImpl<T,Descriptor::d>::n> const& PiNeq, T omega )
{
    const T jSqr = VectorTemplateImpl<T,Descriptor::d>::normSqr(j);
    f[0] = dynamicsTemplatesImpl<T,Descriptor>::bgk_ma2_equilibrium(0, rhoBar, invRho, j, jSqr)
                  + ((T)1-omega) *
                        offEquilibriumTemplatesImpl<T,Descriptor>::fromPiToFneq(0, PiNeq);
    for (plint iPop=1; iPop<=Descriptor::q/2; ++iPop) {
        f[iPop] = dynamicsTemplatesImpl<T,Descriptor>::
                bgk_ma2_equilibrium(iPop, rhoBar, invRho, j, jSqr);
        f[iPop+Descriptor::q/2] = dynamicsTemplatesImpl<T,Descriptor>::
                bgk_ma2_equilibrium(iPop+Descriptor::q/2, rhoBar, invRho, j, jSqr);
        T fNeq = ((T)1-omega) *
                 offEquilibriumTemplatesImpl<T,Descriptor>::fromPiToFneq(iPop, PiNeq);
        f[iPop] += fNeq;
        f[iPop+Descriptor::q/2] += fNeq;
    }
    return jSqr*invRho*invRho;
}

static T bgk_ma2_constRho_collision(Array<T,Descriptor::q>& f, T rhoBar, Array<T,Descriptor::d> const& j, T ratioRho, T omega) {
    T invRho = Descriptor::invRho(rhoBar);
    const T jSqr = VectorTemplateImpl<T,Descriptor::d>::normSqr(j);
    for (plint iPop=0; iPop < Descriptor::q; ++iPop) {
        T feq = dynamicsTemplatesImpl<T,Descriptor>::
                    bgk_ma2_equilibrium(iPop, rhoBar, invRho, j, jSqr );
        f[iPop] =
          ratioRho*feq + Descriptor::t[iPop]*(ratioRho-(T)1) +
          ((T)1-omega)*(f[iPop]-feq);
    }
    return jSqr*invRho*invRho;
}

static T precond_bgk_ma2_equilibrium(plint iPop, T rhoBar, T invRho, Array<T,Descriptor::d> const& j, T jSqr, T invGamma) {
    T c_j = Descriptor::c[iPop][0]*j[0];
    for (int iD=1; iD < Descriptor::d; ++iD) {
       c_j += Descriptor::c[iPop][iD]*j[iD];
    }
    return Descriptor::t[iPop] * (
           rhoBar + Descriptor::invCs2 * c_j +
           invGamma * Descriptor::invCs2/(T)2 * invRho * (
               Descriptor::invCs2 * c_j*c_j - jSqr )
       );
}

static T precond_bgk_ma2_collision(Array<T,Descriptor::q>& f, T rhoBar, Array<T,Descriptor::d> const& j, T omega, T invGamma) {
    T invRho = Descriptor::invRho(rhoBar);
    const T jSqr = VectorTemplateImpl<T,Descriptor::d>::normSqr(j);
    for (plint iPop=0; iPop < Descriptor::q; ++iPop) {
        f[iPop] *= (T)1-omega;
        f[iPop] += omega * dynamicsTemplatesImpl<T,Descriptor>::precond_bgk_ma2_equilibrium (
                                iPop, rhoBar, invRho, j, jSqr, invGamma );
    }
    return jSqr*invRho*invRho;
}

static T bgk_ma4_equilibrium(plint iPop, T rhoBar, T invRho, Array<T,Descriptor::d> const& j, T jSqr, T rhoThetaBar) {
    typedef Descriptor L;
    T c_j = (T)L::c[iPop][0]*j[0];
    for (int iD=1; iD < L::d; ++iD) {
        c_j += (T)L::c[iPop][iD]*j[iD];
    }
    
    T c_u = c_j*invRho;
    T uSqr = jSqr*invRho*invRho;
    T thetaBar = invRho*rhoThetaBar;
    
    return L::t[iPop] * (
        rhoBar + L::invCs2 * c_j +
        L::invCs2/(T)2 * ( invRho * (L::invCs2 * c_j*c_j - jSqr )
                                    + rhoThetaBar * ((T)L::cNormSqr[iPop]-L::cs2*(T)L::d))
        + L::invCs2*L::invCs2*L::invCs2/(T)6 * c_j * (
            c_u*c_u-(T)3*L::cs2*uSqr+(T)3*L::cs2*thetaBar*(L::cNormSqr[iPop]-L::cs2*((T)L::d+(T)2))
        )
        + L::invCs2*L::invCs2*L::invCs2*L::invCs2/(T)24 * (
            c_j*c_u*c_u*c_u - (T)6*L::cs2*uSqr*c_j*c_u + (T)3*L::cs2*L::cs2*jSqr*invRho*uSqr
            + (T)6*L::cs2*rhoThetaBar*(c_u*c_u*((T)L::cNormSqr[iPop]-L::cs2*((T)L::d+(T)4))+L::cs2*uSqr*(L::cs2*((T)L::d+(T)2)-(T)L::cNormSqr[iPop]))
            + (T)3*L::cs2*L::cs2*rhoThetaBar*rhoThetaBar*invRho*(
                (T)L::cNormSqr[iPop]*(T)L::cNormSqr[iPop]-(T)2*L::cs2*(L::d+2)*L::cNormSqr[iPop]+L::cs2*L::cs2*(T)L::d*((T)L::d+(T)2)
            )
        )
    );
}

static T bgk_ma4_collision(Array<T,Descriptor::q>& f, T rhoBar, Array<T,Descriptor::d> const& j, T rhoThetaBar, T omega) {
    T invRho = Descriptor::invRho(rhoBar);
    const T jSqr = VectorTemplateImpl<T,Descriptor::d>::normSqr(j);
    for (plint iPop=0; iPop < Descriptor::q; ++iPop) {
        f[iPop] *= (T)1-omega;
        f[iPop] += omega * dynamicsTemplatesImpl<T,Descriptor>::bgk_ma4_equilibrium (
            iPop, rhoBar, invRho, j, jSqr, rhoThetaBar);
    }
    return jSqr*invRho*invRho;
}

};  // struct dynamicsTemplatesImpl

}  // namespace plb

#include "latticeBoltzmann/dynamicsTemplates2D.h"
#include "latticeBoltzmann/dynamicsTemplates3D.h"

#endif  // DYNAMICS_TEMPLATES_H