/*
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
Second version - intermediate step to parallel version
------------ --------
( imin & imax in ii direction for domain decomposition )
Build & use of code
gcc poisson_scalar_modular_2.c -o poisson_scalar_modular_2
poisson_scalar_modular_2 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 MAX_x 350 /* maximum dimensions */
#define MAX_y 350
#define EPS 1.0e-7 /* tolerance in solution */
#define IMAGE 1 /* image - yes (1) / no (0) */
/* global array & file for solution */
int ibound[MAX_x][MAX_y];
double phi[MAX_x][MAX_y], ff[MAX_x][MAX_y];
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 );
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 );
}
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" ); }
*N_iterations
= atoi(argv
[1]); *imin = 0;
*imax = *ii;
if ( *ii != *jj )
{
fprintf( stderr
, " Code only setup for N_x = N_y \n"); }
if ( *ii > MAX_x || *jj > MAX_y)
{
fprintf( stderr
, " Grid dimensions too large \n"); }
*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;
/* 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_2.dat","w+"))==NULL
) { fprintf(stderr
,"error in file open \n"); }
fprintf( outfile
, " Grid size: %d x %d \n", ii
, jj
);
/* output initial phi field */
fprintf( outfile
, " Initial - PHI field \n"); for ( i=imin ; i<imax ; i++ )
{
for ( j=0 ; j<jj ; j++)
{
fprintf( outfile
, " %f", phi
[i
][j
] ); }
}
for ( i=imin ; i<imax ; i++ )
{
for ( j=0 ; j<jj ; j++)
{
fprintf( outfile
, " %d", ibound
[i
][j
] ); }
}
}
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
] ); }
}
/* 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=0 ; i<ii ; 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;
}