TEXT   27

# 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);
32.     Nx = atoi(argv);
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 = 1;
52.     phi = 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. }