Commit 2487fe8b authored by David Goz's avatar David Goz 😴
Browse files

mpi/omp_comm round-off error spotting

parent ec1c275b
Loading
Loading
Loading
Loading
+1 −1
Original line number Diff line number Diff line
# set the MPI install path

MPI_INSTALL_PATH = /home/gozzilla/software/openmpi/openmpi-5.0.3
MPI_INSTALL_PATH = /home/darkenergy/software/openmpi/openmpi-5.0.3



+79 −61
Original line number Diff line number Diff line
@@ -218,21 +218,7 @@ int main(int argc, char **argv)
     -------------------------------------------------------- */
    
  const int ibeg = NGHOST;
  /* const int iend   = ibeg + NX_GLOB - 1; */
  /* const int nx     = iend - ibeg + 1; */
  /* const int nx_tot = nx + 2 * NGHOST; */
    
  const int jbeg = NGHOST;
  /* const int jend   = jbeg + NY_GLOB - 1; */
  /* const int ny     = jend - jbeg + 1; */
  /* const int ny_tot = ny + 2 * NGHOST; */

  /* if (rank == MASTER) */
  /*   { */
  /*     printf("\n\t Grid indices:"); */
  /*     printf("\n\t\t ibeg, iend = %d, %d; nx_tot = %d"    ,ibeg, iend, nx_tot); */
  /*     printf("\n\t\t jbeg, jend = %d, %d; ny_tot = %d\n\n",jbeg, jend, ny_tot); */
  /*   } */
  
  /* --------------------------------------------------------
     2. Generate grids, allocate memory
@@ -426,7 +412,28 @@ void JacobiAlgorithm(MyData **const restrict Phi,
	  Phi[j][i] = 0.25 * (Phi0[j][i-1] + Phi0[j][i+1] +
			      Phi0[j-1][i] + Phi0[j+1][i]);

	  *error += delta[X] * delta[Y] * fabs(Phi[j][i] - Phi0[j][i]);
	  *error += (delta[X] * delta[Y] * fabs(Phi[j][i] - Phi0[j][i]));
	} /* loop over columns */
    } /* loop over rows */

  return;
}

