fd2d.c.txt

FDTD的并行算法

#include<stdio.h>
#include<stdlib.h>
#include<memory.h>
#include<math.h>
#include <mpi.h>
#include <complex.h>
#include <iostream.h>
#include <time.h>
#include "vectors.h"
#include "buffers.h"
#include "rec.h"
#include "mark.h"

#undef hz

#define c 2.99795637e8
#define eps0 8.854e-12
#define mu0 1.2566306e-6

#define dx 1e-3
#define dy 1e-3
#define ny 100
#define nmax 500

#define w_pulse 15e-12


void main(int argc, char *argv[])
{
	int record_flag=0;
	float t0, dt;
	int i, j, n, lim;
	float	**ey, **ex, **hz, **cbx,
		**dbx, **cby, **dby, **erx, **ery, **erec,
		kx, ky;

	float	**exbufl1, **exbufl2, **exbufr1, **exbufr2,
		**eybufl1, **eybufl2, **eybufr1, **eybufr2;
		
	float	t, exc, e_rec;
	float 	**erecf;
	
	time_t proc_start, proc_finish;
	
	FILE	*fpgnu;
	
	MPI_Comm com = MPI_COMM_WORLD;
    	int size, rank;
    	MPI_Status status;
   
    	MPI_Init(&argc, &argv);
    	MPI_Comm_size(MPI_COMM_WORLD, &size);
    	MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 	
	
	int nxo=2000, nx;

//======================================================================================
	int *nx_size;
	nx_size=(int*)calloc(size,sizeof(int));

	bench_test(size, rank, status, nxo, &nx, nx_size);	
	if(rank==0)
	{
		printf("Benchmark Complete.\n");
	}


//	printf("%d ",nx);
	
	if(rank==0)
	{
	for(i=0;i<=size-1;i++)
	{
		printf("%d ",nx_size[i]);
	}
	printf("\n");
	}
	
//======================================================================================
	proc_start=time(&proc_start);
	
	if(rank==0)
	{
	fpgnu=fopen("pulse.gnu","w");
	fprintf(fpgnu,"se zra[0:1]\n");
	}
		
	dt=1/(sqrt(1/sq(dx)+1/sq(dy))*c);
	t0=4.0*w_pulse;
	kx=(c*dt-dx)/(c*dt+dx);
	ky=(c*dt-dy)/(c*dt+dy);
	
	ex=exa(nx);
	ey=eya(nx);
	hz=hza(nx);
	cbx=cbxa(nx);
	dbx=dbxa(nx);
	cby=cbya(nx);
	dby=dbya(nx);
	erx=erxa(nx);
	ery=erya(nx);
	erec=ereca(nx);

	if(rank==0)
	{
		erecf=erecfa(nxo);
	}
	
	exbufl1=exbufl1a(nx);
	exbufl2=exbufl2a(nx);
	exbufr1=exbufr1a(nx);
	exbufr2=exbufr2a(nx);
	eybufl1=eybufl1a(nx);
	eybufl2=eybufl2a(nx);
	eybufr1=eybufr1a(nx);
	eybufr2=eybufr2a(nx);
	
	for(i=0;i<=nx+1;i++)
	{
		for(j=0;j<=ny+1;j++)
		{
			erx[i][j]=1;
			ery[i][j]=1;	
		}	
	}


	for(i=0;i<=nx+1;i++)
	{
		for(j=0;j<=ny+1;j++)
		{
			cbx[i][j]=dt/(eps0*erx[i][j]*dx);	
			dbx[i][j]=dt/(mu0*dx);	
			cby[i][j]=dt/(eps0*ery[i][j]*dy);	
			dby[i][j]=dt/(mu0*dy);
		}	
	}

	
//================================================================
//                    MAIN FDTD LOOP
//================================================================
	
	for(n=0;n<0.4*nmax;n++)
	{
		for(i=1;i<=nx;i++)
		{
			for(j=0;j<=ny+1;j++)
			{
					ey[i][j]=ey[i][j]-cbx[i][j]*(hz[i][j]-hz[i-1][j]);
					
			}
		}

		for(i=0;i<=nx+1;i++)
		{
			for(j=1;j<=ny;j++)
			{
					ex[i][j]=ex[i][j]+cby[i][j]*(hz[i][j]-hz[i][j-1]);
					
			}
		}
	
//============================================================
//			E TRADE
	if(rank!=0)
	{
		float sendbuf[ny+2];
		for(j=0;j<=ny+1;j++)
		{
			sendbuf[j]=ey[1][j];
		}
		
		MPI_Send(&sendbuf,ny+1,MPI_FLOAT,rank-1,8,MPI_COMM_WORLD);			
	}

	if((rank+1)!=size)
	{
		float recbuf[ny+2];

		MPI_Recv(&recbuf,ny+1,MPI_FLOAT,rank+1,8,MPI_COMM_WORLD,&status);			

		for(j=0;j<=ny+1;j++)
		{
			ey[nx+1][j]=recbuf[j];
		}
		
	}

//============================================================	
	
	
	
		t=n*dt;
		exc=exp((-1*sq(t-t0))/(sq(w_pulse)));
		
		if(rank==(size/2))
		{
			int centre;
			centre=nx/2;
			hz[centre][49]+=0.1*exc;
		}
//===============================================================
//                   APPLY MUR ABC
//===============================================================			
		
		for(i=0;i<=nx+1;i++)
		{
			exbufl1[i][1]=exbufl1[i][0];
			exbufl1[i][0]=ex[i][1];
			exbufr1[i][1]=exbufr1[i][0];
			exbufr1[i][0]=ex[i][ny];
			ex[i][0]=exbufl1[i][1]+kx*(exbufl1[i][0]-exbufl1[i][1]);
			ex[i][ny+1]=exbufr1[i][1]+kx*(exbufr1[i][0]-exbufr1[i][1]);
		}
		
		if(rank==0)
		{
		for(j=0;j<=ny+1;j++)
		{
			eybufl1[1][j]=eybufl1[0][j];
			eybufl1[0][j]=ey[1][j];
			ey[0][j]=eybufl1[1][j]+ky*(eybufl1[0][j]-eybufl1[1][j]);
		}
		}
		if((rank+1)==size)
		{
		for(j=0;j<=ny+1;j++)
		{
			eybufr1[1][j]=eybufr1[0][j];
			eybufr1[0][j]=ey[nx][j];
			ey[nx+1][j]=eybufr1[1][j]+ky*(eybufr1[0][j]-eybufr1[1][j]);
		}
		}

		
		
		
		for(i=1;i<=nx;i++)
		{
			for(j=1;j<=ny;j++)
			{
				erec[i][j]=sqrt(sq(ey[i][j])+sq(ex[i][j]));
			}
		}
			
		
		
	
		for(i=0;i<=nx;i++)
		{
			for(j=0;j<=ny;j++)
			{
				hz[i][j]=hz[i][j]-dbx[i][j]*(ey[i+1][j]-ey[i][j])+
								  dby[i][j]*(ex[i][j+1]-ex[i][j]);
			}
		}
		
//============================================================
//			H TRADE

	if((rank+1)!=size)
	{
		float sendbuf[ny+2];
		for(j=0;j<=ny+1;j++)
		{
			sendbuf[j]=hz[nx][j];
		}
		
		MPI_Send(&sendbuf,ny+1,MPI_FLOAT,rank+1,99,MPI_COMM_WORLD);			
	}
	

	if(rank!=0)
	{
		float recbuf[ny+2];
		MPI_Recv(&recbuf,ny+1,MPI_FLOAT,rank-1,99,MPI_COMM_WORLD,&status);			

		for(j=0;j<=ny+1;j++)
		{
			hz[0][j]=recbuf[j];
		}

	}

//============================================================	

		if((n%20==0)&&(record_flag!=0))
		{
		if(rank==0)
		{
			rec_data(erec,fpgnu,rank,erecf,nxo,status,size,nx,nx_size,n);
		}
		else
		{
		int	s=0;
		s=nx*ny;
		float *sendbuf;
		sendbuf=(float*)calloc(s,sizeof(float));
		for(i=1;i<=nx;i++)
		{
			for(j=1;j<=ny;j++)
			{
				sendbuf[(i-1)+(j-1)*nx]=erec[i][j];
			}
		}
			MPI_Send(&sendbuf,s,MPI_FLOAT,0,22,MPI_COMM_WORLD);			
		}
		}
	
	}
	
	if(rank==0)
	{
		fclose(fpgnu);
		proc_finish=time(&proc_finish);
		printf("\nTotal Time = %d secs\n\n",(int)proc_finish-(int)proc_start);

	}
	
	MPI_Finalize();
	
}