C   19

poisson scalar modular 2d version 2 c

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

Raw Paste


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