Trying to understand the source code

Palabos General
#1

Hi All,

I tried to understand the source code the last days but there are some questions I don’t get answered myself. Any help would be highly appreciated.

First thing I struggled with was the function ComputeRhoU contained in lbhelpers.h:


    static void computeRhoU (
        Cell<T,Lattice> const& cell, T& rho, T u[Lattice<T>::d] )
    {
        rho = T();
        for (int iD=0; iD < Lattice<T>::d; ++iD) {
            u[iD] = T();
        }
        for (int iPop=0; iPop < Lattice<T>::q; ++iPop) {
            rho += cell[iPop];
            for (int iD=0; iD < Lattice<T>::d; ++iD) {
                u[iD] += cell[iPop]*Lattice<T>::c[iPop][iD];
            }
        }
        [b]rho += (T)1;[/b]
        for (int iD=0; iD < Lattice<T>::d; ++iD) {
            u[iD] /= rho;
        }
    }

Why do you add 1 to rho?

Also the substraction of the lattice weight from the equilibrium function ((lbhelpers.h, D2Q9) in a mystery to me:


static T equilibrium( int iPop, T rho, const T u[Lattice<T>::d],
                          const T uSqr )
    {
        T c_u = T();
        for (int iD=0; iD < Lattice<T>::d; ++iD) {
           c_u += Lattice<T>::c[iPop][iD]*u[iD];
        }
        return rho * Lattice<T>::t[iPop] * (
               (T)1 + Lattice<T>::invCs2 * c_u +
               Lattice<T>::invCs2 * Lattice<T>::invCs2/(T)2 * c_u*c_u -
               Lattice<T>::invCs2/(T)2 * uSqr
           ) [b]- Lattice<T>::t[iPop][/b];
    }

In the case of MRT the lattice weights are also added to the current and equilibrium distributions before transforming them into momentum space:


/// Computation of all equilibrium distribution (in momenta space)
    static void computeEquilibrium( T momentaEq[Lattice<T>::q], 
                                 T rho, const T u[Lattice<T>::d],
                                 const T uSqr )
    {
        for (int iPop = 0; iPop < Lattice<T>::q; ++iPop)
        {
            momentaEq[iPop] = T();
            for (int jPop = 0; jPop < Lattice<T>::q; ++jPop)
            {
                momentaEq[iPop] += Lattice<T>::M[iPop][jPop] * 
                        (lbHelpers<T,Lattice>::equilibrium(jPop,rho,u,uSqr) + 
                        [b]Lattice<T>::t[jPop])[/b];
            }
        }
    }
    
    static void computeMomenta(T momenta[Lattice<T>::q], Cell<T,Lattice> &cell)
    {
        for (int iPop = 0; iPop < Lattice<T>::q; ++iPop)
        {
            momenta[iPop] = T();
            for (int jPop = 0; jPop < Lattice<T>::q; ++jPop)
            {
                momenta[iPop] += Lattice<T>::M[iPop][jPop] * 
                        (cell[jPop] [b]+ Lattice<T>::t[jPop][/b]);
            }
        }
    }

Sorry to bother you guys but I really couldn’t resolve that until now. Thanks in advance,

Martin

#2

To all these questions there’s one quick answer:

http://www.lbmethod.org/doku.php?id=models:tricks

#3

You just have to know where to find it! Thanks a lot Jonas, that really helps!

#4

Hi Jonas;
I cannot use the link that you attached. I am new to LBM and there are some parameter from the above code that I could not figure out. It would be great if you could help me to explain. What are t[Pop] and c[iPop][iD]
Thank you for your time