fdtd_xypbc_zpml.cpp

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

									( (Ez[i+1][j][k] - Ez[i][j][k])*inv_dx -
									  (Ex[i][j][k+1] - Ex[i][j][k])*inv_dz );
				}
			}
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Hz field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::calc_Hz()
{
	long i, j, k;
	double mu_r;
	double Gz_r;
	double dtORinv_dx, dtORinv_dy;

	#pragma omp parallel default(shared) private(i,j,k,mu_r,dtORinv_dx,dtORinv_dy)
	{
		dtORinv_dx = dt*inv_dx;
		dtORinv_dy = dt*inv_dy;

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

					Gz_r = Gz[i][j][k];

					Gz[i][j][k] = Gz[i][j][k] +  (Ex[i][j+1][k] - Ex[i][j][k])*dtORinv_dy - 
												 (Ey[i+1][j][k] - Ey[i][j][k])*dtORinv_dx;

					Hz[i][j][k] = Hz[i][j][k] + (K_E6_a[k]*Gz[i][j][k] - K_E6_b[k]*Gz_r)/mu_r;
				}
			}
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Periodic boundary condition for Ex
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::PBC_Ex()
{
	long i, k;
	//j = 0
	#pragma omp parallel for private(i,k) schedule(dynamic,nr_threads)
	for (i = 0; i < nx; i++)
	{	
		for (k = 1; k < nz-1; k++) //PEC at k = 0, k = nz-1
		{
			Ex[i][0][k] =  Ex[i][ny-1][k];
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Periodic boundary condition for Ey
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::PBC_Ey()
{
	long j, k;
	//i = 0
	#pragma omp parallel for private(j,k) schedule(dynamic,nr_threads)
	for (j = 0; j < ny; j++)
	{	
		for (k = 1; k < nz-1; k++) //PEC at k = 0, k = nz-1  
		{
			Ey[0][j][k] =  Ey[nx-1][j][k];
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Periodic boundary condition for Ez
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::PBC_Ez()
{
	long i, j, k;
	// i = 0;
	#pragma omp parallel for private(j,k) schedule(dynamic,nr_threads)
	for (j = 0; j < ny; j++)
	{	
		for (k = 0; k < nz-1; k++)
		{
			Ez[0][j][k] =  Ez[nx-1][j][k];
		}
	}

	//j = 0
	#pragma omp parallel for private(i,k) schedule(dynamic,nr_threads)
	for (i = 0; i < nx; i++)
	{	
		for (k = 0; k < nz-1; k++)
		{
			Ez[i][0][k] =  Ez[i][ny-1][k];
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Periodic boundary condition for Hx
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::PBC_Hx()
{
	long i, k;
	// j = ny - 1;
	#pragma omp parallel for private(i,k) schedule(dynamic,nr_threads)
	for (i = 0; i < nx; i++)
	{	
		for (k = 0; k < nz-1; k++)
		{
			Hx[i][ny-1][k] = Hx[i][0][k];
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Periodic boundary condition for Hy
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::PBC_Hy()
{
	long j, k;
	// i = nx-1;
	#pragma omp parallel for private(j,k) schedule(dynamic,nr_threads)
	for (j = 0; j < ny; j++)
	{	
		for (k = 0; k < nz-1; k++)
		{
			Hy[nx-1][j][k] = Hy[0][j][k];
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Periodic boundary condition for Hz
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::PBC_Hz()
{
	long i, j, k;
	//Impose the periodic boundary condition
	// i = nx-1;
	#pragma omp parallel for private(j,k) schedule(dynamic,nr_threads)
	for (j = 0; j < ny; j++)
	{	
		for (k = 1; k < nz-1; k++) //PEC at k = 0, k = nz-1
		{
			Hz[nx-1][j][k] = Hz[0][j][k];
		}
	}

	//j = ny-1
	#pragma omp parallel for private(i,k) schedule(dynamic,nr_threads)
	for (i = 0; i < nx; i++)
	{	
		for (k = 1; k < nz-1; k++) //PEC at k = 0, k = nz-1
		{
			Hz[i][ny-1][k] = Hz[i][0][k];
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Collect the function calls for the FDTD algorithm
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_xyPBC_zPML::fdtd_marching(long n_max_steep)
{
	double time = 0;
	long iter, it, num_iter_0 = 0;

	#ifndef run_enviroment //IBM AIX
		//load the workspace data 
		char Path_Data[512];
		
		if (load_workspace_data == 1)
		{
			//iter
			strcpy(Path_Data,Path_Load_Workspace_Data);
			strcat(Path_Data,"/Data_Workspace/iter_0.dat");
			if(load_1D_long(&num_iter_0, 1, Path_Data))
			{
				cout << "Error loading iter_0.dat workspace data" << endl;
				return 1;
			} 
			cout << "num_iter_0 = " << num_iter_0 << endl;
			
			if (jel_fourier_transf_vol == 1)
			{
				sprintf(Path_Data,"%s/Ex_fourier_TR_re_%d.dat",Path_Load_Workspace_Data,num_iter_0);
				load_2D(Ex_fourier_TR_re,n_f_Points,n_frec,Path_Data);
				sprintf(Path_Data,"%s/Ex_fourier_TR_im_%d.dat",Path_Load_Workspace_Data,num_iter_0);
				load_2D(Ex_fourier_TR_im,n_f_Points,n_frec,Path_Data);

				sprintf(Path_Data,"%s/Ey_fourier_TR_re_%d.dat",Path_Load_Workspace_Data,num_iter_0);
				load_2D(Ey_fourier_TR_re,n_f_Points,n_frec,Path_Data);
				sprintf(Path_Data,"%s/Ey_fourier_TR_im_%d.dat",Path_Load_Workspace_Data,num_iter_0);
				load_2D(Ey_fourier_TR_im,n_f_Points,n_frec,Path_Data);

				sprintf(Path_Data,"%s/Ez_fourier_TR_re_%d.dat",Path_Load_Workspace_Data,num_iter_0);
				load_2D(Ez_fourier_TR_re,n_f_Points,n_frec,Path_Data);
				sprintf(Path_Data,"%s/Ez_fourier_TR_im_%d.dat",Path_Load_Workspace_Data,num_iter_0);
				load_2D(Ez_fourier_TR_im,n_f_Points,n_frec,Path_Data);

				sprintf(Path_Data,"%s/Hx_fourier_TR_re_%d.dat",Path_Load_Workspace_Data,num_iter_0);
				load_2D(Hx_fourier_TR_re,n_f_Points,n_frec,Path_Data);
				sprintf(Path_Data,"%s/Hx_fourier_TR_im_%d.dat",Path_Load_Workspace_Data,num_iter_0);
				load_2D(Hx_fourier_TR_im,n_f_Points,n_frec,Path_Data);

				sprintf(Path_Data,"%s/Hy_fourier_TR_re_%d.dat",Path_Load_Workspace_Data,num_iter_0);
				load_2D(Hy_fourier_TR_re,n_f_Points,n_frec,Path_Data);
				sprintf(Path_Data,"%s/Hy_fourier_TR_im_%d.dat",Path_Load_Workspace_Data,num_iter_0);
				load_2D(Hy_fourier_TR_im,n_f_Points,n_frec,Path_Data);

				sprintf(Path_Data,"%s/Hz_fourier_TR_re_%d.dat",Path_Load_Workspace_Data,num_iter_0);
				load_2D(Hz_fourier_TR_re,n_f_Points,n_frec,Path_Data);
				sprintf(Path_Data,"%s/Hz_fourier_TR_im_%d.dat",Path_Load_Workspace_Data,num_iter_0);
				load_2D(Hz_fourier_TR_im,n_f_Points,n_frec,Path_Data);
			}

			strcpy(Path_Data,Path_Load_Workspace_Data);
			strcat(Path_Data,"/Data_Workspace");
			cout << "Load workspace data from: " << Path_Data << endl;

			if(load_FE_BH(Path_Data))
			{
				cout << "Error loading FDTD_3D workspace data" << endl;
				return 1;
			}

			if (jel_plane_wave ==1)
			{
				fdtd_1D.Load_FDTD_1D_Workspace(Path_Data);
			}

			time = num_iter_0*dt;
			num_iter_0++;
		}
	#endif

	for (iter = num_iter_0; iter < n_max_steep; iter++)
	{
		time += dt;
		
		//FDTD update for the electric field components
		calc_Ex();
		//Total field scattered field formulation
		if (jel_plane_wave ==1) 
		{
			//1D FDTD updates
			fdtd_1D.calc_Hy();
			TS_Ex();
		}
		calc_Ey();
		calc_Ez();

		//Periodic boundary conditions
		PBC_Ex();
		PBC_Ey();
		PBC_Ez();
		
		//FDTD update for the magnetic field components
		calc_Hx();
		calc_Hy();
		//Total field scattered field formulation
		if (jel_plane_wave ==1) 
		{
			fdtd_1D.calc_Ex(time);
			TS_Hy();
		}
		calc_Hz();

		//Periodic boundary conditions
		PBC_Hx();
		PBC_Hy();
		PBC_Hz();

		//point source
		if (jel_plane_wave == 0)
		{
			ptSource_J(Hz, Par_ptSource, time);
		}
		
		//the followed field components
	    Set_Data_Followed(iter);

		////////////////////////////////////////////////////////////////////
		//Calculate the averaged fields
		////////////////////////////////////////////////////////////////////
		if (aver_field_volume == 1)
		{
			aver_xyz_WignerC(Ex, nx_a_av, nx_b_av, ny_a_av, ny_b_av, nz_a_av, nz_b_av,
				             &Ex_Foll_av[iter], Vol_av);
			aver_xyz_WignerC(Ey, nx_a_av, nx_b_av, ny_a_av, ny_b_av, nz_a_av, nz_b_av,
				             &Ey_Foll_av[iter], Vol_av);
			aver_xyz_WignerC(Ez, nx_a_av, nx_b_av, ny_a_av, ny_b_av, nz_a_av, nz_b_av,
				             &Ez_Foll_av[iter], Vol_av);
			aver_xyz_WignerC(Hx, nx_a_av, nx_b_av, ny_a_av, ny_b_av, nz_a_av, nz_b_av,
				             &Hx_Foll_av[iter], Vol_av);
			aver_xyz_WignerC(Hy, nx_a_av, nx_b_av, ny_a_av, ny_b_av, nz_a_av, nz_b_av,
				             &Hy_Foll_av[iter], Vol_av);
			aver_xyz_WignerC(Hz, nx_a_av, nx_b_av, ny_a_av, ny_b_av, nz_a_av, nz_b_av,
				             &Hz_Foll_av[iter], Vol_av);
		}
		
		////////////////////////////////////////////////////////////////////
		//Calculate the integrals for Bragg diffraction
		////////////////////////////////////////////////////////////////////
		if (jel_aver_Bragg_WignerUC == 1)
		{
			for (it = 0; it <n_Gxy; it++)
			{
				aver_Bragg_WignerC(Ex, kx, ky, Gxy[it][0], Gxy[it][1], &Ex_Br_av_re[it][iter], 
								   &Ex_Br_av_im[it][iter], Surf_Bragg);
				aver_Bragg_WignerC(Ey, kx, ky, Gxy[it][0], Gxy[it][1], &Ey_Br_av_re[it][iter], 
								   &Ey_Br_av_im[it][iter], Surf_Bragg);
				aver_Bragg_WignerC(Ez, kx, ky, Gxy[it][0], Gxy[it][1], &Ez_Br_av_re[it][iter], 
								   &Ez_Br_av_im[it][iter], Surf_Bragg);
				aver_Bragg_WignerC(Hx, kx, ky, Gxy[it][0], Gxy[it][1], &Hx_Br_av_re[it][iter], 
								   &Hx_Br_av_im[it][iter], Surf_Bragg);
				aver_Bragg_WignerC(Hy, kx, ky, Gxy[it][0], Gxy[it][1], &Hy_Br_av_re[it][iter], 
								   &Hy_Br_av_im[it][iter], Surf_Bragg);
				aver_Bragg_WignerC(Hz, kx, ky, Gxy[it][0], Gxy[it][1], &Hz_Br_av_re[it][iter], 
								   &Hz_Br_av_im[it][iter], Surf_Bragg);
			}
		}

		////////////////////////////////////////////////////////////////////
		//Calculate the Fourier transform of the field components in the subvolume
		////////////////////////////////////////////////////////////////////
		if (jel_fourier_transf_vol == 1 && iter >= start_fourier)
		{
			fourier_transf_vol(Ex, nx_a_f, nx_b_f, ny_a_f, ny_b_f, nz_a_f, nz_b_f,
				               Ex_fourier_TR_re, Ex_fourier_TR_im, fourier_omega, 
						       n_frec, time);
			fourier_transf_vol(Ey, nx_a_f, nx_b_f, ny_a_f, ny_b_f, nz_a_f, nz_b_f,
				               Ey_fourier_TR_re, Ey_fourier_TR_im, fourier_omega, 
						       n_frec, time);
			fourier_transf_vol(Ez, nx_a_f, nx_b_f, ny_a_f, ny_b_f, nz_a_f, nz_b_f,
				               Ez_fourier_TR_re, Ez_fourier_TR_im, fourier_omega, 
						       n_frec, time);

			fourier_transf_vol(Hx, nx_a_f, nx_b_f, ny_a_f, ny_b_f, nz_a_f, nz_b_f,
				               Hx_fourier_TR_re, Hx_fourier_TR_im, fourier_omega, 
						       n_frec, time);
			fourier_transf_vol(Hy, nx_a_f, nx_b_f, ny_a_f, ny_b_f, nz_a_f, nz_b_f,
				               Hy_fourier_TR_re, Hy_fourier_TR_im, fourier_omega, 
						       n_frec, time);
			fourier_transf_vol(Hz, nx_a_f, nx_b_f, ny_a_f, ny_b_f, nz_a_f, nz_b_f,
				               Hz_fourier_TR_re, Hz_fourier_TR_im, fourier_omega, 
						       n_frec, time);
		}

		if (jel_Save_Slice == 1 && iter%nr_save == 1 && iter>=saveFROMinst && iter<=saveTOinst)
		{
			//Save_Field(iter);
			Save_FieldSlices(iter);
		}
	
		if (iter%100 == 1)
		{
			cout<< iter << endl;
			#ifndef run_enviroment //IBM AIX
				//check the elapsed time
				gettimeofday(&tv, &tz);
				el_time = (double) tv.tv_sec + (double) tv.tv_usec / 1000000.0 - time_start;
				//save the workspace data
				if (el_time > limit_time) 
				{
					mode_t Mode = S_IRWXU;

					strcpy(Path_Data,Path_Save_Workspace_Data);
					strcat(Path_Data,"/Data_Workspace");
					mkdir(Path_Data,Mode);//for UNIX-AIX
					
					save_FE_BH(iter,time,Path_Data);
					
					if (jel_plane_wave ==1)
					{
						fdtd_1D.Save_FDTD_1D_Workspace(Path_Data);
					}

					return iter;
				}
			#endif
		}
			
	}
	return iter; //the number of executed steeps
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Gaussian plane wave
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::Init_plWave_Gauss(double E0, double t0, double tw, long n_z_TS)
{
	double phi = 0;
	
	double eps_r = Mat[0][0];
	double mu_r  = Mat[0][1];
	
	n_1D = nz - n_z_TS;

	fdtd_1D.Init(n_1D, dt, dz, eps_r, mu_r);
	
	fdtd_1D.Init_Gauss_Source(E0, t0, tw);

	fdtd_1D.Get_Data(Ex_1D,Hy_1D);

	fdtd_1D.Init_nr_THR(nr_threads);
	
	n_TS_z = n_z_TS;

	jel_plane_wave = 1;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init sinusoidal plane wave
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::Init_plWave_Sin(double E0, double omega, long n_z_TS)
{
	double phi = 0;

	double eps_r = Mat[0][0];
	double mu_r  = Mat[0][1];
	
	n_1D = nz - n_z_TS;

	fdtd_1D.Init(n_1D, dt, dz, eps_r, mu_r);
	
	fdtd_1D.Init_Sin_Source(E0, omega, phi);
	
	fdtd_1D.Get_Data(Ex_1D,Hy_1D);
	
	fdtd_1D.Init_nr_THR(nr_threads);

	n_TS_z = n_z_TS;

	jel_plane_wave = 1;

}

///////////////////////////////////////////////////////////////////////////////////////
//Init sinusoidal modulated Gaussian plane wave
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::Init_plWave_GaussSin(double E0, double omega, double t0, 
										    double tw, long n_z_TS)
{
	double phi = 0;
	
	double eps_r = Mat[0][0];
	double mu_r  = Mat[0][1];
	
	n_1D = nz - n_z_TS;

	fdtd_1D.Init(n_1D, dt, dz, eps_r, mu_r);
	
	fdtd_1D.Init_GaussSin(E0, omega, phi, t0, tw);
	
	fdtd_1D.Get_Data(Ex_1D,Hy_1D);

	fdtd_1D.Init_nr_THR(nr_threads);
	
	n_TS_z = n_z_TS;

	jel_plane_wave = 1;
}

///////////////////////////////////////////////////////////////////////////////////////
//Total field scattered field interface in z direction - Ex update
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::TS_Ex()
{
	//Total field scattered field interface in z direction
	long i, j;
	double eps_r, K_1;

	K_1 = dt/(eps_0*dz);

	#pragma omp parallel for default(shared) private(i,j,eps_r) firstprivate(K_1) schedule(dynamic,nr_threads)
	for (i = 0; i<nx; i++) 
	{
		for (j = 0; j<ny; j++) 
		{
			eps_r = Mat[Ind[i][j][n_TS_z]][0];

			Ex[i][j][n_TS_z] += K_1*Hy_1D[1]/eps_r;
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Total field scattered field interface in z direction - Hy update
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::TS_Hy()
{
	long i, j;
	double mu_r, K_1;

	K_1 = dt/(mu_0*dz);

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

			Hy[i][j][n_TS_z-1] += K_1*Ex_1D[2]/mu_r;
		}
	}