Решение задачи Пуассона, страница 2

Eps=0.1;

LoadFirst(Phi, F, BlockSize, cntY, cntZ, Rank, Hx, Hy, Hz);

LoadInitApprox(Phi, BlockSize, cntY, cntZ, Rank, PrCount, FullBlocks);

Flag=0;

Iter=0;

A=0;

temp=(2/(Hx*Hx)+2/(Hy*Hy)+2/(Hz*Hz)+A);

if (Rank==PrCount-1)

{

IniFile= fopen("process.log", "w");

} 

Phi[Coord(BlockSize-1,cntY-1,cntZ)]=0;

MPI_Barrier(MPI_COMM_WORLD);

start=time(0);

do //while (!Flag)

{

Iter++;

CurEps=0;

Layer1[Coord(0,cntY-1,cntZ)]=0;

if (Rank!=0)

{

MPI_Isend(&Phi[Coord(0,0,0)],cntY*cntZ+1,MPI_FLOAT,prev_pr(Rank,PrCount),TheTag,MPI_COMM_WORLD, &rec1);

MPI_Recv(&Layer1[Coord(0,0,0)],cntY*cntZ+1,MPI_FLOAT,prev_pr(Rank,PrCount),TheTag,MPI_COMM_WORLD,&st);

CurEps=Layer1[Coord(0,cntY-1,cntZ)];

for (j=1; j<cntY-1; j++)

for (k=1; k<cntZ-1; k++)

{

temp2=Phi[Coord(0,j,k)];

Phi[Coord(0,j,k)]=(Layer1[Coord(0,j,k)]+Phi[Coord(1,j,k)])/(Hx*Hx);

Phi[Coord(0,j,k)]+=(Phi[Coord(0,j-1,k)]+Phi[Coord(0,j+1,k)])/(Hy*Hy);

Phi[Coord(0,j,k)]+=(Phi[Coord(0,j,k-1)]+Phi[Coord(0,j,k+1)])/(Hz*Hz);

Phi[Coord(0,j,k)]-=F[Coord(0,j,k)];

Phi[Coord(0,j,k)]/=temp;

CurEps+=fabsf(Phi[Coord(0,j,k)]-temp2);

}

}

for (i=1; i<BlockSize-1; i++)

for (j=1; j<cntY-1; j++)

for (k=1; k<cntZ-1; k++)

{

temp2=Phi[Coord(i,j,k)];

Phi[Coord(i,j,k)]=(Phi[Coord(i-1,j,k)]+Phi[Coord(i+1,j,k)])/(Hx*Hx);

Phi[Coord(i,j,k)]+=(Phi[Coord(i,j-1,k)]+Phi[Coord(i,j+1,k)])/(Hy*Hy);

Phi[Coord(i,j,k)]+=(Phi[Coord(i,j,k-1)]+Phi[Coord(i,j,k+1)])/(Hz*Hz);

Phi[Coord(i,j,k)]-=F[Coord(i,j,k)];

Phi[Coord(i,j,k)]/=temp;

CurEps+=fabsf(Phi[Coord(i,j,k)]-temp2);

}

if ((Rank==PrCount-1)&&(CurEps<Eps))

{

tempPhi=Phi[Coord(0,cntY-1,cntZ)];

Phi[Coord(0,cntY-1,cntZ)]=-1;

MPI_Send(&Phi[Coord(0,0,0)],cntY*cntZ+1,MPI_FLOAT,prev_pr(Rank,PrCount),TheTag,MPI_COMM_WORLD);

Phi[Coord(0,cntY-1,cntZ)]=tempPhi;

Flag=1;

}

if (Rank!=PrCount-1)

{

MPI_Recv(&Layer1[Coord(0,0,0)],cntY*cntZ+1,MPI_FLOAT,next_pr(Rank,PrCount),TheTag,MPI_COMM_WORLD,&st);

if (Layer1[Coord(0,cntY-1,cntZ)]<0)

{

tempPhi=Phi[Coord(0,cntY-1,cntZ)];

Phi[Coord(0,cntY-1,cntZ)]=-1;

if (Rank!=0)

MPI_Send(&Phi[Coord(0,0,0)],cntY*cntZ+1,MPI_FLOAT,prev_pr(Rank,PrCount),TheTag,MPI_COMM_WORLD);

Phi[Coord(0,cntY-1,cntZ)]=tempPhi; 

Flag=1;

}

for (j=1; j<cntY-1; j++)

for (k=1; k<cntZ-1; k++)

{

temp2=Phi[Coord(BlockSize-1,j,k)];

Phi[Coord(BlockSize-1,j,k)]=(Layer1[Coord(0,j,k)]+Phi[Coord(BlockSize-2,j,k)])/(Hx*Hx);

Phi[Coord(BlockSize-1,j,k)]+=(Phi[Coord(BlockSize-1,j-1,k)]+Phi[Coord(BlockSize-1,j+1,k)])/(Hy*Hy);

Phi[Coord(BlockSize-1,j,k)]+=(Phi[Coord(BlockSize-1,j,k-1)]+Phi[Coord(BlockSize-1,j,k+1)])/(Hz*Hz);

Phi[Coord(BlockSize-1,j,k)]-=F[Coord(BlockSize-1,j,k)];

Phi[Coord(BlockSize-1,j,k)]/=temp;

CurEps+=fabsf(Phi[Coord(BlockSize-1,j,k)]-temp2);

}

Phi[Coord(BlockSize-1,cntY-1,cntZ)]=CurEps; 

if (Flag==0)

MPI_Isend(&Phi[Coord(BlockSize-1,0,0)],cntY*cntZ+1,MPI_FLOAT,next_pr(Rank,PrCount),TheTag,MPI_COMM_WORLD, &rec2); 

}

if (Iter==100000) Flag=1;

if (Rank==PrCount-1) fprintf(IniFile, "Last process: Iteration %i, CurEps=%f\n", Iter, CurEps);   

} while(!Flag);

