C   18

poisson parallel 1d c

Guest on 1st August 2022 12:23:32 AM

  1. /******************************************************************
  2. *    1-D Gauss Siedel solution to Poisson equation                *
  3. *                                                                 *
  4. *    Syntax                                                       *
  5. *         Poisson  N  N_x                                         *
  6. *            ( Max # iteration, grid in x )                       *
  7. *                                                                 *
  8. *    Compilation                                                  *
  9. *    mpicc poisson_parallel_1d.c NR.c -lm -o poisson_parallel_1d  *
  10. *                                                                 *
  11. *    Run                                                          *
  12. *    mpiexec -np 4 ./poisson_parallel_1d 1000 100                 *
  13. *                                                                 *
  14. *                                             Parallel Version    *
  15. *******************************************************************/
  16.                                 /* Michel Vallieres */
  17.  
  18. #include <stdio.h>
  19. #include <string.h>
  20. #include <stdlib.h>
  21. #include <math.h>
  22. #include <mpi.h>
  23.  
  24.  
  25. #define EPS 1.0e-7             /* tolerable error     */
  26. #define MAX_PROCS 100          /* max # of processors */
  27. #define MAX_GRID_POINTS 1000   /* max grid points     */
  28.  
  29.  
  30. double phi[MAX_GRID_POINTS], source[MAX_GRID_POINTS];
  31. int    ibound[MAX_GRID_POINTS];
  32.  
  33.  
  34. /* Gauss-Siedel in local piece of lattice */
  35. void Gauss_Siedel_step( int imin, int imax, int iter,
  36.                                     double *diff, double h )
  37. {
  38.     double old, chi;
  39.     int i;
  40.                                          /* over-relaxation */
  41.     chi = 1.0;
  42.     if ( iter > 30 ) chi = 1.2;
  43.     if ( iter > 50 ) chi = 1.8;
  44.  
  45.     *diff = 0.0;
  46.                                          /* scan local lattice */
  47.     for ( i=imin; i<=imax ; i++ )
  48.        {
  49.          if ( ibound[i] == 0 )
  50.            { old = phi[i];
  51.              phi[i] = 0.5 * ( phi[i-1] + phi[i+1] - h*h*source[i] );
  52.              phi[i] = (1.0-chi) * old + chi*phi[i];
  53.              if ( fabs( old - phi[i] ) > *diff )
  54.                             *diff = fabs( old - phi[i] );
  55.            }
  56.        }
  57. }
  58.  
  59.  
  60. /* to exchange the ghost rows */
  61. void exchange ( int myid, int numprocs, int imin, int imax )
  62. {
  63.     MPI_Status  recv_status;
  64.     int         top_process, bot_process;
  65.  
  66.     top_process = myid + 1;
  67.     bot_process = myid - 1;
  68.                                         /* odd - even scheme for safety */
  69.                                         /* & efficiency                 */
  70.     if ( myid % 2 != 0 )
  71.                                         /* odd process */
  72.       {
  73.         if ( myid < numprocs-1 )
  74.            MPI_Ssend( &phi[imax], 1, MPI_DOUBLE, top_process, 121,
  75.                                   MPI_COMM_WORLD);
  76.         if ( myid > 0 )
  77.            MPI_Recv( &phi[imin-1], 1, MPI_DOUBLE, bot_process, 122,
  78.                                   MPI_COMM_WORLD, &recv_status);
  79.         if ( myid > 0 )
  80.            MPI_Ssend( &phi[imin], 1, MPI_DOUBLE, bot_process, 123,
  81.                                   MPI_COMM_WORLD);
  82.         if ( myid < numprocs-1 )
  83.            MPI_Recv( &phi[imax+1], 1, MPI_DOUBLE, top_process, 124,
  84.                                   MPI_COMM_WORLD, &recv_status);
  85.       }
  86.     else
  87.                                         /* even process */
  88.       {
  89.         if ( myid > 0 )
  90.            MPI_Recv( &phi[imin-1], 1, MPI_DOUBLE, bot_process, 121,
  91.                                   MPI_COMM_WORLD, &recv_status);
  92.         if ( myid < numprocs-1 )
  93.            MPI_Ssend( &phi[imax], 1, MPI_DOUBLE, top_process, 122,
  94.                                   MPI_COMM_WORLD);
  95.         if ( myid < numprocs-1 )
  96.            MPI_Recv( &phi[imax+1], 1, MPI_DOUBLE, top_process, 123,
  97.                                   MPI_COMM_WORLD, &recv_status);
  98.         if ( myid > 0 )
  99.            MPI_Ssend( &phi[imin], 1, MPI_DOUBLE, bot_process, 124,
  100.                                   MPI_COMM_WORLD);
  101.       }
  102. }
  103.  
  104.  
  105. /* set up local domain - slices */
  106. void decompose_domain ( int myid, int numprocs, int Nx, int *begin_slice,
  107.                        int *end_slice, int *size_slice, int *imin,
  108.                         int *imax )
  109. {
  110.   int ip, j, average, leftover;
  111.   int temp, slice, N_slices;
  112.                                         /* check dimensions */
  113.     if ( numprocs > MAX_PROCS )
  114.       {
  115.         if ( myid == 0 )  fprintf( stderr,
  116.                "Code written for fewer processes than %d \n", MAX_PROCS);
  117.         MPI_Finalize();
  118.         exit(1);
  119.       }
  120.                                         /* partition domain */
  121.     N_slices =  numprocs;
  122.     average = Nx / N_slices;
  123.     leftover = Nx % N_slices;
  124.     temp = -1;
  125.     for ( slice=0 ; slice<N_slices ; slice++ )
  126.     {
  127.       begin_slice[slice] = temp + 1;
  128.       end_slice[slice] = begin_slice[slice] + average - 1;
  129.       if ( leftover >0 )
  130.         {
  131.           end_slice[slice]++;
  132.           leftover--;
  133.         }
  134.       temp = end_slice[slice];
  135.       size_slice[slice] = end_slice[slice] - begin_slice[slice] + 1;
  136.     }
  137.                                         /* echo partitions */
  138.     if( myid == 0 )
  139.       {
  140.         fprintf( stderr, "\n Local grid: \n ");
  141.         for ( ip=0 ; ip<numprocs ; ip++ )
  142.            fprintf( stderr, " ( %d , %d )", begin_slice[ip],
  143.                                     end_slice[ip] );
  144.         fprintf( stderr, "\n\n" );
  145.       }
  146.                                         /* local sub-domain */
  147.                                         /* careful with ghost line */
  148.     *imin = begin_slice[myid];
  149.     *imax = end_slice[myid];
  150. }
  151.  
  152.  
  153. /* combine sub domains in one -- checks & output */
  154. void gather_check_output ( int myid, int numprocs, int *begin_slice,
  155.                               int *end_slice, int *size_slice,
  156.                                   int Nx, int iter, double h )
  157. {
  158.     int         i, ip;
  159.     double      diff, difference, totaldiff;
  160.     double      *final_phi;
  161.     MPI_Status  recv_status;
  162.  
  163.     if( myid == 0 )
  164.        fprintf( stderr,"\n Converged in %d iterations \n", iter );
  165.  
  166.                                         /* exchange the Ghost values */
  167.     exchange ( myid, numprocs, begin_slice[myid], end_slice[myid] );
  168.  
  169.                                         /* check solution */
  170.     diff = 0.0;
  171.     for ( i=begin_slice[myid] ; i<=end_slice[myid] ; i++ )
  172.       {
  173.         if ( ibound[i] == 0 )
  174.           {
  175.             difference = 0.5 * ( phi[i-1] - 2*phi[i] + phi[i+1]
  176.                                    - h*h*source[i] );
  177.             if ( fabs( difference ) > diff )
  178.                            diff = fabs( difference );
  179.           }
  180.       }
  181.                                         /* find max amongst processes */
  182.     MPI_Allreduce ( &diff, &totaldiff, 1, MPI_DOUBLE,
  183.                                             MPI_MAX, MPI_COMM_WORLD );
  184.     if ( myid == 0 )
  185.         fprintf( stderr," Largest error: %e \n\n", totaldiff );
  186.  
  187.  
  188.                                         /* gather pieces of phi */
  189.     if ( myid == 0 )
  190.       {
  191.                                         /* space to gather final phi */
  192.          final_phi = (double *)malloc( Nx*sizeof(double) );
  193.  
  194.                                         /* transfer own sub-domain */
  195.                                         /* piece of phi            */
  196.          for ( i=begin_slice[myid]; i<=end_slice[myid]; i++ )
  197.            {
  198.              final_phi[i] = phi[i];
  199.            }
  200.                                         /* receives pieces of    */
  201.                                         /* phi from sub-domains  */
  202.          for ( ip=1 ; ip<numprocs ; ip++ )
  203.            {
  204.               MPI_Recv( &final_phi[begin_slice[ip]],
  205.                       size_slice[ip], MPI_DOUBLE, ip, 229,
  206.                                   MPI_COMM_WORLD, &recv_status);
  207.            }
  208.                                         /* output results */
  209.          for ( i=0; i<Nx ; i++ )
  210.                    printf( " %f %f \n", i*h, final_phi[i] );
  211.  
  212.                                         /* free up memory space */
  213.          free( (char *) final_phi );
  214.       }
  215.                                         /* nodes send phi to node 0 */
  216.     else
  217.       {
  218.          MPI_Ssend( &phi[begin_slice[myid]], size_slice[myid],
  219.                        MPI_DOUBLE, 0, 229, MPI_COMM_WORLD );
  220.       }
  221.  
  222. }
  223.  
  224.  
  225.  
  226. int main ( int argc, char *argv[] )
  227. {
  228.     int    iter, Nx, i, is, N_iterations;
  229.     double difference, diff, totaldiff, h, chi;
  230.     int    myid, numprocs;
  231.     int    begin_slice[MAX_PROCS], end_slice[MAX_PROCS],
  232.            size_slice[MAX_PROCS], imin, imax;
  233.  
  234.                                         /* join the MPI virtual machine */
  235.     MPI_Init(&argc, &argv);
  236.     MPI_Comm_rank(MPI_COMM_WORLD, &myid);
  237.     MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
  238.  
  239.                                         /* lattice definition & iteration number */
  240.     if ( argc != 3 )
  241.       {
  242.         if ( myid == 0 )
  243.             printf( " Syntax: poisson  N  N_x \n" );
  244.         MPI_Finalize();
  245.         exit(1);
  246.       }
  247.     N_iterations = atoi(argv[1]);
  248.     Nx = atoi(argv[2]);
  249.     h = 1.0/(double)(Nx-1);
  250.  
  251.                                         /* domain decomposition */
  252.     decompose_domain ( myid, numprocs, Nx, begin_slice,
  253.                          end_slice, size_slice, &imin, &imax );
  254.  
  255.                                         /* initial fields */
  256.     for ( i=imin; i<=imax ; i++ )
  257.       {
  258.         ibound[i] = 0;
  259.         phi[i] = 0.0;
  260.         source[i] = 1000 * exp( - 0.1*(i-Nx/2)*(i-Nx/2) );
  261.       }
  262.                                         /* boundary conditions */
  263.     if ( myid == numprocs - 1 )
  264.       {
  265.          ibound[imax] = 1;
  266.          phi[imax] = -1.0;
  267.       }
  268.     if ( myid == 0 )
  269.       {
  270.         ibound[imin] = 1;
  271.         phi[imin] = 0.0;
  272.       }
  273.                                         /* isolated points in grid */
  274.     is = Nx/3;
  275.     if ( is >= imin || is <= imax )
  276.       {
  277.          ibound[is] = 1;
  278.          phi[is] = +1.0;
  279.       }
  280.  
  281.                                         /* iterate the solution */
  282.     for ( iter=0; iter<N_iterations ; iter++)
  283.       {
  284.  
  285.                                         /* exchange the Ghost values */
  286.         exchange ( myid, numprocs, imin, imax );
  287.  
  288.                                        /* one Gauss-Siedel iteration */
  289.         Gauss_Siedel_step( imin, imax, iter, &diff, h );
  290.  
  291.                                        /* maximum of error from */
  292.                                        /* each sub-domain       */
  293.         MPI_Allreduce ( &diff, &totaldiff, 1, MPI_DOUBLE,
  294.                                             MPI_MAX, MPI_COMM_WORLD );
  295.  
  296.                                        /* accurate enough ? */
  297.         if ( totaldiff < EPS ) break;
  298.       }
  299.  
  300.                                        /* gather, check, & output results */
  301.     gather_check_output ( myid, numprocs, begin_slice,
  302.                        end_slice, size_slice,  Nx, iter, h );
  303.  
  304.                                        /* we are done */
  305.     MPI_Finalize();
  306.     exit(1);
  307. }

Raw Paste


Login or Register to edit or fork this paste. It's free.