parallelism problems about the creation of lattice

Hi, I am still new to Palabos. I am running the “permeability” example with a big 3D matrix, 375x375x375. My work station has 32 cores.

When I run as “mpirun -np 6”, the program runs with 6 working cores and everything works fine.

When I run as “mpirun -np 8”, the program still runs with 8 working cores successfully, but it will be stuck during the procedure “The creation of lattice”. It seems it takes forever to create the 3D lattice by using 8 cores. It cannot proceed to the next procedures.

Does anyone know what the problem is? Thanks a lot!


thank you for reporting this potential bug. I would like to ask you a couple of questions to understand more about the problem.

You use your own data I guess. What I mean is that you created your own “dat” file which you give as the first command-line argument to the executable. Is this the case?

If yes, then have you made sure that the “dimensions” of your data and the command-line parameters (375 x 375 x 375) are the same?

When you say that with 6 cores everything works fine, you also mean that the solution also looks fine, or just that the code does not hang?

Thanks a lot!

Thanks for reply.

Yes, I create my own dat file and the size is 375x375x375. When I run 6 cores, everything works fine and the result is good. When I run 8 cores, the program is running with 8 cores, but stuck during the procedure “Creation of the Lattice”. My command input is correct.

In the tutorial permeability folder, I also tried “TwoSpheres” example, whose size is only 64x64x48. Still 6 cores works while 8 cores don’t.

I tested the Tutorial 1.6: 3D Channel flow. 8 cores works perfectly. I can even run it and get the result successfully by using 24 cores.

