#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "header.h"
#include "mpi.h"
/*global variables:
* X...........extent of the finest lattice in x-direction
* Y...........extent of the finest lattice in y-direction
* XL..........extent of the coarsest lattice in x-direction
* YL..........extent of the coarsest lattice in y-direction
* L...........size of the box, set to max(XL-1,YL-1)[later on]
* NP..........# iterations used to solve the laplace equation or precision
* NPHI........# of different potentials set at the start
* NGROUND.....# of different grounded areas
* outdir[]....output directory
* method[]....method specified in the infile [jacobi / red-black / gauss-seidel]
* exerc[].....exercise-label
* *phi0.......pointer to the start-values of the potential
* **ground....grounded areas.=>ground[#grounded areas][coordinate]
* example: ground[0][x1],ground[0][x2],ground[0][y1],ground[0][y2];
* */
int X, Y, XL, YL, L, F, NPHI, NGROUND, p=0;
char outdir[1000], method[100], exerc[100];
pot *phi0;
int **ground;
double NP;
/*
* reads the input parameters and checks them for validity.
*
* variables are explained in readme.in
*
* */
void read_infile(char *in) {
FILE *file;
int n;
check_file(in, "r");
file = fopen(in, "r");
fscanf(file, "%s\n", exerc);
fscanf(file, "%s\n", outdir);
fscanf(file, "%s\n", method);
if(!(strcmp(method, "gauss-seidel")==0 || strcmp(method, "jacobi")==0 || strcmp(method, "red-black")==0)){
printf("%s is not a valid method.\n", method);
abort();
}/*abort if the method is unknown.*/
fscanf(file, "X %d\n", &X);
fscanf(file, "Y %d\n", &Y);
fscanf(file, "%*s %d\n", &F);
fscanf(file, "%*s %d\n", &XL);
fscanf(file, "%*s %d\n", &YL);
fscanf(file, "%*s %lf\n", &NP);
fscanf(file, "%*s %d\n", &NPHI);
phi0 = malloc(NPHI * sizeof(pot));
for(n=0; n<NPHI; n++){
fscanf(file, "%d %d %d %d %lf\n", &phi0[n].pos[0], &phi0[n].pos[1],&phi0[n].pos[2],&phi0[n].pos[3],&phi0[n].val);
if(phi0[n].pos[0]<0 || phi0[n].pos[1]>XL || phi0[n].pos[2]<0 || phi0[n].pos[3]>YL){
printf("X=0..%d Y=0..%d\tpotential outside boundary of the small grid\n", XL-1, YL-1);
abort();
}/*abort if the potential is outside of the boundaries of the coarsest lattice*/
}
fscanf(file, "%*s %d\n", &NGROUND);
if(NGROUND>0 && (strcmp(method, "red-black")==0 || strcmp(method, "gauss-seidel")==0)){
printf("Grounded areas are not implemented for the %s method.\n", method);
abort();
}
ground=malloc(NGROUND*sizeof(int *));
for(n=0; n<NGROUND; n++){
ground[n]=malloc(4*sizeof(int));
fscanf(file, "%d %d %d %d\n", &ground[n][0],&ground[n][1],&ground[n][2],&ground[n][3]);
if(ground[n][0]<0 || ground[n][1]>=XL || ground[n][2]<0 || ground[n][3]>=YL){
printf("X=0..%d Y=0..%d\tgrounded area outside or on the boundary\n", XL, YL);
abort();
}/*grounded areas must be inside of the box, they may not overlap with the boundary*/
}
fscanf(file, "plot %d", &p);
fclose(file);
/*prints the infile in the terminal for checks.*/
printf("%s\n", exerc);
printf("%s\n", outdir);
printf("%s\n", method);
printf("X %d\n", X);
printf("Y %d\n", Y);
printf("F %d\n", F);
printf("XL %d\n", XL);
printf("YL %d\n", YL);
printf("NP %.2e\n", NP);
printf("npotential %d\n", NPHI);
for(n=0; n<NPHI; n++){
printf("phi[%d]=% .2e at %d %d %d %d\n", n, phi0[n].val, phi0[n].pos[0], phi0[n].pos[1], phi0[n].pos[2],phi0[n].pos[3]);
}
printf("nground %d\n", NGROUND);
for(n=0; n<NGROUND; n++){
printf("ground[%d] at %d %d %d %d\n", n, ground[n][0],ground[n][1],ground[n][2],ground[n][3]);
}
}
int main(int arg, char *args[]) {
int i, x, y, Xstep, Ystep, ix, iy, step;
double *phi, *phibuf, *range, x1, x2, y1, y2, buf;
char out[1000], title[1000], xlabel[100], ylabel[100], plot[1000], dif[1000];
static int my_rank;
FILE *file;
plotinfo *pi;
if(arg!=2){
printf("incorrect number of input parameters.\n");
abort();
}
read_infile(args[1]);/*get input data and test it*/
/*initialization of the MPI-process, required to measure processor-time*/
MPI_Init(&arg, &args);
phi=malloc(XL*YL*sizeof(double));
pi=malloc(3*sizeof(plotinfo));
range=malloc(6*sizeof(double));
for(i=0; i<6; i++)
range[i]=0.;
/*initialize potential*/
for(x=0; x<XL; x++){
for(y=0; y<YL; y++){
phi[y*XL+x]=0.;
}
}
for(i=0; i<NPHI; i++){
printf("nphi[%d]\n", i);
for(x=phi0[i].pos[0]; x<=phi0[i].pos[1]; x++){
for(y=phi0[i].pos[2]; y<=phi0[i].pos[3]; y++){
phi[y*XL+x]=phi0[i].val;
printf("%d %d %d %f\n", i, x, y, phi[y*XL+x]);
}
}
}
Xstep=XL;/*lattice extent of the current step*/
Ystep=YL;
step=1;
while(Xstep<=X && Ystep<=Y){
printf("solving lattice: %dx%d\n", Xstep, Ystep);
sprintf(dif, "%s/dif_X%d_Y%d_%s_%s_m%d.dat", outdir, Xstep, Ystep,exerc, method, step);
sprintf(out, "%s/%s_%s.log", outdir, exerc, method);
/*compute the solution of the laplace equation using the functions defined in functions.c*/
if(0==strcmp(method, "jacobi"))
jacobi(Xstep, Ystep, NP, phi, NPHI, phi0, NGROUND, ground, dif, out);
else if(0==strcmp(method, "red-black"))
red_black(Xstep, Ystep, NP, phi, NPHI, phi0, dif, out);
else if(0==strcmp(method, "gauss-seidel"))
gauss_seidel(Xstep, Ystep, NP, phi, NPHI, phi0, dif, out);
/*define the plotranges*/
if(Xstep>Ystep)
L=Xstep-1;
else
L=Ystep-1;
range[0]=0, range[1]=(double)(Xstep-1)/L,range[2]=0, range[3]=(double)(Ystep-1)/L;/*plotranges*/
/*generate the output for the potential*/
sprintf(out, "%s/phi_X%d_Y%d_%s_%s.dat", outdir, Xstep, Ystep,exerc, method);
printf("%s\n", out);
check_file(out, "w");
file=fopen(out, "w");
fprintf(file, "#Plotdata for:\n");
fprintf(file, "#%s\n", exerc);
fprintf(file, "#%s\n", outdir);
fprintf(file, "#%s\n", method);
fprintf(file, "#X %d\n", X);
fprintf(file, "#Y %d\n", Y);
fprintf(file, "#NP %.2e\n", NP);
fprintf(file, "#npotential %d\n", NPHI);
for(i=0; i<NPHI; i++){
fprintf(file, "#phi[%d]=% .2e at %d %d %d %d\n", i, phi0[i].val,
phi0[i].pos[0], phi0[i].pos[1], phi0[i].pos[2],phi0[i].pos[3]);
}
fprintf(file, "#nground %d\n", NGROUND);
for(i=0; i<NGROUND; i++){
fprintf(file, "#ground[%d] at %d %d %d %d\n", i,
ground[i][0],ground[i][1],ground[i][2],ground[i][3]);
}
for(x=0; x<Xstep; x++){
for(y=0; y<Ystep; y++){
fprintf(file, "%4f %4f % .6e\n",(double)x/L, (double)y/L, phi[y * Xstep + x]);
}
fprintf(file, "\n");/*to create "set pm3d map"
each matrix line has to be separated by the next by an additional line.*/
}
fclose(file);
if(p==1){
sprintf(pi[0].color, "red");
sprintf(pi[0].data, "phi_X%d_Y%d_%s_%s.dat", Xstep, Ystep,exerc, method);
sprintf(pi[0].label, "$\\Phi(NP=%.2e)$", NP);
pi[0].lw=1;
pi[0].pt=1;
pi[0].x=1;
pi[0].y=2;
sprintf(xlabel, "x[a]");
sprintf(ylabel, "y[a]");
sprintf(title, "$\\Phi(x,y)$ at NP=%.2e and $m=%d$", NP, step);
sprintf(plot, "potential_X%d_Y%d_%s_%s", Xstep, Ystep,exerc, method);
plot_data_color(outdir, plot, title, xlabel, ylabel, range, Xstep, Ystep, 1, pi);
}
/*generate the output for the electric field*/
sprintf(out, "%s/E_X%d_Y%d_%s_%s.dat", outdir, Xstep, Ystep,exerc, method);
printf("%s\n", out);
check_file(out, "w");
file=fopen(out, "w");
for(x=0; x<Xstep-1; x++){
for(y=0; y<Ystep-1; y++){
/*discretized gradient of phi:*/
x1=(phi[y*Xstep+x]+phi[(y+1)*Xstep+x])/2;
x2=(phi[y*Xstep+x+1]+phi[(y+1)*Xstep+x+1])/2;
y1=(phi[y*Xstep+x]+phi[y*Xstep+x+1])/2;
y2=(phi[(y+1)*Xstep+x]+phi[(y+1)*Xstep+x+1])/2;
buf=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))/2;
fprintf(file, "%f %f %f %f\n", (double)(2*x+1)/2,(double)(2*y+1)/2, -(x2-x1)/2/buf, -(y2-y1)/2/buf);
}
fprintf(file, "\n");
}
fclose(file);
if(p==1){
sprintf(pi[0].color, "black");
sprintf(pi[0].data, "E_X%d_Y%d_%s_%s.dat", Xstep, Ystep,exerc, method);
sprintf(pi[0].label, "$\\vec{E}$");
pi[0].lw=3;
pi[0].x=1;
pi[0].y=2;
pi[0].u=3;
pi[0].v=4;
sprintf(xlabel, "x[a]");
sprintf(ylabel, "y[a]");
sprintf(title, "$\\vec{E}(x,y)=-\\nabla\\Phi$ at $m=%d$", step);
sprintf(plot, "E_X%d_Y%d_%s_%s", Xstep, Ystep,exerc, method);
range[0]=0., range[1]=(double)Xstep, range[2]=0., range[3]=(double)Ystep;
plot_data_vecfield(outdir, plot, title, xlabel, ylabel, range, Xstep, Ystep, 1, pi);
}
/*generate the output for the magnetic field*/
sprintf(out, "%s/B_X%d_Y%d_%s_%s.dat", outdir, Xstep, Ystep,exerc, method);
printf("%s\n", out);
check_file(out, "w");
file=fopen(out, "w");
for(x=0; x<Xstep; x++){
for(y=0; y<Ystep; y++){
if((x>0 && x<Xstep-1) && (y>0 && y<Ystep-1)){
/*discretized \nabla x \phi:*/
x1=(phi[(y+1)*Xstep+x]-phi[(y-1)*Xstep+x])/2.;
y1=-(phi[y*Xstep+x+1]-phi[y*Xstep+x-1])/2.;
buf=sqrt((x1*x1+y1*y1));
fprintf(file, "%f %f %f %f\n", (double)x,(double)y, x1/buf, y1/buf);
}
else{
fprintf(file, "%f %f %f %f\n", (double)x,(double)y, 0., 0.);
}
}
fprintf(file, "\n");
}
fclose(file);
if(p==1){
sprintf(pi[0].color, "black");
sprintf(pi[0].data, "B_X%d_Y%d_%s_%s.dat", Xstep, Ystep,exerc, method);
sprintf(pi[0].label, "$\\vec{E}$");
pi[0].lw=3;
pi[0].x=1;
pi[0].y=2;
pi[0].u=3;
pi[0].v=4;
sprintf(xlabel, "x[a]");
sprintf(ylabel, "y[a]");
sprintf(title, "$\\vec{B}(x,y)=\\nabla\\times A$ at $m=%d$", step);
sprintf(plot, "B_X%d_Y%d_%s_%s", Xstep, Ystep,exerc, method);
range[0]=0., range[1]=(double)Xstep, range[2]=0., range[3]=(double)Ystep;
plot_data_vecfield(outdir, plot, title, xlabel, ylabel, range, Xstep, Ystep, 1, pi);
sprintf(out, "rm %s/*.dat", outdir);
printf("%s\n", out);
system(out);
}
/*allocate memory for a temporary array and store the current potential*/
phibuf=malloc(Xstep*Ystep*sizeof(double));
for(x=0; x<Xstep; x++){
for(y=0; y<Ystep; y++){
phibuf[y*Xstep+x]=phi[y*Xstep+x];
}
}
free(phi);/*free the memory used by the current step's potential*/
Xstep*=F;/*increase the lattice-extent to the next level*/
Ystep*=F;
/*allocate memory for the next step's potential and store the old
* potential, such that one lattice site of the previous step is stored in
* F*F lattice sites at the appropriate coordinates.
* Example: F=2: x-> xx
* xx*/
phi=malloc(Xstep*Ystep*sizeof(double));
for(x=0, ix=0; x<Xstep; x++){
for(y=0, iy=0; y<Ystep; y++){
phi[y*Xstep+x]=phibuf[iy*(int)Xstep/F+ix];
if((y+1)%F==0)
iy++;
}
if((x+1)%F==0)
ix++;
}
free(phibuf);/*free memory used by the temporary potential*/
/*
* as the lattice size is increased the predefined potentials and grounded
* areas must grow by the same factor:
* */
for(i=0; i<NPHI; i++){
for(x=0; x<4; x++){
phi0[i].pos[x]*=F;
}
}
for(i=0; i<NGROUND; i++){
for(x=0; x<4; x++){
ground[i][x]*=F;
}
}
step++;
}/*end loop over lattice-volumes*/
return 0;
}