TEXT   49

Gauss Siedel solution to Poisson equation

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

Raw Paste

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