Palabos程序代码解读

程序代码解读"/>

Palabos程序代码解读

定义部分

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算法原理。