Palabos Multi-direct forcing scheme IBM RBC论文和代码叙述（不涉及FEM方面）

概览

在之前的浸入边界法文章中，我主要介绍的是Inamuro的IBM算法，近期科研进展涉及到提高IBM算法精度，加之Palabos本身自带Multi-direct forcing scheme，便以此作为较好的一个参考。

论文

Bridging the computational gap between mesoscopic and continuum modeling of
red blood cells for fully resolved blood flow

在我之前的博文里，IBM算法就是以这篇文章的公式去解释的。依稀记得当时Palabos仅有这篇文章预发表，软件尚未更新到最新版本（即包含RBC模拟的版本），所以我仅仅只能看到Inamuro IBM的代码。

介绍

让我们再次回到这里的Fimm循环，顾名思义就是我们要多次做direct-forcing，palabos实现的也很简单：

	for (pluint i = 0; i < sp.inamuro_iT; i++)
	{
		actions3.addProcessor(
			new MultiDirectForcingImmersedBoundaryIteration3D<T,
			DESCRIPTOR>(1. / omega, sp.incompressibleModel, sp.wallMesh),
			particleID_act3, rhoBarID_act3, jID_act3,
			group.getBoundingBox());
		actions3.addCommunication(jID_act3, modif::staticVariables);
	}

这里的sp.inamuro_iT是我们设置的参数，这块代码即是实现IBM的步骤。/代码来源：bloodFlowDefoBodies.cpp line 1969.

下面的整体步骤对你的编程具有一定的参考意义。

整体程序概览

2046	// At every iteration we want to advance both the bodies and the fluid from t to t+1
	// Thus we choose the sequence below:
	// 1. Compute density, velocity and stress fields (t)
	// 2. Fext(t) on solid bodies (& colliding neighbors)
	// 3. Solve the bodies with Fext(t) - But do not advance yet at t+1
	// 4. Change the momentum of the fluid using IBM & any other external force (Shan-Chen forcing scheme).
	//    The fluid needs to know the position & velocity of the immersed body at time t and based on this
	//    we are going to advance to t+1. Our schemes here are explicit and so we compute t+1 with the info
	//    that we have at t. This is why we did not advance at the previous step
	// 5. Collide & Stream (FLUID at t+1)
2056	// 6. Advance particles and copy new vels (SOLID at t+1)

/代码来源：bloodFlowDefoBodies.cpp line 2046.

MultiDirectForcingImmersedBoundaryIteration3D

template<typename T, template<typename U> class Descriptor>
MultiDirectForcingImmersedBoundaryIteration3D<T,Descriptor>::MultiDirectForcingImmersedBoundaryIteration3D(
        T tau_, bool incompressibleModel_,
        RawConnectedTriangleMesh<T>* wallMesh_)
    : tau(tau_),
      incompressibleModel(incompressibleModel_),
      wallMesh(wallMesh_)
{ }

/代码来源：coupledSimulators\npFEM\rbcParticles.h line 766.

Interpolation部分

/代码来源：coupledSimulators\npFEM\rbcParticles.h line 799.

