spiral.cpp

对二维的反应扩散方程使作欧拉差分的并行计算

#define MPICH_SKIP_MPICXX
#define m 100
#define n 100
#include <stdio.h>
#include <mpi.h>
#include <malloc.h>


float f(float,float,float,float,float);
float ff(float);

void master(void);
void slave(void);

int main(int argc,char **argv)
{
	int rank, size;
	
	MPI_Init( &argc, &argv );
	MPI_Comm_rank( MPI_COMM_WORLD, &rank );
	MPI_Comm_size( MPI_COMM_WORLD,&size);
	if(size  != 5) MPI_Abort(MPI_COMM_WORLD,1);
	if (rank == 0)
		master(); /*进程0作为主进程*/
	else
		slave(); /*其它进程作为从进程*/
	MPI_Finalize( );
	return 0;
}

void master()
{
	int i,j,size;
	float u[m][n]; //, v[m][n];
	float buffer[m][n];
	MPI_Comm_size(MPI_COMM_WORLD,&size);
	int *counts,*disp;

	disp = (int *)malloc(size*sizeof(int));
	counts = (int *)malloc(size*sizeof(int));
	//int counts[5],disp[5];
	FILE *fp1;
	//FILE *fp2;
	fp1=fopen("u.dat","w");
	//fp2=fopen("v.dat","w");
	//MPI_Comm_rank( MPI_COMM_WORLD, &rank );

	float sendarray=0;
	for(i=0;i<m;i++)
		for(j=0;j<m;j++)
			u[i][j] = 0;
	
	disp[0]=0;
	counts[0]=1;
	for(i=1;i<5;i++)
	{
		disp[i]=(i-1)*2500;
		counts[i]=2500;
	}
	MPI_Gatherv(&sendarray,1,MPI_FLOAT,buffer,counts,disp,MPI_FLOAT,0,MPI_COMM_WORLD);

	for(i=0;i<m;i++)
	{
		for(j=0;j<n;j++)
		{
			fprintf(fp1,"%12.9lf ",buffer[i][j]);
			//fprintf(fp2,"%12.9lf ",v[i][j]);
		}
		fprintf(fp1,"\n");
		//fprintf(fp2,"\n");
	}
	fclose(fp1);
	//fclose(fp2);
}

