nbody_1000.cpp

MPI编程－Nbody问题并行程序（计算1000个天体的移动

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fstream.h>
#include <iostream.h>

/*-------------------构造点的类型-------------------*/
struct Point{
        double x;
        double y;
};

/*--------------------全局变量----------------------*/
double G;               //万有引力常量  
int N;                  //天体的数目
Point* P;               //每个天体的位置
Point* v;               //每个天体的速度
Point* f;               //每个天体受到的力
double* m;              //每个天体的质量
Point size;             //空间的大小
float C;                //一常数
double dt;              //很小的时间周期

/*---------------初始化各天体的情况-----------------*/
void Initial(char* filename){
	ifstream f1(filename);
	if(!f1){
		cout<<"open file fail!"<<endl;
		return;
	}
	f1>>N;         //从文件读入天体的数目      
	P = new Point[N];
	v = new Point[N];
	f = new Point[N];
	m = new double[N];
	f1>>G>>size.x>>size.y>>dt>>C;      //从文件读入G，size，dt
	for(int i=0;i<N;i++){
		f1>>m[i]>>P[i].x>>P[i].y>>v[i].x>>v[i].y;     //从文件读入各天体的质量，位置和速度
	}
	for(int j=0;j<N;j++) { f[j].x = 0.0; f[j].y = 0.0; }   //初始化每个天体的受力为0
	f1.close();
}

/*-----------计算dt时间内每个天体受到的力-----------*/
void calculateForces(int i){
	double distance,magnitude;
	Point direction;
	for(int j=0;j<N;j++){
		if(j!=i){
			distance = (float)sqrt((P[i].x-P[j].x)*(P[i].x-P[j].x)+(P[i].y-P[j].y)*(P[i].y-P[j].y));    //计算两点间距离
			if(distance!=0){
				magnitude = (G*m[i]*m[j])/(distance*distance+C*C);   //计算万有引力
				direction.x = P[j].x-P[i].x;                         //记录力的x方向
				direction.y = P[j].y-P[i].y;                         //记录力的y方向
				f[i].x = f[i].x+magnitude*direction.x/distance;      //把收到的所有万有引力累加
				f[i].y = f[i].y+magnitude*direction.y/distance;
			}
		}
	}
}

/*----------计算dt时间后，天体的位置和速度----------*/
void moveBodies(int i){
	Point dv;
	Point dp;
	Point tempPoint;
	dv.x = f[i].x/m[i]*dt;       dv.y = f[i].y/m[i]*dt;         //计算x，y方向的dv＝a*dt,a=F/m
	dp.x = (v[i].x+dv.x/2)*dt;   dp.y = (v[i].y+dv.y/2)*dt;     //计算x，y偏移量dp＝(v+dv/2)dt
	tempPoint.x=P[i].x+dp.x;     tempPoint.y = P[i].y+dp.y;     //对位置的偏移量进行累加
	if(tempPoint.x<0||tempPoint.x>size.x||tempPoint.y<0||tempPoint.y>size.y){  //判断下一位置是否越界
		P[i].x = P[i].x;     //如果越界，天体保留原来位置，速度方向取反
		P[i].y = P[i].y;
		v[i].x = -v[i].x;
		v[i].y = -v[i].y;
	}else{                   //记录新的位置和方向
		P[i].x = tempPoint.x;
		P[i].y = tempPoint.y;
		v[i].x = v[i].x + dv.x;
		v[i].y = v[i].y + dv.y;
	}
	f[i].x = 0.0; f[i].y = 0.0;
}

/*-------------------输出各天体的情况---------------*/
void End(char* filename){
	ofstream f2(filename);
	if(!f2){
		cout<<"open file fail!"<<endl;
		return;
	}
	f2<<N<<endl<<G<<endl<<size.x<<" "<<size.y<<endl<<dt<<endl<<C<<endl;        //输出n，G，size，dt到结果文件
	for(int i=0;i<N;i++){
		f2<<m[i]<<" "<<P[i].x<<" "<<P[i].y<<" "<<v[i].x<<" "<<v[i].y<<endl;    //输出各天体的新情况到结果文件
	}
	delete[] P;
	delete[] v;
	delete[] f;
	delete[] m;
	f2.close();
}




