Palabos程序代码解读 | particlesInCone

定义部分

Box3D injectionDomain;
plint startParticleIter = 1200; // Start iterating the particles after this iteration.


T particleProbabilityPerCell = 4.e-6; // Per-cell rate of particle injection.
T fluidCompliance = 6.e-5;            // The fluid->particle force is equal to the velocity difference times the fluid-compliance.
Array<T,3> gravity(-3.0e-7, 0., 0.);  // A body force acting on each particle.
T forceAmplitude = 1.e-4;             // Amplitude for the pairwise particle interaction force.
T cutOffLength = 3.5;                 // Cutoff length for particle interaction in lattice units (make sure the interaction force
                                      // is negligible at this distance).
plint commEnvelope = (plint)cutOffLength +1;  // Width of communication envelope needed to cope with neighborhoods within the cutoff length.
                                              // Must be bigger than blockSize/2.


typedef DenseParticleField3D<T,DESCRIPTOR> ParticleFieldT;

MultiParticleField3D<ParticleFieldT>* particles    =0;

Particles之间作用的数据处理器。

template<typename T, template<typename U> class Descriptor>
class ParticleInteractionForce : public BoxProcessingFunctional3D
{
public:
    ParticleInteractionForce(T forceAmplitude_, Array<T,3> gravity_)
        : forceAmplitude(forceAmplitude_),
          gravity(gravity_)
    { }
    // This is the central part of the code, where the particle interaction is defined.
    virtual void processGenericBlocks(Box3D domain, std::vector<AtomicBlock3D*> blocks)
    {
        PLB_PRECONDITION( blocks.size()==1 );
        ParticleField3D<T,Descriptor>& particleField = *dynamic_cast<ParticleField3D<T,Descriptor>*>(blocks[0]);
        Dot3D offset = particleField.getLocation();

        std::vector<Particle3D<T,Descriptor>*> particles;
        std::vector<Particle3D<T,Descriptor>*> neighbors;
        particleField.findParticles(domain, particles);
        // Loop over all particles assigned to this data processor.
        for (pluint iParticle=0; iParticle<particles.size(); ++iParticle) {
            Particle3D<T,Descriptor>* nonTypeParticle = particles[iParticle];
            if (nonTypeParticle->getTag()!=0) {  // Exclude the wall particles.
                VerletParticle3D<T,Descriptor>* particle =
                    dynamic_cast<VerletParticle3D<T,Descriptor>*>(nonTypeParticle);
                // Compute a neighborhood, in integer coordinates, which contains at least
                // all particles inside a radius of cutOffLength.
                Array<T,3> position = particle->getPosition();
                // Convert global particle coordinates to local coordinates for the domain
                // treated by this data processor.
                plint x = util::roundToInt(position[0])-offset.x;
                plint y = util::roundToInt(position[1])-offset.y;
                plint z = util::roundToInt(position[2])-offset.z;
                Box3D neighborhood(x-commEnvelope,x+commEnvelope,
                                   y-commEnvelope,y+commEnvelope, z-commEnvelope,z+commEnvelope);
                neighbors.clear();
                T rCritical = 2.0;
                T rCriticalSqr = util::sqr(rCritical);
                // Use the particle hash to find neighboring particles efficiently.
                particleField.findParticles(neighborhood, neighbors);
                Array<T,3> force((T)0.,(T)0.,(T)0.);
                for (pluint iNeighbor=0; iNeighbor<neighbors.size(); ++iNeighbor) {
                    Array<T,3> r = particle->getPosition() - neighbors[iNeighbor]->getPosition();
                    T rSqr = normSqr(r);
                    if (rSqr<util::sqr(cutOffLength) && rSqr>1.e-6) {
                        // The potential is amplitude*r^{-6} for r>rCritical, and r*amplitude for r<=rCritical.
                        // This is to aovid numerical instability due to huge forces when two
                        // particles come too close (for example in a badly designed initial condition).
                        if (rSqr>rCriticalSqr) {
                            // The potential is r^{-6}, so the force is r^{-7}. An additional r^{-1}
                            // is needed to normalize the vector r.
                            force += r/(rSqr*rSqr*rSqr*rSqr) * forceAmplitude;
                        }
                        else {
                            // Inside a radius of rCritical (in lattice units), the force is constant.
                            T rNorm = std::sqrt(rSqr);
                            force += r/(rNorm*rCritical*rCriticalSqr*rCriticalSqr*rCriticalSqr) * forceAmplitude;
                        }
                    }
                }
                // Particle acceleration according to Newton's law.
                particle->set_a (
                        particle->get_a() + (gravity+force) * particle->get_invRho() );
            }
        }
    }
    virtual ParticleInteractionForce<T,Descriptor>* clone() const
    {
        return new ParticleInteractionForce<T,Descriptor>(*this);
    }
    virtual void getTypeOfModification(std::vector<modif::ModifT>& modified) const {
        modified[0] = modif::dynamicVariables;
    }
private:
    T forceAmplitude;
    Array<T,3> gravity;
};

