air_pml_lsj.cpp

3d-fdtd with upm ABC

# include<math.h>
# include<iostream>
# include<iomanip>

using namespace std;

#define  C            299792458
#define  Pi           3.1415926535897932384626433832795
#define  epsilon0     8.85418e-12
#define  mu0          1.256637e-6


double ***Alloc(int nx,int ny,int nz,double val)
{
	double ***p;
	int i,j,k;
	p=new double **[nx];
	for(i=0;i<nx;i++){
		p[i]=new double *[ny];
		for(j=0;j<ny;j++){
			p[i][j]=new double[nz];
			for(k=0;k<nz;k++){
				p[i][j][k]=val;
			}
		}
	}
	return p;
}
double *Alloc(int nx,double val)
{
	if (nx<=0) {
		return NULL;
	}
	double *p;
	int i;
	p=new double [nx];
	for(i=0;i<nx;i++){
		p[i]=val;
	}
	return p;
}
int main ()
{
	int		I, J, K,It,Jt,Kt, N;
	int		i, j, k, is,js,ks,n;
	double	max_frequency; 	  
	double	dx,delta_t;									
	
	//高斯脉冲的变量
	double   puls_gauss;							
	//PML边界条件的变量
	double	refmax;                                         //最大反射系数      
	
	double  temp,depth_pml;                                                
	double  ca,cb,da,db,re,rm;
	double  sigmam_air,rhomam_air,eta0;
	
	int    n_pml=8;                                  //完全匹配层网格数
	eta0=sqrt(mu0/epsilon0);
	refmax=1e-5;    
	max_frequency=20e9;           
	double srcT=1/(2*max_frequency);
	double srcT0=3*srcT;
	//-------空气网格数目----------	 
	I=50;                           
	J=20;   
	K=100;
	//-------加入PML层后在空间所画网格数目----------
	It=I+2*n_pml;          
	Jt=J+2*n_pml;
	Kt=K+2*n_pml;
	//激励源的位置
	is=It/2;
	js=Jt/2;                        
	ks=Kt/2;
	//
	dx=(C/max_frequency)/15;//0.2e-3;                      //空间步长
	delta_t=0.8*dx/(sqrt(3)*C);               //时间步长
	//------空气对应的电场场计算相关系数----------------
	ca=1;
	cb=delta_t/(epsilon0*dx); 
	da=1;
	db=delta_t/(mu0*dx);

	//allocate PML coefficients
	double  *sig=Alloc(n_pml,0);
	double  *rho=Alloc(n_pml,0);
	double  *ca0=Alloc(n_pml,0);
	double  *cb0=Alloc(n_pml,0);
	double  *da0=Alloc(n_pml,0);
	double  *db0=Alloc(n_pml,0);

	double  *sigf=Alloc(n_pml,0);
	double  *rhof=Alloc(n_pml,0);
	double  *sigh=Alloc(n_pml,0);
	double  *rhoh=Alloc(n_pml,0);
	double  *ca0f=Alloc(n_pml,0);
	double  *ca0h=Alloc(n_pml,0);
	double  *cb0f=Alloc(n_pml,0);
	double  *cb0h=Alloc(n_pml,0);
	double  *da0f=Alloc(n_pml,0);
	double  *da0h=Alloc(n_pml,0);
	double  *db0f=Alloc(n_pml,0);
	double  *db0h=Alloc(n_pml,0);


	//------申请电磁场及其相关计算系数的动态数组-------------          
	
	double ***Ex=	Alloc(It,Jt+1,Kt+1,0);
	double ***Exy=	Alloc(It,Jt+1,Kt+1,0);
	double ***Exz=	Alloc(It,Jt+1,Kt+1,0);
	double ***caExy=Alloc(It,Jt+1,Kt+1,ca);
	double ***cbExy=Alloc(It,Jt+1,Kt+1,cb);
	double ***caExz=Alloc(It,Jt+1,Kt+1,ca);
	double ***cbExz=Alloc(It,Jt+1,Kt+1,cb);
	//
	double ***Ey=	Alloc(It+1,Jt,Kt+1,0);
	double ***Eyx=	Alloc(It+1,Jt,Kt+1,0);
	double ***Eyz=	Alloc(It+1,Jt,Kt+1,0);
	double ***caEyx=Alloc(It+1,Jt,Kt+1,ca);
	double ***cbEyx=Alloc(It+1,Jt,Kt+1,cb);
	double ***caEyz=Alloc(It+1,Jt,Kt+1,ca);
	double ***cbEyz=Alloc(It+1,Jt,Kt+1,cb);
	//
	double ***Ez=	Alloc(It+1,Jt+1,Kt,0);
	double ***Ezx=	Alloc(It+1,Jt+1,Kt,0);
	double ***Ezy=	Alloc(It+1,Jt+1,Kt,0);
	double ***caEzx=Alloc(It+1,Jt+1,Kt,ca);
	double ***cbEzx=Alloc(It+1,Jt+1,Kt,cb);
	double ***caEzy=Alloc(It+1,Jt+1,Kt,ca);
	double ***cbEzy=Alloc(It+1,Jt+1,Kt,cb);
	//
	double ***Hx=	Alloc(It+1,Jt,Kt,0);
	double ***Hxy=	Alloc(It+1,Jt,Kt,0);
	double ***Hxz=	Alloc(It+1,Jt,Kt,0);
	double ***daHxy=Alloc(It+1,Jt,Kt,da);
	double ***daHxz=Alloc(It+1,Jt,Kt,da);
	double ***dbHxy=Alloc(It+1,Jt,Kt,db);
	double ***dbHxz=Alloc(It+1,Jt,Kt,db);
	//
	double ***Hy=	Alloc(It,Jt+1,Kt,0);
	double ***Hyx=	Alloc(It,Jt+1,Kt,0);
	double ***Hyz=	Alloc(It,Jt+1,Kt,0);
	double ***daHyx=Alloc(It,Jt+1,Kt,da);
	double ***daHyz=Alloc(It,Jt+1,Kt,da);
	double ***dbHyx=Alloc(It,Jt+1,Kt,db);
	double ***dbHyz=Alloc(It,Jt+1,Kt,db);
	//
	double ***Hz=	Alloc(It,Jt,Kt+1,0);
	double ***Hzx=	Alloc(It,Jt,Kt+1,0);
	double ***Hzy=	Alloc(It,Jt,Kt+1,0);
	double ***daHzx=Alloc(It,Jt,Kt+1,da);
	double ***daHzy=Alloc(It,Jt,Kt+1,da);
	double ***dbHzx=Alloc(It,Jt,Kt+1,db);
	double ***dbHzy=Alloc(It,Jt,Kt+1,db);
	//


	//------------修正PML区域电磁场计算参数---------------
	int m;
	double retmp,rmtmp;
	m=3;
	depth_pml=dx*n_pml;
	sigmam_air=-(m+1)*log(refmax)/(2*eta0*depth_pml);  //-(m+1)ln[R(0)]/(2*d*eta)  3<=m<=4        //空气对应PML的最大电导
	rhomam_air=sigmam_air*eta0*eta0;                               //空气对应PML的最大磁阻
	retmp=delta_t/epsilon0;
	rmtmp=delta_t/mu0;
	for(n=0;n<n_pml;n++) 
	{
		 sig[n]=sigmam_air*pow((n+0.5)/(n_pml+0.5),2);
		 rho[n]=rhomam_air*pow((n+1)/(n_pml+0.5),2);
		 re=sig[n]*delta_t/epsilon0;
		 rm=rho[n]*delta_t/mu0;
		 ca0[n]=exp(-re);                                             //空气对应的PML电磁场计算参数
		 cb0[n]=(1-ca0[n])/(sig[n]*dx);
		 da0[n]=exp(-rm);
		 db0[n]=(1-da0[n])/(rho[n]*dx);

	}

//----------------PML区域电场相关参数------------------------
//-------------left and right PML---------------
	//Eyx
	 for(j=0;j<Jt;j++)
	 {
		 for(k=0;k<=Kt;k++)
		 {
			 n=n_pml-1;
			 for(i=1;i<=n_pml;i++)
			 {
				 caEyx[i][j][k]=ca0[n];
				 cbEyx[i][j][k]=cb0[n];
				 n--;
			 }
			 n=0;
			 for(i=n_pml+I;i<It;i++)
			 {
				 caEyx[i][j][k]=ca0[n];
				 cbEyx[i][j][k]=cb0[n];       
				 n++;
			 }
		 }
	 }
	//Ezx
	 for(j=0;j<=Jt;j++)
	 {
		 for(k=0;k<Kt;k++)
		 {
			 n=n_pml-1;
			 for(i=1;i<=n_pml;i++)
			 {
				 caEzx[i][j][k]=ca0[n];
				 cbEzx[i][j][k]=cb0[n];
				 n--;
			 }
			 n=0;
			 for(i=n_pml+I;i<It;i++)
			 {
				 caEzx[i][j][k]=ca0[n];
				 cbEzx[i][j][k]=cb0[n];       
				 n++;
			 }
		 }
	 }
//------------front and back PML----------------
	 //Exy
	for(i=0;i<It;i++)
	{
		for(k=0;k<=Kt;k++)
		{
			n=n_pml-1;
			for(j=1;j<=n_pml;j++)
			{
				caExy[i][j][k]=ca0[n];
				cbExy[i][j][k]=cb0[n];   
				n--;
			}
			n=0;
			for(j=n_pml+J;j<Jt;j++)
			{
				caExy[i][j][k]=ca0[n];
				cbExy[i][j][k]=cb0[n];   
				n++;
			}
		}
	}
	 //Ezy
	for(i=0;i<=It;i++)
	{
		for(k=0;k<Kt;k++)
		{
			n=n_pml-1;
			for(j=1;j<=n_pml;j++)
			{
				caEzy[i][j][k]=ca0[n];
				cbEzy[i][j][k]=cb0[n];   
				n--;
			}
			n=0;
			for(j=n_pml+J;j<Jt;j++)
			{
				caEzy[i][j][k]=ca0[n];
				cbEzy[i][j][k]=cb0[n];   
				n++;
			}
		}
	}

//------------down and up PML----------------
	//Exz
	for(i=0;i<It;i++)
	{
		for(j=0;j<=Jt;j++)
		{
			n=n_pml-1;
			for(k=1;k<=n_pml;k++)
			{
				caExz[i][j][k]=ca0[n];
				cbExz[i][j][k]=cb0[n];
				n--;
			}       
			n=0;
			for(k=n_pml+K;k<Kt;k++)
			{
				caExz[i][j][k]=ca0[n];
				cbExz[i][j][k]=cb0[n];
				n++;
			}
		}
	}
	//Eyz
	for(i=0;i<=It;i++)
	{
		for(j=0;j<Jt;j++)
		{
			n=n_pml-1;
			for(k=1;k<=n_pml;k++)
			{
				caEyz[i][j][k]=ca0[n];
				cbEyz[i][j][k]=cb0[n];
				n--;
			}       
			n=0;
			for(k=n_pml+K;k<Kt;k++)
			{
				caEyz[i][j][k]=ca0[n];
				cbEyz[i][j][k]=cb0[n];
				n++;
			}
		}
	}


 //----------------PML区域磁场相关参数设置------------------------ 
   //-------------left and right PML---------------
	//Hyx
	 for(j=0;j<=Jt;j++)
	 {
		 for(k=0;k<Kt;k++)
		 {
			 n=n_pml-1;
			 for(i=0;i<n_pml;i++)
			 {
				 daHyx[i][j][k]=da0[n];
				 dbHyx[i][j][k]=db0[n];
				 n--;
			 }
			 n=0;
			 for(i=n_pml+I;i<It;i++)
			 {
				 daHyx[i][j][k]=da0[n];
				 dbHyx[i][j][k]=db0[n];     
				 n++;
			 }
		 }
	 }
	//Hzx
	 for(j=0;j<Jt;j++)
	 {
		 for(k=0;k<=Kt;k++)
		 {
			 n=n_pml-1;
			 for(i=0;i<n_pml;i++)
			 {
				 daHzx[i][j][k]=da0[n];
				 dbHzx[i][j][k]=db0[n];
				 n--;
			 }
			 n=0;
			 for(i=n_pml+I;i<It;i++)
			 {
				 daHzx[i][j][k]=da0[n];
				 dbHzx[i][j][k]=db0[n];     
				 n++;
			 }
		 }
	 }
//------------Front and back PML----------------
	 //Hxy
	for(i=0;i<=It;i++)
	{
		for(k=0;k<Kt;k++)
		{
			n=n_pml-1;
			for(j=0;j<n_pml;j++)
			{
				daHxy[i][j][k]=da0[n];
				dbHxy[i][j][k]=db0[n]; 
				n--;
			}
			n=0;
			for(j=n_pml+J;j<Jt;j++)
			{
				daHxy[i][j][k]=da0[n];
				dbHxy[i][j][k]=db0[n]; 
				n++;
			}
		}
	}
	 //Hzy
	for(i=0;i<It;i++)
	{
		for(k=0;k<=Kt;k++)
		{
			n=n_pml-1;
			for(j=0;j<n_pml;j++)
			{
				daHzy[i][j][k]=da0[n];
				dbHzy[i][j][k]=db0[n]; 
				n--;
			}
			n=0;
			for(j=n_pml+J;j<Jt;j++)
			{
				daHzy[i][j][k]=da0[n];
				dbHzy[i][j][k]=db0[n]; 
				n++;
			}
		}
	}
//------------down and up PML-----------------------
	//Hxz
	for(i=0;i<=It;i++)
	{
		for(j=0;j<Jt;j++)
		{
			n=n_pml-1;
			for(k=0;k<n_pml;k++)
			{
				daHxz[i][j][k]=da0[n];
				dbHxz[i][j][k]=db0[n];
				n--;
			}       
			n=0;
			for(k=n_pml+K;k<Kt;k++)
			{
				daHxz[i][j][k]=da0[n];
				dbHxz[i][j][k]=db0[n];
				n++;
			}
		}
	}
	//Hyz
	for(i=0;i<It;i++)
	{
		for(j=0;j<=Jt;j++)
		{
			n=n_pml-1;
			for(k=0;k<n_pml;k++)
			{
				daHyz[i][j][k]=da0[n];
				dbHyz[i][j][k]=db0[n];
				n--;
			}       
			n=0;
			for(k=n_pml+K;k<Kt;k++)
			{
				daHyz[i][j][k]=da0[n];
				dbHyz[i][j][k]=db0[n];
				n++;
			}
		}
	}

	//--------------------计算传播时间和高斯脉冲生成时间-------------------
	
	N=600;
	//----------------主循环--------------------
	cout<<"主循环"<<endl;
	cout<<"完成进度:"<<setw(4)<<setw(4)<<N<<endl;
	putchar(8);putchar(8);putchar(8);putchar(8);putchar(8);
	putchar(8);putchar(8);putchar(8);putchar(8);

	double t;
	for(n=1;n<=N;n++)
	{
		//----------------计算高斯脉冲-----------------------
		t=n*delta_t;
		temp= (t-srcT0)/srcT;
		puls_gauss=exp(-temp*temp);

		cout<<"step="<<n<<endl;
		//------------------电场的循环迭代---------------------
		for(i=0;i<It;i++)
		{
			for(j=1;j<Jt;j++)
			{
				for(k=1;k<Kt;k++)
				{
					Exy[i][j][k]=caExy[i][j][k]*Exy[i][j][k]+cbExy[i][j][k]*(Hz[i][j][k]-Hz[i][j-1][k]);
					Exz[i][j][k]=caExz[i][j][k]*Exz[i][j][k]+cbExz[i][j][k]*(-Hy[i][j][k]+Hy[i][j][k-1]);
					Ex[i][j][k]=Exy[i][j][k]+Exz[i][j][k];
				}
			}
		}
		//boundary
		for(k=0;k<=Kt;k+=Kt){
			for(i=0;i<It;i++){
				for(j=0;j<=Jt;j++){
					Ex[i][j][k]=0;
				}
			}
		}
		for(j=0;j<=Jt;j+=Jt){
			for(i=0;i<It;i++){
				for(k=0;k<=Kt;k++){
					Ex[i][j][k]=0;
				}
			}
		}
		////
		for(i=1;i<It;i++)
		{
			for(j=0;j<Jt;j++)
			{
				for(k=1;k<Kt;k++)
				{
					Eyx[i][j][k]=caEyx[i][j][k]*Eyx[i][j][k]+cbEyx[i][j][k]*(-Hz[i][j][k]+Hz[i-1][j][k]);
					Eyz[i][j][k]=caEyz[i][j][k]*Eyz[i][j][k]+cbEyz[i][j][k]*(Hx[i][j][k]-Hx[i][j][k-1]);
					Ey[i][j][k]=Eyx[i][j][k]+Eyz[i][j][k];
				}
			}
		}		
		//boundary
		for(k=0;k<=Kt;k+=Kt){
			for(i=0;i<=It;i++){
				for(j=0;j<Jt;j++){
					Ey[i][j][k]=0;
				}
			}
		}
		for(i=0;i<=It;i+=It){
			for(j=0;j<Jt;j+=Jt){
				for(k=0;k<=Kt;k++){
					Ey[i][j][k]=0;
				}
			}
		}
		////
		for(i=1;i<It;i++)
		{
			for(j=1;j<Jt;j++)
			{
				for(k=0;k<Kt;k++)
				{
					Ezx[i][j][k]=caEzx[i][j][k]*Ezx[i][j][k]+cbEzx[i][j][k]*(Hy[i][j][k]-Hy[i-1][j][k]);
					Ezy[i][j][k]=caEzy[i][j][k]*Ezy[i][j][k]+cbEzy[i][j][k]*(-Hx[i][j][k]+Hx[i][j-1][k]);
					Ez[i][j][k]=Ezx[i][j][k]+Ezy[i][j][k];
				}
			}
		}		
		//boundary
		for(j=0;j<=Jt;j+=Jt){
			for(i=0;i<=It;i++){
				for(k=0;k<Kt;k++){
					Ez[i][j][k]=0;
				}
			}
		}
		for(i=0;i<=It;i+=It){
			for(j=0;j<=Jt;j+=Jt){
				for(k=0;k<Kt;k++){
					Ez[i][j][k]=0;
				}
			}
		}
		////
		//---------------------------加入高斯脉冲-----------------		
		Ex[is][js][ks]=delta_t*puls_gauss/epsilon0+Ex[is][js][ks];

		//-----------------------------循环迭代H-----------------------------
		//hx
		for(i=0;i<=It;i++){
			for(j=0;j<Jt;j++){
				for(k=0;k<Kt;k++){
					Hxy[i][j][k]=daHxy[i][j][k]*Hxy[i][j][k]-dbHxy[i][j][k]*(Ez[i][j+1][k]-Ez[i][j][k]);
					Hxz[i][j][k]=daHxz[i][j][k]*Hxz[i][j][k]+dbHxz[i][j][k]*(Ey[i][j][k+1]-Ey[i][j][k]);
					Hx[i][j][k]=Hxy[i][j][k]+Hxz[i][j][k];
				}
			}
		}
		//hy
		for(i=0;i<It;i++){
			for(j=0;j<=Jt;j++){
				for(k=0;k<Kt;k++){
					Hyx[i][j][k]=daHyx[i][j][k]*Hyx[i][j][k]+dbHyx[i][j][k]*(Ez[i+1][j][k]-Ez[i][j][k]);
					Hyz[i][j][k]=dbHyz[i][j][k]*Hyz[i][j][k]-dbHyz[i][j][k]*(Ex[i][j][k+1]-Ex[i][j][k]);
					Hy[i][j][k]=Hyx[i][j][k]+Hyz[i][j][k];
				}
			}
		}
		//hz
		for(i=0;i<It;i++){
			for(j=0;j<Jt;j++){
				for(k=0;k<=Kt;k++){
					Hzx[i][j][k]=daHzx[i][j][k]*Hzx[i][j][k]-dbHzx[i][j][k]*(Ey[i+1][j][k]-Ey[i][j][k]);
					Hzy[i][j][k]=daHzy[i][j][k]*Hzy[i][j][k]+dbHzy[i][j][k]*(Ex[i][j+1][k]-Ex[i][j][k]);
					Hz[i][j][k]=Hzx[i][j][k]+Hzy[i][j][k];				
				}
			}
		}
	}
	
	//-----------------------------------结束-----------------------------------
	return 0;
}