并行源程序.txt

计算并行运算的代码,用mpI实现并行运算,

/*并行源程序*/
#include <mpi.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <conio.h>

#define pi  3.1415926
#define c   2.998e-1
#define E1  8.85e-3
#define E2  2.58*E1
#define dt  0.1374642e-2
#define f0  10.5
#define dr  0.713809e-3
#define t0  0.15
#define t   0.05
#define df  1/(218*dt)
#define u   1.25e+3
#define n0  0
#define N   258
#define rn  218
#define dc  pi/90
#define ld  c/9.5
#define ita 120*pi
#define z0  35*dr
#define s   0.577
#define X1  80
#define Y   100
#define Z1  75
#define X2  20
#define Z2  100

float ex1[X1+1][Y+1][Z1+1],ey1[X1+1][Y+1][Z1+1],ez1[X1+1][Y+1][Z1+1];
float hx1[X1+1][Y+1][Z1+1],hy1[X1+1][Y+1][Z1+1],hz1[X1+1][Y+1][Z1+1];
float ex2[X2+1][Y+1][Z2+1],ey2[X2+1][Y+1][Z2+1],ez2[X2+1][Y+1][Z2+1];
float hx2[X2+1][Y+1][Z2+1],hy2[X2+1][Y+1][Z2+1],hz2[X2+1][Y+1][Z2+1];

float S_hx1[X2+1][Y+1],S_hy1[X2+1][Y+1],S_hz1[X2+1][Y+1];
float R_hx2[X2+1][Y+1],R_hy2[X2+1][Y+1],R_hz2[X2+1][Y+1];
float R_ex1[X2+1][Y+1],R_ey1[X2+1][Y+1],R_ez1[X2+1][Y+1];
float S_ex2[X2+1][Y+1],S_ey2[X2+1][Y+1],S_ez2[X2+1][Y+1];

float S_MY00_ex2[X2+1][Z2+1],S_MY00_ez2[X2+1][Z2+1];
float S_MY01_ex2[X2+1][Z2+1],S_MY01_ez2[X2+1][Z2+1];
float S_MYY0_ex2[X2+1][Z2+1],S_MYY0_ez2[X2+1][Z2+1];
float S_MYY1_ex2[X2+1][Z2+1],S_MYY1_ez2[X2+1][Z2+1];
float R_MY00_ex2[X2+1][Z2+1],R_MY00_ez2[X2+1][Z2+1];
float R_MY01_ex2[X2+1][Z2+1],R_MY01_ez2[X2+1][Z2+1];
float R_MYY0_ex2[X2+1][Z2+1],R_MYY0_ez2[X2+1][Z2+1];
float R_MYY1_ex2[X2+1][Z2+1],R_MYY1_ez2[X2+1][Z2+1];

float S_MZ00_ex2[X2+1][Y+1],S_MZ00_ey2[X2+1][Y+1];
float S_MZ01_ex2[X2+1][Y+1],S_MZ01_ey2[X2+1][Y+1];
float R_MZ00_ex2[X2+1][Y+1],R_MZ00_ey2[X2+1][Y+1];
float R_MZ01_ex2[X2+1][Y+1],R_MZ01_ey2[X2+1][Y+1];

float S_LX00_ex2[X2-1],S_LX00_ey2[X2-1],S_LX00_ez2[X2-1];
float S_LX01_ex2[X2-1],S_LX01_ey2[X2-1],S_LX01_ez2[X2-1];
float R_LX00_ex2[X2-1],R_LX00_ey2[X2-1],R_LX00_ez2[X2-1];
float R_LX01_ex2[X2-1],R_LX01_ey2[X2-1],R_LX01_ez2[X2-1];

float S_MY20_ex2[X2-1][Z2],S_MY20_ez2[X2-1][Z2];
float S_MY30_ex2[X2-1][Z2],S_MY30_ez2[X2-1][Z2];
float R_MY20_ex2[X2-1][Z2],R_MY20_ez2[X2-1][Z2];
float R_MY30_ex2[X2-1][Z2],R_MY30_ez2[X2-1][Z2];

float S_MZ20_ex2[X2-1][Y-1],S_MZ20_ey2[X2-1][Y-1];
float R_MZ20_ex2[X2-1][Y-1],R_MZ20_ey2[X2-1][Y-1];

float S_LX20_ex2[X2-1],S_LX20_ey2[X2-1],S_LX20_ez2[X2-1];
float S_LX21_ex2[X2-1],S_LX21_ey2[X2-1],S_LX21_ez2[X2-1];
float R_LX20_ex2[X2-1],R_LX20_ey2[X2-1],R_LX20_ez2[X2-1];
float R_LX21_ex2[X2-1],R_LX21_ey2[X2-1],R_LX21_ez2[X2-1];

