fdtd_xypbc_zpml.cpp

利用c++语言写的三维FDTD

#include "stdafx.h"
#include "FDTD_xyPBC_zPML.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CFDTD_xyPBC_zPML::CFDTD_xyPBC_zPML(void)
{
	//The real and imaginary field components
	Ex = NULL; Ey = NULL; Fz = NULL; Ez = NULL; 
	Hx = NULL; Hy = NULL; Gz = NULL; Hz = NULL;
	
	//Coefficients containing the PML boundary parameters
	K_E1_a = NULL; K_E1_b = NULL; K_E3_a = NULL; K_E3_b = NULL; 
	K_E4_a = NULL; K_E4_b = NULL; K_E6_a = NULL; K_E6_b = NULL;

	//the incident field for total field scattered field formulation
	Ex_1D = NULL; Hy_1D = NULL;

	//Contains the followed field components
	Ex_Foll = NULL;	 Ey_Foll = NULL;	 Ez_Foll = NULL;
	Hx_Foll = NULL;	 Hy_Foll = NULL;	 Hz_Foll = NULL;

	//Contains the followed field components
	Ex_Foll_av = NULL;	 Ey_Foll_av = NULL;	 Ez_Foll_av = NULL;
	Hx_Foll_av = NULL;	 Hy_Foll_av = NULL;	 Hz_Foll_av = NULL;

	/////////////////////////
	//Fourier transform in a subvolume
	/////////////////////////
	Ex_fourier_TR_re = NULL; Ex_fourier_TR_im = NULL; 
	Ey_fourier_TR_re = NULL; Ey_fourier_TR_im = NULL; 
	Ez_fourier_TR_re = NULL; Ez_fourier_TR_im = NULL; 
	Hx_fourier_TR_re = NULL; Hx_fourier_TR_im = NULL; 
	Hy_fourier_TR_re = NULL; Hy_fourier_TR_im = NULL; 
	Hz_fourier_TR_re = NULL; Hz_fourier_TR_im = NULL; 
	fourier_omega = NULL;

	//Bragg diffraction
	Ex_Br_av_re = NULL; Ey_Br_av_re = NULL; Ez_Br_av_re = NULL;
	Hx_Br_av_re = NULL; Hy_Br_av_re = NULL; Hz_Br_av_re = NULL;
	Ex_Br_av_im = NULL; Ey_Br_av_im = NULL; Ez_Br_av_im = NULL;
	Hx_Br_av_im = NULL; Hy_Br_av_im = NULL; Hz_Br_av_im = NULL;
	Gxy = NULL;

	Mask_Wigner = NULL;

	//to save the field components
	path_name_Ex = NULL; 
	path_name_Ey = NULL; 
	path_name_Ez = NULL;
	path_name_Hx = NULL;
	path_name_Hy = NULL;
	path_name_Hz = NULL;

	path_name_Ex_1D = NULL;

	#ifndef run_enviroment
		Path_Load_Workspace_Data = NULL;
		Path_Save_Workspace_Data = NULL;
	#endif

	nr_threads = 1;
}

