fdtd_3d_xypbc_zpml.cpp

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

//The plane wave propagates in z direction parpendicular to the surface of the photonic crystal
//In this code in the Bloch condition kx and ky are always zero
//The fields are real

#include "stdafx.h"

#include "FDTD_xyPBC_zPML.h"
#include "read_input_file.h"

int main(int argc, char* argv[])
{
	char *file_name = NULL;
	file_name = (char *) calloc(512,sizeof(char));
	if (!file_name)
	{
		cout << "Memory allocation problem - *file_name" << endl;
		return 1;
	}
	
	if (argc > 1)
	{
		strcpy(file_name,argv[1]);
	}
	else
	{
		strcpy(file_name,"Input_Data.txt");
	}

	/////////////////////////////////////////////////////////////////////////////
	//to measure the ellapsed time
	/////////////////////////////////////////////////////////////////////////////
	#ifdef run_enviroment //WIN
		clock_t start, finish;
		double  duration;

		start = clock();
	#else //IBM AIX
		double time_start, time_end, el_time;
		struct timeval tv;
		struct timezone tz;

		gettimeofday(&tv, &tz);
		time_start = (double) tv.tv_sec + (double) tv.tv_usec / 1000000.0;
	#endif

	///////////////////////////////////////////////////////
	//initialization to read from file the input data
	///////////////////////////////////////////////////////
	Input_Data *Inp_D = NULL;
	Inp_D = (Input_Data *) calloc(1,sizeof(Input_Data));

	Inp_D->Path_Load_Workspace_Data = NULL;
	Inp_D->path_name_CoordPointSource = NULL;
	Inp_D->path_name_Index = NULL;
	Inp_D->path_name_MatParam = NULL;
	Inp_D->path_name_ParamPointSource = NULL;
	Inp_D->path_name_Gxy = NULL;
	Inp_D->path_name_Wigner_Cell = NULL;

	if (readFile_fdtd_data(Inp_D,file_name))
	{
		cout << "Error reading input data: " << file_name << endl;
		return 1;
	}
	else
	{
		cout << "Input file: " << file_name << " loaded succesfull" << endl;
	}
	///////////////////////////////////////////////////////

	/////////////////////////////////////////////////////////////////////////////
	//set up the OpenMP enviroment
	/////////////////////////////////////////////////////////////////////////////
	int num_threads = 1; 
	#ifdef _OPENMP
		omp_set_dynamic(1);
		num_threads = Inp_D->num_threads;
		omp_set_num_threads(num_threads);
	#endif

	double pi = 180.0*atan(1.0)/45.0;
	//permittivity of free space 
	double eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]
	//permeability of free space 
    double mu_0  = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6;  // [H/m]
	//the speed of the light in the vacuum
	double speed_c = 1.0/sqrt(eps_0*mu_0); // [m/s]

	//the parameters of the photonic crystal
		
	long n_x = Inp_D->nx; //total number of cells in x direction
	long n_y = Inp_D->ny; //total number of cells in y direction
	long n_z = Inp_D->nz;  //total number of cells in z direction
		
	long n_PML = Inp_D->n_PML; //[cells] the size of PML

	double dx = Inp_D->dx;
	double dy = Inp_D->dy;
	double dz = Inp_D->dz;
	
	double dt = ( 1.0/( speed_c*sqrt(1.0/(dx*dx)+1.0/(dy*dy)+1.0/(dz*dz)) ) )*0.9;
	cout << "dt = " << dt << endl;

    long num_iter = Inp_D->num_iter; //the number of FDTD iterations

	long nx_a_av, nx_b_av, ny_a_av, ny_b_av, nz_a_av, nz_b_av;
	//long Vol_av;
	if (Inp_D->aver_field_volume == 1)
	{
		nx_a_av = Inp_D->nx_a_av; nx_b_av = Inp_D->nx_b_av; 
		ny_a_av = Inp_D->ny_a_av; ny_b_av = Inp_D->ny_b_av; 
		nz_a_av = Inp_D->nz_a_av; nz_b_av = Inp_D->nz_b_av; 
	}

	/////////////////////////
	//Fourier transform in a subvolume
	/////////////////////////
	long nx_a_f, nx_b_f, ny_a_f, ny_b_f, nz_a_f, nz_b_f;
	double frec_1, frec_2;
	long n_frec, start_fourier, n_f_Points;
	if (Inp_D->fourier_transf_vol == 1) 
	{
		//compute the Fourier transform in this volume
		nx_a_f = Inp_D->nx_a_f; nx_b_f = Inp_D->nx_b_f; 
		ny_a_f = Inp_D->ny_a_f; ny_b_f = Inp_D->ny_b_f; 
		nz_a_f = Inp_D->nz_a_f; nz_b_f = Inp_D->nz_b_f; 

		n_f_Points = (nx_b_f - nx_a_f + 1)*(ny_b_f - ny_a_f + 1)*(nz_b_f - nz_a_f + 1);

		//between these frequencies
		frec_1 = Inp_D->frec_1;
		frec_2 = Inp_D->frec_2;
		n_frec = Inp_D->n_frec; //nr of frequencies
		start_fourier = Inp_D->start_fourier;
	}

	//the excitation
	long source_Type = Inp_D->source_type; // 1- Gauss; 2- Sin; 3- Gauss-Sin

	long jel_plane_wave = Inp_D->jel_plane_wave; //the type of the excitation 0 - point  source


	long n_Coord;
	long **Coord_ptSource = NULL;
	double **Param_Source = NULL;
	long **Mask = NULL;
	Mask = Init_Matrix_2D<long>(n_x,n_y);

	double E0, tw, t0, omega, phase;
	long n_TS, n_TS_z;
	
	long i;

	switch (jel_plane_wave)
	{
		case 0:
			{
				n_Coord = 1;

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

				/*for (i = -10; i<10; i++)
				{
					for (k = -3; k<3; k++)
					{
						Coord_ptSource[ii][0] = n_x/2+i; 
						Coord_ptSource[ii][1] = n_PML+5; 
						Coord_ptSource[ii][2] = n_z/2+k;
						ii++;
					}
				}*/
				Coord_ptSource[0][0] = 70;
				Coord_ptSource[0][1] = 5;
				Coord_ptSource[0][2] = n_z/2;
		

				if (source_Type == 1 || source_Type == 2)
				{
					Param_Source = Init_Matrix_2D<double>(n_Coord,3);
				}
				else
				{
					Param_Source = Init_Matrix_2D<double>(n_Coord,5);
				}
				if (!Param_Source)
				{
					cout << "Memory allocation problem - **Param_Source" << endl;
					return 1;
				}
				
				E0 = 100;
				tw = 350*dt;
				t0 = 4*tw;
				omega = 2*pi*0.25*speed_c;///a_lat;

				switch (source_Type)
				{
					case 1:
						for (i = 0; i<n_Coord; i++)
						{
							//Gaussian source
							Param_Source[i][0] = E0; //J0
							Param_Source[i][1] = tw; //tw
							Param_Source[i][2] = t0; //t0							
						}
						break;
					case 2:
						for (i = 0; i<n_Coord; i++)
						{
							//Sinusoidal source
							Param_Source[i][0] = E0;    //J0
							Param_Source[i][1] = omega; //omega
							Param_Source[i][2] = -pi/2; //phase
						}							
						break;
					case 3:
						for (i = 0; i<n_Coord; i++)
						{
							//Sinusoidal modulated Gaussian source
							Param_Source[i][0] = E0; //J0
							Param_Source[i][1] = tw; //tw
							Param_Source[i][2] = t0; //t0
							Param_Source[i][3] = omega; //omega
							Param_Source[i][4] = 0; //phase
						}
						break;
				}
			}		
			break;			   
		case 1: //plane wave excitation with total field - scattered field formulation
			
			E0 = Inp_D->X0; //amplituide
			switch ( Inp_D->source_type)
			{
				case 1: //Gaussian source
					tw = Inp_D->tw;
					t0 = Inp_D->t0;
					break;
				case 2: //Sin source
					omega      = Inp_D->omega;
					phase      = Inp_D->phase;
					break;
				case 3://Gaussian-Sin source
					tw    = Inp_D->tw;
					t0    = Inp_D->t0;
					omega = Inp_D->omega;
					phase = Inp_D->phase;
					break;
			}
			
			n_TS = Inp_D->n_TS; //the size of the scattered field zone

			n_TS_z = n_PML + n_TS;
	       
			break;
	}


	long n_x_a = 0;
	long n_x_b = n_x;
	long n_y_a = 0;
	long n_y_b = n_y;
	long n_z_a = 0; //n_PML;
	long n_z_b = n_z;// - n_PML;

	/////////////////////////////////////////////////////////////////////////////
	//Initialization of the geometry
	/////////////////////////////////////////////////////////////////////////////
	long ***Index = NULL;
	Index = Init_Matrix_3D<long>(n_x,n_y,n_z);
	if (!Index)
	{
		cout << "Memory allocation problem - ***Index" << endl;
		return 1;
	}

	//Reads from file the elements of the Geometry matrix 
    switch ( load_3D_Geom_long(Index, n_x, n_y, n_z, Inp_D->path_name_Index) )
	{
		case 0:
				printf("Geom file loaded succesfull\n");
				break;
		case 1:
				printf("faild to open the Geom file \n");
				return 2;
		case 2:
				printf("wrong Geom file content\n");
				return 3;
	}

	/////////////////////////////////////////////////////////////////////////////
	double **Materials = NULL;
	long n_Mat = Inp_D->n_Mat;
	Materials = Init_Matrix_2D<double>(n_Mat,2);
	if (!Materials)
	{
		//cout << "Memory allocation problem - **Mat" << endl;
		printf("Memory allocation problem - **Mat\n");
		return -1;
	}

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

	/////////////////////////////////////////////////////////////////////////////
	//Initializations of the FDTD algorithm
	/////////////////////////////////////////////////////////////////////////////
	CFDTD_xyPBC_zPML FDTD_3D;

	FDTD_3D.Init_nr_THR(num_threads);

	if(!FDTD_3D.Init(Index, n_x, n_y, n_z, Materials, n_Mat, n_PML, dt, dx, dy, dz))
	{
		cout << "Memory allocation problem - Init - FDTD" << endl;
		return 1;
	}

	FDTD_3D.Init_PML_Par(Inp_D->PML_eps_r_z_1, Inp_D->PML_mu_r_z_1,
		                 Inp_D->PML_eps_r_z_2, Inp_D->PML_mu_r_z_2);

	if (jel_plane_wave==0)
	{
		FDTD_3D.Init_ptSource(Coord_ptSource, n_Coord, Param_Source, source_Type);
	}
	else
	{
		switch (source_Type)
		{
			case 1:
				FDTD_3D.Init_plWave_Gauss(E0,t0,tw,n_TS_z);
				break;
			case 2:
				FDTD_3D.Init_plWave_Sin(E0, omega, n_TS_z);
				break;
			case 3:
				FDTD_3D.Init_plWave_GaussSin(E0, omega, t0, tw, n_TS_z);
				break;
		}
	}

	/////////////////////////////////////////////////////////////////////////////
	//Initializations for data save
	/////////////////////////////////////////////////////////////////////////////

	//Create the directory to save output data
	char *Path = NULL, *Path_Data = NULL;
	Path =(char *) calloc(512,sizeof(char));
	Path_Data =(char *) calloc(512,sizeof(char));
	if (!Path || !Path_Data)
	{
		cout << "Memory allocation problem - *Path or *Path_Data" << endl;
		return 1;
	}

	#ifdef run_enviroment 
		//for WIN
		GetCurrentDirectory(512, Path); 
		strcpy(Path_Data,Path);
		strcat(Path_Data,"//data");
		CreateDirectory(Path_Data,0);
	#else
		//for UNIX-AIX
		strcpy(Path,Inp_D->path_save_data);
		strcpy(Path_Data,Path);
		mode_t Mode = S_IRWXU;
		int cont = 1;
		while (mkdir(Path_Data, Mode)) //create the main directory
		{
			if(errno == EEXIST) //check if the directory name exists 
			{
				sprintf(Path_Data, "%s_%i", Path, cont);
			}
			else
			{
				cout << "Directory creation problem (Data)" << endl;
				return -1;
			}
			cont++;
		}
	#endif

	long n_x_yz, n_y_xz, n_z_xy;
	long saveFROMinst, saveTOinst;
	long nr_Save;
	
	int save_field_slices = Inp_D->save_field; //to save data - Ex, Ey, Ez, Hx, Hy, Hz slices

	if (save_field_slices == 1)
	{
		#ifdef run_enviroment 
			//for WIN
			if (jel_plane_wave == 1)
			{
				strcpy(Path,Path_Data);
				strcat(Path,"/data_Ex_1D");
				CreateDirectory(Path,0);//for WIN
			}

			strcpy(Path,Path_Data);
			strcat(Path,"/data_Ex");
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_Ey");
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_Ez");
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_Hx");
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_Hy");
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_Hz");
			CreateDirectory(Path,0);//for WIN
		#else
			if (jel_plane_wave == 1)
			{
				strcpy(Path,Path_Data);
				strcat(Path,"/data_Ex_1D");
				mkdir(Path,Mode);//for UNIX-AIX
			}

			//for UNIX-AIX
			strcpy(Path,Path_Data);
			strcat(Path,"/data_Ex");
			mkdir(Path,Mode);//for UNIX-AIX
									
			strcpy(Path,Path_Data);
			strcat(Path,"/data_Ey");
			mkdir(Path,Mode);//for UNIX-AIX
						
			strcpy(Path,Path_Data);
			strcat(Path,"/data_Ez");
			mkdir(Path,Mode);//for UNIX-AIX
						
			strcpy(Path,Path_Data);
			strcat(Path,"/data_Hx");
			mkdir(Path,Mode);//for UNIX-AIX