/*
    Over-relaxation solution to Poisson equation

    Syntax
         Poisson  N  N_x  N_y
            ( Max # iteration, grid in x, grid in y )

    Note: only written for equal x and y meshes -- grid(x) = grid(y)

    IBOUND:  positions on grid corresponding to boundaries

    Modular code - simple main code - logical subroutines

    Third version - intermediate step to parallel version
                    ------------         --------

    ( imin & imax in ii direction for domain decomposition )
    ( better handling of arrays )


    Build & use of code

gcc poisson_scalar_modular_3.c NR_utilities.c -o  poisson_scalar_modular_3

poisson_scalar_modular_3  N  N_x  N_y | make_image -s N_x N_y | convert - - | xv -
 
*/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

#define EPS 1.0e-7                      /* tolerance in solution */

#define IMAGE 1                         /* image - yes (1) - no (0) */

int    **ibound;                        /* global array & file for solution */
double **phi, **ff;
FILE   *outfile;
                                        /* functions prototypes */
void set_grid ( int argc, char *argv[], int *N_iterations, 
                int *imin, int *imax, int *ii, int *jj, double *h );
void initial_and_boundary ( int imin, int imax, int ii, int jj );
void output_results ( int imin, int imax, int ii, int jj, double h );
void Gauss_Siedel_step( int imin, int imax, int ii, int jj, 
                        int iter, double *diff, double h );

                                        /* numerical routines prototypes */
                                        /* from nrutil.h */
double **dmatrix(long nrl, long nrh, long ncl, long nch);
int **imatrix(long nrl, long nrh, long ncl, long nch);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);



int main ( int argc, char *argv[] )
{
    double h;
    int    imin, imax, ii, jj, i, j;
    double diff, aux, chi;
    int    N_iterations, iter;

                                        /* lattice definition & 
                                                       iteration number */
    set_grid ( argc, argv, &N_iterations, &imin, &imax, &ii, &jj, &h );

                                        /* initial fields & 
                                                       boundary conditions */
    initial_and_boundary ( imin, imax, ii, jj );


                                        /* iterate the solution */
    for ( iter=1 ; iter<N_iterations ; iter++ )
      {
	                                /* Gauss-Siedel iteration */
        Gauss_Siedel_step( imin, imax, ii, jj, iter, &diff, h );
                                        /* solution accurate enough ? */
        if ( diff < EPS ) break;
      }

                                        /* output the phi field */
    output_results ( imin, imax, ii, jj, h );

                                        /* free space allocated to */
                                        /* matrices -- good habit  */
    free_dmatrix( phi, imin, imax, 0, jj );
    free_dmatrix( ff,  imin, imax, 0, jj );
    free_imatrix( ibound,  imin, imax, 0, jj );

}


void set_grid ( int argc, char *argv[], int *N_iterations, 
                int *imin, int *imax, int *ii, int *jj, double *h )
{
                                        /* lattice definition & 
                                                        iteration number */
    if ( argc != 4 ) 
      {
        printf( " Syntax: poisson  N  N_x  N_y   ( where N_x = N_y ) \n" );
        exit(1);
      }
    *N_iterations = atoi(argv[1]);
    *ii = atoi(argv[2]);
    *imin = 0;
    *imax = *ii;
    *jj = atoi(argv[3]);
    if ( *ii != *jj ) 
      {
        fprintf( stderr, " Code only setup for N_x = N_y \n");
        exit(1);
      }
    *h = 1.0/(double)(*ii-1);
}