返回一个圆形区域。

class CircularInjection {
public:
    CircularInjection ( T radius_, Array<T,3> center_ )
        : radius(radius_),
          center(center_)
    { }
    bool operator()(Array<T,3> const& pos) const {
        return (util::sqr(pos[1]-center[1]) + util::sqr(pos[2]-center[2]))
               < util::sqr(radius);
    }
private:
    T radius;
    Array<T,3> center;
};

建立区域，使其与流场相等。

MultiBlockManagement3D const& latticeManagement(fluidLattice->getMultiBlockManagement());
    MultiBlockManagement3D particleManagement (
            latticeManagement.getSparseBlockStructure(),
            latticeManagement.getThreadAttribution().clone(),
            commEnvelope, latticeManagement.getRefinementLevel() );

    // Allocation of the data for the particles. This only allocates the particle-hash
    // (a grid on which the particles are stored). The particles themselves are injected
    // later on.
    particles = new MultiParticleField3D<ParticleFieldT> (
        particleManagement, defaultMultiBlockPolicy3D().getCombinedStatistics() );

数据处理器的集成。

 // Add a particle on each vertex of the cone wall. These particles don't move, but
    // they interact with the injected particles, so that these rebounce from the wall.
    addWallParticlesGeneric<T,DESCRIPTOR,ParticleFieldT>(*particles, *triangleBd);

    // In the following the data processors for the equations of motion and the
    // interaction terms are manually added to the particle field. In other situations,
    // data processors are added implicitly, such as for example the data processors
    // for the implementation of a boundary condition in the function call
    // boundaryCondition->insert() above.

    // Prepare the arguments to be provided to the data processors: the particle field
    // for Verlet integration, and particles+fluid for the coupling terms.
    std::vector<MultiBlock3D*> particleArg;
    particleArg.push_back(particles);

    std::vector<MultiBlock3D*> particleFluidArg;
    particleFluidArg.push_back(particles);
    particleFluidArg.push_back(fluidLattice);

    // Calls the function advance() on each particle which, for Verlet particles,
    // updates the position according to v(t+0.5)
    integrateProcessingFunctional (
            new AdvanceParticlesEveryWhereFunctional3D<T,DESCRIPTOR>(-1.0),
            fluidLattice->getBoundingBox(), particleArg, 0);
    // Compute fluid->particle force (which includes friction).
    integrateProcessingFunctional (
            new FluidToParticleCoupling3D<T,DESCRIPTOR>(1.),
            fluidLattice->getBoundingBox(), particleFluidArg, 1 );
    // Compute particle->particle and gravity forces, with help of the data processor
    // ParticleInteractionForce defined above.
    integrateProcessingFunctional (
            new ParticleInteractionForce<T,DESCRIPTOR>(forceAmplitude, gravity),
            fluidLattice->getBoundingBox(), particleArg, 1 );
    // Integrate the velocity according to the Verlet algorithm.
    integrateProcessingFunctional (
            new VerletUpdateVelocitySelective3D<T,DESCRIPTOR>(new util::SelectLargerEqualInt(1)),
            fluidLattice->getBoundingBox(), particleArg, 1 );

    // Definition of a domain from which particles will be injected in the flow field.
    //   The specific domain is close to the inlet of the tube.
    injectionDomain = Box3D(fluidLattice->getBoundingBox());
    injectionDomain.x0 = inletCenter[0] +3;
    injectionDomain.x1 = inletCenter[0] +5;

    // Definition of an absorbtion domains for the particles, one at the bottom and one at the top.
    Box3D absorbtionDomain1(fluidLattice->getBoundingBox());
    absorbtionDomain1.x0 = outletCenter[0] -2;
    absorbtionDomain1.x1 = outletCenter[0] -2;
    Box3D absorbtionDomain2(fluidLattice->getBoundingBox());
    absorbtionDomain2.x0 = inletCenter[0] +1;
    absorbtionDomain2.x1 = inletCenter[0] +1;

    // Functional which absorbs the particles in the specified domains, to avoid that they leave the
    // computational domain.
    integrateProcessingFunctional (
            new AbsorbParticlesFunctionalSelective3D<T,DESCRIPTOR>(new util::SelectLargerEqualInt(1)),
            absorbtionDomain1, particleArg, 0 );
    integrateProcessingFunctional (
            new AbsorbParticlesFunctionalSelective3D<T,DESCRIPTOR>(new util::SelectLargerEqualInt(1)),
            absorbtionDomain2, particleArg, 0 );

    // Execute all data processors once to start the simulation off with well-defined initial values.
    particles->executeInternalProcessors();

