TEXT   17

# poisson scalar modular 2d c

Guest on 1st August 2022 12:27:37 AM

1. /*
2.     Over-relaxation solution to Poisson equation
3.
4.     Syntax
5.          Poisson  N  N_x  N_y
6.             ( Max # iteration, grid in x, grid in y )
7.
8.     Note: only written for equal x and y meshes -- N_x = N_y
9.
10.     IBOUND:  positions on grid corresponding to boundaries
11.
12.     Modular code - simple main code - logical subroutines
13.
14.
15.     Build & use of code
16.
17. gcc poisson_scalar_modular_2d.c -o  poisson_scalar_modular_2d
18.
19. poisson_scalar_modular_2  N  N_x  N_y | make_image -s N_x N_y | convert - - | xv -
20.
21. */
22.
23. #include <stdio.h>
24. #include <string.h>
25. #include <stdlib.h>
26. #include <math.h>
27.
28. #define MAX_x 350                       /* maximum dimensions   */
29. #define MAX_y 350
30.
31. #define EPS 1.0e-7                      /* tolerance in solution */
32.
33. #define IMAGE 1                         /* image - yes (1) - no (0) */
34.
35.                                         /* global array & file for solution */
36. int   ibound[MAX_x][MAX_y];
37. double phi[MAX_x][MAX_y], ff[MAX_x][MAX_y];
38. FILE   *outfile;
39.
40.                                         /* functions prototypes */
41. void set_grid ( int argc, char *argv[], int *N_iterations,
42.                 int *ii, int *jj, double *h );
43. void initial_and_boundary ( int ii, int jj );
44. void output_results ( int ii, int jj, double h );
45. void Gauss_Siedel_step( int ii, int jj, int iter, double *diff, double h );
46.
47.
48. int main ( int argc, char *argv[] )
49. {
50.     double h;
51.     int    ii, jj, i, j;
52.     double diff, aux, chi;
53.     int    N_iterations, iter;
54.
55.                                         /* lattice definition & iteration number */
56.     set_grid ( argc, argv, &N_iterations, &ii, &jj, &h );
57.                                         /* initial fields & boundary conditions */
58.     initial_and_boundary ( ii, jj );
59.                                         /* iterate the solution */
60.     for ( iter=1 ; iter<N_iterations ; iter++ )
61.       {
62.                                         /* Gauss-Siedel iteration */
63.         Gauss_Siedel_step( ii, jj, iter, &diff, h );
64.                                         /* solution accurate enough ? */
65.         if ( diff < EPS ) break;
66.       }
67.                                         /* output the phi field */
68.     output_results ( ii, jj, h );
69.
70. }
71.
72.
73. void set_grid ( int argc, char *argv[], int *N_iterations,
74.                 int *ii, int *jj, double *h )
75. {
76.                                         /* lattice definition & iteration number */
77.     if ( argc != 4 )
78.       {
79.         printf( " Syntax: poisson  N  N_x  N_y   ( where N_x = N_y ) \n" );
80.         exit(1);
81.       }
82.     *N_iterations = atoi(argv[1]);
83.     *ii = atoi(argv[2]);
84.     *jj = atoi(argv[3]);
85.     if ( *ii != *jj )
86.       {
87.         fprintf( stderr, " Code only setup for N_x = N_y \n");
88.         exit(1);
89.       }
90.     if ( *ii > MAX_x || *jj > MAX_y)
91.       {
92.         fprintf( stderr, " Grid dimensions too large \n");
93.         exit(1);
94.       }
95.     *h = 1.0/(double)(*ii-1);
96. }
97.
98.
99. void initial_and_boundary ( int ii, int jj )
100. {
101.                                         /* initial phi values & boundary conditions */
102.     int i, j, i_s;
103.                                         /* initial fields */
104.     for ( i=0 ; i<ii ; i++ )
105.       {
106.         for ( j=0 ; j<jj ; j++)
107.           {
108.             ff[i][j] = 0.0;
109.             ibound[i][j] = 0;
110.             phi[i][j] = 0.0;
111.           }
112.       }
113.                                         /* Boundary Conditions */
114.     for ( j=0 ; j<jj ; j++)
115.       {
116.         ibound[0][j] = 1;
117.         phi[0][j] = 2.0;
118.         ibound[ii-1][j] = 1;
119.         phi[ii-1][j] = 2.0;
120.       }
121.     for ( i=0 ; i<ii ; i++ )
122.       {
123.         ibound[i][0] = 1;
124.         phi[i][0] = 2.0;
125.         ibound[i][jj-1] = 1;
126.         phi[i][jj-1] = 2.0;
127.       }
128.                                         /* single value in center of square */
129.     i_s = 4*ii/5;
130.     ibound[i_s][3*jj/4] = 1;
131.     phi[i_s][3*jj/4] = -25.0;
132.     i_s = ii/3;
133.     ibound[i_s][jj/3] = 1;
134.     phi[i_s][jj/3] = -10.0;
135.                                         /* open output file */
136.     if ((outfile=fopen("Field_1.dat","w+"))==NULL)
137.         { fprintf(stderr,"error in file open \n");
138.           exit(1);
139.         }
140.
141.     fprintf( outfile, " Grid size: %d x %d \n", ii, jj );
142.
143.                                         /* output initial phi field */
144.     fprintf( outfile, " Initial - PHI field \n");
145.     for ( i=0 ; i<ii ; i++ )
146.       {
147.       for ( j=0 ; j<jj ; j++)
148.         {
149.             fprintf( outfile, " %f", phi[i][j] );
150.         }
151.         fprintf( outfile, "\n");
152.       }
153.
154.     fprintf( outfile, " Boundary  \n");
155.     for ( i=0 ; i<ii ; i++ )
156.       {
157.       for ( j=0 ; j<jj ; j++)
158.         {
159.             fprintf( outfile, " %d", ibound[i][j] );
160.         }
161.         fprintf( outfile, "\n");
162.       }
163. }
164.
165.
166. void output_results( int ii, int jj, double h )
167.                                         /* output results */
168. {
169.     int    i, j;
170.     double aux, phi_min, phi_max, residual, largest_residual;
171.
172.                                         /* output results */
173.     fprintf( outfile, " Final - PHI field \n");
174.     for ( i=0 ; i<ii ; i++ )
175.       {
176.       for ( j=0 ; j<jj ; j++)
177.         {
178.             fprintf( outfile, " %f", phi[i][j] );
179.         }
180.         fprintf( outfile, "\n");
181.       }
182.                                         /* check if solution is */
183.                                         /* really a solution */
184.     fprintf( outfile, " Check - ( LHS - RHS ) \n");
185.     largest_residual = 0.0;
186.     for ( i=0 ; i<ii ; i++ )
187.       {
188.       for ( j=0 ; j<jj ; j++)
189.         {
190.                                         /* min & max of phi as an aside */
191.           if ( phi[i][j] < phi_min ) phi_min = phi[i][j];
192.           if ( phi[i][j] > phi_max ) phi_max = phi[i][j];
193.
194.           if ( ibound[i][j] == 0 )
195.             {
196.             residual =  phi[i+1][j] + phi[i-1][j] +
197.                              phi[i][j+1] + phi[i][j-1] -
198.                                4.0 * phi[i][j] - h * h * ff[i][j];
199.             }
200.           else
201.             {
202.             residual = 0.0;
203.             }
204.           if( fabs( residual ) > largest_residual )
205.                           largest_residual = fabs( residual );
206.         }
207.       }
208.       fprintf( outfile,"largest residual : %e \n", largest_residual );
209.
210.                                         /* output phi for image */
211.     if ( IMAGE == 1 )
212.     {
213.       aux = 255.0 / ( phi_max - phi_min );
214.       for ( i=0 ; i<ii ; i++ )
215.         {
216.           for ( j=0 ; j<jj ; j++)
217.                printf( " %f ", aux * ( phi[i][j] - phi_min ) );
218.         }
219.     }
220. }
221.
222.
223. void Gauss_Siedel_step( int ii, int jj, int iter,
224.            double *diff, double h )
225.                                         /* one Gauss-Siedel iteration */
226. {
227.       int    i, j;
228.       double difference, aux, chi;
229.
230.                                         /* over-relaxation control */
231.       chi = 1.0;
232.       if ( iter > 30 ) chi = 1.2;
233.       if ( iter > 50 ) chi = 1.8;
234.                                         /* Gauss-Siedel update */
235.       difference = 0.0;
236.       for ( i=0 ; i<ii ; i++ )
237.         {
238.         for ( j=0 ; j<jj ; j++)
239.           {
240.           if ( ibound[i][j] == 0 )
241.             {
242.               aux = phi[i][j];
243.               phi[i][j] = 0.25 * ( phi[i+1][j] + phi[i-1][j] +
244.                              phi[i][j+1] +phi[i][j-1] - h * h * ff[i][j] );
245.               phi[i][j] = (1-chi)*aux + chi*phi[i][j];
246.               if ( difference <  fabs( aux - phi[i][j] ) )
247.                   difference = fabs( aux - phi[i][j] );
248.             }
249.           }
250.         }
251.       fprintf( stderr, " iter: %d - error %e \n", iter, difference );
252.       *diff=difference;
253. }