fdtd_xypbc_zpml.cpp

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

	if (save_3D_xz(Ez,n_x_a,n_x_b,ny_xz,n_z_a,n_z_b,iter,path_name_Ez))
		return -1;
	if (save_3D_yz(Ez,nx_yz,n_y_a,n_y_b,n_z_a,n_z_b,iter,path_name_Ez))
		return -1;

	if (save_3D_xy(Hx,n_x_a,n_x_b,n_y_a,n_y_b,nz_xy,iter,path_name_Hx))
		return -1;
	if (save_3D_xz(Hx,n_x_a,n_x_b,ny_xz,n_z_a,n_z_b,iter,path_name_Hx))
		return -1;
	if (save_3D_yz(Hx,nx_yz,n_y_a,n_y_b,n_z_a,n_z_b,iter,path_name_Hx))
		return -1;

	if (save_3D_xy(Hy,n_x_a,n_x_b,n_y_a,n_y_b,nz_xy,iter,path_name_Hy))
		return -1;
	if (save_3D_xz(Hy,n_x_a,n_x_b,ny_xz,n_z_a,n_z_b,iter,path_name_Hy))
		return -1;
	if (save_3D_yz(Hy,nx_yz,n_y_a,n_y_b,n_z_a,n_z_b,iter,path_name_Hy))
		return -1;

	if (save_3D_xy(Hz,n_x_a,n_x_b,n_y_a,n_y_b,nz_xy,iter,path_name_Hz))
		return -1;
	if (save_3D_xz(Hz,n_x_a,n_x_b,ny_xz,n_z_a,n_z_b,iter,path_name_Hz))
		return -1;
	if (save_3D_yz(Hz,nx_yz,n_y_a,n_y_b,n_z_a,n_z_b,iter,path_name_Hz))
		return -1;

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Average in x, y and z 
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_xyPBC_zPML::aver_xyz_WignerC(double ***X, long nx_a_av, long nx_b_av, long ny_a_av, long ny_b_av,
			                           long nz_a_av, long nz_b_av, double *X_av, double Vol)
{
	long i, j, k;
	double x_aver = 0.0;

	#pragma omp parallel for default(shared) private(i,j,k) \
	                         reduction(+: x_aver) schedule(dynamic,nr_threads)
	for (i = nx_a_av; i <= nx_b_av; i++)
	{
		for (j = ny_a_av; j <= ny_b_av; j++)
		{
			for (k = nz_a_av; k <= nz_b_av; k++)
			{
				if (Mask_Wigner[i][j][k-nz_a_av] > 0)
				{
					x_aver += X[i][j][k];
	  			}
			}
		}
	}

	*X_av = x_aver/Vol;

	return 0;	
}

