0

I have a mesh file and I did a partitioning of it using METIS(in 4 parts/processes).METIS provided me with the partitioning file of the mesh(gave me a file with the number of process where each element of the mesh belongs to).My job now is to input these information to my parallel code.I tried to do it by letting each process to have access to the same mesh file and read the data that it wants based on partitioning file.

#include <iostream>
#include <fstream>
#include <sstream>
#include "mpi.h"
using namespace std;
//each process stores the partitioning
int* PartitionFile(){
 ifstream file ("epart.txt");
 int NE=14;
 int part,i=0;
 int *partition=new int[14];
 while(file>>part){
   partition[i]=part;
   i++;
 }
 file.close();
 return partition;
}
int FindSizeOfLocalElements(int *epart,int rank){
  int size=0;
  for (int i=0;i<14;i++){
   if(epart[i]==rank){
    size+=1;
   }
  }
  return size;
}
//stores the elements of each process
int * LocalPartition(int* epart,int size,int rank){

  int *localPart=new int[size];
  int j=0;
 for(int i=0;i<14;i++){
    if (epart[i]==rank){
        localPart[j]=i+1;//+1 because elements start from 1(global numbering)
        j+=1;
    }
  }
return localPart;
}
int **ElementConnectivityMeshFile(int* localPart,int size){
  ifstream file ("mesh.txt");
  int node1,node2,node3;
  int elem=1;
  int i=0;
  int **elemConn=new int*[size];
  for(int j=0;j<size;j++){
    elemConn[j]=new int[3];//each element has 3 nodes.Here elements has local numbering.Global numbering is stored in localPart
 }
  while(file>>node1>>node2>>node3){
    if (elem==localPart[i]){
        elemConn[i][0]=node1;
        elemConn[i][1]=node2;
        elemConn[i][2]=node3;
        i+=1;
    }
    elem+=1;
    if(elem>14){break;}
  }
  file.close();
  return elemConn;
}




int main(){
  MPI_Init(NULL, NULL);
  int numProc;
  MPI_Comm_size(MPI_COMM_WORLD, &numProc);
  int rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  int *epart=PartitionFile();
  int size=FindSizeOfLocalElements(epart,rank);
  int *elem=LocalPartition(epart,size,rank);
  int **elemConn=ElementConnectivityMeshFile(elem,size);
  MPI_Finalize();
  return 0;
}

This part of code gives me the desired results,however I want to know how efficient is letting MPI processes read the same file,by using c++ standard functions, and if that can affect scalability and performance.For this demostration i used a mesh of 14 elements and 4 processes.

mesh file
1 3 2
2 3 4
3 5 4
4 5 6
5 7 6
8 7 5
3 8 5
9 7 8
9 8 3
1 9 3
10 9 1
11 10 1
11 1 12
12 1 2

