fdtd_3d_complex.cpp

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

///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_3D_COMPLEX::Init_ptSource(long **Coord_ptSource, long n_Coord, double **Param_Source, int source_Type,
									int Pt_Source_Ex, int Pt_Source_Ey, int Pt_Source_Ez,
	                                int Pt_Source_Hx, int Pt_Source_Hy, int Pt_Source_Hz)
{
	long nr_param;
	switch (source_Type)
	{
		case 1:
			nr_param = 3;
			Param_Source_re = Init_Matrix_2D<double>(n_Coord,3);
			Param_Source_im = Init_Matrix_2D<double>(n_Coord,3);
			break;
		case 2:
			nr_param = 4;
			Param_Source_re = Init_Matrix_2D<double>(n_Coord,4);
			Param_Source_im = Init_Matrix_2D<double>(n_Coord,4);
			break;
		case 3:
			nr_param = 5;
			Param_Source_re = Init_Matrix_2D<double>(n_Coord,5);
			Param_Source_im = Init_Matrix_2D<double>(n_Coord,5);
			break;
	}
	if (!Param_Source_re || !Param_Source_im)
	{
		cout << "Memory allocation problem - **Param_Source" << endl;
		return 1;
	}

	long i,j;
	for (i = 0; i<n_Coord; i++)
	{
		for (j = 0; j<nr_param; j++)
		{
			Param_Source_re[i][j] = Param_Source[i][j]; 
			Param_Source_im[i][j] = Param_Source[i][j]; 
		}
	}

	//the phase difference of the real and imaginary part is pi/2
	if(source_Type == 2)
	{
		for (i = 0; i<n_Coord; i++)
		{
			Param_Source_im[i][2] -= pi/2.0; //phase
		}
	}

	if(source_Type == 3)
	{
		for (i = 0; i<n_Coord; i++)
		{
			Param_Source_im[i][4] -= pi/2.0; //phase
		}
	}

	FDTD_3D_re.Init_ptSource(Coord_ptSource, n_Coord, Param_Source_re, source_Type,
		                     Pt_Source_Ex,  Pt_Source_Ey,  Pt_Source_Ez,
	                         Pt_Source_Hx,  Pt_Source_Hy,  Pt_Source_Hz);
	FDTD_3D_im.Init_ptSource(Coord_ptSource, n_Coord, Param_Source_im, source_Type,
		                     Pt_Source_Ex,  Pt_Source_Ey,  Pt_Source_Ez,
	                         Pt_Source_Hx,  Pt_Source_Hy,  Pt_Source_Hz);

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init for the Total Field - Scattered Field formulation
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_3D_COMPLEX::Init_TS(double Teta, double Phi, double Gamma, long N_TS,
				              long nPML_1D, int source_type, double X0, double t0,
				              double tw, double omega, double phase, double PML_eps_r_1D,
				              double PML_mu_r_1D, int *Jel_TS_planes)
{
	double phase_re = phase;
	double phase_im = phase - pi/2.0;

	FDTD_3D_re.Init_TS(Teta, Phi, Gamma, N_TS, nPML_1D, source_type, X0, t0,
				       tw, omega, phase_re, PML_eps_r_1D, PML_mu_r_1D, 
					   Jel_TS_planes);

	FDTD_3D_im.Init_TS(Teta, Phi, Gamma, N_TS, nPML_1D, source_type, X0, t0,
				       tw, omega, phase_im, PML_eps_r_1D, PML_mu_r_1D, 
					   Jel_TS_planes);
	
	return 0;
}


///////////////////////////////////////////////////////////////////////////////////////
//Init the data save -- write in file slices from the computational volume
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_3D_COMPLEX::Init_Save_FieldSlices(long nxa, long nxb, long nya, long nyb, long nza, long nzb, 
		                                    long nz_xy, long ny_xz, long nx_yz, long nr_Save,  long saveFROMinst,
							                long saveTOinst, char *path_data)
{
	char path[512];
	strcpy(path,path_data);
	strcat(path,"/data_EH_re");
	if (FDTD_3D_re.Init_Save_FieldSlices(nxa, nxb, nya, nyb, nza, nzb, 
				 				         nz_xy, ny_xz, nx_yz, nr_Save,  saveFROMinst,
							             saveTOinst, path))
	{
		cout << "Memory allocation problem - Init_Save_FieldSlices_re - FDTD" << endl;
		return 1;
	}

	strcpy(path,path_data);
	strcat(path,"/data_EH_im");
	if (FDTD_3D_im.Init_Save_FieldSlices(nxa, nxb, nya, nyb, nza, nzb, 
				 				         nz_xy, ny_xz, nx_yz, nr_Save,  saveFROMinst,
							             saveTOinst, path))
	{
		cout << "Memory allocation problem - Init_Save_FieldSlices_im - FDTD" << endl;
		return 1;
	}

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init the followed point
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_3D_COMPLEX::Init_Followed(long **Ind_Followed, long n_Ind_F, long num_iter)
{

	if ( FDTD_3D_re.Init_Followed(Ind_Followed, n_Ind_F, num_iter) )
	{
		cout << "Memory allocation problem - Init_Followed_re" << endl;
		return 1;
	}

	if ( FDTD_3D_im.Init_Followed(Ind_Followed, n_Ind_F, num_iter) )
	{
		cout << "Memory allocation problem - Init_Followed_im" << endl;
		return 1;
	}

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Get the pointers to the followed data
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_COMPLEX::Get_Data_Followed(double **&hx_foll_re, double **&hy_foll_re, double **&hz_foll_re, 
										 double **&ex_foll_re, double **&ey_foll_re, double **&ez_foll_re,
										 double **&hx_foll_im, double **&hy_foll_im, double **&hz_foll_im, 
										 double **&ex_foll_im, double **&ey_foll_im, double **&ez_foll_im)
{
	FDTD_3D_re.Get_Data_Followed(hx_foll_re, hy_foll_re, hz_foll_re, 
							     ex_foll_re, ey_foll_re, ez_foll_re);
	FDTD_3D_im.Get_Data_Followed(hx_foll_im, hy_foll_im, hz_foll_im, 
							     ex_foll_im, ey_foll_im, ez_foll_im);
}

///////////////////////////////////////////////////////////////////////////////////////
//The fdtd marching
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_3D_COMPLEX::fdtd_marching(long n_max_steep)
{
	FDTD_3D_re.Get_Electric_Fields(Ex_re,Ey_re,Ez_re);
	FDTD_3D_im.Get_Electric_Fields(Ex_im,Ey_im,Ez_im);

	long iter, num_iter_0 = 0;
	double time = 0;
	
	#ifndef run_enviroment //IBM AIX
		//load the workspace data 
		char Path_Data[512];
		cout << "load_workspace_data = " << load_workspace_data << endl;
		if (load_workspace_data == 1)
		{
			cout << "Path = " << Path_Load_Workspace_Data << endl;
			//iter
			strcpy(Path_Data,Path_Load_Workspace_Data);
			strcat(Path_Data,"/Data_Re/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;
			num_iter_0++;
			time = num_iter_0*dt;

			strcpy(Path_Data,Path_Load_Workspace_Data);
			strcat(Path_Data,"/Data_Re");
			if(FDTD_3D_re.load_FGE_BH(Path_Data))
			{
				cout << "Error loading FDTD_3D_re workspace data" << endl;
				return 1;
			}


			strcpy(Path_Data,Path_Load_Workspace_Data);
			strcat(Path_Data,"/Data_Im");
			if(FDTD_3D_im.load_FGE_BH(Path_Data))
			{
				cout << "Error loading FDTD_3D_im workspace data" << endl;
				return 1;
			}

		}
	#endif

	for (iter = num_iter_0; iter < n_max_steep; iter++)
	{
		time += dtDIV2;

		//update the real part of the electric field
		FDTD_3D_re.calc_E(time);
		//update the imaginary part of the electric field
		FDTD_3D_im.calc_E(time);

		//Bloch type boundary conditions in x direction
		if (nPML_x_1 == 1 && nPML_x_2 == 1)
		{
			bc_Bloch_Ey_iDir(Ey_re, Ey_im, cos_kxORa, sin_kxORa);
			bc_Bloch_Ez_iDir(Ez_re, Ez_im, cos_kxORa, sin_kxORa);
		}
		
		//Bloch type boundary conditions in y direction
		if (nPML_y_1 == 1 && nPML_y_2 == 1)
		{
			bc_Bloch_Ex_jDir(Ex_re, Ex_im, cos_kyORb, sin_kyORb);
			bc_Bloch_Ez_jDir(Ez_re, Ez_im, cos_kyORb, sin_kyORb);
		}
		
		//Bloch type boundary conditions in z direction
		if (nPML_z_1 == 1 && nPML_z_2 == 1)
		{
			bc_Bloch_Ex_kDir(Ex_re, Ex_im, cos_kzORc, sin_kzORc);
			bc_Bloch_Ey_kDir(Ey_re, Ey_im, cos_kzORc, sin_kzORc);
		}

		time += dtDIV2;

		//update the real part of the magnetic field
		FDTD_3D_re.calc_H(time,iter);
		//update the imaginary part of the magnetic field
		FDTD_3D_im.calc_H(time,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;

					//re
					strcpy(Path_Data,Path_Save_Workspace_Data);
					sprintf(Path_Data, "%s/%i", Path_Data, myrank);
					if ( mkdir(Path_Data,Mode) )
					{
						if ( errno != EEXIST ) //the directory was already created
						{
							cout << "Directory creation problem errno : " << errno << endl;
						}
					}
					strcat(Path_Data,"/Data_Re");
					cout << Path_Data << endl;
					if ( mkdir(Path_Data,Mode) )//for UNIX-AIX
					{
						cout << "Directory creation problem errno : " << errno << endl;
					}
					FDTD_3D_re.save_FGE_BH(iter,time,Path_Data);

					//im
					strcpy(Path_Data,Path_Save_Workspace_Data);
					sprintf(Path_Data, "%s/%i", Path_Data, myrank);
					if ( mkdir(Path_Data,Mode) )
					{
						if ( errno != EEXIST ) //the directory was already created
						{
							cout << "Directory creation problem errno : " << errno << endl;
						}
					}
					strcat(Path_Data,"/Data_Im");
					cout << Path_Data << endl;
					if ( mkdir(Path_Data,Mode) )//for UNIX-AIX
					{
						cout << "Directory creation problem errno : " << errno << endl;
					}
					FDTD_3D_im.save_FGE_BH(iter,time,Path_Data);

					return iter;
				}
			#endif
		}

	}

	return iter;
}

///////////////////////////////////////////////////////////////////////////////////////
//Reset the field components
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_COMPLEX::Reset_Field_Comp()
{
	FDTD_3D_re.Reset_Field_Comp();
	FDTD_3D_im.Reset_Field_Comp();
}

///////////////////////////////////////////////////////////////////////////////////////
//Init for the Bloch BC
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_COMPLEX::Init_KxORa_KyORb_KzORc(double kx, double ky, double kz)
{
	double lat_aORkx = (nx-1)*dx*kx;
	double lat_bORky = (ny-1)*dy*ky;
	double lat_cORkz = (nz-1)*dz*kz;

	//the wave vectors
	sin_kxORa = sin(lat_aORkx);
	cos_kxORa = cos(lat_aORkx);
	sin_kyORb = sin(lat_bORky);
	cos_kyORb = cos(lat_bORky);
	sin_kzORc = sin(lat_cORkz);
	cos_kzORc = cos(lat_cORkz);
}