Palabos源码：collideAndStream

几乎每个算例里都会有这一行代码，它背后的源代码如下：

/** This operation is more efficient than a successive application of
 * collide(int,int,int,int,int,int) and stream(int,int,int,int,int,int),
 * because memory is traversed only once instead of twice.
 */
template<typename T, template<typename U> class Descriptor>
void BlockLattice3D<T,Descriptor>::collideAndStream(Box3D domain) {
    // Make sure domain is contained within current lattice
    PLB_PRECONDITION( contained(domain, this->getBoundingBox()) );

    global::profiler().start("collStream");
    global::profiler().increment("collStreamCells", domain.nCells());

    static const plint vicinity = Descriptor<T>::vicinity;

    // First, do the collision on cells within a boundary envelope of width
    // equal to the range of the lattice vectors (e.g. 1 for D3Q19)
    collide(Box3D(domain.x0,domain.x0+vicinity-1,
                  domain.y0,domain.y1,
                  domain.z0,domain.z1) );
    collide(Box3D(domain.x1-vicinity+1,domain.x1,
                  domain.y0,domain.y1,
                  domain.z0,domain.z1) );
    collide(Box3D(domain.x0+vicinity,domain.x1-vicinity,
                  domain.y0,domain.y0+vicinity-1,
                  domain.z0,domain.z1) );
    collide(Box3D(domain.x0+vicinity,domain.x1-vicinity,
                  domain.y1-vicinity+1,domain.y1,
                  domain.z0,domain.z1) );
    collide(Box3D(domain.x0+vicinity,domain.x1-vicinity,
                  domain.y0+vicinity,domain.y1-vicinity,
                  domain.z0,domain.z0+vicinity-1) );
    collide(Box3D(domain.x0+vicinity,domain.x1-vicinity,
                  domain.y0+vicinity,domain.y1-vicinity,
                  domain.z1-vicinity+1,domain.z1) );

    // Then, do the efficient collideAndStream algorithm in the bulk,
    // excluding the envelope (this is efficient because there is no
    // if-then-else statement within the loop, given that the boundary
    // region is excluded)
    bulkCollideAndStream(Box3D(domain.x0+vicinity,domain.x1-vicinity,
                               domain.y0+vicinity,domain.y1-vicinity,
                               domain.z0+vicinity,domain.z1-vicinity) );

    // Finally, do streaming in the boundary envelope to conclude the
    // collision-stream cycle
    boundaryStream(domain, Box3D(domain.x0,domain.x0+vicinity-1,
                                 domain.y0,domain.y1,
                                 domain.z0,domain.z1) );
    boundaryStream(domain, Box3D(domain.x1-vicinity+1,domain.x1,
                                 domain.y0,domain.y1,
                                 domain.z0,domain.z1) );
    boundaryStream(domain, Box3D(domain.x0+vicinity,domain.x1-vicinity,
                                 domain.y0,domain.y0+vicinity-1,
                                 domain.z0,domain.z1) );
    boundaryStream(domain, Box3D(domain.x0+vicinity,domain.x1-vicinity,
                                 domain.y1-vicinity+1,domain.y1,
                                 domain.z0,domain.z1) );
    boundaryStream(domain, Box3D(domain.x0+vicinity,domain.x1-vicinity,
                                 domain.y0+vicinity,domain.y1-vicinity,
                                 domain.z0,domain.z0+vicinity-1) );
    boundaryStream(domain, Box3D(domain.x0+vicinity,domain.x1-vicinity,
                                 domain.y0+vicinity,domain.y1-vicinity,
                                 domain.z1-vicinity+1,domain.z1) );
    global::profiler().stop("collStream");
}

/** At the end of this method, finalizeIteration() and
 * executeInternalProcessors() are automatically invoked.
 * \sa collideAndStream(int,int,int,int,int,int) */
template<typename T, template<typename U> class Descriptor>
void BlockLattice3D<T,Descriptor>::collideAndStream() {
    collideAndStream(this->getBoundingBox());

    implementPeriodicity();

    this->executeInternalProcessors();
    this->evaluateStatistics();
    this->incrementTime();
}

collide

template<typename T, template<typename U> class Descriptor>
void BlockLattice3D<T,Descriptor>::collide(Box3D domain) {
    // Make sure domain is contained within current lattice
    PLB_PRECONDITION( contained(domain, this->getBoundingBox()) );

    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) {
                grid[iX][iY][iZ].collide(this->getInternalStatistics());
                grid[iX][iY][iZ].revert();
            }
        }
    }
}

/** \sa collide(int,int,int,int) */
template<typename T, template<typename U> class Descriptor>
void BlockLattice3D<T,Descriptor>::collide() {
    collide(this->getBoundingBox());
}