epart file
2
2
0
0
0
1
0
1
1
3
3
3
2
2
spyros
  • 127
  • 5
  • 2
    It very much depends on your filesystem, but having the same file open by large numbers of processes can cause issues with contention. Regardless of this, the internal network bandwidth is normally much higher than IO to disk so, as a first step, reading on a single process and broadcasting the result is better than everyone reading. However, this will only work here for reading epart as I see that processes do not store all of the mesh (and the code probably wouldn't be scalable in memory if they did). The best solution would be to use MPI IO but you'd need to convert flles to binary format. – David Henty May 25 '20 at 08:35
  • So,MPI IO can let me read "random" elements from the mesh file?I'm asking because I am not familiar with MPI IO – spyros May 25 '20 at 12:25
  • Yes - you would use MPI_Type_create_indexed_block as the filetype - https://www.rookiehpc.com/mpi/docs/mpi_type_create_indexed_block.php... I'll see if I can knock up an example, The only restriction compared to normal usage of MPI derived types is that the displacements have to be non-decreasing, i.e. you can't jump around in the file and must pick them off in order which sometimes requires you to do a sort first. – David Henty May 25 '20 at 15:15

1 Answers1

1

I think this program illustrates the basic approach using MPI-IO with binary files:

#include <stdio.h>
#include <mpi.h>

#define NELEM 14
#define NVERT  3

int meshfile[NELEM][NVERT] =
  { { 1,  3,  2},
    { 2,  3,  4},
    { 3,  5,  4},
    { 4,  5,  6},
    { 5,  7,  6},
    { 8,  7,  5},
    { 3,  8,  5},
    { 9,  7,  8},
    { 9,  8,  3},
    { 1,  9,  3},
    {10,  9,  1},
    {11, 10,  1},
    {11,  1, 12},
    {12,  1,  2}, };

int partfile[NELEM] = {2, 2, 0, 0, 0, 1, 0, 1, 1, 3, 3, 3, 2, 2};

int main(void)
{
  int i;

  int part[NELEM];
  int mesh[NELEM][NVERT];

  /* Should really malloc smaller mesh based on local size */

  FILE *fp;

  int rank, size;
  MPI_Comm comm;
  MPI_Status status;
  MPI_File fh;
  MPI_Datatype filetype;
  int disp[NELEM];
  int nelemlocal;

  /* Should really malloc smaller displ based on nelemlocal */

  comm = MPI_COMM_WORLD;

  MPI_Init(NULL, NULL);

  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);

  if (rank == 0)
    {
      printf("Running on %d processes\n", size);

      // data files should already exist but create them here so we
      // have a self-contained program

      fp = fopen("mesh.dat", "w");
      fwrite(meshfile, sizeof(int), NELEM*NVERT, fp);
      fclose(fp);

      fp = fopen("part.dat", "w");
      fwrite(partfile, sizeof(int), NELEM, fp);
      fclose(fp);
    }

  // could read on rank 0 and broadcast, but using MPI-IO then
  // "readall" should take an efficient collective approach
  // every rank read the whole partition file

  MPI_File_open(comm, "part.dat", MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);
  MPI_File_set_view(fh, 0, MPI_INT, MPI_INT, "native", MPI_INFO_NULL);
  MPI_File_read_all(fh, part, NELEM, MPI_INT, &status);
  MPI_File_close(&fh);

  nelemlocal = 0;

  // pick out local elements and record displacements

  for (i=0; i < NELEM; i++)
    {
      if (part[i] == rank)
        {
          disp[nelemlocal] = i*NVERT;
          nelemlocal += 1;
        }
    }

  printf("on rank %d, nelemlocal = %d\n", rank, nelemlocal);

  // create the MPI datatype to use as the filetype, which is
  // effectively a mask that selects only the elements for this rank

  MPI_Type_create_indexed_block(nelemlocal, NVERT, disp, MPI_INT, &filetype);
  MPI_Type_commit(&filetype);

  MPI_File_open(comm, "mesh.dat", MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);

  // set the file view appropriate to this rank
  MPI_File_set_view(fh, 0, MPI_INT, filetype, "native", MPI_INFO_NULL);

  // each rank only reads its own set of elements based on file view
  MPI_File_read_all(fh, mesh, nelemlocal*NVERT, MPI_INT, &status);

  MPI_File_close(&fh);

  // check we got the correct data

  for (i=0; i < nelemlocal; i++)
    {
      printf("on rank %d, mesh[%d] = %d, %d, %d\n",
         rank, i, mesh[i][0], mesh[i][1], mesh[i][2]);
    }

  MPI_Finalize();
}

and it seems to give the right answer;

dsh@laptop$ mpicc -o metisio metisio.c
dsh@laptop$ mpirun -n 4 ./metisio | sort
on rank 0, mesh[0] = 3, 5, 4
on rank 0, mesh[1] = 4, 5, 6
on rank 0, mesh[2] = 5, 7, 6
on rank 0, mesh[3] = 3, 8, 5
on rank 0, nelemlocal = 4
on rank 1, mesh[0] = 8, 7, 5
on rank 1, mesh[1] = 9, 7, 8
on rank 1, mesh[2] = 9, 8, 3
on rank 1, nelemlocal = 3
on rank 2, mesh[0] = 1, 3, 2
on rank 2, mesh[1] = 2, 3, 4
on rank 2, mesh[2] = 11, 1, 12
on rank 2, mesh[3] = 12, 1, 2
on rank 2, nelemlocal = 4
on rank 3, mesh[0] = 1, 9, 3
on rank 3, mesh[1] = 10, 9, 1
on rank 3, mesh[2] = 11, 10, 1
on rank 3, nelemlocal = 3
Running on 4 processes
David Henty
  • 1,694
  • 9
  • 11
  • In the code I noticed that if I use a dynamic 2D array for mesh i get a segmant fault.MPI IO needs to store file data in contigeous memory blocks? – spyros May 26 '20 at 14:39