// Michel Vallieres  //

                          // Fall 2010         //
                          // Poisson equation  // 

// use                              //
//   ./p | ./plot_image.py -s 64 64 //


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

#define  real double

#define  N 64                // grid size -- adjust to need

#define  SIZE 10.0           // domain size - arbitrary
#define  PHI0  5.0           // phi on rhs of domain
#define  MAXIT 50000         // maximum # iterations
#define  EPS   1e-8          // required accuracy new-old field

#define  DEBUG 0

// gaussian function (source component)
real gauss( real x, real y, real X, real Y )
{
  return exp( -( (x-X)*(x-X) + (y-Y)*(y-Y) ) );
}


// RHS of Poisson eq - sum of 5 gaussians
real source( real x, real y )
{
  double value;
  value =  - 6.0*gauss(x,y,5.0,5.0) + 2.0*gauss(x,y,2.0,5.0)
           + 3.0*gauss(x,y,5.0,8.0) + 4.0*gauss(x,y,8.0,5.0) 
           + 5.0*gauss(x,y,5.0,2.0);
  return value;
}


// domain setup - a square
// Dirichlet boundary conditions
// 0 left, PHI0 on right, linear interpolation top-bottom
void domain_bc_initial( real *delta, real phi[N][N], real rho[N][N] )
{
  int i, j;
  real x, y, dxdy;

  dxdy = SIZE / (real)( N - 1 );
  *delta = dxdy;

  for ( i=0; i<N; i++ )
    for ( j=0; j<N ; j++ )
      {
	x = j * dxdy;;
	y = i * dxdy;
	rho[i][j] = source( x, y );
      }

   for ( j=0; j<N ; j++ )
     {
       	x = j * dxdy;;
        phi[0][j] = x*PHI0/SIZE;
        phi[N-1][j] = phi[0][j];
     }
  for ( i=0; i<N; i++ )
    {
      phi[i][0] = 0.0;
      phi[i][N-1] = PHI0;
    }
  for ( i=1; i<N-1; i++ )
    for ( j=1; j<N-1 ; j++ )
      phi[i][j] = phi[0][j];
}


// print or plot (pipe in plot_image.py)
void print_function( real phi[N][N], int plot )
{
  int i, j;

  if ( plot==1 )
    {
      for ( i=0; i<N; i++ )
	for ( j=0; j<N ; j++ )
	  printf( "%f\n", phi[i][j] );
    }
  else
    {
  for ( i=0; i<N; i++ )
    {
      for ( j=0; j<N ; j++ )
	printf( "%f ", phi[i][j] );
      printf( "\n" );
    }
  printf( "\n" );
    }
}


// find solution via Gauss-Siedel iterations
// SOR (mix of new/old field) to improve convergence
int Gauss_Siedel_SOR( real delta, real phi[N][N], real rho[N][N] )
{
  int it, i, j;
  real diff, old, new;
  real SORfactor;

  SORfactor = 0.5;                        
  for ( it=1; it<=MAXIT; it++ )               // iterations
    {
      if ( it > 50 )   SORfactor = 1.0;       // arbitrary SOR
      if ( it > 500 )   SORfactor = 1.2;      // convergence factor
      if ( it > 1000 )   SORfactor = 1.4;
      if ( it > 1500 )   SORfactor = 1.6;
      diff = 0.0;                             // accuracy check
      for ( i=1; i<N-1; i++ )
	for ( j=1; j<N-1 ; j++ )
	  {
            old = phi[i][j];                  // finite difference
	    new = (  phi[i+1][j] + phi[i-1][j] +
                     phi[i][j+1] + phi[i][j-1] - 
                     rho[i][j] * delta * delta ) / 4.0;
	                                      // SOR mix of fields
            phi[i][j] = (1.0-SORfactor)*old + SORfactor*new;
	    if ( it%10 == 0 )                 // diff old-new fields
	      if ( fabs(phi[i][j]-old) > diff ) 
		diff = fabs(phi[i][j]-old);
	  }
      if ( it%10 == 0 )                       // convergence check
	  if ( it%10 == 0 && diff < EPS ) break;
    }
  return it;
}


// check solution by back substitution & finite difference
// to evaluate the lhs
int check_solution( real delta, real phi[N][N], real rho[N][N] )
{
  int it, i, j;
  real largest_error;
  real lhs, rhs, diff;

  largest_error = 0.0;
  for ( i=1; i<N-1; i++ )
    for ( j=1; j<N-1 ; j++ )
      {
	lhs = ( phi[i+1][j] + phi[i-1][j] +
                phi[i][j+1] + phi[i][j-1] - 
		4.0*phi[i][j] ) / (delta*delta);
	rhs = rho[i][j];
	diff = fabs( lhs - rhs );
	if ( diff > largest_error )
	  largest_error = diff;
       }
  fprintf( stderr, " Largest error: %g \n\n", largest_error );
}


// main code - overall logic
int main()
{
  real phi[N][N], rho[N][N];
  real delta;
  int iterations;

  domain_bc_initial( &delta, phi, rho );

  if ( DEBUG == 1 ) print_function( phi, 0 );

  if ( DEBUG == 1 ) print_function( rho, 0 );
  iterations = Gauss_Siedel_SOR( delta, phi, rho );

  print_function( phi, 1 );

  fprintf( stderr, "\n No of iterations:  %d\n\n", iterations );

  check_solution( delta, phi, rho );

}