void slave()
{
	FILE *fpp;
	fpp=fopen("test.dat","aw");
	int rbuf;
	int x,y;
	int rank,size;
	int  i, j, ss;
	int i_first, i_last;
	float a = 0.84f,b=0.07f,eps=0.04f,tao=0.02f,h=0.5f;
	float init=0,fint=14.0,t;
    MPI_Status status;
	MPI_Comm_rank( MPI_COMM_WORLD, &rank );
	MPI_Comm_size( MPI_COMM_WORLD,&size);
	float u[(m/4)+2][n],v[(m/4)+2][n];
	float uu[(m/4)+2][n],vv[(m/4)+2][n],pxx[(m/4)+2][n];
	float send[25][100];

	i_first=1;
	i_last=m/4;
	if(rank==1) i_first++;
	if(rank==4) i_last--;
	//initial values
	for(i=0;i<(m/4)+2;i++)
		for(j=0;j<n;j++)
		{
			u[i][j]=0;
			v[i][j]=0;
			uu[i][j]=0;
			vv[i][j]=0;
			pxx[i][j]=0;
		}

	if(rank==2)
	{
		//row 20,21
		for(i=m/4-5;i<=m/4-4;i++)
			for(j=0;j<n/2;j++)
				u[i][j]=0.9f;
		//row 22,23
		for(i=m/4-3;i<=m/4-2;i++)
			for(j=0;j<n/2;j++)
			{
				u[i][j]=0.7f;
				v[i][j]=0.7f;
			}
		//row 24,25
		for(i=m/4-1;i<=m/4;i++)
			for(j=0;j<n/2;j++)
				v[i][j]=0.9f;
	}
	t=init;
	ss = -1;
	for(;;)
	{
		//send down unless I'm at the bottom,
		//then receive from up
		if(rank<size-1)
		{
			MPI_Sendrecv(u[m/4],n,MPI_FLOAT,rank+1,1001,\
				u[m/4+1],n,MPI_FLOAT,rank+1,1003,MPI_COMM_WORLD,&status);
			MPI_Sendrecv(v[m/4],n,MPI_FLOAT,rank+1,1002, \
				v[m/4+1],n,MPI_FLOAT,rank+1,1004,MPI_COMM_WORLD, &status);
		}
		//send up unless I'm at the top
		if(rank>1)
		{
			MPI_Sendrecv(u[1],n,MPI_FLOAT,rank-1,1003, \
				u[0],n,MPI_FLOAT,rank-1,1001,MPI_COMM_WORLD, &status);
			MPI_Sendrecv(v[1],n,MPI_FLOAT,rank-1,1004, \
				v[0],n,MPI_FLOAT,rank-1,1002,MPI_COMM_WORLD, &status);
		}
		
		//calculate ppx
		//center
		for(i=i_first;i<=i_last;i++)
			for(j=1;j<n-1;j++)
				pxx[i][j]=u[i+1][j]+u[i-1][j]+u[i][j+1]+u[i][j-1]-4*u[i][j];
		//left and right boundary
		for(i=i_first;i<=i_last;i++)
		{
			pxx[i][0]=u[i+1][0]+u[i-1][0]+2*u[i][1]-4*u[i][0];
			pxx[i][n-1]=u[i+1][n-1]+u[i-1][n-1]+2*u[i][n-2]-4*u[i][n-1];
		}
		//up and down boundary
		if(rank==1)
		{
			for(j=1;j<n-1;j++)
				pxx[1][j]=u[1][j+1]+u[1][j-1]+2*u[2][j]-4*u[1][j];
			pxx[1][0]=2*u[2][0]+2*u[1][1]-4*u[1][0];
			pxx[1][n-1]=2*u[2][n-1]+2*u[1][n-2]-4*u[1][n-1];
		}
		if(rank==4)
		{
			for(j=1;j<n-1;j++)
				pxx[m/4][j]=u[m/4][j+1]+u[m/4][j-1]+2*u[m/4-1][j]-4*u[m/4][j];
			pxx[m/4][0]=2*u[m/4-1][0]+2*u[m/4][1]-4*u[m/4][0];
			pxx[m/4][n-1]=2*u[m/4-1][n-1]+2*u[m/4][n-2]-4*u[m/4][n-1];
		}

		//calcuate the function
		for(i=1;i<=m/4;i++)
		{
			for(j=0;j<n;j++)
			{

				uu[i][j]=u[i][j]+tao*(f(a,b,eps,u[i][j],\
					     v[i][j])+pxx[i][j]/(h*h));
				vv[i][j]=v[i][j]+tao*(ff(u[i][j])-v[i][j]);

			}			
		}
		//count the cycle
		ss=ss+1;
		if(ss>=500)
		{
			ss=0;
			printf("%f",t);
		}
		t=t+tao;
		if(t>=fint)
		{			
			for(i=1;i<=m/4;i++)
				for(j=0;j<n;j++)
				{	
					send[i-1][j]=u[i][j];
				//	buffer2[j]=v[i][j];
				}
				fprintf(fpp,"process:%d\n",rank);
			for(i=0;i<m/4;i++)
			{	for(j=0;j<4;j++)
					fprintf(fpp,"%12.9f",send[i][j]);
				fprintf(fpp,"\n");
			}


			MPI_Gatherv(send,2500,MPI_FLOAT,&rbuf,&x,&y,MPI_FLOAT,0,MPI_COMM_WORLD);
			break;
		}
		for(i=1;i<=m/4;i++)
			for(j=0;j<n;j++)
			{
				u[i][j]=uu[i][j];
				v[i][j]=vv[i][j];
			}
	}
	fclose(fpp);
}

float f(float aa,float bb,float epsl,float u, float v)
{
	return u*(1-u)*(u-(v+bb)/aa)/epsl;
}

float ff(float u)
{
	float temp;
	if(u < 1/3.0)
		temp = 0.0;
	if(u > 1)
		temp = 1.0;
	if(u >= 1/3.0 && u <= 1.0) 
		temp =1-6.75*u*(u-1)*(u-1);
	return temp;
}