Simulation hanging before entering main loop

[gaturpin]

Hey, I'm having an issue where my simulation hangs during the setup, before entering the main loop when using refinement = static. The sim doesn't crash or print any errors, it just keeps running without making it to the main loop.
I don't run into any problems with refinement = none however.

I don't believe this is a bug with Athena as I can run the provided pgens with static mesh refinement without issue (e.g. slotted_cylinder). I also don't think this a problem with the cluster I am using as I see the same behaviour on another cluster.

Some more info:
Compilation/config command:
python configure.py --prob pgen --eos adiabatic -mpi -hdf5 --hdf5_path=/share/apps/centos7/hdf5-parallel/1.8.18/intel --mpiccmd h5pcc --cflag="-diag-disable 3180"

Setup loop:

void MeshBlock::ProblemGenerator(ParameterInput *pin) {
  Real v1(0.0), v2(0.0), v3(0.0);
  Real x_cell(0.0), y_cell(0.0), z_cell(0.0);
  Real r_cell[2], theta_cell[2], ax_cell[2], ay_cell[2];
  ax_cell[0] = ax_cell[1] = ay_cell[0] = ay_cell[1] = 0.0;
  Real ax_ind(0.0), ay_ind(0.0);
  Real disc_mass;
  //  Initialize density and momenta
  for (int k=ks; k<=ke; ++k) {
    for (int j=js; j<=je; ++j) {
      for (int i=is; i<=ie; ++i) {
        GetCartCoord(pcoord,x_cell,y_cell,z_cell,i,j,k);
        phydro->u(IDN,k,j,i) = DenProfile(x_cell,y_cell,z_cell); 
        VelProfile(x_cell,y_cell,z_cell,v1,v2,v3);

        phydro->u(IM1,k,j,i) = phydro->u(IDN,k,j,i)*v1; // Store momentum
        phydro->u(IM2,k,j,i) = phydro->u(IDN,k,j,i)*v2;
        phydro->u(IM3,k,j,i) = phydro->u(IDN,k,j,i)*v3;

        if (NON_BAROTROPIC_EOS) {
          Real cs_squared = SoundSpeedSQR(x_cell,y_cell,z_cell);
          phydro->u(IEN,k,j,i) = cs_squared*phydro->u(IDN,k,j,i)/(gamma_gas - 1.0); // Store energy
          phydro->u(IEN,k,j,i) += 0.5*(SQR(phydro->u(IM1,k,j,i))+SQR(phydro->u(IM2,k,j,i))
                                       + SQR(phydro->u(IM3,k,j,i)))/phydro->u(IDN,k,j,i);
        }
        disc_mass += (x1size/nx1)*(x2size/nx2)*(x3size/nx3)*phydro->u(IDN,k,j,i);
      }
    }
  }

#ifdef MPI_PARALLEL
  {
    if (Globals::my_rank == 0) {
      MPI_Reduce(&disc_mass, &total_disc_mass, 1, MPI_ATHENA_REAL, MPI_SUM, 0,
                 MPI_COMM_WORLD);

    } 
    else {
      MPI_Reduce(&disc_mass, &total_disc_mass, 1, MPI_ATHENA_REAL, MPI_SUM, 0,
                 MPI_COMM_WORLD);

    }
  }  
  MPI_Bcast(&total_disc_mass,1,MPI_ATHENA_REAL,0,MPI_COMM_WORLD);

  MPI_Barrier(MPI_COMM_WORLD);
#endif

  /*if (Globals::my_rank == 0){
    std::cout << Globals::my_rank << " After MPI Barrier " << std::endl;
  }*/

  return;
}

Boundary conditions are of the form:

void DiskInnerX1(MeshBlock *pmb,Coordinates *pco, AthenaArray<Real> &prim, FaceField &b,
                 Real time, Real dt,
                 int il, int iu, int jl, int ju, int kl, int ku, int ngh) {
  Real x(0.0), y(0.0), z(0.0);
  Real v1(0.0), v2(0.0), v3(0.0);


  //std::cout << "Boundary 1 Inner" << std::endl;


  for (int k=kl; k<=ku; ++k) {
    for (int j=jl; j<=ju; ++j) {
      for (int i=1; i<=ngh; ++i) {

        GetCartCoord(pco,x,y,z,il-i,j,k);
        prim(IDN,k,j,il-i) = DenProfile(x,y,z);
        VelProfile(x,y,z,v1,v2,v3);

        prim(IM1,k,j,il-i) = v1;
        prim(IM2,k,j,il-i) = v2;
        prim(IM3,k,j,il-i) = v3;

        if (NON_BAROTROPIC_EOS)
          prim(IEN,k,j,il-i) = SoundSpeedSQR(x, y, z)*prim(IDN,k,j,il-i);
      }
    }
  }
}

During the setup, in my pgen, I manage to get print statements out of:
MeshBlock::ProblemGenerator - Each rank appears to get through without issue as I get print statements after the MPI_Barrier()
Boundary - DiskInnerX1
Boundary - DiskInnerX2
Boundary - DiskOuterX2

But no other boundary functions
I have tried other boundary conditions (periodic rather than user etc.) and I still get hanging sims.

Does anyone have any ideas?

Cheers!
George

Version info

• Athena++ version: 21.0
• Compiler and version: icc (ICC) 17.0.1
• Operating system: CentOS Linux 7

[c-white]

Adding new MPI commands is somewhat dangerous. In this case, there may well be a problem with the fact that Athena++ has

1 MPI rank <-> 1 Mesh <-> 1 or more MeshBlocks

ProblemGenerator() is called by each MeshBlock on the Mesh in turn. So if you have more MeshBlocks than ranks, the global communication is almost certainly not doing what you expect. If different Meshes have different numbers of MeshBlocks, then this code will hang, as different ranks will call the MPI routines different numbers of times.

The code might work correctly if you have exactly 1 MeshBlock per rank. Otherwise, the communication pattern needs to be

1. Have each MeshBlock set itself up.
2. Have the Mesh aggregate results from all its local MeshBlocks.
3. Have the Meshes communicate with MPI.
4. Have the Mesh distribute whatever data is needed to each local MeshBlock.

[gaturpin]

Ah yes that is the case - matching the number of ranks to the number of MeshBlocks solved this.
Thank you!