fdtd_3d_xypbc_zpml.cpp

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

						
			strcpy(Path,Path_Data);
			strcat(Path,"/data_Hy");
			mkdir(Path,Mode);//for UNIX-AIX
						
			strcpy(Path,Path_Data);
			strcat(Path,"/data_Hz");
			mkdir(Path,Mode);//for UNIX-AIX
		#endif

		//middle planes of the computational space - for data save
		n_x_yz = (long) (n_x/2);
		n_y_xz = (long) (n_y/2); 
		n_z_xy = (long) (n_z/2);	

		saveFROMinst =  1; 
		saveTOinst   = 10000;

		nr_Save = 10;

		if (!FDTD_3D.Init_Save_FieldSlices(n_x_a, n_x_b, n_y_a, n_y_b, n_z_a, n_z_b, 
				 						   n_z_xy, n_y_xz, n_x_yz, saveFROMinst, saveTOinst,
										   nr_Save, Path_Data))
		{
			cout << "Memory allocation problem - Init_Save_FieldSlices - FDTD" << endl;
			return 1;
		}
	}

	//The followed field components
	long n_Ind_F = 60;
	long **Ind_F = NULL;

	Ind_F = Init_Matrix_2D<long>(n_Ind_F,3);
	if(!Ind_F)
	{
		cout << "Memory allocation problem - Ind_F" << endl;
		return -1;
	}

	long step_x = (n_x-1)/19;
	for (i = 0; i<20; i++)
	{
		Ind_F[i][0] = 1 + i*step_x; Ind_F[i][1] = n_y/2; Ind_F[i][2] = n_z/2;
	}
	long step_y = (n_y-1)/19;
	for (i = 20; i<40; i++)
	{
		Ind_F[i][0] = n_x/2; Ind_F[i][1] = 1 + (i-20)*step_y; Ind_F[i][2] = n_z/2;
	}
	long step_z = (n_z-2*n_PML-2*n_TS_z-4)/19;
	for (i = 40; i<60; i++)
	{
		Ind_F[i][0] = n_x/2; Ind_F[i][1] = n_y/2; Ind_F[i][2] = n_PML + n_TS_z + 1 + (i-40)*step_z;
	}

	//Save Ind_F
	strcpy(file_name,Path_Data);
	strcat(file_name,"/Ind_Foll");
	long x1 =0, y1 = 0, y2 = 3, it = 0;
	if ( save_2D_long(Ind_F,x1,n_Ind_F,y1,y2,it,file_name) )
		return 2;
		
	FDTD_3D.Init_Followed(Ind_F, n_Ind_F, num_iter);

    /////////////////////////////////////////////////////////////////////////////
	//Save the followed data
	/////////////////////////////////////////////////////////////////////////////
	double **hx_foll = NULL, **hy_foll = NULL, **hz_foll = NULL;
	double **ex_foll = NULL, **ey_foll = NULL, **ez_foll = NULL; 

	FDTD_3D.Get_Data_Followed(hx_foll, hy_foll, hz_foll, 
							  ex_foll, ey_foll, ez_foll);

	long ***Mask_Wigner = NULL;
	//load the geometry file of the Wigner cell
	if (Inp_D->aver_field_volume == 1 || Inp_D->jel_aver_Bragg_WignerUC == 1)
	{
		Mask_Wigner = Init_Matrix_3D<long>(n_x,n_y,1);
		if (!Mask_Wigner)
		{
			cout << "Memory allocation problem - **Mask_Wigner" << endl;
			return 1;
		}

		switch ( load_3D_Geom_long(Mask_Wigner, n_x, n_y, 1, Inp_D->path_name_Wigner_Cell) )
		{
			case 0:
					printf("Mask_Wigner file loaded succesfull\n");
					break;
			case 1:
					printf("faild to open the Mask_Wigner file \n");
					return 2;
			case 2:
					printf("wrong Mask_Wigner file content\n");
					return 3;
		}
		FDTD_3D.Init_Mask_Wigner(Mask_Wigner,n_x,n_y,1);
	}

	/////////////////////////////////////////////////////////////////////////////
	//Compute the averaged field components
	/////////////////////////////////////////////////////////////////////////////
	double *ex_foll_av = NULL, *ey_foll_av = NULL, *ez_foll_av = NULL;
	double *hx_foll_av = NULL, *hy_foll_av = NULL, *hz_foll_av = NULL;
	if (Inp_D->aver_field_volume == 1)
	{
		FDTD_3D.Init_Aver_xyz(nx_a_av, nx_b_av, ny_a_av, 
			                  ny_b_av, nz_a_av, nz_b_av, num_iter);

		FDTD_3D.Get_Data_Aver(hx_foll_av, hy_foll_av, hz_foll_av, 
							  ex_foll_av, ey_foll_av, ez_foll_av);
	}

	/////////////////////////////////////////////////////////////////////////////
	//Bragg diffraction -- square primitive cell
	/////////////////////////////////////////////////////////////////////////////
	double **ex_avBr_re = NULL, **ey_avBr_re = NULL, **ez_avBr_re = NULL;
	double **hx_avBr_re = NULL, **hy_avBr_re = NULL, **hz_avBr_re = NULL;
	double **ex_avBr_im = NULL, **ey_avBr_im = NULL, **ez_avBr_im = NULL;
	double **hx_avBr_im = NULL, **hy_avBr_im = NULL, **hz_avBr_im = NULL;
	double **Gxy = NULL;
	long n_Gxy = Inp_D->n_Gxy;

	if (Inp_D->jel_aver_Bragg_WignerUC == 1)
	{
		Gxy = Init_Matrix_2D<double>(n_Gxy,2);
		if (!Gxy)
		{
			cout << "Memory allocation problem - **Gxy" << endl;
			return 1;
		}

		switch (load_2D(Gxy, n_Gxy, 2, Inp_D->path_name_Gxy))
		{
			case 0:
					printf("Gxy file loaded succesfull\n");
					break;
			case 1:
					printf("faild to open the Gxy file \n");
					return 2;
			case 2:
					printf("wrong Gxy file content\n");
					return 3;
			case 3:
					printf("wrong Gxy file dimension\n");
					return 3;
		}

		FDTD_3D.Init_aver_Bragg_WignerC(Inp_D->nz_Bragg, 0.0, 0.0, 
			                            Gxy, n_Gxy, num_iter);

		FDTD_3D.Get_Data_Bragg_Aver(ex_avBr_re, ex_avBr_im, ey_avBr_re, ey_avBr_im, 
			                        ez_avBr_re, ez_avBr_im, hx_avBr_re, hx_avBr_im, 
							        hy_avBr_re, hy_avBr_im, hz_avBr_re, hz_avBr_im);
	}
	
	/////////////////////////////////////////////////////////////////////////////
	//Fourier transform in a subvolume
	/////////////////////////////////////////////////////////////////////////////
	double **Ex_fourier_TR_re = NULL, **Ex_fourier_TR_im = NULL; 
	double **Ey_fourier_TR_re = NULL, **Ey_fourier_TR_im = NULL; 
	double **Ez_fourier_TR_re = NULL, **Ez_fourier_TR_im = NULL; 
	double **Hx_fourier_TR_re = NULL, **Hx_fourier_TR_im = NULL; 
	double **Hy_fourier_TR_re = NULL, **Hy_fourier_TR_im = NULL; 
	double **Hz_fourier_TR_re = NULL, **Hz_fourier_TR_im = NULL; 
	double *fourier_omega = NULL;

	if (Inp_D->fourier_transf_vol == 1)
	{
		FDTD_3D.Init_Fourier_Transf(nx_a_f, nx_b_f, ny_a_f, ny_b_f, nz_a_f, nz_b_f,
			                        frec_1, frec_2, n_frec, start_fourier);

		FDTD_3D.Get_Data_Fourier(Ex_fourier_TR_re, Ex_fourier_TR_im, Ey_fourier_TR_re, Ey_fourier_TR_im,
					             Ez_fourier_TR_re, Ez_fourier_TR_im, Hx_fourier_TR_re, Hx_fourier_TR_im,
								 Hy_fourier_TR_re, Hy_fourier_TR_im, Hz_fourier_TR_re, Hz_fourier_TR_im,
								 fourier_omega);
	}

	/////////////////////////////////////////////////////////////////////////////
	//start of the main FDTD loop
	/////////////////////////////////////////////////////////////////////////////
	#ifndef run_enviroment //IBM AIX 
		FDTD_3D.Init_Mesure_Ellapsed_Time(time_start, Inp_D->limit_time, Inp_D->load_workspace_data,
			                              Inp_D->Path_Load_Workspace_Data, Path_Data);
	#endif

	int iter = FDTD_3D.fdtd_marching(num_iter); //returns the number of FDTD steps

	/////////////////////////////////////////////////////////////////////////////
	//end of the main FDTD loop
	/////////////////////////////////////////////////////////////////////////////
	if (iter > 1) //there is no error
	{
		//hx_foll
		strcpy(file_name,Path_Data);
		strcat(file_name,"/Hx_foll");
		if ( save_2D(hx_foll,x1,n_Ind_F,y1,iter,iter,file_name) )
		//if ( !save_2D_binary(hx_foll,n_Ind_F,iter,iter,file_name) )
			return 2;

		//hy_foll
		strcpy(file_name,Path_Data);
		strcat(file_name,"/Hy_foll");
		if ( save_2D(hy_foll,x1,n_Ind_F,y1,iter,iter,file_name) )
			return 2;

		//hz_foll
		strcpy(file_name,Path_Data);
		strcat(file_name,"/Hz_foll");
		if ( save_2D(hz_foll,x1,n_Ind_F,y1,iter,iter,file_name) )
			return 2;

		//ex_foll
		strcpy(file_name,Path_Data);
		strcat(file_name,"/Ex_foll");
		if ( save_2D(ex_foll,x1,n_Ind_F,y1,iter,iter,file_name) )
			return 2;

		//ey_foll
		strcpy(file_name,Path_Data);
		strcat(file_name,"/Ey_foll");
		if ( save_2D(ey_foll,x1,n_Ind_F,y1,iter,iter,file_name) )
			return 2;

		//ez_foll
		strcpy(file_name,Path_Data);
		strcat(file_name,"/Ez_foll");
		if ( save_2D(ez_foll,x1,n_Ind_F,y1,iter,iter,file_name) )
			return 2;

		if (Inp_D->aver_field_volume == 1)
		{
			//Hx_foll_av
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hx_foll_av");
			if ( save_1D(hx_foll_av,0,iter,iter,file_name) )
				return 2;

			//Hy_foll_av
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hy_foll_av");
			if ( save_1D(hy_foll_av,0,iter,iter,file_name) )
				return 2;

			//Hz_foll_av
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hz_foll_av");
			if ( save_1D(hz_foll_av,0,iter,iter,file_name) )
				return 2;

			//Ex_foll_av
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ex_foll_av");
			if ( save_1D(ex_foll_av,0,iter,iter,file_name) )
				return 2;

			//Ey_foll_av
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ey_foll_av");
			if ( save_1D(ey_foll_av,0,iter,iter,file_name) )
				return 2;

			//Ez_foll_av
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ez_foll_av");
			if ( save_1D(ez_foll_av,0,iter,iter,file_name) )
				return 2;
		}

		//save -- Bragg diffraction
		if (Inp_D->jel_aver_Bragg_WignerUC == 1)
		{
			//ex_avBr_re
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ex_avBr_re");
			if ( save_2D(ex_avBr_re,0,n_Gxy,0,iter,iter,file_name) )
				return 2;
			//ex_avBr_im
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ex_avBr_im");
			if ( save_2D(ex_avBr_im,0,n_Gxy,0,iter,iter,file_name) )
				return 2;
			//ey_avBr_re
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ey_avBr_re");
			if ( save_2D(ey_avBr_re,0,n_Gxy,0,iter,iter,file_name) )
				return 2;
			//ey_avBr_im
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ey_avBr_im");
			if ( save_2D(ey_avBr_im,0,n_Gxy,0,iter,iter,file_name) )
				return 2;
			//ez_avBr_re
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ez_avBr_re");
			if ( save_2D(ez_avBr_re,0,n_Gxy,0,iter,iter,file_name) )
				return 2;
			//ez_avBr_im
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ez_avBr_im");
			if ( save_2D(ez_avBr_im,0,n_Gxy,0,iter,iter,file_name) )
				return 2;

			//hx_avBr_re
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hx_avBr_re");
			if ( save_2D(hx_avBr_re,0,n_Gxy,0,iter,iter,file_name) )
				return 2;
			//hx_avBr_im
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hx_avBr_im");
			if ( save_2D(hx_avBr_im,0,n_Gxy,0,iter,iter,file_name) )
				return 2;
			//hy_avBr_re
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hy_avBr_re");
			if ( save_2D(hy_avBr_re,0,n_Gxy,0,iter,iter,file_name) )
				return 2;
			//hy_avBr_im
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hy_avBr_im");
			if ( save_2D(hy_avBr_im,0,n_Gxy,0,iter,iter,file_name) )
				return 2;
			//hz_avBr_re
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hz_avBr_re");
			if ( save_2D(hz_avBr_re,0,n_Gxy,0,iter,iter,file_name) )
				return 2;
			//hz_avBr_im
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hz_avBr_im");
			if ( save_2D(hz_avBr_im,0,n_Gxy,0,iter,iter,file_name) )
				return 2;
		}

		if (Inp_D->fourier_transf_vol == 1)
		{
			//Ex_fourier_TR_re
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ex_fourier_TR_re");
			if ( save_2D(Ex_fourier_TR_re,0,n_f_Points,0,n_frec,iter,file_name) )
				return 2;
			//Ex_fourier_TR_im
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ex_fourier_TR_im");
			if ( save_2D(Ex_fourier_TR_im,0,n_f_Points,0,n_frec,iter,file_name) )
				return 2;

			//Ey_fourier_TR_re
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ey_fourier_TR_re");
			if ( save_2D(Ey_fourier_TR_re,0,n_f_Points,0,n_frec,iter,file_name) )
				return 2;
			//Ey_fourier_TR_im
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ey_fourier_TR_im");
			if ( save_2D(Ey_fourier_TR_im,0,n_f_Points,0,n_frec,iter,file_name) )
				return 2;

			//Ez_fourier_TR_re
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ez_fourier_TR_re");
			if ( save_2D(Ez_fourier_TR_re,0,n_f_Points,0,n_frec,iter,file_name) )
				return 2;
			//Ez_fourier_TR_im
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Ez_fourier_TR_im");
			if ( save_2D(Ez_fourier_TR_im,0,n_f_Points,0,n_frec,iter,file_name) )
				return 2;

			//Hx_fourier_TR_re
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hx_fourier_TR_re");
			if ( save_2D(Hx_fourier_TR_re,0,n_f_Points,0,n_frec,iter,file_name) )
				return 2;
			//Hx_fourier_TR_im
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hx_fourier_TR_im");
			if ( save_2D(Hx_fourier_TR_im,0,n_f_Points,0,n_frec,iter,file_name) )
				return 2;

			//Hy_fourier_TR_re
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hy_fourier_TR_re");
			if ( save_2D(Hy_fourier_TR_re,0,n_f_Points,0,n_frec,iter,file_name) )
				return 2;
			//Hy_fourier_TR_im
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hy_fourier_TR_im");
			if ( save_2D(Hy_fourier_TR_im,0,n_f_Points,0,n_frec,iter,file_name) )
				return 2;

			//Hz_fourier_TR_re
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hz_fourier_TR_re");
			if ( save_2D(Hz_fourier_TR_re,0,n_f_Points,0,n_frec,iter,file_name) )
				return 2;
			//Hz_fourier_TR_im
			strcpy(file_name,Path_Data);
			strcat(file_name,"/Hz_fourier_TR_im");
			if ( save_2D(Hz_fourier_TR_im,0,n_f_Points,0,n_frec,iter,file_name) )
				return 2;

			//fourier_omega
			strcpy(file_name,Path_Data);
			strcat(file_name,"/fourier_omega");
			if ( save_1D(fourier_omega,0,n_frec,iter,file_name) )
				return 2;
		}
	}
	else
	{
		cout << "Error in fdtd run" << endl;
	}

	
	/////////////////////////////////////////////////////////////////////////////
	//Free the allocated memory
	/////////////////////////////////////////////////////////////////////////////
	FDTD_3D.Free_Mem(); //Free the memory

	//Free the memory allocated in the main program
	if (Path)
	{
		free(Path);
		Path = NULL;
	}
	if (Path_Data)
	{
		free(Path_Data);
		Path_Data = NULL;
	}
	if (file_name)
	{
		free(file_name);
		file_name = NULL;
	}
	if (Ind_F)
			Ind_F = Free_Matrix_2D<long>(Ind_F);
	if (Materials)
			Materials = Free_Matrix_2D<double>(Materials);
	if (Index)
			Index = Free_Matrix_3D<long>(Index,n_x);
	if (Mask_Wigner)
		    Mask_Wigner = Free_Matrix_3D<long>(Mask_Wigner,n_x);

	//the elapsed time
	#ifdef run_enviroment //WIN
		finish = clock();
		duration = (double)(finish - start) / CLOCKS_PER_SEC;
		cout<< "duration: " << duration << " seconds" <<endl;
	#else //IBM AIX
		gettimeofday(&tv, &tz);
		time_end = (double) tv.tv_sec + (double) tv.tv_usec / 1000000.0;
		cout << "Ellapsed time = " << time_end-time_start << endl;
	#endif
	
	cout << "success" << endl;
	
	return 0;
}