CFDTD_xyPBC_zPML::~CFDTD_xyPBC_zPML(void)
{
	Free_Mem();
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize the number of threads
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::Init_nr_THR(int nr_Threads)
{
	nr_threads = nr_Threads;
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize the parameters of the class, allocates the memory for matrices
///////////////////////////////////////////////////////////////////////////////////////
bool CFDTD_xyPBC_zPML::Init(long ***&Index, long n_x, long n_y, long n_z, double **&Mater, 
								 long nMat, long n_PML, double d_t, double d_x, double d_y, 
								 double d_z)
{   
	pi = 180.0*atan(1.0)/45.0;
	//permittivity of free space 
	eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]
	//permeability of free space 
    mu_0  = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6;  // [H/m]

	TwoORpi = 2.0*pi;

	Ind = Index; //contains the indices of the material type. Dimensions [nx][ny][nz]

	//dimensions of the computational space
    nx = n_x;
	ny = n_y; 
	nz = n_z;

	//time increment
	dt = d_t;

	//elementary cell sizes
	dx = d_x;
	dy = d_y;
	dz = d_z;
	inv_dx = 1.0/dx;
	inv_dy = 1.0/dy;
	inv_dz = 1.0/dz;

	//coordinates of the PML region
	n_z_l = n_PML;
	n_z_h = nz - n_PML;

	NORMA = 1.0;

	//contains the permittivity and conductivity of each material from the geometry
	Mat = Mater;
	//number of materials presents in the investigated geometry
	n_Mat = nMat;

	//Allocate memory for matrices, initializing them with zero
	//Field components
	Ex = Init_Matrix_3D<double>(n_x,n_y,n_z);
	if(!Ex)
	{
		Free_Mem();
		return false;
	}

	Ey = Init_Matrix_3D<double>(n_x,n_y,n_z);
	if(!Ey)
	{
		Free_Mem();
		return false;
	}

	Fz = Init_Matrix_3D<double>(n_x,n_y,n_z-1);
	if(!Fz)
	{
		Free_Mem();
		return false;
	}

	Ez = Init_Matrix_3D<double>(n_x,n_y,n_z-1);
	if(!Ez)
	{
		Free_Mem();
		return false;
	}

	Hx = Init_Matrix_3D<double>(n_x,n_y,n_z-1);
	if(!Hx)
	{
		Free_Mem();
		return false;
	}

	Hy = Init_Matrix_3D<double>(n_x,n_y,n_z-1);
	if(!Hy)
	{
		Free_Mem();
		return false;
	}

	Gz = Init_Matrix_3D<double>(n_x,n_y,n_z);
	if(!Gz)
	{
		Free_Mem();
		return false;
	}

	Hz = Init_Matrix_3D<double>(n_x,n_y,n_z);
	if(!Hz)
	{
		Free_Mem();
		return false;
	}

	//Coefficients containing the PML boundary parameters
	K_E1_a = (double *)calloc(nz,sizeof(double)); 
	if(!K_E1_a)
	{
		Free_Mem();
		return false;
	}

	K_E1_b = (double *)calloc(nz,sizeof(double)); 
	if(!K_E1_b)
	{
		Free_Mem();
		return false;
	}

	K_E3_a = (double *)calloc(nz-1,sizeof(double)); 
	if(!K_E3_a)
	{
		Free_Mem();
		return false;
	}

	K_E3_b = (double *)calloc(nz-1,sizeof(double)); 
	if(!K_E3_b)
	{
		Free_Mem();
		return false;
	}

	K_E4_a = (double *)calloc(nz-1,sizeof(double)); 
	if(!K_E4_a)
	{
		Free_Mem();
		return false;
	}

	K_E4_b = (double *)calloc(nz-1,sizeof(double)); 
	if(!K_E4_b)
	{
		Free_Mem();
		return false;
	}

	K_E6_a = (double *)calloc(nz,sizeof(double)); 
	if(!K_E6_a)
	{
		Free_Mem();
		return false;
	}

	K_E6_b = (double *)calloc(nz,sizeof(double)); 
	if(!K_E6_b)
	{
		Free_Mem();
		return false;
	}

	return true;
}

///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::Init_PML_Par(double eps_r_a, double mu_r_a,
									double eps_r_b, double mu_r_b)
{
	long k; 
	//Outside of the PML region
	for (k = 0; k < nz-1; k++)
	{
		K_E1_a[k] = 1.0;
		K_E1_b[k] = dt/eps_0;

		K_E3_a[k] = 1.0/eps_0;
		K_E3_b[k] = 1.0/eps_0;

		K_E4_a[k] = 1.0;
		K_E4_b[k] = dt/mu_0;

		K_E6_a[k] = 1.0/mu_0;
		K_E6_b[k] = 1.0/mu_0;
	}
	K_E1_a[nz-1] = 1.0;
	K_E1_b[nz-1] = dt/eps_0;

	K_E6_a[nz-1] = 1.0/mu_0;
	K_E6_b[nz-1] = 1.0/mu_0;

	//PML parameters
	double ka_max = 1.0;
	double exponent = 4.0;
	double R_err = 1.0e-16;
	
	double eta = sqrt(mu_0*mu_r_a/eps_0/eps_r_a);

	double sigma_z, ka_z;
	long n_PML = n_z_l;
	double sigma_max= -(exponent + 1.0)*log(R_err)/(2.0*eta*n_PML*dz);
	
	//1
	for (k = 0; k < n_PML; k++)
	{
		sigma_z         = sigma_max*pow( (n_PML - k)/((double) n_PML) ,exponent);
		ka_z            = 1.0 + (ka_max - 1.0)*pow( (n_PML - k)/((double) n_PML) ,exponent);

		K_E1_a[k]       = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
		
		K_E1_b[k]       = 2.0*dt/(2.0*eps_0*ka_z + sigma_z*dt);
		
		K_E6_a[k]       = (2.0*eps_0*ka_z + sigma_z*dt)/(2.0*eps_0*mu_0);
		
		K_E6_b[k]       = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*mu_0);
		
		sigma_z         = sigma_max*pow( (n_PML - k - 0.5)/n_PML ,exponent);
		ka_z            = 1.0 + (ka_max - 1.0)*pow( (n_PML - k - 0.5)/n_PML ,exponent);

        K_E3_a[k]       = (2.0*eps_0*ka_z + sigma_z*dt)/(2.0*eps_0*eps_0);
		
		K_E3_b[k]       = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*eps_0);
		
		K_E4_a[k]       = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
		
		K_E4_b[k]       = 2.0*eps_0*dt/(2.0*eps_0*ka_z + sigma_z*dt)/mu_0;
	}

	eta = sqrt(mu_0*mu_r_b/eps_0/eps_r_b);
	sigma_max= -(exponent + 1.0)*log(R_err)/(2.0*eta*n_PML*dz);
    //2
	for (k = 0; k < n_PML; k++)
	{
		sigma_z         = sigma_max*pow( (n_PML - k)/((double) n_PML) ,exponent);
		ka_z            = 1.0 + (ka_max - 1.0)*pow( (n_PML - k)/((double) n_PML) ,exponent);

		K_E1_a[nz-k-1]  = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
		
		K_E1_b[nz-k-1]  = 2.0*dt/(2.0*eps_0*ka_z + sigma_z*dt);
		
		K_E6_a[nz-k-1]  = (2.0*eps_0*ka_z + sigma_z*dt)/(2.0*eps_0*mu_0);
		
		K_E6_b[nz-k-1]  = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*mu_0);
		
		sigma_z         = sigma_max*pow( (n_PML - k - 0.5)/n_PML ,exponent);
		ka_z            = 1.0 + (ka_max - 1.0)*pow( (n_PML - k - 0.5)/n_PML ,exponent);

        K_E3_a[nz-k-2]  = (2.0*eps_0*ka_z + sigma_z*dt)/(2.0*eps_0*eps_0);
		
		K_E3_b[nz-k-2]  = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*eps_0);
		
		K_E4_a[nz-k-2]  = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
		
		K_E4_b[nz-k-2]  = 2.0*eps_0*dt/(2.0*eps_0*ka_z + sigma_z*dt)/mu_0;
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Ex field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::calc_Ex()
{
	long i, j, k;
	double eps_r;

	#pragma omp parallel default(shared) private(i,j,k,eps_r)
	{
		#pragma omp for schedule(dynamic,nr_threads)
		for (i = 0; i < nx; i++)
		{	
			for (j = 1; j < ny; j++)
			{	
				for (k = 1; k < nz-1; k++)//PEC at k = 0, k = nz-1
				{
					eps_r = Mat[Ind[i][j][k]][0];
					
					Ex[i][j][k] = K_E1_a[k]*Ex[i][j][k] + K_E1_b[k]/eps_r*
									( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
										(Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );
				}
			}
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Ey field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::calc_Ey()
{
	long i, j, k;
	double eps_r;

	#pragma omp parallel default(shared) private(i,j,k,eps_r)
	{
		#pragma omp for schedule(dynamic,nr_threads)
		for (i = 1; i < nx; i++)
		{	
			for (j = 0; j < ny; j++)
			{	
				for (k = 1; k < nz-1; k++)//PEC at k = 0, k = nz-1
				{
					eps_r = Mat[Ind[i][j][k]][0];
					
					Ey[i][j][k] = K_E1_a[k]*Ey[i][j][k] + K_E1_b[k]/eps_r*
									( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
									(Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );
				}
			}
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Ez field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::calc_Ez()
{
	long i, j, k;
	double eps_r;
	double Fz_r;
	double dtORinv_dx, dtORinv_dy, kkx, kky;
	
	#pragma omp parallel default(shared) private(i,j,k,eps_r,Fz_r,kkx,kky)
	{
		dtORinv_dx = dt*inv_dx;
	    dtORinv_dy = dt*inv_dy;

		//1 - First PML region
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (i = 1; i < nx; i++)
		{	
			for (j = 1; j < ny; j++)
			{	
				for (k = 0; k < n_z_l; k++)
				{
					eps_r = Mat[Ind[i][j][k]][0];

					Fz_r = Fz[i][j][k];

					Fz[i][j][k] += (Hy[i][j][k] - Hy[i-1][j][k])*dtORinv_dx -
								(Hx[i][j][k] - Hx[i][j-1][k])*dtORinv_dy;
		
					Ez[i][j][k] += (K_E3_a[k]*Fz[i][j][k] - K_E3_b[k]*Fz_r)/eps_r;
				}
			}
		}
	
		//2 - Outside of the PML
		kkx = dtORinv_dx/eps_0;
		kky = dtORinv_dy/eps_0;

		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (i = 1; i < nx; i++)
		{	
			for (j = 1; j < ny; j++)
			{	
				for (k = n_z_l; k < n_z_h-1; k++)
				{
					eps_r = Mat[Ind[i][j][k]][0];

					Ez[i][j][k] += ( (Hy[i][j][k] - Hy[i-1][j][k])*kkx -
									 (Hx[i][j][k] - Hx[i][j-1][k])*kky)/eps_r;
				}
			}
		}

		//3 - Second PML region
		#pragma omp for schedule(dynamic,nr_threads)
		for (i = 1; i < nx; i++)
		{	
			for (j = 1; j < ny; j++)
			{	
				for (k = n_z_h-1; k < nz-1; k++)
				{
					
					eps_r = Mat[Ind[i][j][k]][0];

					Fz_r = Fz[i][j][k];

					Fz[i][j][k] += (Hy[i][j][k] - Hy[i-1][j][k])*dtORinv_dx -
								   (Hx[i][j][k] - Hx[i][j-1][k])*dtORinv_dy;

					Ez[i][j][k] += (K_E3_a[k]*Fz[i][j][k] - K_E3_b[k]*Fz_r)/eps_r;
				}
			}
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Hx field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::calc_Hx()
{
	long i, j, k;
	double mu_r;

	#pragma omp parallel default(shared) private(i,j,k,mu_r)
	{
		#pragma omp for schedule(dynamic,nr_threads)
		for (i = 0; i < nx; i++)
		{	
			for (j = 0; j < ny-1; j++)
			{	
				for (k = 0; k < nz-1; k++)
				{
					mu_r = Mat[Ind[i][j][k]][1];

					Hx[i][j][k] = K_E4_a[k]*Hx[i][j][k] + K_E4_b[k]/mu_r*
									 ( (Ey[i][j][k+1] - Ey[i][j][k])*inv_dz -
									   (Ez[i][j+1][k] - Ez[i][j][k])*inv_dy);
				}
			}
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Hy field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::calc_Hy()
{
	long i, j, k;
	double mu_r;

	#pragma omp parallel default(shared) private(i,j,k,mu_r)
	{
		#pragma omp for schedule(dynamic,nr_threads)
		for (i = 0; i < nx-1; i++)
		{	
			for (j = 0; j < ny; j++)
			{	
				for (k = 0; k < nz-1; k++)
				{
					mu_r = Mat[Ind[i][j][k]][1];

					Hy[i][j][k] = K_E4_a[k]*Hy[i][j][k] + K_E4_b[k]/mu_r*