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);
35.     ii = atoi(argv);
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 = 1;
55.     phi = 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. }