src/poisson_2d/P_2d.c

   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 functions
  13
  14
  15     Build & use of code
  16
  17 gcc poisson_2d.c -lm -o  poisson_2d
  18
  19 ./poisson_2d  N  N_x  N_y | plot_image.py -s N_x N_y -t "2d Poisson"
  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 to record solution */
  36 int   ibound[MAX_x][MAX_y];
  37 double phi[MAX_x][MAX_y], ff[MAX_x][MAX_y];
  38
  39
  40 /* lattice definition & iteration number */
  41 void set_grid ( int argc, char *argv[], int *N_iterations,
  42                 int *NI, int *NJ, double *h )
  43 {
  44     if ( argc != 4 )
  45       {
  46         printf( " Syntax: poisson  N  N_x  N_y   ( where N_x = N_y ) \n" );
  47         exit(1);
  48       }
  49     *N_iterations = atoi(argv[1]);
  50     *NI = atoi(argv[2]);
  51     *NJ = atoi(argv[3]);
  52
  53     if ( *NI != *NJ )
  54       {
  55         fprintf( stderr, " Code only setup for N_x = N_y \n");
  56         exit(1);
  57       }
  58
  59     if ( *NI > MAX_x || *NJ > MAX_y)
  60       {
  61         fprintf( stderr, " Grid dimensions too large \n");
  62         exit(1);
  63       }
  64
  65     *h = 1.0/(double)(*NI-1);
  66 }
  67
  68
  69 /* initial phi values, source & boundary conditions */
  70 void initial_source_and_boundary ( int NI, int NJ )
  71 {
  72     int i, j, i_s;
  73     double r2;
  74                                         /* initial fields & source*/
  75     for ( i=0 ; i<NI ; i++ )
  76       {
  77         for ( j=0 ; j<NJ ; j++)
  78           {
  79             r2 = (i-NI/5)*(i-NI/5) + (j-3*NJ/4)*(j-3*NJ/4);
  80             ff[i][j] = 7000 * exp( - 0.01 * r2 );
  81             ibound[i][j] = 0;
  82             phi[i][j] = 0.0;
  83           }
  84       }
  85                                         /* Boundary Conditions */
  86     for ( j=0 ; j<NJ ; j++)
  87       {
  88         ibound[0][j] = 1;
  89         phi[0][j] = 2.0;
  90         ibound[NI-1][j] = 1;
  91         phi[NI-1][j] = 2.0;
  92       }
  93     for ( i=0 ; i<NI ; i++ )
  94       {
  95         ibound[i][0] = 1;
  96         phi[i][0] = 2.0;
  97         ibound[i][NJ-1] = 1;
  98         phi[i][NJ-1] = 2.0;
  99       }
 100                                         /* single values somewhere in square */
 101     i_s = 4*NI/5;
 102     ibound[i_s][3*NJ/4] = 1;
 103     phi[i_s][3*NJ/4] = -25.0;
 104     i_s = NI/3;
 105     ibound[i_s][NJ/3] = 1;
 106     phi[i_s][NJ/3] = -10.0;
 107
 108
 109                                         /* echo input */
 110     fprintf( stderr, " Grid size: %d x %d \n", NI, NJ );
 111
 112 }
 113
 114
 115 /* check solution & output field data for image */
 116 void output_results( int NI, int NJ, double h )
 117 {
 118     int    i, j;
 119     double aux, phi_min, phi_max, residual, largest_residual;
 120
 121                                         /* check if solution is */
 122                                         /* really a solution */
 123     fprintf( stderr, " Check - ( LHS - RHS ) \n");
 124     largest_residual = 0.0;
 125     for ( i=0 ; i<NI ; i++ )
 126       {
 127       for ( j=0 ; j<NJ ; j++)
 128         {
 129                                         /* min & max of phi as an aside */
 130           if ( phi[i][j] < phi_min ) phi_min = phi[i][j];
 131           if ( phi[i][j] > phi_max ) phi_max = phi[i][j];
 132
 133           if ( ibound[i][j] == 0 )
 134             {
 135             residual =  phi[i+1][j] + phi[i-1][j] +
 136                              phi[i][j+1] + phi[i][j-1] -
 137                                4.0 * phi[i][j] - h * h * ff[i][j];
 138             }
 139           else
 140             {
 141             residual = 0.0;
 142             }
 143           if( fabs( residual ) > largest_residual )
 144                           largest_residual = fabs( residual );
 145         }
 146       }
 147
 148       fprintf( stderr, "largest residual : %e \n", largest_residual );
 149
 150                                         /* output phi for image */
 151     if ( IMAGE == 1 )
 152     {
 153       aux = 255.0 / ( phi_max - phi_min );
 154       for ( i=0 ; i<NI ; i++ )
 155         {
 156           for ( j=0 ; j<NJ ; j++)
 157                printf( " %f \n", aux * ( phi[i][j] - phi_min ) );
 158         }
 159     }
 160 }
 161
 162
 163 /* one Gauss-Siedel iteration */
 164 void Gauss_Siedel_step( int NI, int NJ, int iter,
 165            double *diff, double h )
 166 {
 167       int    i, j;
 168       double difference, aux, chi;
 169
 170                                         /* over-relaxation control */
 171       chi = 1.0;
 172       if ( iter > 30 ) chi = 1.2;
 173       if ( iter > 50 ) chi = 1.8;
 174                                         /* Gauss-Siedel update */
 175       difference = 0.0;
 176       for ( i=0 ; i<NI ; i++ )
 177         {
 178         for ( j=0 ; j<NJ ; j++)
 179           {
 180           if ( ibound[i][j] == 0 )
 181             {
 182               aux = phi[i][j];
 183               phi[i][j] = 0.25 * ( phi[i+1][j] + phi[i-1][j] +
 184                              phi[i][j+1] +phi[i][j-1] - h * h * ff[i][j] );
 185               phi[i][j] = (1-chi)*aux + chi*phi[i][j];
 186               if ( difference <  fabs( aux - phi[i][j] ) )
 187                   difference = fabs( aux - phi[i][j] );
 188             }
 189           }
 190         }
 191       fprintf( stderr, " iter: %d - error %e \n", iter, difference );
 192       *diff=difference;
 193 }
 194
 195
 196
 197
 198 int main ( int argc, char *argv[] )
 199 {
 200     double h;
 201     int    NI, NJ, i, j;
 202     double diff, aux, chi;
 203     int    N_iterations, iter;
 204
 205                                         /* lattice definition & */
 206                                         /* iteration number */
 207     set_grid ( argc, argv, &N_iterations, &NI, &NJ, &h );
 208                                         /* initial fields, source */
 209                                         /* & boundary conditions */
 210     initial_source_and_boundary ( NI, NJ );
 211                                         /* iterate the solution */
 212     for ( iter=1 ; iter<N_iterations ; iter++ )
 213       {
 214                                         /* Gauss-Siedel iteration */
 215         Gauss_Siedel_step( NI, NJ, iter, &diff, h );
 216                                         /* solution accurate enough ? */
 217         if ( diff < EPS ) break;
 218       }
 219                                         /* check & output the phi field */
 220     output_results ( NI, NJ, h );
 221
 222 }
 223