Data processor for non-local operation

Palabos General
#1

Dear all,

I need to exchange the densities and velocities in Plane1 (0, iY, iZ) and Plane2 (iY, 0, iZ), and then re-initialize the off-lattice at every time step.

Firstly, I wrote a loop and perform initializeAtEquilibrium at every cell in these two planes. It works, but slow.

Then, I used data processor with the following code

[code=“cpp”]
template<typename T, template class Descriptor>
class BC
{
public:
BC(MultiBlockLattice3D<T,Descriptor>& lattice_)
: lattice(lattice_)
{ }
void operator()(plint iX, plint iY, plint iZ, T& rho, Array<T,3>& u) const
{
rho = lattice.get(iY,iX,iZ).computeDensity();
lattice.get(iY,iX,iZ).computeVelocity(u);
}
private:
MultiBlockLattice3D<T,Descriptor>& lattice;
};



and 
[code="cpp"]
initializeAtEquilibrium (*lattice, Plane1, BC<T,DESCRIPTOR>(*lattice) );

The code can be compiled successfully but will lead to MPI error during simulation. I think this is because operation can’t be non-local by more than one cell. Am I right? Then, how can I fix this problem? Is there any other method to achieve my purpose?

Thanks.

yong

#2

Dear Yong,

The reason why your code is slow and fails in parallel is that it violates a principle of parallelism in Palabos. The code contained in the operator() is executed in parallel: the loop over iX, iY, iZ is split up and distributed to several processors. This code is then no longer allowed to access the interface of a MultiBlockLattice3D, because when it does, Palabos doesn’t know how to parallelize the operation.

There exists no predefined class in Palabos that does exactly what you need. It is however straightforward to write your own. Here is the code:

[code=“cpp”]
template<typename T, template class Descriptor>
class SetToLocalEquilibrium3D : public BoxProcessingFunctional3D_L<T,Descriptor>
{
virtual void process(Box3D domain, BlockLattice3D<T,Descriptor>& lattice) {
typedef Descriptor D;
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) {
Cell<T,Descriptor>& cell = lattice.get(iX,iY,iZ);
T rhoBar; Array<T,3> j; Array<T,D::q> fEq;
cell.getDynamics().computeRhoBarJ(cell, rhoBar, j);
T jSqr = normSqr(j);
cell.getDynamics().computeEquilibria(fEq, rhoBar, j, jSqr);
for (plint iPop=0; iPop<D::q; ++iPop) {
cell[iPop] = fEq[iPop];
}
}
}
}
}
void getTypeOfModification(std::vectormodif::ModifT& modified) const {
modified[0] = modif::staticVariables;
}
virtual SetToLocalEquilibrium3D<T,Descriptor>* clone() const {
return new SetToLocalEquilibrium3D<T,Descriptor>(*this);
}
};



which you call by writing

[code="cpp"]
applyProcessingFunctional(new SetToLocalEquilibrium3D<T,DESCRIPTOR>, domain, lattice);

Cheers,
Jonas

#3

Dear Jonas,

Thanks very much for your reply.

Yes, I noticed that after I read some stuff about data processor and wrote another code with BoxProcessingFunctional3D_LL and applyProcessingFunctional. I used iniCellAtEquilibrium(cell, density, velocity) after I exchanged the density and velocity. The code can be compiled and run without error. However, the results (density and velocity fields) are incorrect and the same as simulation without this boundary condition. In other words, data processor doesn’t work.

Moreover, I considered your suggestion and used cell.getDynamics().computeEquilibria. There is no change.

The code I have tried are:

