Wave Equation Implementation (Version 2)

Wave Equation Implementation (Version 2)

Take a look at the code below and note the MPI functions. Click on the name of each MPI routine to read a short description of that routine's purpose and syntax.

/*-------------------------------------------------------------------
 *   This program implements the concurrent wave equation described  
 *   in Chapter 5 of Fox et al., 1988, Solving Problems on Concurrent
 *   Processors, vol 1.   
 *
 *   A vibrating string is decomposed into points.  Each process is  
 *   responsible for updating the amplitude of a number of points over
 *   time.
 *       
 *   At each iteration, each process exchanges boundary points with
 *   nearest neighbors.  This version uses low level sends and receives
 *   to exchange boundary points.
 *
 *   Revised:  11/27/97  by Saleh Elmohamed
 *------------------------------------------------------------------- */

#include 
#include 
#include 
#include 
#include 

#define MASTER 0 
#define MAXPOINTS 1000     /* the recommended max number of points */
#define MAXSTEPS 10000     /* the recommended max number of steps  */
#define PI 3.14159265

int E_OUT1 = 30;
int E_OUT2 = 40;

void init_master(void);
void init_workers(void);
void init_line(void);
void update (void);
void output_master(void);
void output_workers(void);

MPI_Status status;
MPI_Request request;
int procid,                /* task ID */
    nproc,                 /* number of processes  */
    nsteps,                /* number of time steps */
    tpoints,               /* total points along string */
    npoints,               /* number of points handled by this processor */
    first,                 /* index of 1st point handled by this processor */
    rcode;                 /* generic return code */
double values[MAXPOINTS+2],  /* values at time t */
       oldval[MAXPOINTS+2],  /* values at time (t-dt) */
       newval[MAXPOINTS+2];  /* values at time (t+dt) */

double begin_time, stop_time;

/* time_t  t1, t2; */
/*------------------------------------------------------------
 *      Master obtains input values from user
 *------------------------------------------------------------ */

void init_master(void)
{
  char tchar[8];
  int buffer[2];

  /*==============================================*/
  /* Set number of points, number of iterations   */
  /* These numbers can be set by the user. It is  */
  /* recommended that tpoints not to exceed 1000  */
  /* and nsteps not to exceed 10000.              */
  /*==============================================*/

  tpoints = 100;
  nsteps = 1000;

#if 0
  while ((tpoints < nproc) || (tpoints > MAXPOINTS))
    {
      printf("Enter number of points along vibrating string\n");
      scanf("%s", tchar);
      tpoints = atoi(tchar);
      if ((tpoints < nproc) || (tpoints > MAXPOINTS))
	printf("enter value between %d and %d\n", nproc, MAXPOINTS);
    }
  while ((nsteps < 1) || (nsteps > MAXSTEPS))
    {
      printf("Enter number of time steps\n");
      scanf("%s", tchar);
      nsteps = atoi(tchar);
      if ((nsteps < 1) || (nsteps > MAXSTEPS))
	printf("enter value between 1 and %d\n", MAXSTEPS);
    }
#endif

  printf("%d: points = %d, steps = %d\n", procid, tpoints, nsteps);

  /* then broadcast time advance parameter, total points, time steps */
  buffer[0] = tpoints;
  buffer[1] = nsteps;
  MPI_Bcast(&buffer, 2, 
    MPI_INT, MASTER, MPI_COMM_WORLD);
}

/*------------------------------------------------------------
 *    Workers receive input values from master
 *------------------------------------------------------------*/
  
void init_workers(void)
{
  int buffer[2];

  /* receive total points, time steps */
  MPI_Bcast(&buffer, 2, 
    MPI_INT, MASTER, MPI_COMM_WORLD);
  tpoints = buffer[0];
  nsteps = buffer[1];
}

/*---------------------------------------------------------
 *   All processes initialize points on line
 *--------------------------------------------------------- */

void init_line(void)
{
  int nmin, nleft, npts, i, j, k;
  double x, fac;

  /* calculate initial values based on sine curve */
  nmin = tpoints/nproc;
  nleft = tpoints%nproc;
  fac = 2.0 * PI;
  for (i = 0, k = 0; i < nproc; i++)
    {
      npts = (i < nleft) ? nmin + 1 : nmin;
      if (procid == i)
	{
	  first = k + 1;
	  npoints = npts;
	  printf ("%d: first = %d, npoints = %d\n", procid, first, npts);
	  for (j = 1; j <= npts; j++, k++)
	    {
	      x = (double)k/(double)(tpoints - 1);
	      values[j] = sin (fac * x);
	    }  
	}
      else k += npts;
    }
  for (i = 1; i <= npoints; i++) 
    oldval[i] = values[i];
}

/*-----------------------------------------------------------
 *    Calculate new values using wave equation
 *-----------------------------------------------------------*/
  
void do_math(int i)
{
  const double dtime = 0.3;
  const double c = 1.0;
  const double dx = 1.0; 
  double tau, sqtau;

  tau = (c * dtime / dx);
  sqtau = tau * tau;
  newval[i] = (2.0 * values[i]) - oldval[i]  
    + (sqtau * (values[i-1] - (2.0 * values[i]) + values[i+1]));
}

