fdtd_2d_tm_pml.cpp

2维时域有限差分模拟

// fdtd_2D_TM_PML.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "fdtd_2d_TM.h"

int main(int argc, char* argv[])
{
	#ifdef senddata
		Engine *ep;
		ep = engOpen(NULL);
	#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]

	int n_x = 5500; //total number of cells in x direction
	int n_y = 3000; //total number of cells in y direction
	int n_PML = 15;
		
	double a_lat = 1;

	double dx = 0.0075e-6;
	double dy = 0.0075e-6;
	
	double dt = ( 1.0/speed_c/sqrt( 1.0/(dx*dx) + 1.0/(dy*dy) ) )*0.9;

	int num_iter = 40000; //the number of FDTD iterations

	//the excitation
	int jel_plane_wave = 1; // 0 - point source, 1 -plane wave
	int source_Type = 2; // 1- Gauss; 2- Sin; 3- Gauss-Sin

	int n_Coord;
	int **Coord_ptSource = NULL;
	double **Param_Source = NULL;
	
	double E0, tw, t0, omega, phi;
	int n_TS_xa, n_TS_xb, n_TS_ya, n_TS_yb;
	int n_TS = 5;

	int i;

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

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

			Coord_ptSource[0][0] = n_x/2; 
			Coord_ptSource[0][1] = n_y/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 = 40*dt;
			t0 = 4*tw;
			omega = 2*pi*0.5*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 - scattered field formulation
			
			//Gaussian pulse
			E0 = 100;   //[V/m]
			tw = 60*dt; //[s]
			t0 = 4*tw;  //[s]
			phi = 0*pi/180; //the direction of propagation
			//Sin
			omega = 2.975752870946056e+015;//2.0*pi*0.5*speed_c/a_lat; //[rad/s]]
			//the position of the total field-scattered field region 
			n_TS_xa = n_PML + n_TS;
	        n_TS_xb = n_x - n_PML - n_TS;
	        n_TS_ya = n_PML + n_TS;
			n_TS_yb = n_y - n_PML  - n_TS;
	        break;
	}

				
	int n_x_a = 0;
	int n_x_b = n_x;
	int n_y_a = 0;
	int n_y_b = n_y;

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

	//Reads from file the elements of the Geometry matrix 
	char path_name[512] = "Geom_Lens_5500_3000.geom";
	switch ( load_2D_Geom_int(Index, n_x, n_y, path_name) )
	{
		case 0:
			    printf("Geom file loaded successfully \n");
			    break;
		case 1:
				printf("faild to open the Geom file \n");
				return 2;
		case 2:
				printf("wrong file content\n");
				return 3;
	    case 3:
				printf("wrong Geom file dimension\n");
				return 3;
		/*default:
		{
			char name_Index[] = "Index";
			int a1 = 0, b1 =0, c1 = 0;
			if ( !save_2D_int( Index,a1,n_x,b1,n_y,c1,name_Index) )
						return false;
		}*/
	}
	/////////////////////////////////////////////////////////////////////////////

	double **Materials = NULL;
	int n_Mat = 2;
	Materials = Init_Matrix_2D<double>(n_Mat,2);
	if (!Materials)
	{
		cout << "Memory allocation problem - **Materials" << endl;
		return 1;
	}

	Materials[0][0] = 1.0;  //eps_r_1;
	Materials[0][1] = 1.0;  //mu_r_1

	Materials[1][0] = 4.0; //eps_r_2;
	Materials[1][1] = 1.0;  //mu_r_2;


	/////////////////////////////////////////////////////////////////////////////
	//Initializations of the FDTD algorithm
	/////////////////////////////////////////////////////////////////////////////
	CFDTD_2D_TM FDTD_2D;
	#ifdef senddata
		FDTD_2D.set_MATLABengine(ep);
	#endif

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

	FDTD_2D.Init_PML_Param();

	if (jel_plane_wave == 0)
	{
		FDTD_2D.Init_ptSource(Coord_ptSource, n_Coord, Param_Source, source_Type);
	}
	else
	{
		switch (source_Type)
		{
			case 1:
				if (FDTD_2D.Init_plWave_Gauss(E0, t0, tw, phi, n_TS_xa, n_TS_xb,
						                      n_TS_ya, n_TS_yb))
				{
					cout << "memory allocation problem" << endl;
					return 1;
				}
				break;
			case 2:
				if (FDTD_2D.Init_plWave_Sin(E0, omega, phi, n_TS_xa, n_TS_xb,
						                    n_TS_ya, n_TS_yb))
				{
					cout << "memory allocation problem" << endl;
					return 1;
				}
				break;
			case 3:
				if (FDTD_2D.Init_plWave_GaussSin(E0, omega, t0, tw, phi, n_TS_xa,
						                         n_TS_xb, n_TS_ya, n_TS_yb))
				{
					cout << "memory allocation problem" << endl;
					return 1;
				}
				break;
		}
	}
	
	/////////////////////////////////////////////////////////////////////////////
	//Initializations for data save
	/////////////////////////////////////////////////////////////////////////////

	//Create the directory to save output data
	char *Path = NULL;
	Path =(char *) calloc(512,sizeof(char));
	if (!Path)
	{
		cout << "Memory allocation problem - *Path" << endl;
		return 1;
	}
	GetCurrentDirectory(512, Path); 
	strcat(Path,"\\data");
	CreateDirectory(Path,0);
		
	char *Path_Ez_Hxy = NULL;
	int nr_Save;
	
	int save_field = 1; // to save data

	if (save_field == 1)
	{
		Path_Ez_Hxy =(char *) calloc(512,sizeof(char));
		if (!Path_Ez_Hxy)
		{
			cout << "Memory allocation problem - *Path_Ez_Hxy" << endl;
			return 1;
		}
		strcat(Path_Ez_Hxy,Path);
		strcat(Path_Ez_Hxy,"\\data_EH");
		CreateDirectory(Path_Ez_Hxy,0);

		nr_Save = 10;

		if (!FDTD_2D.Init_Save_Field(n_x_a, n_x_b, n_y_a, n_y_b, nr_Save, Path_Ez_Hxy))
		{
			cout << "Memory allocation problem - Init_Save_Field - FDTD" << endl;
			return 1;
		}
	}

	//The followed field components
	int n_Ind_F = 4;
	int **Ind_F = NULL;
	Ind_F = Init_Matrix_2D<int>(n_Ind_F,2);
	if(!Ind_F)
	{
		cout << "Memory allocation problem - Ind_F" << endl;
		return 1;
	}

	if (jel_plane_wave ==1)
	{
		Ind_F[0][0]  = n_TS_xa + 10;    Ind_F[0][1]  = n_y/2;
		Ind_F[1][0]  = n_TS_xb - 10; 	Ind_F[1][1]  = n_y/2; 
		Ind_F[2][0]  = n_x/2;           Ind_F[2][1]  = n_TS_ya + 5;
		Ind_F[3][0]  = n_x/2;           Ind_F[3][1]  = n_TS_yb - 5; 
	}
	else
	{
		Ind_F[0][0]  = n_PML + 10;          Ind_F[0][1]  = n_y/2;
		Ind_F[1][0]  = n_x - n_PML - 10; 	Ind_F[1][1]  = n_y/2; 
		Ind_F[2][0]  = n_x/2;               Ind_F[2][1]  = n_PML + 5;
		Ind_F[3][0]  = n_x/2;               Ind_F[3][1]  = n_y - n_PML - 5; 
	}

	//Save Ind_F
	char *file_name = NULL;
	file_name = (char *) calloc(512,sizeof(char));
	if (!file_name)
	{
		cout << "Memory allocation problem - *file_name" << endl;
		return 1;
	}
	strcat(file_name,Path);
	strcat(file_name,"/Ind_Foll");
	int x1 =0, y1 = 0, y2 = 2, it = 0;
	if ( !save_2D_int(Ind_F,x1,n_Ind_F,y1,y2,it,file_name) )
		return 2;
		
	FDTD_2D.Init_Followed(Ind_F, n_Ind_F, num_iter);

    /////////////////////////////////////////////////////////////////////////////
	//Save the followed data
	/////////////////////////////////////////////////////////////////////////////
	double **hx_foll = NULL, **hy_foll = NULL, **ez_foll = NULL;
	
	FDTD_2D.Get_Data_Followed(hx_foll, hy_foll, ez_foll);

	/////////////////////////////////////////////////////////////////////////////
	//to measure the ellapsed time
	/////////////////////////////////////////////////////////////////////////////
	clock_t start, finish;
	double  duration;

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

	/////////////////////////////////////////////////////////////////////////////
	//end of the main FDTD loop
	/////////////////////////////////////////////////////////////////////////////

	finish = clock();
	duration = (double)(finish - start) / CLOCKS_PER_SEC;
	cout<< "duration: " << duration << " seconds" <<endl;

	//hx_foll
	strcpy(file_name,Path);
	strcat(file_name,"/Hx_foll");
	if ( !save_2D(hx_foll,x1,n_Ind_F,y1,num_iter,iter,file_name) )
	//if ( !save_2D_binary(hx_foll,n_Ind_F,num_iter,iter,file_name) )
		return 2;

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

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

	/////////////////////////////////////////////////////////////////////////////
	//Free the allocated memory
	/////////////////////////////////////////////////////////////////////////////
	FDTD_2D.Free_Mem(); //Free the memory allocated in the FDTD_2D object

	//Free the memory allocated in the main program
	if (Path)
	{
		free(Path);
		Path = NULL;
	}
	if (Path_Ez_Hxy)
	{
		free(Path_Ez_Hxy);
		Path_Ez_Hxy = NULL;
	}
	if (file_name)
	{
		free(file_name);
		file_name = NULL;
	}
	if (Ind_F)
			Ind_F = Free_Matrix_2D<int>(Ind_F);
	if (Materials)
			Materials = Free_Matrix_2D<double>(Materials);
	if (Index)
			Index = Free_Matrix_2D<int>(Index);
	if (Param_Source)
		    Param_Source = Free_Matrix_2D<double>(Param_Source);
	
	cout << "success" << endl;

	return 0;
}