Dear All,
I am a new Palabos user. I am doing single phase simulations for random pack porous media. I need to calculate the porosity and tortuosity. I write two data processors. For the porosity calculations, I just count the voxel whose values are 0 (0, 1 and 2 represents the pore voxel, grain boundary and inactive voxels, respectively). For the tortuosity calculations, I sum the velocity vector in x direction.
//For porosity calculation
template
class Porosity: public BoxProcessingFunctional3D_S
{
public :
Porosity(plint& pcount_) : pcount(pcount_)
{ };
virtual void process (Box3D domain, ScalarField3D<T>& mediadomain)
{
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) {
if (mediadomain(iX,iY,iZ)==0) {
pcount=pcount+1;
}
}
}
}
};
virtual Porosity<T>* clone() const
{
return new Porosity<T>(*this);
};
private :
plint pcount;
};
//For tortuosity calculation
template<typename T, template class Descriptor>
class VelocityTortuosity : public BoxProcessingFunctional3D_L<T,Descriptor>
{
public :
VelocityTortuosity (T& sumvelocity_) : sumvelocity(sumvelocity_)
{ };
virtual void process (Box3D domain, Blocklattice3D<T,Descriptor>& lattice)
{
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) {
Array<T,3> U;
lattice.get(iX,iY,iZ).computeVelocity(U);
sumvelocity = sumvelocity + sqrt(U[0]U[0]);
}
}
}
};
virtual VelocityTortuosity<T, Descriptor> clone() const
{
return new VelocityTortuosity<T, Descriptor>(*this);
};
private :
T sumvelocity;
};
Before callcing the processors, I defined the variables,
plint pcount = 0;
T sumvelocity = 0.0;
Then I call the processors.
applyProcessingFunctional(new Porosity(sumvelocity),
lattice.getBoundingBox(), mediadomain);
applyProcessingFunctional(new VelocityTortuosity<T,DESCRIPTOR>(sumvelocity),
lattice.getBoundingBox(),lattice);
After that I calculate the porosity and tortuosity as;
T Porosity= (T) pcount/T (nxnynz);
T averageUnorm = computeAverage( *computeVelocityNorm (lattice, Box3D(1, nx-2, 1, ny-2, 1, nz-2)) );
T tortuosity = averageUnorm * averageUnorm * nx * nx * ny * ny * nz * nz/(sumvelocity * sumvelocity);
When I try to compile, I have recieved an error. I cannot figure out what is the problem with the data processors. Actually, I try to write looking the example files. I would appreciate if you could inform me about what I doing wrong. Thanks.
Best Regards,
Asya