#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "mpi.h"
#include "header.h"
#ifndef FUNCTIONS_C
#define FUNCTIONS_C
/*
* compilation of the used functions to solve the laplace equation.
* variables are named consistently with the input parameters, thus
* further no explanations are given considering these.
* */
/*
* this function performs a check whether
* the directory 'char *in' is available or not.
* */
void check_file(char *in, char *tags){
FILE *file;
if(0==strcmp(tags, "r")){
if(NULL==(file=fopen(in, "r"))){
printf("failed to open %s for reading.\n", in);
abort();
}
fclose(file);
}
else if(0==strcmp(tags, "w")){
if(NULL==(file=fopen(in, "w"))){
printf("failed to open %s for writing.\n", in);
abort();
}
fclose(file);
}
else if(0==strcmp(tags, "rw")){
if(NULL==(file=fopen(in, "r"))){
if(NULL==(file=fopen(in, "w"))){
printf("failed to open %s for writing.\n", in);
abort();
}
}
else{
printf("%s already exists.\n", in);
abort();
}
}
else if(0==strcmp(tags, "a")){
if(NULL==(file=fopen(in, "a"))){
printf("failed to open %s for appending\n", in);
abort();
}
fclose(file);
}
}
int map(int ix, int *x, int *y, int X, int Y){
int i=0,j=0,k=-1;
while(i<X && ix!=k){
j=0;
while(j<Y && ix!=k){
k=j*X+i;
j++;
}
i++;
}
if(k==ix){
*(x)=i-1;
*(y)=j-1;
return 0;
}
else{
printf("int map(%d,...,X=%d, Y=%d) failed to find coordinates\n", ix, X, Y);
return 1;
}
}
/*
* check if (ix, iy) is set to a potential
* */
int check_potential(int ix, int iy, int NPHI, pot *phi0){
int n;
for(n=0; n<NPHI; n++){
if(ix>=phi0[n].pos[0] && ix<=phi0[n].pos[1] && iy>=phi0[n].pos[2] && iy<=phi0[n].pos[3])
return 1;
}
return 0;
}
/*
* check if (ix, iy) is inside of the grounded area
*/
int check_ground(int ix, int iy, int NGROUND, int **ground){
int n;
for(n=0; n<NGROUND; n++){
if(ix>ground[n][0] && ix <ground[n][1] && iy>ground[n][2] && iy<ground[n][3])
return 1;
}
return 0;
}
/*
* average over grounded lattice sites. separate grounded areas are averaged independently.
*
* NOTE: only the boundary of the grounded sites is computed void jacobi() thus only these lattice sites
* may contribute to the average. the interior of the grounded areas is set to this average.
* */
void ground_it(int NGROUND, int **ground, int X, double *phi){
int j, x, y, N=0;
double buf;
for(j=0; j<NGROUND; j++){
for(N=0, buf=0., x=ground[j][0], y=ground[j][2]; x<=ground[j][1]; x++, N++)
buf+=phi[y*X+x];
for(x=ground[j][0], y=ground[j][3]; x<=ground[j][1]; x++, N++)
buf+=phi[y*X+x];
for(y=ground[j][2]+1, x=ground[j][0]; y<=ground[j][3]-1; y++, N++)
buf+=phi[y*X+x];
for(y=ground[j][2]+1, x=ground[j][1]; y<=ground[j][3]-1; y++, N++)
buf+=phi[y*X+x];
buf/=(double)N;
for(x=ground[j][0]; x<=ground[j][1]; x++){
for(y=ground[j][2]; y<=ground[j][3]; y++){
phi[y*X+x]=buf;
}
}
}
}
/*
* Implementation of the jacobi algorithm to solve a 2d-laplace equation.
* char *dif is the output directory for the precision and iteration-number.
* char *logf is the output directory for the total required number of iterations,
* the lattice volume of the solution, the total time for the solution and its actual
* precision.
* */
void jacobi(int X, int Y, double NP, double *phi, int NPHI,
pot *phi0, int NGROUND, int **ground, char *dif, char *logf){
double *dum, pres=1.+NP, buf=0.;
int x, y, i, ix;
double t1, t2;
FILE *file;
/*MPI_Wtime() returns the current time of the process in seconds with double-precision. meaning
* that this function allows for a much more precise estimate for the required time than any function
* defined in time.h where the highest precision is in seconds. For example consider a lattice of the size 8x8,
* the solution of this lattice is computed in the order of 10e-4 seconds.
* */
t1 = MPI_Wtime();
/*
* double *dum is used to store the new solution of the laplace equation.
* */
dum=malloc(X*Y*sizeof(double));
for(ix=0; ix<X*Y; ix++)
dum[ix]=phi[ix];
check_file(dif, "w");
file=fopen(dif, "w");
i=0;
/*
* loop over either the precision (pres) if NP<1.
* or the iteration number (i) if NP>1.
* */
while((i<(int)NP && NP>1.) || (pres>NP && NP<1.)){
for(x=1; x<X-1; x++){/*loop over the volume of the lattice, excluding the boundary*/
for(y=1; y<Y-1; y++){
ix=y*X+x;
/*
* check if (x,y) is inside of the grounded area or if it is
* predefined by a potential.
* */
if(check_potential(x, y, NPHI, phi0)==0 && check_ground(x, y, NGROUND, ground)==0){
dum[ix]=(phi[(y-1)*X+x]+phi[(y+1)*X+x]+phi[y*X+(x-1)]+phi[y*X+(x+1)])/4.;
/*dum[ix]=(1.-r)*phi[ix]+r*0.25*(phi[(y-1)*X+x]+phi[(y+1)*X+x]+phi[y*X+(x-1)]+phi[y*X+(x+1)]);*/
}
}
}/*end loop over lattice-volume
average over the grounded area
*/
ground_it(NGROUND, ground, X, dum);
buf=0.;/*buffer for the precision*/
for(x=1; x<X-1; x++){
for(y=1; y<Y-1; y++){
ix=y*X+x;
if(fabs(phi[ix]-dum[ix])>buf)
buf=fabs(phi[ix]-dum[ix]);/*find the maximum absolute change from the previous step to this one*/
phi[ix]=dum[ix];/*copy the new solution*/
}
}
fprintf(file, "%d %.6e\n", i, buf);/*output for the current precision and iteration-number*/
pres=buf;
i++;
}
fclose(file);
/*
* logfile output: lattice extent, method, volume, iteration-number, precision and total time.
* */
check_file(logf, "a");
file=fopen(logf, "a");
fprintf(file, "###################################################\n");
fprintf(file, "#X = %d\t Y = %d\t jacobi\n", X, Y);
fprintf(file, "#\tX*Y\tN\tP\tR\tTime\n");
t2 = MPI_Wtime();
fprintf(file, "%8d %8d %.2e %.4e\n", X*Y, i-1, pres, t2-t1);
fclose(file);
free(dum);/*free memory required for the temporary array*/
}
/**
* Implementation of the red-black algorithm for the solution of the 2d-laplace equation.
* Since it is largely the same as the jacobi-algorithm only changes to the algorithm are commented.
*
* NOTE: grounded areas are not computed with this method.
*/
void red_black(int X, int Y, double NP, double *phi, int NPHI,
pot *phi0, char *dif, char *logf){
double *dum, pres=1.+NP, buf=0.;
int x, y, i,ix, eo;
double t1, t2;
FILE *file;
t1 = MPI_Wtime();
dum=malloc(X*Y*sizeof(double));
for(ix=0; ix<X*Y; ix++)
dum[ix]=phi[ix];
check_file(dif, "w");
file=fopen(dif, "w");
i=0;
while((i<(int)NP && NP>1.) || (pres>NP && NP<1.)){
/*loop over even and odd lattice sites, which are ordered like:
* eoeoeo
* oeoeoe
* eoeoeo
* oeoeoe
* */
for(eo=0; eo<2; eo++){
for(x=1; x<X-1; x++){
for(y=1; y<Y-1; y++){
if((x+y)%2==eo && check_potential(x, y, NPHI, phi0)==0){/*compute either even or odd lattice sites*/
dum[y*X+x]=(phi[(y-1)*X+x]+phi[(y+1)*X+x]+phi[y*X+(x-1)]+phi[y*X+(x+1)])/4.;
}
}
}
/*partial update of either even or odd lattice sites.*/
for(x=1, buf=0.; x<X-1; x++){
for(y=1; y<Y-1; y++){
if((x+y)%2==eo){
ix=y*X+x;
if(fabs(phi[ix]-dum[ix])>buf)
buf=fabs(phi[ix]-dum[ix]);
phi[ix]=dum[ix];
}
}
}
}
fprintf(file, "%d %.6e\n", i, buf);
pres=buf;
i++;
}
fclose(file);
check_file(logf, "a");
file=fopen(logf, "a");
fprintf(file, "###################################################\n");
fprintf(file, "#X = %d\t Y = %d\t red-black\n", X, Y);
fprintf(file, "#\tX*Y\tN\tP\tR\tTime\n");
t2 = MPI_Wtime();
fprintf(file, "%8d %8d %.2e %.4e\n", X*Y, i-1, pres, t2-t1);
fclose(file);
free(dum);
}
/**
* Implementation of the gauss-seidel algorithm for the solution of the 2d-laplace equation.
* Since it is largely the same as the jacobi-algorithm only changes to the algorithm are commented.
*
* NOTE: grounded areas are not computed with this method.
*/
void gauss_seidel(int X, int Y, double NP, double *phi, int NPHI,
pot *phi0, char *dif, char *logf){
double pres=1.+NP, buf=0., buf2=0.;
int x, y, i, ix;
double t1, t2;
FILE *file;
t1 = MPI_Wtime();
check_file(dif, "w");
file=fopen(dif, "w");
i=0;
while((i<(int)NP && NP>1.) || (pres>NP && NP<1.)){
for(y=1, buf2=0.; y<Y-1; y++){
for(x=1; x<X-1; x++){
ix=y*X+x;
if(check_potential(x, y, NPHI, phi0)==0){
/*the buffer is only needed to compute the current precision*/
buf=(phi[(y-1)*X+x]+phi[(y+1)*X+x]+phi[y*X+(x-1)]+phi[y*X+(x+1)])/4.;
if(fabs(phi[ix]-buf)>buf2){
buf2=fabs(phi[ix]-buf);
}
phi[ix]=buf;/*update of the new step within the same loop as the actual computation.*/
}
}
}
pres=buf2;
fprintf(file, "%d %.6e\n", i, pres);
i++;
}
fclose(file);
check_file(logf, "a");
file=fopen(logf, "a");
fprintf(file, "###################################################\n");
fprintf(file, "#X = %d\t Y = %d\t gauss-seidel\n", X, Y);
fprintf(file, "#\tX*Y\tN\tP\tR\tTime\n");
t2 = MPI_Wtime();
fprintf(file, "%8d %8d %.2e %.4e\n", X*Y, i-1, pres, t2-t1);
fclose(file);
}
#endif