/*------------------------------------------------------------------
 * All processes update their points a specified number of times  
 *------------------------------------------------------------------*/
  
void update(void)
{
  int i, j, msglen, left, right;

  /* update values for each point along string */
  msglen = sizeof (double);
  for (i = 1; i <= nsteps; i++)
    {
      left = procid -1;
      if(left < 0) 
	left = MPI_PROC_NULL;
      right = procid + 1;
      if(right == nproc ) 
	right = MPI_PROC_NULL;
      /* shift data to the left (non-circular) */
      MPI_Sendrecv(&values[1], 
       1, MPI_DOUBLE, left, 0,&values[npoints+1], 1, MPI_DOUBLE,
       right, 0, MPI_COMM_WORLD, &status);
      
      /* shift data to the right (non-circular) */
      MPI_Sendrecv(&values[npoints], 
       1, MPI_DOUBLE, right, 0, &values[0], 1, MPI_DOUBLE, left, 0, 
       MPI_COMM_WORLD, &status);

      /* update points along line */
      for (j = 1; j <= npoints; j++)
	{
	  /* global endpoints */
	  if ((first + j - 1 == 1) || (first + j - 1 == tpoints))
	    newval[j] = 0.0;
	  else
	    do_math(j);
	}
      for (j = 1; j <= npoints; j++)
	{
	  oldval[j] = values[j];
	  values[j] = newval[j];
	}
    }
}

/*-------------------------------------------------------------
 * Master receives results from workers and prints
 *------------------------------------------------------------- */

void output_master(void)  
{
  int i, start, npts, tpts, buffer[2];  
  double results[MAXPOINTS];
  
  /* store worker's results in results array */
  for (i = 1; i < nproc; i++)
    {
     /* receive number of points and first point */
     MPI_Recv(&buffer, 2, 
       MPI_INT, MPI_ANY_SOURCE, E_OUT1, MPI_COMM_WORLD,
	       &status);

      start = buffer[0];
      npts = buffer[1];

      /* receive results */
      MPI_Recv(&results[start-1], 
        npts, MPI_DOUBLE, MPI_ANY_SOURCE,
               E_OUT2, MPI_COMM_WORLD, &status);
    }

  /* store master's results in results array */
  for (i = first; i < first + npoints; i++)
    results[i-1] = values[i];

  tpts = (tpoints < 10) ? tpoints: 10;  
  printf("The first %d points (just for validation):\n", tpts);
  for (i = 0; i < tpts; i++) 
    printf("%4.2f  ", results[i]);
  printf("\n");
}

/*--------------------------------------------------------------
 *      Workers send the updated values to the master
 *--------------------------------------------------------------*/
  
void output_workers(void)
{
  int buffer[2];

  /* send first point and number of points handled to master */
  buffer[0] = first;
  buffer[1] = npoints;
  MPI_Isend(&buffer, 2, 
    MPI_INT, MASTER, E_OUT1, MPI_COMM_WORLD, &request);
  MPI_Wait(&request, &status);
  /* send results to master */
  MPI_Isend(&values[1], 
    npoints, MPI_DOUBLE, MASTER, E_OUT2, MPI_COMM_WORLD,
	    &request);
  MPI_Wait(&request, &status);
}

/*------------------------------------------------------
 *      Main program
 *-----------------------------------------------------*/

main(int argc, char **argv)
{
  /* learn number of tasks and task ID */
  rcode = MPI_Init(&argc, &argv);
  begin_time = MPI_Wtime();
  MPI_Comm_rank(MPI_COMM_WORLD,
    &procid); 
  MPI_Comm_size(MPI_COMM_WORLD,&nproc);
  if (rcode >= 0) 
    printf ("Processor %d: Wave Program running\n", procid);

  /* get program parameters and initialize wave values */
  if (procid == MASTER) 
    init_master();
  else 
    init_workers();

  init_line();

  /*  t1 = time(NULL); */
  /* update values along the line for nstep time steps */
  update();

  /* collect results and print */
  if (procid == MASTER) 
    output_master();
  else 
    output_workers();

  stop_time = (MPI_Wtime() - begin_time);
  printf("... it took about %6.3f seconds\n", stop_time);

  /*  t2 = time(NULL); */

  /*
  printf("The whole thing took about %4.2f seconds\n", 
         difftime(t2, t1));
  */
 MPI_Finalize();
}


A 4-processor run: Processor 1: Wave Program running Processor 0: Wave Program running Processor 2: Wave Program running Processor 3: Wave Program running 0: points = 100, steps = 1000 0: first = 1, npoints = 25 1: first = 26, npoints = 25 2: first = 51, npoints = 25 3: first = 76, npoints = 25 ... it took about 1.085 seconds ... it took about 1.085 seconds ... it took about 1.085 seconds The first 10 points (just for validation): 0.00 0.06 0.12 0.19 0.25 0.31 0.36 0.42 0.48 0.53 ... it took about 1.089 seconds