MPI_Barrier(MPI_COMM_WORLD);

start=time(0)-start;

for (h=0; h<PrCount; h++)

{

MPI_Barrier(MPI_COMM_WORLD);

if (h==Rank)

{

if (h!=0)

ResFile = fopen("result.txt", "a");

else ResFile = fopen("result.txt", "w");

for (i=0; i<BlockSize; i++)

for (j=0; j<cntY; j++)

for (k=0; k<cntZ; k++)

fprintf(ResFile, "Phi[%i,%i,%i]=%f\n", i+Rank*BlockSize,j,k,Phi[Coord(i,j,k)]);

fclose(ResFile);     

}     

} 

MPI_Barrier(MPI_COMM_WORLD); 

printf("process %i make %i iterations\n ", Rank, Iter);   

if (Rank==PrCount-1)

{

fprintf(IniFile,"Calculation time: %d", start);

fclose(IniFile);

} 

MPI_Finalize();

return 0;          

}

int next_pr(int cur, int size)

{

return (cur+size+1)%size;

}

int prev_pr(int cur, int size)

{

return (cur+size-1)%size;

}

int my_init(float *Hx, float *Hy, float *Hz, int *cX, int *cY, int *cZ)

{

FILE *IniFile;

float dX, dY, dZ, tmp1,tmp2,tmp3;

IniFile= fopen("init.dat", "r");

fscanf(IniFile, "%f %f %f", &dX, &dY, &dZ);

fscanf(IniFile, "%i %i %i", cX, cY, cZ);

tmp1=*cX; tmp2=*cY; tmp3=*cZ;

*Hx=dX/(tmp1-1);

*Hy=dY/(tmp2-1);

*Hz=dZ/(tmp3-1);

fclose(IniFile);

return 0;

}

int Coord(int x, int y, int z)

{

return (x*(MaxY*MaxZ)+y*MaxZ+z);

}

int LoadFirst(float *Phi, float *F, int BlockSize, int cntY, int cntZ, int Rank, float Hx, float Hy, float Hz)

{

int i,j,k;

float CurX, CurY, CurZ;

float StartX;

StartX=(Rank*BlockSize)*Hx;

for (i=0; i<BlockSize; i++)

for (j=0; j<cntY; j++)

for (k=0; k<cntZ; k++)

{

CurX=StartX+i*Hx; CurY=j*Hy; CurZ=k*Hz;

Phi[Coord(i,j,k)]=CurX*CurX*CurX+CurY*CurY*CurY+CurZ*CurZ*CurZ;

F[Coord(i,j,k)]=6*CurX+6*CurY+6*CurZ;

}

return 0;

}

int LoadInitApprox(float *Phi, int BlockSize, int cntY,int  cntZ,int  Rank,int PrCount,int  FullBlocks)

{

int i,j,k;

float StI, FinI;

float StartX;

if (Rank==0) StI=1; else StI=0;

if (Rank==PrCount-1) FinI=BlockSize-1; else FinI=BlockSize;

for (i=StI; i<FinI; i++)

for (j=1; j<cntY-1; j++)

for (k=1; k<cntZ-1; k++)

{

Phi[Coord(i,j,k)]=1;

}

return 0;  

}