float xey[Y+1][Z1+1][10],xez[Y+1][Z1+1][10];
float yexb[X1+1][Z1+1][10],yezb[X1+1][Z1+1][10];
float yexs[X2+1][Z2+2][10],yezs[X2+1][Z2+2][10];
float zex[X1+3][Y+1][10],zey[X1+3][Y+1][10];
float yey1[Y+1][30][14],yez1[Y+1][30][14];
float yex2[Y+1][30][14],yey2[Y+1][30][14];
float zey1[Z1+1][30][14],zez1[Z1+1][30][14];
float zex2[Z1+1][30][14],zez2[Z1+1][30][14];
float xex1[X1+1][30][14],xey1[X1+1][30][14];
float xex2[X1+1][30][14],xez2[X1+1][30][14];
float A,B,C;
float tt0[3][16],tt[8][16];
float rex[X2+1][Y+1][rn+1],rey[X2+1][Y+1][rn+1];
float rhx[X2+1][Y+1][rn+1],rhy[X2+1][Y+1][rn+1];
float fex[X2+1][Y+1][2],fey[X2+1][Y+1][2];
float fhx[X2+1][Y+1][2],fhy[X2+1][Y+1][2];
float Er1[X2+1][Y+1][91],Ei1[X2+1][Y+1][91];
float Er2[X2+1][Y+1][91],Ei2[X2+1][Y+1][91];
float Er3[X2+1][Y+1][91],Ei3[X2+1][Y+1][91];
float Er4[X2+1][Y+1][91],Ei4[X2+1][Y+1][91];
float mex[X2+1][Y+1],mey[X2+1][Y+1];
float mhx[X2+1][Y+1],mhy[X2+1][Y+1];
float se1[91][2],se2[91][2];
float Me1[91],Me2[91];
float se3[91][2],se4[91][2];
float Me3[91],Me4[91];
float Vs[rn+1],Is[rn+1];
float Vsw[2],Isw[2],pin[2],Mp;
float G1[91],G2[91],G3[91],G4[91];
int c7[8]={1,7,21,35,35,21,7,1};

double starttime0,endtime0,starttime,endtime;
double starttime1,endtime1;
MPI_Status status;
int my_rank,process_num;


void main(int argc, char **argv)
{ 
  FILE  *fp14,*fp15,*fp16,*fp17,*fp18;
 int m,i,j,k,n;
  float y;

  void bigarea_h();
  void smallarea_h();
  void bigarea_e();
  void smallarea_e();

  void Send_h1();
  void Recv_h2();
  void Recv_e1();
  void Send_e2();
  
  void Send_MY();
  void Send_MZ();
  void Send_LX();
  void Recv_MY();
  void Recv_MZ();
  void Recv_LX();
  
  void Send_MY1();
  void Send_MZ1();
  void Send_LX1();
  void Recv_MY1();
  void Recv_MZ1();
  void Recv_LX1();
  
  void init_liao();
  void murx(int g);
  void mury1(int g);
  void mury2(int g);
  void murz(int g,int xs,int xf,int yf);
  void liaoy(int xx,int zz);
  void liaoz(int xx,int yy);
  void liaox(int xs,int xf,int yy,int zz);
  void fourier(float e[X2+1][Y+1][rn+1]);
  void fourier1(float e[rn+1]);
  void mo(float e[X2+1][Y+1][2]);
  void sss();
  void sss1();
  void putdata();

  printf("hello\n");
  MPI_Init(&argc,&argv);
  MPI_Comm_size(MPI_COMM_WORLD, &process_num);
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
  starttime=MPI_Wtime(); 

  if(my_rank==0) {
    init_liao();          
  }
  
  for(n=1;n<=N;n++)
  { 
   if(n==1||n==50||n==100||n==150||n==200||n==N) {
	 starttime0=MPI_Wtime(); 
	  }
	  y=(((n*dt-t0)/t)*((n*dt-t0)/t));
    ex1[30][50][13]=1000*exp((-1.0*y))*sin(2*pi*f0*(n-n0)*dt);
    
   if(my_rank==0){
    bigarea_h();

    Send_h1();
  }
  if(my_rank==1){
  smallarea_h();
  Recv_h2();
}
   if(my_rank==0){
   if(n<=rn)
    {  Vs[n]=ex1[30][50][13];
       Is[n]=fabs(hz1[30][51][13])+fabs(hz1[30][50][13])+fabs(hy1[30][50][14])+fabs(hy1[30][50][13]); 
      }
   bigarea_e();
  }
  if(my_rank==1){
   smallarea_e();
	}
	 
  
  	//接收1进程传过来的用于计算mury2,murz,liaox的数据
  	if(my_rank==0){
   Recv_MY();
   Recv_MZ();  
   Recv_LX(); 
  }
 if(my_rank==0){
    murx(0); 
    murx(X1);
    mury1(0);
    mury2(0);
    mury2(Y);
	mury1(Y);

    murz(Z1,1,X1-1,Y-1);
    murz(0,1,29,Y-1);
    murz(0,51,X1-1,Y-1);
    murz(Z2,1,X2-1,Y-1);
    liaoy(0,0);
    liaoy(0,Z1);
    liaoy(X1,0);
    liaoy(X1,Z1);
    liaoz(0,0);
    liaoz(0,Y);
    liaoz(X1,0);
    liaoz(X1,Y);
    liaox(0,X1,0,Z1);
    liaox(0,X1,Y,Z1);
    liaox(0,29,0,0);
    liaox(0,29,Y,0);
    liaox(51,X1,0,0);
    liaox(51,X1,Y,0);
    liaox(1,X2-1,0,Z2);
    liaox(1,X2-1,Y,Z2);
 
    Send_MY1(); 
 	  Send_MZ1(); 
    Send_LX1();
 
	  Recv_e1(); 
 
 
} 
/////////////////////////////////////////////////////////////////////////////////////////
   if(my_rank==1) {    

//将数据传给0进程计算mury2,murz,liaox

 	Send_MY();
 	Send_MZ();
	Send_LX();

 //接收0进程mury2传过来的数据

 	Recv_MY1();
 	Recv_MZ1();
 	Recv_LX1();

  Send_e2();

  } 
 if ((n>=36)&&(n<=36+rn) )
 	  {m=n-36;
      for (i=30;i<=50;i++)
 	    for (j=0;j<=Y;j++)
        { rex[i-30][j][m]=ex1[i][j][35];
          rey[i-30][j][m]=ey1[i][j][35];
          rhx[i-30][j][m]=hx1[i][j][35];
          rhy[i-30][j][m]=hy1[i][j][35]; 
              }  
     }
   if(n==1||n==50||n==100||n==150||n==200||n==N){

	   endtime0=MPI_Wtime(); 
	   printf("rank%d:the %d time running time is:%f\n",my_rank,n,endtime0-starttime0);
   }

MPI_Barrier(MPI_COMM_WORLD);
printf("end[%d]\n",n) ;
 }

 if(my_rank==0) {          
 fourier(rex);
 fourier(rey);
 fourier(rhx);
 fourier(rhy);
 mo(fex);
 mo(fey);
 mo(fhx);
 mo(fhy);
 sss();
 sss1();
 fourier1(Vs);
 fourier1(Is);
 
 putdata();
}
endtime=MPI_Wtime();
printf("the whole running time is %f\n",endtime-starttime); 
MPI_Finalize();
printf("rank%d:the %d time running time is:%f\n",my_rank,n,endtime0-starttime0);
}




