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. }

Raw Paste


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