void initial_and_boundary ( int imin, int imax, int ii, int jj )
{
                                        /* initial phi values & 
                                                       boundary conditions */
    int i, j, i_s;

                                        /* make space for the field matrices */
    phi = dmatrix( imin, imax, 0, jj );
    ff = dmatrix( imin, imax, 0, jj );
    ibound = imatrix( imin, imax, 0, jj );

                                        /* initial fields */
    for ( i=imin ; i<imax ; i++ )
      {
        for ( j=0 ; j<jj ; j++)
	  {
            ff[i][j] = 0.0;
            ibound[i][j] = 0;
            phi[i][j] = 0.0;
	  } 
      }
                                        /* Boundary Conditions */
    for ( j=0 ; j<jj ; j++)
      {
        ibound[imin][j] = 1;
        phi[imin][j] = 2.0;
        ibound[imax-1][j] = 1;
        phi[imax-1][j] = 2.0;
      }
    for ( i=imin ; i<imax ; i++ )
      {
        ibound[i][0] = 1;
        phi[i][0] = 2.0;
        ibound[i][jj-1] = 1;
        phi[i][jj-1] = 2.0;
      }
                                     /* single value in center of square */
    i_s = 4*ii/5;
    ibound[i_s][3*jj/4] = 1;
    phi[i_s][3*jj/4] = -25.0;
    i_s = ii/3;
    ibound[i_s][jj/3] = 1;
    phi[i_s][jj/3] = -10.0;

                                        /* open output file */
    if ((outfile=fopen("Field_3.dat","w+"))==NULL)
	{ fprintf(stderr,"error in file open \n");
	  exit(1);
	}

    fprintf( outfile, " Grid size: %d x %d \n", ii, jj );

                                        /* output initial phi field */
    fprintf( outfile, " Initial - field \n");
    for ( i=imin ; i<imax ; i++ )
      {
      for ( j=0 ; j<jj ; j++)
        {
            fprintf( outfile, " %f", phi[i][j] );
        }
        fprintf( outfile, "\n");
      }

    fprintf( outfile, " Boundary  \n");
    for ( i=imin ; i<imax ; i++ )
      {
      for ( j=0 ; j<jj ; j++)
        {
            fprintf( outfile, " %d", ibound[i][j] );
        }
        fprintf( outfile, "\n");
      }
}


void output_results( int imin, int imax, int ii, int jj, double h )
                                        /* output results */
{
    int    i, j;
    double aux, phi_min, phi_max, residual, largest_residual;

                             /* output results */
    fprintf( outfile, " Final - PHI field \n");
    for ( i=imin ; i<imax ; i++ )
      {
      for ( j=0 ; j<jj ; j++)
        {
            fprintf( outfile, " %f", phi[i][j] );
        }
        fprintf( outfile, "\n");
      }
                                        /* check if solution is */
                                        /* really a solution */
    fprintf( outfile, " Check - ( LHS - RHS ) \n");
    largest_residual = 0.0;
    for ( i=imin ; i<imax ; i++ )
      {
      for ( j=0 ; j<jj ; j++)
        {
                                        /* min & max of phi as an aside */
          if ( phi[i][j] < phi_min ) phi_min = phi[i][j];
          if ( phi[i][j] > phi_max ) phi_max = phi[i][j];

          if ( ibound[i][j] == 0 )
	    {
            residual =  phi[i+1][j] + phi[i-1][j] + 
                             phi[i][j+1] + phi[i][j-1] -
                               4.0 * phi[i][j] - h * h * ff[i][j];
	    }
          else
	    {
            residual = 0.0;
	    }
          if( fabs( residual ) > largest_residual ) 
                          largest_residual = fabs( residual );
        }
      }
      fprintf( outfile,"largest residual : %e \n", largest_residual );

                                        /* output phi for image */
    if ( IMAGE == 1 )
      {
        aux = 255.0 / ( phi_max - phi_min );
        for ( i=imin ; i<imax ; i++ )
          {
             for ( j=0 ; j<jj ; j++)
                  printf( " %f ", aux * ( phi[i][j] - phi_min ) );
          }
      }
}


void Gauss_Siedel_step( int imin, int imax, int ii, int jj, int iter, 
           double *diff, double h )
                                        /* one Gauss-Siedel iteration */
{
      int    i, j;
      double difference, aux, chi;

                                        /* over-relaxation control */
      chi = 1.0;
      if ( iter > 30 ) chi = 1.2;
      if ( iter > 50 ) chi = 1.8;
	                                /* Gauss-Siedel update */
      difference = 0.0;
      for ( i=imin ; i<imax ; i++ )
        {
        for ( j=0 ; j<jj ; j++)
	  {
          if ( ibound[i][j] == 0 )
	    {
              aux = phi[i][j];
              phi[i][j] = 0.25 * ( phi[i+1][j] + phi[i-1][j] + 
                             phi[i][j+1] +phi[i][j-1] - h * h * ff[i][j] );
              phi[i][j] = (1-chi)*aux + chi*phi[i][j];
              if ( difference <  fabs( aux - phi[i][j] ) )
                  difference = fabs( aux - phi[i][j] );    
	    }
	  }
        }
      fprintf( stderr, " iter: %d - error %e \n", iter, difference );
      *diff=difference;
}