Palabos源码学习 | 浸入边界法 | 论文公式对应源码讲解

Conbined multi-direct forcing method and IBM算法概览

截图源自论文：Bridging the computational gap between mesoscopic and continuum modeling of
red blood cells for fully resolved blood flow

其算法是IBM的一个版本，Multi-direct forcing method的IBM，论文标题：Combined multi-direct forcing and immersed boundary method for simulating flows with moving particles

前言

在前面一篇文章中我已经说明了一些关于IBM在程序中的调用和用法，本篇主要是单独讲一下论文上公式对应代码的部分，故独立出来。但不再重复代码结构功能的解释，比较初学的同学可以去先看前面那篇文章。

主要目标是梳理算法思路，找到对应的代码公式关系，方便优化算法的同学在Palabos基础上修改程序。

Palabos源码

位置：src\offLattice\immersedWalls3D.h和.hh。

以InamuroIteration3D为例。

> /* ******** InamuroIteration3D ************************************ */
> 
> template<typename T, class VelFunction>
> InamuroIteration3D<T,VelFunction>::InamuroIteration3D(VelFunction
> velFunction_, T tau_, bool incompressibleModel_)
>     : velFunction(velFunction_),
>       tau(tau_),
>       incompressibleModel(incompressibleModel_) { }
> 
> template<typename T, class VelFunction> void
> InamuroIteration3D<T,VelFunction>::processGenericBlocks (
>         Box3D domain, std::vector<AtomicBlock3D*> blocks ) {
>     PLB_PRECONDITION( blocks.size()==3 );
>     ScalarField3D<T>* rhoBar = dynamic_cast<ScalarField3D<T>*>(blocks[0]);
>     TensorField3D<T,3>* j = dynamic_cast<TensorField3D<T,3>*>(blocks[1]);
>     AtomicContainerBlock3D* container = dynamic_cast<AtomicContainerBlock3D*>(blocks[2]);
>     PLB_ASSERT( rhoBar );
>     PLB_ASSERT( j );
>     PLB_ASSERT( container );
>     Dot3D location = rhoBar->getLocation();
>     Dot3D ofsJ = computeRelativeDisplacement(*rhoBar, *j);
>     ImmersedWallData3D<T>* wallData = 
>         dynamic_cast<ImmersedWallData3D<T>*>( container->getData() );
>     PLB_ASSERT(wallData);
> 
>     std::vector< Array<T,3> > const& vertices = wallData->vertices;
>     std::vector<T> const& areas = wallData->areas;
>     PLB_ASSERT( vertices.size()==areas.size() );
>     std::vector<Array<T,3> > deltaG(vertices.size());
>     std::vector<Array<T,3> >& g = wallData->g;
>     PLB_ASSERT( vertices.size()==g.size() );
> 
>     // In this iteration, the force is computed for every vertex.
>     if (incompressibleModel) {
>         for (pluint i=0; i<vertices.size(); ++i) {
>             Array<T,3> const& vertex = vertices[i];
>             Array<plint,3> intPos((plint) vertex[0] - location.x, (plint) vertex[1] - location.y, (plint) vertex[2] - location.z);
>             const Array<plint,2> xLim((vertex[0] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
>             const Array<plint,2> yLim((vertex[1] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
>             const Array<plint,2> zLim((vertex[2] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
>             const Array<T,3> fracPos(util::frac(vertex[0]), util::frac(vertex[1]), util::frac(vertex[2]));
>             Array<T,3> averageJ; averageJ.resetToZero();
>             // Use the weighting function to compute the average momentum
>             // and the average density on the surface vertex.
>             // x   x . x   x
>             for (plint dx = xLim[0]; dx <= xLim[1]; dx++) {
>                 for (plint dy = yLim[0]; dy <= yLim[1]; dy++) {
>                     for (plint dz = zLim[0]; dz <= zLim[1]; dz++) {
>                         Array<plint,3> pos(intPos+Array<plint,3>(dx,dy,dz));
>                         Array<T,3> nextJ = j->get(pos[0]+ofsJ.x, pos[1]+ofsJ.y, pos[2]+ofsJ.z);
>                         Array<T,3> r((T)dx-fracPos[0],(T)dy-fracPos[1],(T)dz-fracPos[2]);
>                         T W = inamuroDeltaFunction<T>().W(r);
>                         averageJ += W*nextJ;//此处为论文的Step3，Lagrangian点通过Dirac函数插值得到速度
>                     }
>                 }
>             }
>             //averageJ += (T)0.5*g[i];
>             Array<T,3> wallVelocity = velFunction(vertex);//这里应当对应step3，更新Lagrangian点速度，但这个Palabos算例的固体点速度是直接给定的，具体我不在这里展开。
>             deltaG[i] = areas[i]*(wallVelocity-averageJ);//此处为论文的Step4，Lagrangian点Xk的体积力计算
>             g[i] += deltaG[i];//g会被留下来，供有需要的时候提取，注意这里不仅仅只是赋值，而是+=，累计的效果。
>         }
>     } else { // Compressible model.
>         for (pluint i=0; i<vertices.size(); ++i) {
>             Array<T,3> const& vertex = vertices[i];
>             Array<plint,3> intPos((plint) vertex[0] - location.x, (plint) vertex[1] - location.y, (plint) vertex[2] - location.z);
>             const Array<plint,2> xLim((vertex[0] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
>             const Array<plint,2> yLim((vertex[1] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
>             const Array<plint,2> zLim((vertex[2] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
>             const Array<T,3> fracPos(util::frac(vertex[0]), util::frac(vertex[1]), util::frac(vertex[2]));
>             Array<T,3> averageJ; averageJ.resetToZero();
>             T averageRhoBar = T();
>             // Use the weighting function to compute the average momentum
>             // and the average density on the surface vertex.
>             // x   x . x   x
>             for (plint dx = xLim[0]; dx <= xLim[1]; dx++) {
>                 for (plint dy = yLim[0]; dy <= yLim[1]; dy++) {
>                     for (plint dz = zLim[0]; dz <= zLim[1]; dz++) {
>                         Array<plint,3> pos(intPos+Array<plint,3>(dx,dy,dz));
>                         T nextRhoBar = rhoBar->get(pos[0], pos[1], pos[2]);
>                         Array<T,3> nextJ = j->get(pos[0]+ofsJ.x, pos[1]+ofsJ.y, pos[2]+ofsJ.z);
>                         Array<T,3> r((T)dx-fracPos[0],(T)dy-fracPos[1],(T)dz-fracPos[2]);
>                         T W = inamuroDeltaFunction<T>().W(r);
>                         averageJ += W*nextJ;
>                         averageRhoBar += W*nextRhoBar;
>                     }
>                 }
>             }
>             //averageJ += (T)0.5*g[i];
>             Array<T,3> wallVelocity = velFunction(vertex);
>             deltaG[i] = areas[i]*((averageRhoBar+(T)1.)*wallVelocity-averageJ);
>             //g[i] += deltaG[i];
>             g[i] += deltaG[i]/((T)1.0+averageRhoBar);
>         }
>     }
>     //上面的是压缩模型的部分，这里不作解释。
>     //下面的循环即表示force spreading的部分。
>     // In this iteration, the force is applied from every vertex to the grid nodes.
>     for (pluint i=0; i<vertices.size(); ++i) {
>         Array<T,3> const& vertex = vertices[i];
>         Array<plint,3> intPos((plint) vertex[0] - location.x, (plint) vertex[1] - location.y, (plint) vertex[2] - location.z);
>         const Array<plint,2> xLim((vertex[0] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
>         const Array<plint,2> yLim((vertex[1] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
>         const Array<plint,2> zLim((vertex[2] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
>         const Array<T,3> fracPos(util::frac(vertex[0]), util::frac(vertex[1]), util::frac(vertex[2]));
>         for (plint dx = xLim[0]; dx <= xLim[1]; dx++) {
>             for (plint dy = yLim[0]; dy <= yLim[1]; dy++) {
>                 for (plint dz = zLim[0]; dz <= zLim[1]; dz++) {
>                     Array<plint,3> pos(intPos+Array<plint,3>(dx,dy,dz));
>                     Array<T,3> nextJ = j->get(pos[0]+ofsJ.x, pos[1]+ofsJ.y, pos[2]+ofsJ.z);
>                     Array<T,3> r((T)dx-fracPos[0],(T)dy-fracPos[1],(T)dz-fracPos[2]);
>                     T W = inamuroDeltaFunction<T>().W(r);
>                     nextJ += tau*W*deltaG[i]; //Corrected velocity//此处代表论文的Algorithm2.3
>                     j->get(pos[0]+ofsJ.x, pos[1]+ofsJ.y, pos[2]+ofsJ.z) = nextJ;//此处是2.3的3，注意△t在这个操作中约掉了。
>                 }
>             }
>         }
>     } }

代码中pos[] 是Eulerian点的坐标,但加上ofsJ.xyz便代表速度velocity。

这个方法是Palabos中通过computeRelativeDisplacement实现的，这个功能在User guide和算例boussinesqThermal中都有介绍，此处不过多介绍。