is there any bug in the “permeability.cpp”? I suspect the problem comes from reading the data:

  for (plint iX=0; iX<nx-1; ++iX) {
        string fname = createFileName("slice_", iX, 4)+"_truc.dat";
        pcout << "Reading slice " << fname;
        plb_ifstream geometryFile(fNameIn.c_str());
        if (!geometryFile.is_open()) {
                pcout << "Error: could not open geometry file " << fNameIn << endl;
                return -1;
        geometryFile >> slice;
        copy(slice, slice.getBoundingBox(), geometry, Box3D(iX,iX,0,ny-1,0,nz-1));


Problem solved! The input codes in “permeability.cpp” do have problems for parallelism. We should change it to the following format:

plb_ifstream ifile(“energy.dat”);
if (ifile.is_open()) {
plint nx, ny;
ifile >> nx >> ny;
// Broadcast the input data to all processors; otherwise, the behavior
// of the program is undefined.
global::mpi().bCast(&nx, 1);
global::mpi().bCast(&ny, 1);
MultiScalarField2D energy(nx,ny);
ifile >> energy;


let me thank you once again for this bug report. There was indeed a bug in the permeability.cpp tutorial code. It was about reading the porous medium geometry input file. The following corrected code will be included in the next Palabos release:

/* This file is part of the Palabos library.
 * Copyright (C) 2011-2016 FlowKit Sarl
 * Route d'Oron 2
 * 1010 Lausanne, Switzerland
 * E-mail contact:
 * The most recent release of Palabos can be downloaded at 
 * <>
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 * The library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * GNU Affero General Public License for more details.
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <>.

/* Main author: Wim Degruyter */

#include "palabos3D.h"
#include "palabos3D.hh"

#include <vector>
#include <cmath>
#include <cstdlib>

using namespace plb;

typedef double T;
#define DESCRIPTOR descriptors::D3Q19Descriptor

// This function object returns a zero velocity, and a pressure which decreases
//   linearly in x-direction. It is used to initialize the particle populations.
class PressureGradient {
    PressureGradient(T deltaP_, plint nx_) : deltaP(deltaP_), nx(nx_)
    { }
    void operator() (plint iX, plint iY, plint iZ, T& density, Array<T,3>& velocity) const
        density = (T)1 - deltaP*DESCRIPTOR<T>::invCs2 / (T)(nx-1) * (T)iX;

    T deltaP;
    plint nx;

void readGeometry(std::string fNameIn, std::string fNameOut, MultiScalarField3D<int>& geometry)
    const plint nx = geometry.getNx();
    const plint ny = geometry.getNy();
    const plint nz = geometry.getNz();

    Box3D sliceBox(0,0, 0,ny-1, 0,nz-1);
    std::auto_ptr<MultiScalarField3D<int> > slice = generateMultiScalarField<int>(geometry, sliceBox);
    plb_ifstream geometryFile(fNameIn.c_str());
    for (plint iX=0; iX<nx-1; ++iX) {
        if (!geometryFile.is_open()) {
            pcout << "Error: could not open geometry file " << fNameIn << std::endl;
        geometryFile >> *slice;
        copy(*slice, slice->getBoundingBox(), geometry, Box3D(iX,iX, 0,ny-1, 0,nz-1));

        VtkImageOutput3D<T> vtkOut("porousMedium", 1.0);
        vtkOut.writeData<float>(*copyConvert<int,T>(geometry, geometry.getBoundingBox()), "tag", 1.0);

        std::auto_ptr<MultiScalarField3D<T> > floatTags = copyConvert<int,T>(geometry, geometry.getBoundingBox());
        std::vector<T> isoLevels;
        typedef TriangleSet<T>::Triangle Triangle;
        std::vector<Triangle> triangles;
        Box3D domain = floatTags->getBoundingBox().enlarge(-1);
        isoSurfaceMarchingCube(triangles, *floatTags, isoLevels, domain);
        TriangleSet<T> set(triangles);
        std::string outDir = fNameOut + "/";
        set.writeBinarySTL(outDir + "porousMedium.stl");

void porousMediaSetup(MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
        OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition,
        MultiScalarField3D<int>& geometry, T deltaP)
    const plint nx = lattice.getNx();
    const plint ny = lattice.getNy();
    const plint nz = lattice.getNz();

    pcout << "Definition of inlet/outlet." << std::endl;
    Box3D inlet (0,0, 1,ny-2, 1,nz-2);
    boundaryCondition->addPressureBoundary0N(inlet, lattice);
    setBoundaryDensity(lattice, inlet, (T) 1.);

    Box3D outlet(nx-1,nx-1, 1,ny-2, 1,nz-2);
    boundaryCondition->addPressureBoundary0P(outlet, lattice);
    setBoundaryDensity(lattice, outlet, (T) 1. - deltaP*DESCRIPTOR<T>::invCs2);

    pcout << "Definition of the geometry." << std::endl;
    // Where "geometry" evaluates to 1, use bounce-back.
    defineDynamics(lattice, geometry, new BounceBack<T,DESCRIPTOR>(), 1);
    // Where "geometry" evaluates to 2, use no-dynamics (which does nothing).
    defineDynamics(lattice, geometry, new NoDynamics<T,DESCRIPTOR>(), 2);

    pcout << "Initilization of rho and u." << std::endl;
    initializeAtEquilibrium( lattice, lattice.getBoundingBox(), PressureGradient(deltaP, nx) );

    delete boundaryCondition;

void writeGifs(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)
    const plint nx = lattice.getNx();
    const plint ny = lattice.getNy();
    const plint nz = lattice.getNz();

    const plint imSize = 600;
    ImageWriter<T> imageWriter("leeloo");

    // Write velocity-norm at x=0.
    imageWriter.writeScaledGif(createFileName("ux_inlet", iter, 6),
            *computeVelocityNorm(lattice, Box3D(0,0, 0,ny-1, 0,nz-1)),
            imSize, imSize );

    // Write velocity-norm at x=nx/2.
    imageWriter.writeScaledGif(createFileName("ux_half", iter, 6),
            *computeVelocityNorm(lattice, Box3D(nx/2,nx/2, 0,ny-1, 0,nz-1)),
            imSize, imSize );

void writeVTK(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)
    VtkImageOutput3D<T> vtkOut(createFileName("vtk", iter, 6), 1.);
    vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", 1.);
    vtkOut.writeData<3,float>(*computeVelocity(lattice), "velocity", 1.);

T computePermeability(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, T nu, T deltaP, Box3D domain )
    pcout << "Computing the permeability." << std::endl;

    // Compute only the x-direction of the velocity (direction of the flow).
    plint xComponent = 0;
    plint nx = lattice.getNx();

    T meanU = computeAverage(*computeVelocityComponent(lattice, domain, xComponent));

    pcout << "Average velocity     = " << meanU                         << std::endl;
    pcout << "Lattice viscosity nu = " << nu                            << std::endl;
    pcout << "Grad P               = " << deltaP/(T)(nx-1)              << std::endl;
    pcout << "Permeability         = " << nu*meanU / (deltaP/(T)(nx-1)) << std::endl;

    return meanU;

int main(int argc, char **argv)
    plbInit(&argc, &argv);

    if (argc!=7) {
        pcout << "Error missing some input parameter\n";
        pcout << "The structure is :\n";
        pcout << "1. Input file name.\n";
        pcout << "2. Output directory name.\n";
        pcout << "3. number of cells in X direction.\n";
        pcout << "4. number of cells in Y direction.\n";
        pcout << "5. number of cells in Z direction.\n";
        pcout << "6. Delta P .\n";
        pcout << "Example: " << argv[0] << " twoSpheres.dat tmp/ 48 64 64 0.00005\n";
        exit (EXIT_FAILURE);
    std::string fNameIn  = argv[1];
    std::string fNameOut = argv[2];

    const plint nx = atoi(argv[3]);
    const plint ny = atoi(argv[4]);
    const plint nz = atoi(argv[5]);
    const T deltaP = atof(argv[6]);


    const T omega = 1.0;
    const T nu    = ((T)1/omega- (T)0.5)/DESCRIPTOR<T>::invCs2;

    pcout << "Creation of the lattice." << std::endl;
    MultiBlockLattice3D<T,DESCRIPTOR> lattice(nx,ny,nz, new BGKdynamics<T,DESCRIPTOR>(omega));
    // Switch off periodicity.

    pcout << "Reading the geometry file." << std::endl;
    MultiScalarField3D<int> geometry(nx,ny,nz);
    readGeometry(fNameIn, fNameOut, geometry);

    pcout << "nu = " << nu << std::endl;
    pcout << "deltaP = " << deltaP << std::endl;
    pcout << "omega = " << omega << std::endl;
    pcout << "nx = " << lattice.getNx() << std::endl;
    pcout << "ny = " << lattice.getNy() << std::endl;
    pcout << "nz = " << lattice.getNz() << std::endl;

    porousMediaSetup(lattice, createLocalBoundaryCondition3D<T,DESCRIPTOR>(), geometry, deltaP);

    // The value-tracer is used to stop the simulation once is has converged.
    // 1st parameter:velocity
    // 2nd parameter:size
    // 3rd parameters:threshold
    // 1st and second parameters ae used for the length of the time average (size/velocity)
    util::ValueTracer<T> converge(1.0,1000.0,1.0e-4);

    pcout << "Simulation begins" << std::endl;
    plint iT=0;

    const plint maxT = 30000;
    for (;iT<maxT; ++iT) {
        if (iT % 20 == 0) {
            pcout << "Iteration " << iT << std::endl;
        if (iT % 500 == 0 && iT>0) {


        if (converge.hasConverged()) {

    pcout << "End of simulation at iteration " << iT << std::endl;

    pcout << "Permeability:" << std::endl << std::endl;
    computePermeability(lattice, nu, deltaP, lattice.getBoundingBox());
    pcout << std::endl;

    pcout << "Writing VTK file ..." << std::endl << std::endl;
    pcout << "Finished!" << std::endl << std::endl;

    return 0;


Hey Dimitris,
I manually updated by permeability.cpp file with the information you provided above, however when I compile and run it in parallel it only seems to print out the first layer of the lattice (x=1) and then solve with rest of the domain empty. I’ve had this with multiple .dats that have run successfully on the old permeability code just not in parallel.
Any help would be greatly appreciated,

Dear Matt,

I really apologize for this, but the code above still had a bug! I edited and corrected my post above (which was created on October 28, 2015 05:10PM). Now the code should be fine and you should be able to run it in parallel with no problem. In any case, please let me know if everything is OK or if your problem persists.

Again, I apologize for the inconvenience.


Hey Dimitris,
The modification you made to the code seems to have worked, it solved as far as I could tell correctly. I also really liked that you added the extraction of geometry in stl and vti format. Also I’d like to compliment you on the incredibly fast response, I didn’t think you’d post the solution to the problem with in twelve hours!
Thanks so much,

Hey Dimitris,

When I use the permeability.cpp in 300X300X300 model, the following error informatio came out. I really appreciate your help.

Creation of the lattice.
terminate called after throwing an instance of ‘std::bad_alloc’
what(): std::bad_alloc
Aborted (core dumped)

Hey Matthew Yu,

When I use the permeability.cpp in 200X200X300 model, the following error informatio came out. I really appreciate your help.

Creation of the lattice.
terminate called after throwing an instance of ‘std::bad_alloc’
what(): std::bad_alloc

But when I use the permeability.cpp in 200X200X200 model, everything is OK.