在之前的推文中理应包含上图，现在贴出来供读者参考。

  if (incompressibleModel) {
        for (pluint iParticle=0; iParticle<particles.size(); ++iParticle) {//running through all particles
            SurfaceParticleWithAreaBase3D<T,Descriptor>* particle =
                dynamic_cast<SurfaceParticleWithAreaBase3D<T,Descriptor>*>(particles[iParticle]);
            PLB_ASSERT(particle);

            Array<T,3> particlePos = particle->getPosition();
            Array<plint,3> intPos((plint) particlePos[0] - location.x,
                                  (plint) particlePos[1] - location.y,
                                  (plint) particlePos[2] - location.z);
            const Array<plint,2> xLim((particlePos[0] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
            const Array<plint,2> yLim((particlePos[1] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
            const Array<plint,2> zLim((particlePos[2] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
            const Array<T,3> fracPos(util::frac(particlePos[0]), util::frac(particlePos[1]), util::frac(particlePos[2]));
            Array<T,3> averageJ; averageJ.resetToZero();//reset to zero in every interation
            // 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));//pos will be local lattices
                        Array<T,3> nextJ = j->get(pos[0]+ofsJ.x, pos[1]+ofsJ.y, pos[2]+ofsJ.z);//get velocity of local lattices from j's
                        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;//Interpolation
                    }
                }
            }

            Array<T,3> const& wallVelocity = particle->getVelocity();
            T area = particle->getArea(wallMesh);//dx*dx(not very sure)

            deltaG[iParticle] = area * (wallVelocity - averageJ);
        }

Force Spreading部分

/代码来源：coupledSimulators\npFEM\rbcParticles.h line 870.

 for (pluint iParticle=0; iParticle<particles.size(); ++iParticle) {
        SurfaceParticleWithAreaBase3D<T,Descriptor>* particle =
            dynamic_cast<SurfaceParticleWithAreaBase3D<T,Descriptor>*>(particles[iParticle]);
        PLB_ASSERT(particle);

        Array<T,3> particlePos = particle->getPosition();
        Array<plint,3> intPos((plint) particlePos[0] - location.x,
                              (plint) particlePos[1] - location.y,
                              (plint) particlePos[2] - location.z);
        const Array<plint,2> xLim((particlePos[0] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
        const Array<plint,2> yLim((particlePos[1] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
        const Array<plint,2> zLim((particlePos[2] < (T) 0 ? Array<plint,2>(-2, 1) : Array<plint,2>(-1, 2)));
        const Array<T,3> fracPos(util::frac(particlePos[0]), util::frac(particlePos[1]), util::frac(particlePos[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[iParticle];
                    j->get(pos[0]+ofsJ.x, pos[1]+ofsJ.y, pos[2]+ofsJ.z) = nextJ;//force spread and get new velocity
                }
            }
        }
    }
}

Initial Estimation

在进行Multi-direct forcing之前，我们需要计算此处的大Uk，这个大Uk是否solid solver计算出来的，而后面的小u星号（大Xk表示Lagrangian点）则通过Interpolation。

大Uk的由来

在Algorithm1中，提到了计算固体点的速度v。

v的计算

xcm是固体点的中心，center of mass
vcm是其速度
L angular momentum
I inertia matrix
w 角速度

我觉得下文代码非常清晰，无需过多解释。

代码

/代码来源：coupledSimulators\npFEM\solver.cpp line 953.

SHAPEOP_INLINE void Solver::PBD_Damping()
{
    velocities_ = (points_ - oldPoints_) / delta_;//get velocity from v=s/t

    Vector3 Xcm(0., 0., 0.), Vcm(0., 0., 0.);
    Matrix3X tmpX = (M_ * points_.transpose()).transpose();
    Matrix3X tmpV = (M_ * velocities_.transpose()).transpose();
    for (int i = 0; i < static_cast<int>(tmpX.cols()); ++i) {
        Xcm += tmpX.col(i);
        Vcm += tmpV.col(i);
    }
    Scalar TotalMass = (MatrixXX(M_)).trace();
    Xcm /= TotalMass;//center of mass
    Vcm /= TotalMass;//velocity of center of mass

    Xcm_ = Xcm;

    Matrix3X r = points_;//get locations of particles
    r.colwise() -= Xcm;//get the distance from center to particle

    Vector3 L(0., 0., 0.);
    for (int i = 0; i < static_cast<int>(tmpX.cols()); ++i) {
        L += (r.col(i)).cross(tmpV.col(i));//get angular momentum
    }

    Matrix33 Icm = Matrix33::Zero();
    for (int i = 0; i < static_cast<int>(tmpX.cols()); ++i) {
        // Cross product converted to matrix multiplication
        Matrix33 Rxi = Matrix33::Zero();
        Rxi(0, 1) = -r.col(i)[2];
        Rxi(0, 2) = r.col(i)[1];
        Rxi(1, 0) = r.col(i)[2];
        Rxi(1, 2) = -r.col(i)[0];
        Rxi(2, 0) = -r.col(i)[1];
        Rxi(2, 1) = r.col(i)[0];

        Icm += M_.coeff(i, i) * Rxi * Rxi.transpose();
    }
    Matrix33 Icm_inv = Icm.inverse();

    Vector3 Ocm = Icm_inv * L;//angular velocity

    for (int i = 0; i < static_cast<int>(velocities_.cols()); ++i)
    {
        // See PBD Muller
        Vector3 dvi = Vcm + Ocm.cross(r.col(i)) - velocities_.col(i);
        velocities_.col(i) += Calpha_ * dvi;//in the paper, damping coefficient it is 0.6
    }

#ifndef NPFEM_SA
    // Lattice Units
    Scalar C_vel = externalSolverSpatialResolution_ / externalSolverTimeResolution_;

    Matrix3X lattice_velocities = velocities_ / C_vel;

    size_t velocityLimit = 0;
    for (int i = 0; i < static_cast<int>(velocities_.cols()); ++i)
    {
        if (lattice_velocities.col(i).norm() >= max_u_lb_)
        {
            lattice_velocities.col(i).normalize();
            lattice_velocities.col(i) *= max_u_lb_fin_;

            velocities_.col(i) = lattice_velocities.col(i) * C_vel;

            ++velocityLimit;
        }
    }

#ifdef ENABLE_LOGS
    if (velocityLimit > 0)
    {
        global::logfile_nonparallel("instability.log")
            .flushEntry("At Palabos iT=" + util::val2str(Palabos_iT_) + " | ShapeOp bodyID=" + util::val2str(bodyID_) +
                " | VELOCITY LIMIT ACTIVATED, on " + util::val2str(velocityLimit) + " vertices");
    }
#endif // ENABLE_LOGS

#endif // !NPFEM_SA

    points_ = oldPoints_ + velocities_ * delta_;//new positions
}