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

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