/*--------------------------主函数------------------*/
int main(int argc,char* argv[]){
	int myid,numprocs;
    double starttime,endtime;               //记录运行的开始时间和结束时间

    MPI_Init(&argc,&argv);                  //启动MPI计算   
    MPI_Comm_size(MPI_COMM_WORLD,&numprocs);//确定进程数
    MPI_Comm_rank(MPI_COMM_WORLD,&myid);    //确定自己的进程标识
    starttime=MPI_Wtime();                  //记录当前时间
    Initial("sample_input.in");             //读取数据
	double* PP = new double[2*N];           //记录天体的x，y坐标
	double* vv = new double[2*N];           //记录天体的x，y方向
	int d=N/numprocs;                       //记录每个进程分配到的天体的数目
	for(int time=0; time<1000;time+=1){     //进行循环计算1000次
		for(int i=myid*d;i<(myid+1)*d;i++){ //每个进程计算d个天体
			calculateForces(i);             //计算第i个天体的力
			moveBodies(i);                  //计算第i个天体新的位置和力
			PP[2*i]=P[i].x;                 //记录天体的坐标，方向
			PP[2*i+1]=P[i].y;
			vv[2*i]=v[i].x;
			vv[2*i+1]=v[i].y;
		}
		if(myid==(numprocs-1)){             //如果是最后一个进程，把剩下的没分配的天体进行计算
			for(int j=(myid+1)*d;j<N;j++){
				calculateForces(j);
				moveBodies(j);
				PP[2*j]=P[j].x;
				PP[2*j+1]=P[j].y;
				vv[2*j]=v[j].x;
				vv[2*j+1]=v[j].y;
			}
		}
		MPI_Barrier(MPI_COMM_WORLD);        //进行各进程同步
		for(int j=0;j<numprocs;j++){        //各进程依次广播自己计算的位置
			MPI_Barrier(MPI_COMM_WORLD);    
			if (j!=(numprocs-1))	MPI_Bcast(&PP[j*2*d],2*d,MPI_DOUBLE,j,MPI_COMM_WORLD);    
			else MPI_Bcast(&PP[j*2*d],(N-j*d)*2,MPI_DOUBLE,j,MPI_COMM_WORLD);
			MPI_Barrier(MPI_COMM_WORLD);    //广播完各进程同步
			for(int k=j*d;k<(j+1)*d;k++){   //收到广播信息后，其他进程更新自己的值
				P[k].x=PP[2*k];
				P[k].y=PP[2*k+1];
			}
			MPI_Barrier(MPI_COMM_WORLD);    //更新完成后同步
		}
	}
	MPI_Barrier(MPI_COMM_WORLD);            //计算完1000后同步
	for(int j=0;j<numprocs;j++){            //各进程依次广播自己的速度
		MPI_Barrier(MPI_COMM_WORLD);        //发送速度前同步
		if (j!=(numprocs-1))	MPI_Bcast(&vv[j*2*d],2*d,MPI_DOUBLE,j,MPI_COMM_WORLD);   
		else MPI_Bcast(&vv[j*2*d],(N-j*d)*2,MPI_DOUBLE,j,MPI_COMM_WORLD);
		MPI_Barrier(MPI_COMM_WORLD);        //广播完同步
		for(int k=j*d;k<(j+1)*d;k++){       //收到广播信息后，其他进程更新自己的值
			v[k].x=vv[2*k];
			v[k].y=vv[2*k+1];
		}
	}
	MPI_Barrier(MPI_COMM_WORLD);            //计算完成后同步
	if(myid==0)     End("result1000.data"); //如果是0进程，输出运行结果到文件中
	endtime=MPI_Wtime();                    //计算每个进程整个程序的运行时间
	if(myid==0) cout<<"总的运行时间："<<endtime-starttime<<"s"<<endl;     //输出每个进程的运行时间 
	delete[] PP;
	delete[] vv; 
	MPI_Finalize();                         //结束MPI运算 
	return 0;
}