具体调用，定义particles模板。

// Create a template for the particles to be injected.
    VerletParticle3D<T,DESCRIPTOR> *particleTemplate=0;
    particleTemplate = new VerletParticle3D<T,DESCRIPTOR>(1, Array<T,3>(0.,0.,0.));
    particleTemplate->setFluidCompliance(fluidCompliance);

if (i>=startParticleIter) {
            std::vector<MultiBlock3D*> particleArg;
            particleArg.push_back(particles);
            // Progressively increase the injection rate.
            T probability = particleProbabilityPerCell;
            //T probability = particleProbabilityPerCell*(1.+(T)i/6000.);
            Particle3D<T,DESCRIPTOR>* particle = particleTemplate->clone();
            // Attribute a color to the first particles by changing their tags (the tags
            // are written in the VTK file).
            if (i<=startParticleIter+20000) {
                particle->setTag(i-startParticleIter);
            }
            applyProcessingFunctional (
                    new AnalyticalInjectRandomParticlesFunctional3D<T,DESCRIPTOR,CircularInjection> (
                        particle, probability, CircularInjection(radius/2.0, inletCenter) ),
                    injectionDomain, particleArg );
        }

 if (i>=startParticleIter) {
            // Iterate the particles..
            particles->executeInternalProcessors();
        }

小总结

定义时选定为DenseParticleField3D，随后的AdvanceParticlesEveryWhereFunctional3D会执行particleField.advanceParticles(particleField.getBoundingBox(), cutOffValue);，找到Dense Particles对应的源码如下：

template<typename T, template<typename U> class Descriptor>
void DenseParticleField3D<T,Descriptor>::advanceParticles(Box3D domain, T cutOffSpeedSqr) {
    Box3D finalDomain;
    if( intersect(domain, particleGrid.getBoundingBox(), finalDomain) )
    {
        std::vector<std::pair<Dot3D,Particle3D<T,Descriptor>*> > nextCellParticles;
        for (plint iX=finalDomain.x0; iX<=finalDomain.x1; ++iX) {
            for (plint iY=finalDomain.y0; iY<=finalDomain.y1; ++iY) {
                for (plint iZ=finalDomain.z0; iZ<=finalDomain.z1; ++iZ) {
                    std::vector<Particle3D<T,Descriptor>*>& particles = particleGrid.get(iX,iY,iZ);
                    std::vector<Particle3D<T,Descriptor>*> newLocalParticles;
                    for (pluint iParticle=0; iParticle<particles.size(); ++iParticle) {
                        Particle3D<T,Descriptor>* particle = particles[iParticle];
                        Array<T,3> oldPos( particle->getPosition() );
                        particle->advance();
                        if (cutOffSpeedSqr>=T() && normSqr(oldPos-particle->getPosition()) < cutOffSpeedSqr)
                        {
                            delete particle;
                        }
                        else {
                            plint newX, newY, newZ;
                            this->computeGridPosition(particle->getPosition(), newX, newY, newZ);
                            if (newX==iX && newY==iY && newZ==iZ) {
                                newLocalParticles.push_back(particle);
                            }
                            else {
                                if (contained(newX,newY,newZ, finalDomain)) {
                                    nextCellParticles.push_back(std::make_pair(Dot3D(newX,newY,newZ),particle));
                                }
                                else {
                                    delete particle;
                                }
                            }
                        }
                    }
                    newLocalParticles.swap(particles);
                }
            }
        }
        for (pluint i=0; i<nextCellParticles.size(); ++i) {
            plint newX = nextCellParticles[i].first.x;
            plint newY = nextCellParticles[i].first.y;
            plint newZ = nextCellParticles[i].first.z;
            particleGrid.get(newX,newY,newZ).push_back(nextCellParticles[i].second);
        }
    }
}

遍历指定区域lattice格点，根据区域lattice坐标得到particles格点坐标，遍历该坐标下所有particles，利用advance获得新位置，接下来判断particle新位置是否合理，如果新位置就是该格点则不变，如果位置超出则被删除，如果在指定区域内，则利用nextCellParticles存储新位置，它的构成是<Dot3D,Particle3D<T,Descriptor>*> ，first是一个Dot3D，second是Particles3D。

最后的操作也很直观，即新particle区域位置设定一个particle。particleGrid.get(newX,newY,newZ).push_back(nextCellParticles[i].second);

后续是VerletUpdateVelocitySelective3D：

