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

文章目录

getDynamics()会去哪？
- 套娃一层
- 套娃二层
- 套娃三层
- 套娃四层
其他地方的computeEquilibrium

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

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

getDynamics()会去哪？

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

套娃一层

h文件

   /// Compute equilibrium distribution functionvirtual 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 stepvirtual void collide(Cell<T,Descriptor>& cell,BlockStatistics& statistics_);/// Implementation of the collision step, with imposed macroscopic variablesvirtual void collideExternal(Cell<T,Descriptor>& cell, T rhoBar,Array<T,Descriptor<T>::d> const& j, T thetaBar, BlockStatistics& stat);/// Compute equilibrium distribution functionvirtual 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);
}