void get_error(MyData **const restrict Phi,
	       MyData **const restrict Phi0,
	       const MyData  *restrict delta,
	       const int               jbeg,
	       const int               jend,
	       const int               ibeg,
	       const int               iend,
	       MyData *const restrict  error)
{
  *error = 0.0;
  for (int j=jbeg ; j<=jend ; j++)
    {
      for (int i=ibeg ; i<=iend ; i++)
	{
	  *error += (delta[X] * delta[Y] * fabs(Phi[j][i] - Phi0[j][i]));
	} /* loop over columns */
    } /* loop over rows */
    
@@ -448,21 +455,25 @@ void Jacobi_Communication(MyData **const restrict Phi,
  const int data_row_size = ThisTask->domain.dim[Y];
  
  /* First task: issue the communication */
  MPI_Request request[8];
  MPI_Request request[4];

  MyData **const restrict buffer = Phi0;
  
  /* 4 nonblocking MPI_Irecv */
  MPI_Irecv(&buffer[ThisTask->domain.local_start[X] - 1][ThisTask->domain.local_start[Y]    ], data_row_size, MPI_MyDatatype, ThisTask->nbrbottom, 0, ThisTask->comm2d, &request[0]);
  MPI_Irecv(&buffer[ThisTask->domain.local_end[X]   + 1][ThisTask->domain.local_start[Y]    ], data_row_size, MPI_MyDatatype, ThisTask->nbrtop,    1, ThisTask->comm2d, &request[1]);
  MPI_Irecv(&buffer[ThisTask->domain.local_start[X]    ][ThisTask->domain.local_start[Y] - 1], 1,             column,         ThisTask->nbrleft,   2, ThisTask->comm2d, &request[2]);
  MPI_Irecv(&buffer[ThisTask->domain.local_start[X]    ][ThisTask->domain.local_end[Y]   + 1], 1,             column,         ThisTask->nbrright,  3, ThisTask->comm2d, &request[3]);
  MPI_Isendrecv(&buffer[ThisTask->domain.local_end[X]      ][ThisTask->domain.local_start[Y]    ], data_row_size, MPI_MyDatatype, ThisTask->nbrtop,    0,
		&buffer[ThisTask->domain.local_start[X] - 1][ThisTask->domain.local_start[Y]    ], data_row_size, MPI_MyDatatype, ThisTask->nbrbottom, 0,
		ThisTask->comm2d, &request[0]);

  MPI_Isendrecv(&buffer[ThisTask->domain.local_start[X]    ][ThisTask->domain.local_start[Y]    ], data_row_size, MPI_MyDatatype, ThisTask->nbrbottom, 1,
		&buffer[ThisTask->domain.local_end[X]   + 1][ThisTask->domain.local_start[Y]    ], data_row_size, MPI_MyDatatype, ThisTask->nbrtop,    1,
		ThisTask->comm2d, &request[1]);

  /* 4 nonblocking MPI_Isend */
  MPI_Isend(&buffer[ThisTask->domain.local_end[X]  ][ThisTask->domain.local_start[Y]], data_row_size, MPI_MyDatatype, ThisTask->nbrtop,    0, ThisTask->comm2d, &request[4]);
  MPI_Isend(&buffer[ThisTask->domain.local_start[X]][ThisTask->domain.local_start[Y]], data_row_size, MPI_MyDatatype, ThisTask->nbrbottom, 1, ThisTask->comm2d, &request[5]);
  MPI_Isend(&buffer[ThisTask->domain.local_start[X]][ThisTask->domain.local_end[Y]  ], 1,             column,         ThisTask->nbrright,  2, ThisTask->comm2d, &request[6]);
  MPI_Isend(&buffer[ThisTask->domain.local_start[X]][ThisTask->domain.local_start[Y]], 1,             column,         ThisTask->nbrleft,   3, ThisTask->comm2d, &request[7]);
  MPI_Isendrecv(&buffer[ThisTask->domain.local_start[X]    ][ThisTask->domain.local_end[Y]      ], 1,             column,         ThisTask->nbrright,  2,
		&buffer[ThisTask->domain.local_start[X]    ][ThisTask->domain.local_start[Y] - 1], 1,             column,         ThisTask->nbrleft,   2,
		ThisTask->comm2d, &request[2]);
  
  MPI_Isendrecv(&buffer[ThisTask->domain.local_start[X]    ][ThisTask->domain.local_start[Y]    ], 1,             column,         ThisTask->nbrleft,   3,
		&buffer[ThisTask->domain.local_start[X]    ][ThisTask->domain.local_end[Y]   + 1], 1,             column,         ThisTask->nbrright,  3,
		ThisTask->comm2d, &request[3]);

  /**************************************** computation ****************************************/
  /* perform the computation with the local data, (i.e. ghost cells are not required) */
@@ -473,39 +484,46 @@ void Jacobi_Communication(MyData **const restrict Phi,
  const int ibeg = ThisTask->domain.local_start[Y] + 1;
  const int iend = ThisTask->domain.local_end[Y]   - 1;
  
  /* wait the data on the boundaries */
  MPI_Waitall(8, request, MPI_STATUSES_IGNORE);
  *error = 0.0;
  JacobiAlgorithm(Phi, Phi0, delta, jbeg, jend, ibeg, iend, error);

  /*********************************************************************************************/
  
  /* wait the data on the boundaries */
  MPI_Waitall(4, request, MPI_STATUSES_IGNORE);
  
  /*  nbrbottom */
  JacobiAlgorithm(Phi, Phi0, delta,
		  ThisTask->domain.local_start[X], ThisTask->domain.local_start[X],
		  ThisTask->domain.local_start[Y], ThisTask->domain.local_end[Y],
		  error);

  /* nbrtop */
  JacobiAlgorithm(Phi, Phi0, delta,
		  ThisTask->domain.local_end[X],   ThisTask->domain.local_end[X],
		  ThisTask->domain.local_start[Y], ThisTask->domain.local_end[Y],
		  error);

  /* nbrleft */
  JacobiAlgorithm(Phi, Phi0, delta,
		  ThisTask->domain.local_start[X], ThisTask->domain.local_end[X],
		  ThisTask->domain.local_start[Y], ThisTask->domain.local_start[Y],
		  error);

  /* nbrright */
  JacobiAlgorithm(Phi, Phi0, delta,
		  ThisTask->domain.local_start[X], ThisTask->domain.local_end[X],
		  ThisTask->domain.local_end[Y],   ThisTask->domain.local_end[Y],
		  error);

  /* Round-off error fixing ??? */
  {
    *error = 0.0;
  /* JacobiAlgorithm(Phi, Phi0, delta, jbeg, jend, ibeg, iend, error); */
  JacobiAlgorithm(Phi, Phi0, delta, jbeg-1, jend+1, ibeg-1, iend+1, error);

  /* /\* wait the data on the boundaries *\/ */
  /* MPI_Waitall(8, request, MPI_STATUSES_IGNORE); */
  
  /* /\*  nbrbottom *\/ */
  /* JacobiAlgorithm(Phi, Phi0, delta, */
  /* 		  ThisTask->domain.local_start[X], ThisTask->domain.local_start[X], */
  /* 		  ThisTask->domain.local_start[Y], ThisTask->domain.local_end[Y], */
  /* 		  error); */

  /* /\* nbrtop *\/ */
  /* JacobiAlgorithm(Phi, Phi0, delta, */
  /* 		  ThisTask->domain.local_end[X],   ThisTask->domain.local_end[X], */
  /* 		  ThisTask->domain.local_start[Y], ThisTask->domain.local_end[Y], */
  /* 		  error); */

  /* /\* nbrleft *\/ */
  /* JacobiAlgorithm(Phi, Phi0, delta, */
  /* 		  ThisTask->domain.local_start[X], ThisTask->domain.local_end[X], */
  /* 		  ThisTask->domain.local_start[Y], ThisTask->domain.local_start[Y], */
  /* 		  error); */

  /* /\* nbrright *\/ */
  /* JacobiAlgorithm(Phi, Phi0, delta, */
  /* 		  ThisTask->domain.local_start[X], ThisTask->domain.local_end[X], */
  /* 		  ThisTask->domain.local_end[Y],   ThisTask->domain.local_end[Y], */
  /* 		  error); */
    get_error(Phi, Phi0, delta,
	      ThisTask->domain.local_start[X], ThisTask->domain.local_end[X],
	      ThisTask->domain.local_start[Y], ThisTask->domain.local_end[Y],
	      error);
  }
  
  MPI_Type_free(&column);

+1 −1

File changed.

Contains only whitespace changes.