template<typename T, template<typename U> class Descriptor>
void VerletUpdateVelocitySelective3D<T,Descriptor>::processGenericBlocks (
        Box3D domain, std::vector<AtomicBlock3D*> blocks )
{
    PLB_PRECONDITION( blocks.size()==1 );
    ParticleField3D<T,Descriptor>& particleField =
        *dynamic_cast<ParticleField3D<T,Descriptor>*>(blocks[0]);

    std::vector<Particle3D<T,Descriptor>*> found;
    particleField.findParticles(domain, found);

    for (pluint iParticle=0; iParticle<found.size(); ++iParticle) {
        Particle3D<T,Descriptor>* nonTypeParticle = found[iParticle];
        if ((*tags)(nonTypeParticle->getTag())) {
            VerletParticle3D<T,Descriptor>* particle =
                dynamic_cast<VerletParticle3D<T,Descriptor>*>(nonTypeParticle);
            PLB_ASSERT( particle );

            Array<T,3> a(particle->get_a());
            particle->set_v(particle->get_vHalfTime() + (T)0.5*a);
            PLB_ASSERT( norm(particle->get_v()) < 1. );
        }
    }
}

其主要作用：

Array<T,3> a(particle->get_a());
particle->set_v(particle->get_vHalfTime() + (T)0.5*a);

这部分作用与particle3D.hh如下代码几乎一样，比起整个advance，实际上略去了更新坐标的内容，毕竟在这个算例里坐标更新是根据AdvanceParticlesEveryWhereFunctional3D来进行的。

template<typename T, template<typename U> class Descriptor>
void VerletParticle3D<T,Descriptor>::advance() {
    vHalfTime = v + (T)0.5*a;
    this->getPosition() += vHalfTime;
}

而这里的a，在particle3D.hh里有如下调用：

Array<T,3> force( (fluidVelocity-this->get_v()) *fluidCompliance);
    this->set_a(force / this->get_rho());

可以发现这里根据速度，fluidCompliance，rho一起计算的数值。

在particlesInCone算例注释里提到advance是根据v(t+0.5)来更新位置的，我想就是这里的部分吧，在算例自带的数据处理器也有set_a部分，为：

 particle->set_a (
particle->get_a() + (gravity+force) * particle->get_invRho() );

因为是rho*g的形式，我想它应该是作为体积力产生作用的，我最近没有找到particles部分对应的参考文献，所以具体理论暂时不清楚。

这个算例主要的不同点是它自己写了一个particles之间相互作用的数据处理器。

    integrateProcessingFunctional (
            new ParticleInteractionForce<T,DESCRIPTOR>(forceAmplitude, gravity),
            fluidLattice->getBoundingBox(), particleArg, 1 );
    // Integrate the velocity according to the Verlet algorithm.

毕竟是需要包含重力影响，所以就没有使用src里数据处理器的坐标更新方式。

附录：关于a的计算

通过velocity和rhobarJ都可以获得a。

template<typename T, template<typename U> class Descriptor>
void VerletParticle3D<T,Descriptor>::velocityToParticle(TensorField3D<T,3>& velocityField, T scaling)
{
    Array<T,3> position(this->getPosition());
    std::vector<Dot3D> pos(8);
    std::vector<T> weights(8);
    linearInterpolationCoefficients(velocityField, position, pos, weights);
    Array<T,3> fluidVelocity;
    fluidVelocity.resetToZero();
    for (plint iCell=0; iCell<8; ++iCell) {
        fluidVelocity += weights[iCell]*velocityField.get(pos[iCell].x,pos[iCell].y,pos[iCell].z)*scaling;
    }

    Array<T,3> force( (fluidVelocity-this->get_v()) *fluidCompliance);
    this->set_a(force / this->get_rho());
}

template<typename T, template<typename U> class Descriptor>
void VerletParticle3D<T,Descriptor>::velocityToParticle(NTensorField3D<T>& velocityField, T scaling)
{
    Array<T,3> position(this->getPosition());
    std::vector<Dot3D> pos(8);
    std::vector<T> weights(8);
    linearInterpolationCoefficients(velocityField, position, pos, weights);
    Array<T,3> fluidVelocity;
    fluidVelocity.resetToZero();
    for (plint iCell=0; iCell<8; ++iCell) {
        T* data = velocityField.get(pos[iCell].x,pos[iCell].y,pos[iCell].z);
        fluidVelocity[0] += weights[iCell]*scaling * data[0];
        fluidVelocity[1] += weights[iCell]*scaling * data[1];
        fluidVelocity[2] += weights[iCell]*scaling * data[2];
    }

    Array<T,3> force( (fluidVelocity-this->get_v()) *fluidCompliance);
    this->set_a(force / this->get_rho());
}

它的具体过程是先通过流场速度插值计算出流速，如果是velocity方式则接着乘以一个系数scaling，最后乘以fluidCompliance，将其设定为a。

至于为什么要这样做，大家可以自行搜索verlet算法原理。