///////////////////////////////////////////////////////////////////////////////////////
//The Bragg reflection -- Wigner cell
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_xyPBC_zPML::aver_Bragg_WignerC(double ***X, double kx, double ky, 
										 double Gx, double Gy, double *X_avBr_re, 
										 double *X_avBr_im, double Surf)
{
	long i, j;
	double x_aver_re = 0.0, x_aver_im = 0.0;

	double ct_x = (kx+Gx)*dx;
	double ct_y = (ky+Gy)*dy;

	#pragma omp parallel for default(shared) private(i,j) firstprivate(ct_x, ct_y)\
	                         reduction(+: x_aver_re, x_aver_im) schedule(dynamic,nr_threads)
	for (i = 0; i < nx; i++)
	{
		for (j = 0; j < ny; j++)
		{
			if (Mask_Wigner[i][j][0] > 0)
			{
	  			x_aver_re += X[i][j][nz_Bragg]*cos( ct_x*i + ct_y*j);
				x_aver_im += X[i][j][nz_Bragg]*sin( ct_x*i + ct_y*j);
			}
		}
	}

	*X_avBr_re = x_aver_re/Surf;
	*X_avBr_im = x_aver_im/Surf;

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Fourier transform of a field component in the specified volume
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_xyPBC_zPML::fourier_transf_vol(double ***U, long nx_a_f, long nx_b_f, long ny_a_f, long ny_b_f, 
				                         long nz_a_f, long nz_b_f, double **FrT_re, double **FrT_im, 
				                         double *FrT_om, long n_FrT_om, double time)
{
	long i, j, k, ff, pp;
	double coef_1, coef_2;
		 
	for (ff = 0; ff < n_FrT_om; ff++)
	{
		coef_1 = cos(FrT_om[ff]*time);
		coef_2 = sin(FrT_om[ff]*time);
		pp = 0;
 		for (i = nx_a_f; i <= nx_b_f; i++)
		{
			for (j = ny_a_f; j <= ny_b_f; j++)
			{
	  			for (k = nz_a_f; k <= nz_b_f; k++)
				{
					FrT_re[pp][ff] += U[i][j][k]*coef_1; 
					FrT_im[pp][ff] += U[i][j][k]*coef_2;
				}
				pp++;
	  		}
		}
	}

	return 0;	
}

#ifndef run_enviroment //IBM AIX
///////////////////////////////////////////////////////////////////////////////////////
//Measure ellapsed time
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_xyPBC_zPML::Init_Mesure_Ellapsed_Time(double TimeStart, double LimitTime, int LoadWorkspaceData, 
												char *path_load, char *path_save)
{
	time_start = TimeStart;
	limit_time = LimitTime;

	load_workspace_data = LoadWorkspaceData;

	Path_Load_Workspace_Data =(char *) calloc(512,sizeof(char));
	Path_Save_Workspace_Data =(char *) calloc(512,sizeof(char));
	if (!Path_Load_Workspace_Data || !Path_Save_Workspace_Data)
	{
		cout << "memory allocation problem in FDTD_3D_Complex.Init_Mesure_Ellapsed_Time - Path" << endl;
		return 1;
	}
	
	strcpy(Path_Load_Workspace_Data,path_load);
	strcpy(Path_Save_Workspace_Data,path_save);
    	
	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Save workspace data in binary file to be able to continue the computations
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_xyPBC_zPML::save_FE_BH(long iter, double time, char *Path)
{
	mode_t Mode = S_IRWXU;
	char *file_name = NULL;
	file_name =(char *) calloc(512,sizeof(char));
	if (!file_name)
	{
		cout << "memory allocation problem in save_FE_BH" << endl;
		return 1;
	}
	
	//iter
	strcpy(file_name,Path);
	strcat(file_name,"/iter");
	if ( save_1D_long(&iter,0,1,0,file_name) )
		return 2;

	//time
	strcpy(file_name,Path);
	strcat(file_name,"/time");
	if ( save_1D(&time,0,1,0,file_name) )
		return 2;

	
	//save Ex	
	strcpy(file_name,Path);
	strcat(file_name,"/Ex");
	save_3D_binary(Ex,nx,ny,nz,0,file_name);
	
	//save Ey
	strcpy(file_name,Path);
	strcat(file_name,"/Ey");
	save_3D_binary(Ey,nx,ny,nz,0,file_name);

	//save Ez
	strcpy(file_name,Path);
	strcat(file_name,"/Fz");
	save_3D_binary(Fz,nx,ny,nz-1,0,file_name);
	
	strcpy(file_name,Path);
	strcat(file_name,"/Ez");
	save_3D_binary(Ez,nx,ny,nz-1,0,file_name);
	
	//save Hx
	strcpy(file_name,Path);
	strcat(file_name,"/Hx");
	save_3D_binary(Hx,nx,ny,nz-1,0,file_name);

	//save Hy
	strcpy(file_name,Path);
	strcat(file_name,"/Hy");
	save_3D_binary(Hy,nx,ny,nz-1,0,file_name);
	
	//save Hz
	strcpy(file_name,Path);
	strcat(file_name,"/Gz");
	save_3D_binary(Gz,nx,ny,nz,0,file_name);

	strcpy(file_name,Path);
	strcat(file_name,"/Hz");
	save_3D_binary(Hz,nx,ny,nz,0,file_name);

	
	if(file_name)
	{
		free(file_name);
		file_name = NULL;
	}

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Load workspace data from binary file to be able to continue the computations
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_xyPBC_zPML::load_FE_BH(char *Path)
{
	char *file_name = NULL;
	file_name =(char *) calloc(512,sizeof(char));
	if (!file_name)
	{
		cout << "memory allocation problem in load_FGE_BH" << endl;
		return 1;
	}

	//load Ex
	sprintf(file_name,"%s/Ex_%isz%isz%i_0.dat",Path,nx,ny,nz);
	if (load_3D_binary(Ex,nx,ny,nz,file_name))
	{
		cout << "error reading Ex: " << file_name << endl;
		return 1;
	}

	//load Ey
	sprintf(file_name,"%s/Ey_%isz%isz%i_0.dat",Path,nx,ny,nz);
	if (load_3D_binary(Ey,nx,ny,nz,file_name))
	{
		cout << "error reading Ey: " << file_name << endl;
		return 1;
	}

	//load Fz, Ez
	sprintf(file_name,"%s/Fz_%isz%isz%i_0.dat",Path,nx,ny,nz-1);
	if (load_3D_binary(Fz,nx,ny,nz-1,file_name))
	{
		cout << "error reading Fz: " << file_name << endl;
		return 1;
	}

	sprintf(file_name,"%s/Ez_%isz%isz%i_0.dat",Path,nx,ny,nz-1);
	if (load_3D_binary(Ez,nx,ny,nz-1,file_name))
	{
		cout << "error reading Ez: " << file_name << endl;
		return 1;
	}

	//load Hx
	sprintf(file_name,"%s/Hx_%isz%isz%i_0.dat",Path,nx,ny,nz-1);
	if (load_3D_binary(Hx,nx,ny,nz-1,file_name))
	{
		cout << "error reading Hx: " << file_name << endl;
		return 1;
	}

	//load Hy
	sprintf(file_name,"%s/Hy_%isz%isz%i_0.dat",Path,nx,ny,nz-1);
	if (load_3D_binary(Hy,nx,ny,nz-1,file_name))
	{
		cout << "error reading Hy: " << file_name << endl;
		return 1;
	}

	//load Gz, Hz
	sprintf(file_name,"%s/Gz_%isz%isz%i_0.dat",Path,nx,ny,nz);
	if (load_3D_binary(Gz,nx,ny,nz,file_name))
	{
		cout << "error reading Gz: " << file_name << endl;
		return 1;
	}

	sprintf(file_name,"%s/Hz_%isz%isz%i_0.dat",Path,nx,ny,nz);
	if (load_3D_binary(Hz,nx,ny,nz,file_name))
	{
		cout << "error reading Hz: " << file_name << endl;
		return 1;
	}
	
	if(file_name)
	{
		free(file_name);
		file_name = NULL;
	}
	
	return 0;
}
#endif

///////////////////////////////////////////////////////////////////////////////////////
//Free the allocated memory
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::Free_Mem()
{
	//The field components
	if(Ex)
		Ex = Free_Matrix_3D<double>(Ex,nx);
	if(Ey)
		Ey = Free_Matrix_3D<double>(Ey,nx);
	if(Fz)
		Fz = Free_Matrix_3D<double>(Fz,nx);
	if(Ez)
		Ez = Free_Matrix_3D<double>(Ez,nx);
	if(Hx)
		Hx = Free_Matrix_3D<double>(Hx,nx);
	if(Hy)
		Hy = Free_Matrix_3D<double>(Hy,nx);
	if(Gz)
		Gz = Free_Matrix_3D<double>(Gz,nx);
	if(Hz)
		Hz = Free_Matrix_3D<double>(Hz,nx);

	//The PML parameters
	if (K_E1_a)
	{
		free(K_E1_a);
		K_E1_a = NULL;
	}

	if (K_E1_b)
	{
		free(K_E1_b);
		K_E1_b = NULL;
	}

	if (K_E3_a)
	{
		free(K_E3_a);
		K_E3_a = NULL;
	}

	if (K_E3_b)
	{
		free(K_E3_b);
		K_E3_b = NULL;
	}

	if (K_E4_a)
	{
		free(K_E4_a);
		K_E4_a = NULL;
	}

	if (K_E4_b)
	{
		free(K_E4_b);
		K_E4_b = NULL;
	}

	if (K_E6_a)
	{
		free(K_E6_a);
		K_E6_a = NULL;
	}

	if (K_E6_b)
	{
		free(K_E6_b);
		K_E6_b = NULL;
	}

	//Save data
	if (path_name_Ex_1D)
	{
		free(path_name_Ex_1D);
		path_name_Ex_1D = NULL;
	}

	if (path_name_Ex)
	{
		free(path_name_Ex);
		path_name_Ex = NULL;
	}

	if (path_name_Ey)
	{
		free(path_name_Ey);
		path_name_Ey = NULL;
	}

	if (path_name_Ez)
	{
		free(path_name_Ez);
		path_name_Ez = NULL;
	}

	if (path_name_Hx)
	{
		free(path_name_Hx);
		path_name_Hx = NULL;
	}

	if (path_name_Hy)
	{
		free(path_name_Hy);
		path_name_Hy = NULL;
	}

	if (path_name_Hz)
	{
		free(path_name_Hz);
		path_name_Hz = NULL;
	}

	//The followed field components
	if(Hx_Foll)
		Hx_Foll = Free_Matrix_2D<double>(Hx_Foll);
	if(Hy_Foll)
		Hy_Foll = Free_Matrix_2D<double>(Hy_Foll);
	if(Hz_Foll)
		Hz_Foll = Free_Matrix_2D<double>(Hz_Foll);
	if(Ex_Foll)
		Ex_Foll = Free_Matrix_2D<double>(Ex_Foll);
	if(Ey_Foll)
		Ey_Foll = Free_Matrix_2D<double>(Ey_Foll);
	if(Ez_Foll)
		Ez_Foll = Free_Matrix_2D<double>(Ez_Foll);

	//The averaged field components
	if (Ex_Foll_av) 
	{
		free(Ex_Foll_av);
		Ex_Foll_av = NULL;
	} 
	if (Ey_Foll_av) 
	{
		free(Ey_Foll_av);
		Ey_Foll_av = NULL;
	}
	if (Ez_Foll_av)
	{
		free(Ez_Foll_av);
		Ez_Foll_av = NULL;
	}
	if (Hx_Foll_av) 
	{
		free(Hx_Foll_av);
		Hx_Foll_av = NULL;
	}
	if (Hy_Foll_av) 
	{
		free(Hy_Foll_av);
		Hy_Foll_av = NULL;
	} 
	if (Hz_Foll_av)
	{
		free(Hz_Foll_av);
		Hz_Foll_av = NULL;
	}

	//Bragg diffraction
	if (Ex_Br_av_re)
		Ex_Br_av_re = Free_Matrix_2D<double>(Ex_Br_av_re);
	if (Ex_Br_av_im)
		Ex_Br_av_im = Free_Matrix_2D<double>(Ex_Br_av_im);
	if (Ey_Br_av_re)
		Ey_Br_av_re = Free_Matrix_2D<double>(Ey_Br_av_re);
	if (Ey_Br_av_im)
		Ey_Br_av_im = Free_Matrix_2D<double>(Ey_Br_av_im);
	if (Ez_Br_av_re)
		Ez_Br_av_re = Free_Matrix_2D<double>(Ez_Br_av_re);
	if (Ez_Br_av_im)
		Ez_Br_av_im = Free_Matrix_2D<double>(Ez_Br_av_im);
	if (Hx_Br_av_re)
		Hx_Br_av_re = Free_Matrix_2D<double>(Hx_Br_av_re);
	if (Hx_Br_av_im)
		Hx_Br_av_im = Free_Matrix_2D<double>(Hx_Br_av_im);
	if (Hy_Br_av_re)
		Hy_Br_av_re = Free_Matrix_2D<double>(Hy_Br_av_re);
	if (Hy_Br_av_im)
		Hy_Br_av_im = Free_Matrix_2D<double>(Hy_Br_av_im);
	if (Hz_Br_av_re)
		Hz_Br_av_re = Free_Matrix_2D<double>(Hz_Br_av_re);
	if (Hz_Br_av_im)
		Hz_Br_av_im = Free_Matrix_2D<double>(Hz_Br_av_im);
	
	//In case of IBM AIX -- load and save workspace data
	#ifndef run_enviroment
		if (Path_Load_Workspace_Data)
		{
			free(Path_Load_Workspace_Data);
			Path_Load_Workspace_Data = NULL;
		}

		if (Path_Save_Workspace_Data)
		{
			free(Path_Save_Workspace_Data);
			Path_Save_Workspace_Data = NULL;
		}
	#endif
}