C   45

Red-Black implementation

Guest on 27th April 2022 01:43:10 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. *    Red-Black implementation                       *
  9. *                                                   *
  10. *****************************************************/
  11.                                 /* Michel Vallieres */
  12.  
  13. #include <stdio.h>
  14. #include <string.h>
  15. #include <stdlib.h>
  16. #include <math.h>
  17.  
  18. #define EPS 1.0e-7
  19.  
  20. int main ( int argc, char *argv[] )
  21. {
  22.     int    iter, ii, i, N_iterations;
  23.     double *phi, *source;
  24.     int    *ibound;
  25.     double difference, diff, old, h, chi;
  26.     int    scan;
  27.  
  28.                                /* lattice definition & iteration number */
  29.     if ( argc != 3 )
  30.       {
  31.         printf( " Syntax: poisson  N  N_x \n" );
  32.         exit(1);
  33.       }
  34.     N_iterations = atoi(argv[1]);
  35.     ii = atoi(argv[2]);
  36.     h = 1.0/(double)(ii-1);
  37.  
  38.                                /* make space for arrays */
  39.     phi = (double *)malloc( ii*sizeof(double) );
  40.     source = (double *)malloc( ii*sizeof(double) );
  41.     ibound = (int *)malloc( ii*sizeof(int) );
  42.  
  43.                                /* initial fields */
  44.                                /* no source      */
  45.     for ( i=0; i<ii ; i++ )
  46.       {
  47.         ibound[i] = 0;
  48.         phi[i] = 0.0;
  49.         source[i] = 1000 * exp( - 0.1*(i-ii/2)*(i-ii/2) );
  50.       }
  51.                                /* boundary conditions */
  52.     ibound[ii-1] = 1;
  53.     phi[ii-1] = -1.0;
  54.     ibound[0] = 1;
  55.     phi[0] = 0.0;
  56.                                /* isolated points in grid */
  57.     ibound[ii/3] = 1;
  58.     phi[ii/3] = +1.0;
  59.  
  60.                                /* iterate the solution */
  61.     for ( iter=0; iter<N_iterations ; iter++)
  62.       {
  63.                                /* over-relaxation */
  64.         chi = 1.0;
  65.         if ( iter > 30 ) chi = 1.2;
  66.         if ( iter > 50 ) chi = 1.8;
  67.         diff = 0.0;
  68.                                /* Red-Black scanning mode */
  69.         for ( scan=0; scan<2 ; scan++ )
  70.           {
  71.                                /* scan lattice */
  72.             for ( i=scan; i<ii ; i=i+2 )
  73.               {
  74.                 if ( ibound[i] == 0 )
  75.                   { old = phi[i];
  76.                     phi[i] = 0.5 * ( phi[i-1] + phi[i+1] - h*h*source[i] );
  77.                     phi[i] = (1.0-chi) * old + chi*phi[i];
  78.                     if ( fabs( old - phi[i] ) > diff )
  79.                               diff = fabs( old - phi[i] );
  80.                   }
  81.               }
  82.           }
  83.                              /* accurate enough ? */
  84.         if ( diff < EPS ) break;
  85.       }
  86.  
  87.     fprintf( stderr,"\n Converged in %d iterations \n", iter );
  88.  
  89.                              /* check solution */
  90.     diff = 0.0;
  91.     for ( i=0; i<ii ; i++ )
  92.       {
  93.         if ( ibound[i] == 0 )
  94.           {
  95.             difference = 0.5 * ( phi[i-1] - 2*phi[i] + phi[i+1]
  96.                                    - h*h*source[i] );
  97.             if ( fabs( difference ) > diff )
  98.                            diff = fabs( difference );
  99.           }
  100.       }
  101.     fprintf( stderr," Largest error: %e \n\n", diff );
  102.  
  103.                              /* output results */
  104.     for ( i=0; i<ii ; i++ )
  105.        printf( " %f %f \n", i*h, phi[i] );
  106.  
  107.  
  108.                              /* free allocated memory */
  109.     free( (char *) phi );
  110.     free( (char *) source);
  111.     free( (char *) ibound);
  112. }

Raw Paste

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