GuoExternalForceBGKdynamics

template<typename T, template<typename U> class Descriptor>
void GuoExternalForceBGKdynamics<T,Descriptor>::collide (
        Cell<T,Descriptor>& cell,
        BlockStatistics& statistics )
{
    T rhoBar = this->computeRhoBar(cell);
    Array<T,Descriptor<T>::d> u, j;
    this->computeVelocity(cell, u);
    T rho = Descriptor<T>::fullRho(rhoBar);
    for (plint iD = 0; iD < Descriptor<T>::d; ++iD)
    {
        j[iD] = rho * u[iD];
    }
    
    T uSqr = dynamicsTemplates<T,Descriptor>::bgk_ma2_collision(cell, rhoBar, j, this->getOmega());
    externalForceTemplates<T,Descriptor>::addGuoForce(cell, u, this->getOmega(), (T)1);
    
    if (cell.takesStatistics()) {
        gatherStatistics(statistics, rhoBar, uSqr);
    }
}

ShanChenExternalForceBGKdynamics

template<typename T, template<typename U> class Descriptor>
void ShanChenExternalForceBGKdynamics<T,Descriptor>::collide (
        Cell<T,Descriptor>& cell,
        BlockStatistics& statistics )
{
    T rhoBar;
    Array<T,Descriptor<T>::d> j;
    momentTemplates<T,Descriptor>::get_rhoBar_j(cell, rhoBar, j);
    
    T invOmega = 1./this->getOmega();
    Array<T,Descriptor<T>::d> force;
    force.from_cArray(cell.getExternal(Descriptor<T>::ExternalField::forceBeginsAt));
    Array<T,Descriptor<T>::d> jCorrected(j+invOmega*Descriptor<T>::fullRho(rhoBar) * force);
    
    dynamicsTemplates<T,Descriptor>::bgk_ma2_collision(cell, rhoBar, jCorrected, this->getOmega());
    if (cell.takesStatistics()) {
        T uSqr = T();
        T invRho = Descriptor<T>::invRho(rhoBar);
        for (int iD=0; iD<Descriptor<T>::d; ++iD) {
            uSqr += util::sqr(j[iD]*invRho + 0.5*force[iD]);
        }
        gatherStatistics(statistics, rhoBar, uSqr);
    }
}

bulkCollideAndStream

/** This method is fast, but it is erroneous when applied to boundary
 * cells.
 * \sa collideAndStream(int,int,int,int,int,int)
 * \sa collideAndStream()
 */
template<typename T, template<typename U> class Descriptor>
void BlockLattice3D<T,Descriptor>::bulkCollideAndStream(Box3D domain) {
    // Make sure domain is contained within current lattice
    PLB_PRECONDITION( contained(domain, this->getBoundingBox()) );

    // if (Descriptor<T>::q==15 || Descriptor<T>::q==19) {
    if (Descriptor<T>::q==19) {
        // On nearest-neighbor lattice, use the cache-efficient
        //   version of collidAndStream.
        blockwiseBulkCollideAndStream(domain);
    }
    else {
        // Otherwise, use the straightforward implementation.
        //   Note that at some point, we should implement the cache-efficient
        //   version for extended lattices as well.
        linearBulkCollideAndStream(domain);
    }
}

boundaryStream

/** This method is slower than bulkStream(int,int,int,int), because it must
 * be verified which distribution functions are to be kept from leaving
 * the domain.
 * \sa stream(int,int,int,int)
 * \sa stream()
 */
template<typename T, template<typename U> class Descriptor>
void BlockLattice3D<T,Descriptor>::boundaryStream(Box3D bound, Box3D domain) {
    // Make sure bound is contained within current lattice
    PLB_PRECONDITION( contained(bound, this->getBoundingBox()) );
    // Make sure domain is contained within bound
    PLB_PRECONDITION( contained(domain, bound) );

    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) {
                for (plint iPop=1; iPop<=Descriptor<T>::q/2; ++iPop) {
                    plint nextX = iX + Descriptor<T>::c[iPop][0];
                    plint nextY = iY + Descriptor<T>::c[iPop][1];
                    plint nextZ = iZ + Descriptor<T>::c[iPop][2];
                    if ( nextX>=bound.x0 && nextX<=bound.x1 &&
                         nextY>=bound.y0 && nextY<=bound.y1 &&
                         nextZ>=bound.z0 && nextZ<=bound.z1 )
                    {
                        std::swap(grid[iX][iY][iZ][iPop+Descriptor<T>::q/2],
                                  grid[nextX][nextY][nextZ][iPop]);
                    }
                }
            }
        }
    }
}