/* 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 #include #include #include #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 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 30 ) chi = 1.2; if ( iter > 50 ) chi = 1.8; /* Gauss-Siedel update */ difference = 0.0; for ( i=imin ; i