[code=“cpp”]
template<typename T, template class Descriptor>
class PeriodicBCInletProcessor3D : public BoxProcessingFunctional3D_LL<T,Descriptor,T,Descriptor>
{
virtual void process(Box3D domain, BlockLattice3D<T,Descriptor>& lattice,
BlockLattice3D<T,Descriptor>& latticeOutlet1)
{
//typedef Descriptor D;
Dot3D absoluteOffset = lattice.getLocation();
for (plint iX=domain.x0; iX<=domain.x1; ++iX)
{
plint absoluteX = absoluteOffset.x + iX;
for (plint iY=domain.y0; iY<=domain.y1; ++iY)
{
plint absoluteY = absoluteOffset.y + iY;
for (plint iZ=domain.z0; iZ<=domain.z1; ++iZ)
{
plint absoluteZ = absoluteOffset.z + iZ;
Cell<T, Descriptor>& cell = lattice.get(iX,iY,iZ);
Cell<T, Descriptor>& cell1 = latticeOutlet1.get(absoluteY,absoluteX+1,absoluteZ);

		  T density = cell1.computeDensity(); 
                      Array<T, 3> velocity;
		  cell1.computeVelocity(velocity); 
		  cell.defineDensity(density);
		  cell.defineVelocity(velocity);
                      iniCellAtEquilibrium(cell, density, velocity);					
				
                      //T rhoBar;
		  //Array<T,3> j;
		  //Array<T,D::q> fEq;
                      //cell1.getDynamics().computeRhoBarJ(cell, rhoBar, j);
                      //T jSqr = normSqr(j);
                      //cell.getDynamics().computeEquilibria(fEq, rhoBar, j, jSqr);
                      //for (plint iPop=0; iPop<D::q; ++iPop) 
		  //{
                      //    cell[iPop] = fEq[iPop];
                      //}
                 }
          }
    }
}
virtual PeriodicBCInletProcessor3D<T,Descriptor>* clone() const
{
    return new PeriodicBCInletProcessor3D<T,Descriptor>(*this);
}
virtual void getModificationPattern(std::vector<bool>& isWritten) const
{
    isWritten[0] = true;
    isWritten[1] = false;
}		
virtual void getTypeOfModification(std::vector<modif::ModifT>& modified) const 
{
    modified[0] = modif::staticVariables;
}
virtual BlockDomain::DomainT appliesTo() const 
{
    return BlockDomain::bulkAndEnvelope;
}

};



and call it by
[code="cpp"]
void setPeriodicBC (MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
		            MultiBlockLattice3D<T,DESCRIPTOR>& latticeInlet1,
		            MultiBlockLattice3D<T,DESCRIPTOR>& latticeOutlet1,		
		            Box3D inletPlane, Box3D outletPlane)
{
	    applyProcessingFunctional(
				new PeriodicBCInletProcessor3D<T,DESCRIPTOR>, 
		                inletPlane, lattice, latticeOutlet1);
	    //applyProcessingFunctional(
		//		new PeriodicBCOutletProcessor3D<T,DESCRIPTOR>, 
		//                outletPlane, lattice, latticeInlet1);		
}

Thanks.

#4

Hi Jonas,

I just find out that BoxProcessingFunctional3D_LL doesn’t work. BoxProcessingFunctional3D_L works but can’t finish what I need as I have to exchange data on different atomic-lattices. Is there any method to fix this problem?

Thanks.

yong

#5

Dear Yong,

A BoxProcessingFunctional3D_L is intended to act on a single lattice, while a BoxProcessingFunctional3D_LL acts at the same time on two lattices, for example to transfer information from one lattice to the other. As examples that are potentially helpful to you, a fully working instance of the BoxProcessingFunctional3D_LL can be found in palabos/src/multiPhysics/boussinesqThermalProcessor3D.h, where it is used to couple a fluid with the temperature field, and other examples are in palabos/src/dataProcessors/dataAnalysisFunctional3D.h, where the data processors are used to transfer data from one lattice to another.

Hope this helps!
Cheers,
Jonas

#6

Hi Jonas,

Yes, I knew these two examples. I also found out the bug in my previous code with BoxProcessingFunctional3D_LL. As for applyProcessingFunctional( functional, DOMAIN, LATTICE1, LATTICE2), seems the domain corresponding the intersection of LATTICE1 and LATTICE2 should be a subset of DOMAIN.

My purpose is to write a data processor for a special periodic BC. It’s different from the previous one in Palabos as the periodic boundaries in my case are not parallel. I have to exchange the data at (iX, iY, iZ) in the inlet plane with the one at (iY, iX, iZ) in the outlet plane.

I planed to use LATTICE1 as the full block, and LATTICE2 a sub-block (such as my outlet plane) therein. If (iX, iY, iZ) is a cell in an atomic-block of LATTICE1, (iY, iX, iZ) is then out of this atomic-block and should be in LATTICE2. Again, code can be compiled but still have mpi error during simulation. Is there any data processor to transfer data from one atomic-block to another one?

Thanks.

Yong