Palabos源码：computeEquilibrium(iPop, rhoBar, j, jSqr)的过程

之前我提到过这个computeEquilibrium可以用来修改cell内部属性，
它的调用平平无奇

cell[iPop]=cell.getDynamics().computeEquilibrium(iPop, rhoBar, j, jSqr);

getDynamics()会去哪？

cell.getDynamics()后面可以接很多函数，而定义这些函数的地方为src/core/dynamics.h和.hh，这个源码值得单独开一篇博文介绍

套娃一层

h文件

   /// Compute equilibrium distribution function
    virtual T computeEquilibrium(plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
                                 T jSqr, T thetaBar=T()) const =0;

hh文件

template<typename T, template<typename U> class Descriptor>
T CompositeDynamics<T,Descriptor>::computeEquilibrium (
        plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j, T jSqr, T thetaBar) const
{
    return baseDynamics -> computeEquilibrium(iPop, rhoBar, j, jSqr, thetaBar);
}

可以看到其他的dynamics里，对于定义为NoDynamics和BounceBack等地方，它会返回T()。

template<typename T, template<typename U> class Descriptor>
T NoDynamics<T,Descriptor>::computeEquilibrium(plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
                                            T jSqr, T thetaBar) const
{
    return T();
}

template<typename T, template<typename U> class Descriptor>
T BounceBack<T,Descriptor>::computeEquilibrium(plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
                                               T jSqr, T thetaBar) const
{
    return T();
}

附近也有computeEquilibria

template<typename T, template<typename U> class Descriptor>
void Dynamics<T,Descriptor>::computeEquilibria (
        Array<T,Descriptor<T>::q>& fEq,  T rhoBar, Array<T,Descriptor<T>::d> const& j, T jSqr, T thetaBar ) const
{
    for (int iPop=0; iPop<Descriptor<T>::q; ++iPop) {
        fEq[iPop] = computeEquilibrium(iPop, rhoBar, j, jSqr, thetaBar);
    }
}

套娃二层

但现在仅仅停留在下面这行代码，我们仍未知道 computeEquilibrium的内部代码。

return baseDynamics -> computeEquilibrium(iPop, rhoBar, j, jSqr, thetaBar);

经过一番寻找，在src/basicDynamics/isoThermalDynamics.h和hh文件中，找到了如下：
以class BGKdynamics内的 computeEquilibrium为例，我们找到了如下代码：

template<typename T, template<typename U> class Descriptor>
T BGKdynamics<T,Descriptor>::computeEquilibrium(plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
                                                T jSqr, T thetaBar) const
{
    T invRho = Descriptor<T>::invRho(rhoBar);
    return dynamicsTemplates<T,Descriptor>::bgk_ma2_equilibrium(iPop, rhoBar, invRho, j, jSqr);
}

套娃三层

接下来我们找一找bgk_ma2_equilibrium，让我们来看看位于src/latticeBoltzmann的dynamicTemplate.h
在struct dynamicsTemplates里面，显然它与Descriptor也有关

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);
}

套娃四层

接下来看src/latticeBoltzmann的dynamicTemplate3D.h
检索一下D3Q19，就能定位到如下代码区域
此处的c_j是microscopic flow velocity vector点乘方向向量，看 是不是很符合书上的feq公式呢？

static T bgk_ma2_equilibrium(plint iPop, T rhoBar, T invRho, Array<T,3> const& j, T jSqr ) {
    T c_j = D::c[iPop][0]*j[0] + D::c[iPop][1]*j[1] + D::c[iPop][2]*j[2];
    return D::t[iPop] * ( rhoBar + (T)3.*c_j + invRho*((T)4.5*c_j*c_j - (T)1.5*jSqr) );
}

其他地方的computeEquilibrium

以GuoExternalForceBGKdynamics为例
（basicDynamics/externalForceDynamics.h)

暂时不清楚这里的函数会被何时调用，猜测是在算例的迭代过程时内部运算。
virtual T computeEquilibrium(plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j, T jSqr, T thetaBar=T()) const;

/// Implementation of O(Ma^2) BGK dynamics with an external force (Guo approach)
template<typename T, template<typename U> class Descriptor>
class GuoExternalForceBGKdynamics : public ExternalForceDynamics<T,Descriptor> {
public:
/* *************** Construction / Destruction ************************ */
    GuoExternalForceBGKdynamics(T omega_);
    GuoExternalForceBGKdynamics(HierarchicUnserializer& unserializer);

    /// Clone the object on its dynamic type.
    virtual GuoExternalForceBGKdynamics<T,Descriptor>* clone() const;

    /// Return a unique ID for this class.
    virtual int getId() const;

/* *************** Collision and Equilibrium ************************* */

    /// Implementation of the collision step
    virtual void collide(Cell<T,Descriptor>& cell,
                         BlockStatistics& statistics_);

    /// Implementation of the collision step, with imposed macroscopic variables
    virtual void collideExternal(Cell<T,Descriptor>& cell, T rhoBar,
                         Array<T,Descriptor<T>::d> const& j, T thetaBar, BlockStatistics& stat);

    /// Compute equilibrium distribution function
    virtual T computeEquilibrium(plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
                                 T jSqr, T thetaBar=T()) const;
private:
    static int id;
};

basicDynamics/externalForceDynamics.hh

template<typename T, template<typename U> class Descriptor>
T GuoExternalForceBGKdynamics<T,Descriptor>::computeEquilibrium (
        plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
        T jSqr, T thetaBar) const
{
    T invRho = Descriptor<T>::invRho(rhoBar);
    return dynamicsTemplates<T,Descriptor>::bgk_ma2_equilibrium(iPop, rhoBar, invRho, j, jSqr);
}