//////////////////////////////////////////////////////////////////////////////
//初始化数组tt,为liao边界条件作准备
 void init_liao() {
 	int l,m0,m1,m2,m3,i,j,j1;
  float tmp;
  tt[1][1]=(2-s)*(1-s)/2;
  tt[1][2]=s*(2-s);
  tt[1][3]=s*(s-1)/2;
  for (l=2;l<=7;l++)
   { m2=2*l+1;
     m1=l-1;
     for (j=1;j<=m2;j++)
      {if (j<(2*l-1)) tt0[1][j]=tt[m1][j];
       else tt0[1][j]=0; }
     m3=2*l;
     tt0[2][1]=0;
     for (j=2;j<=m3;j++)
      { j1=j-1;
        tt0[2][j]=tt[m1][j1]; }
     tt0[2][m2]=0;
     tt0[3][1]=0;
     tt0[3][2]=0;
     for (j=3;j<=m2;j++)
      { j1=j-2;
        tt0[3][j]=tt[m1][j1]; }
     m0=2*l+1;
     for (i=1;i<=m0;i++)
      { tt[l][i]=0;
        for (j=1;j<=3;j++)
           { tmp=tt[1][j]*tt0[j][i];
             tt[l][i]=tt[l][i]+tmp; }    }   }
   A=(c*dt-dr)/(c*dt+dr);
   B=(2*dr)/(c*dt+dr);
   C=pow((c*dt),2)/(2*dr*(c*dt+dr));
    }
///////////////////////////////////////////////////////////////////////////    
//计算大区域磁场
 void bigarea_h(){
  int i,j,k,i1,j1,k1;
  float x1,y1,z1;
    for(i=0;i<=X1-1;i++)
     for(j=0;j<=Y-1;j++)
      for(k=0;k<=Z1-1;k++)
       {i1=i+1;
        j1=j+1;
        k1=k+1;
        x1=ey1[i][j][k1]-ey1[i][j][k]+ez1[i][j][k]-ez1[i][j1][k];
        hx1[i][j][k]=hx1[i][j][k]+dt/(u*dr)*x1;
        y1=ez1[i1][j][k]-ez1[i][j][k]+ex1[i][j][k]-ex1[i][j][k1];
	     hy1[i][j][k]=hy1[i][j][k]+dt/(u*dr)*y1;
        z1=ex1[i][j1][k]-ex1[i][j][k]+ey1[i][j][k]-ey1[i1][j][k];
        hz1[i][j][k]=hz1[i][j][k]+dt/(u*dr)*z1;    }
      }
//////////////////////////////////////////////////////////////////////////////
//计算小区域磁场
void smallarea_h(){
	int i,j,k,i1,j1,k1;
  float x1,y1,z1;
    for(i=0;i<=X2-1;i++)
     for(j=0;j<=Y-1;j++)
      for(k=1;k<=Z2;k++)
       {i1=i+1;
        j1=j+1;
        k1=k-1;
        x1=ey2[i][j][k1]-ey2[i][j][k]+ez2[i][j][k]-ez2[i][j1][k];
        hx2[i][j][k]=hx2[i][j][k]+dt/(u*dr)*x1;