Palabos源码：自由面模型的体积力是如何实现的

源码位置：src/multiphysics/bodyforce3D

template<typename T, template<typename U> class Descriptor>
AddConstForceToMomentum3D<T,Descriptor>::AddConstForceToMomentum3D(Array<T,3> const& force_)
    : force(force_)
{ }

可以看出体积力是通过一个AddConstForceToMomentum3D数据处理器实现的，它的具体代码为：

for (plint iX = domain.x0; iX <= domain.x1; iX++) {
        for (plint iY = domain.y0; iY <= domain.y1; iY++) {
            for (plint iZ = domain.z0; iZ <= domain.z1; iZ++) {
                // Momentum correction.
                Cell<T,Descriptor>& cell = lattice->get(iX, iY, iZ);
                if (cell.getDynamics().hasMoments()) {
                    T dynamicOmega = cell.getDynamics().getDynamicParameter(dynamicParams::dynamicOmega, cell);
                    T tau = 0.0;
                    if (!util::isZero(dynamicOmega)) {
                        tau = (T) 1 / dynamicOmega;
                    } else {
                        tau = (T) 1 / cell.getDynamics().getOmega();
                    }
                    T rho = Descriptor<T>::fullRho(rhoBar->get(iX+ofsR.x, iY+ofsR.y, iZ+ofsR.z));
                    j->get(iX+ofsJ.x, iY+ofsJ.y, iZ+ofsJ.z) += rho * tau * force;
                }
            }
        }

它判断了lattice的cell是否为流体格点，如果是则计算dynamicOmega和tau，提取出密度rho（rhobar和rho不一样），根据rho * tau * force对速度进行改变。在自由面模型的算例里，这个体积力可以是重力。

其他

AddCustomForceToMomentum3D：可根据函数施加自定义的体积力。

AddForceToMomentum3D：从blocks.size() == 4来看，这个应该是用力场导入体积力。

FreeSurfaceAddConstForceToMomentum3D：判断freeSurfaceFlag::isWet进行体积力施加。

FreeSurfaceAddForceToMomentum3D：同之前。

ComputeRotatingFrameForce3D，FreeSurfaceComputeRotatingFrameForce3D等

略