fdtd_1d_eh_pml_loss.cpp

fdtd 3D xyzPML MPI OpenMP

#include "FDTD_1D_EH_PML_LOSS.h"
#include "run_enviroment.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CFDTD_1D_EH_PML_LOSS::CFDTD_1D_EH_PML_LOSS(void)
{
	E_1D   = NULL;
	//F_1 = NULL;
	F_2 = NULL;
	H_1D   = NULL;

	K_F_a = NULL;
	K_F_b = NULL;
	K_H_a = NULL;
	K_H_b = NULL;

	pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679;
	//permittivity of free space 
	eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]
	//permeability of free space 
    mu_0  = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6;  // [H_1D/m]

	num_threads = 1;
}

CFDTD_1D_EH_PML_LOSS::~CFDTD_1D_EH_PML_LOSS(void)
{
	Free_Mem();
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize the nr of threads
///////////////////////////////////////////////////////////////////////////////////////
#ifdef _OPENMP
	void CFDTD_1D_EH_PML_LOSS::Init_nr_THR(int nr_threads)
	{
		num_threads = nr_threads;
	}
#endif

///////////////////////////////////////////////////////////////////////////////////////
//Initialize the class parameters
///////////////////////////////////////////////////////////////////////////////////////
bool CFDTD_1D_EH_PML_LOSS::Init(long N, long n_pml, double d_t, double D, double Eps_r, 
								double Mu_r, double Sigma)
{
	n = N;
	dt = d_t;
	d = D;
	inv_d = 1.0/d;
	n_PML = n_pml;

	n_TS_a = n_PML+10;

	//some usefull constants
	n_PML_MIN_1 = n_PML - 1;
	n_MIN_n_PML_MIN_1 = n - n_PML - 1;
	n_MIN_1 = n - 1;
	n_MIN_n_PML = n - n_PML;

	eps_r = Eps_r;
	mu_r = Mu_r;
	sigma = Sigma;

	beta = sqrt(eps_0*eps_r*mu_0*mu_r);
	eta  = sqrt( mu_0*mu_r/(eps_0*eps_r) );

	E_1D = (double *)calloc(n,sizeof(double));
	if (!E_1D)
	{
		return false;
	}

	/*F_1 = (double *)calloc(n_PML,sizeof(double));
	if (!F_1)
	{
		Free_Mem();
		return false;
	}*/

	F_2 = (double *)calloc(n_PML,sizeof(double));
	if (!F_2)
	{
		Free_Mem();
		return false;
	}

	K_F_a = (double *)calloc(n_PML,sizeof(double));
	if (!K_F_a)
	{
		Free_Mem();
		return false;
	}

	K_F_b = (double *)calloc(n_PML,sizeof(double));
	if (!K_F_b)
	{
		Free_Mem();
		return false;
	}

	H_1D = (double *)calloc(n_MIN_1,sizeof(double));
	if (!H_1D)
	{
		Free_Mem();
		return false;
	}

	K_H_a = (double *)calloc(n_PML,sizeof(double));
	if (!K_H_a)
	{
		Free_Mem();
		return false;
	}

 	K_H_b = (double *)calloc(n_PML,sizeof(double));
	if (!K_H_b)
	{
		Free_Mem();
		return false;
	}

	K_E_a = ( 2.0*eps_r*eps_0 - sigma*dt )/( 2.0*eps_r*eps_0 + sigma*dt );
	K_E_b =  2.0*dt/( 2.0*eps_r*eps_0 + sigma*dt )*inv_d;

	K_H = dt*inv_d/(mu_0*mu_r);
	
	return true;
}

///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_PML_LOSS::Init_PML_Param()
{
	long i;

	//PML parameters
	double ka_max = 2.0;
	double exponent = 4.0;

	double R_err = 1.0e-16;
	double eta = sqrt( mu_0*mu_r/(eps_0*eps_r) );

	double sigma_PML, sigma_max, ka_PML;
		
	sigma_max= -(exponent+1.0)*log(R_err)/(2.0*eta*n_PML*d);

	for (i = 0; i < n_PML; i++)
	{
		sigma_PML     = sigma_max*pow( (n_PML - i)/((double) n_PML) ,exponent);
		ka_PML        = 1.0 + (ka_max - 1.0)*pow( (n_PML-i)/((double) n_PML) ,exponent);

		K_F_a[i]      = (2.0*eps_0*ka_PML - sigma_PML*dt)/(2.0*eps_0*ka_PML + sigma_PML*dt);
		
		K_F_b[i]      = 2.0*eps_0/(2.0*eps_0*ka_PML + sigma_PML*dt);
		
		sigma_PML     = sigma_max*pow( (n_PML - i - 0.5)/n_PML ,exponent);
		ka_PML        = 1.0 + (ka_max - 1.0)*pow( (n_PML - i - 0.5)/n_PML ,exponent);

        K_H_a[i]      = (2.0*eps_0*ka_PML - sigma_PML*dt)/(2.0*eps_0*ka_PML + sigma_PML*dt);
		
		K_H_b[i]      = 2.0*eps_0*dt/( (2.0*eps_0*ka_PML + sigma_PML*dt)*d*mu_0*mu_r);
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize a sinusoidal source
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_PML_LOSS::Init_Sin_Source(double E_0, double om, double Phase)
{
	E0 = E_0;
	omega = om;
	phase = Phase;

	alfa = 1.0*omega/(2.0*pi);

	Source_Type = 1;
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize  a Gaussian source
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_PML_LOSS::Init_Gauss_Source(double E_0, double t_0, double t_w)
{
	E0 = E_0;
	t0 = t_0;
	tw = t_w;
	
	Source_Type = 2;
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize a sinusoidal modulated Gaussian pulse
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_PML_LOSS::Init_GaussSin_Source(double E_0, double om, double Phase, double t_0, double t_w)
{
	E0 = E_0;
	omega = om;
	phase = Phase;
	t0 = t_0;
	tw = t_w;

	Source_Type = 3;
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate E_1D field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_PML_LOSS::Calc_E(double time)
{
	long i, ii, ia;
	double Fz_r;

	/*for (i = 1; i < n_PML; i++)
	{
		Fz_r = F_1[i];
		
		F_1[i] = K_F_a[i]*F_1[i] + K_F_b[i]*( H_1D[i-1] - H_1D[i]);

		E_1D[i] = K_E_a*E_1D[i] + K_E_b*( F_1[i] - Fz_r );
	}*/

	#pragma omp parallel for shared(H_1D,K_E_a,K_E_b,E_1D,n_MIN_n_PML) private(i) num_threads(num_threads)
	for (i = 1; i < n_MIN_n_PML; i++)
	{
		E_1D[i] = K_E_a*E_1D[i] + K_E_b*( H_1D[i-1] - H_1D[i]);
	}

	ia = 0;
	ii = n_PML_MIN_1;
	for (i = n_MIN_n_PML; i < n_MIN_1; i++)
	{
		Fz_r = F_2[ia];
		
		F_2[ia] = K_F_a[ii]*F_2[ia] + K_F_b[ii]*( H_1D[i-1] - H_1D[i]);

		E_1D[i] = K_E_a*E_1D[i] + K_E_b*( F_2[ia] - Fz_r );
		
		ia++;
		ii--;
	}

	switch (Source_Type)
	{
		case 1:
			E_1D[0] = (1.0-exp(-alfa*time))*E0*cos(omega*time + phase);
			break;
		case 2:
			E_1D[0] = E0*exp( - (time - t0)*(time - t0)/(tw*tw) );
			break;
		case 3:
			E_1D[0] = E0*cos(omega*(time-t0) + phase)*exp( - (time - t0)*(time - t0)/(tw*tw) );
			break;
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate H_1D field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_PML_LOSS::Calc_H()
{
	long i, ii;

	/*for (i = 0; i < n_PML; i++)
	{
		H_1D[i] = K_H_a[i]*H_1D[i] + K_H_b[i]*( E_1D[i] - E_1D[i+1]);
	}*/

	#pragma omp parallel for shared(H_1D,K_H,E_1D,n_MIN_n_PML_MIN_1) private(i) num_threads(num_threads)
	for (i = 0; i < n_MIN_n_PML_MIN_1; i++)
	{
		H_1D[i] += K_H*( E_1D[i] - E_1D[i+1]);
	}

	ii = n_PML_MIN_1;
	for (i = n_MIN_n_PML_MIN_1; i < n_MIN_1; i++)
	{
		H_1D[i] = K_H_a[ii]*H_1D[i] + K_H_b[ii]*( E_1D[i] - E_1D[i+1] );
		ii--;
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Get the address of E_1D and H_1D
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_PML_LOSS::Get_Data(double *&e, double *&h)
{
	e = E_1D;
	h = H_1D;
}

///////////////////////////////////////////////////////////////////////////////////////
//Save out the workspace data - the possibility to continue computations
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_1D_EH_PML_LOSS::Save_FDTD_1D_Workspace(char *path)
{	
	char *path_name = NULL, *file_name = NULL;
	path_name =(char *) calloc(256,sizeof(char));
	file_name =(char *) calloc(256,sizeof(char));
	if (!path_name && !file_name)
	{
		if (path_name)
			free(path_name);
		if (file_name)
			free(file_name);
		return 1;
	}
	strcpy(path_name,path);
	strcat(path_name,"/data_fdtd_1D");

	//make directory to save F G E_1D
	#ifdef run_enviroment 
		CreateDirectory(path_name,0);
	#else
		mode_t Mode = S_IRWXU;
		mkdir(path_name,Mode);//for UNIX-AIX
	#endif

	//save F_2
	strcpy(file_name,path_name);
	strcat(file_name,"/F_1D");
	if (save_1D_binary(F_2, n_PML, 0, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}

	//save E_1D
	strcpy(file_name,path_name);
	strcat(file_name,"/E_1D");
	if (save_1D_binary(E_1D, n, 0, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}

	//save H_1D
	strcpy(file_name,path_name);
	strcat(file_name,"/H_1D");
	if(save_1D_binary(H_1D, n-1, 0, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}


	free(path_name);
	free(file_name);

	return 0;
}

int CFDTD_1D_EH_PML_LOSS::Load_FDTD_1D_Workspace(char *path)
{
	char *file_name = NULL;
	file_name =(char *) calloc(256,sizeof(char));
	if (!file_name)
	{
		return 1;
	}
	
	//load F2
	strcpy(file_name,path);
	strcat(file_name,"/F_1D_0.dat");
	if (load_1D_binary(F_2, n_PML, file_name))
	{
		free(file_name);
		return 2; //faild to load data
	}

	//load E_1D
	strcpy(file_name,path);
	strcat(file_name,"/E_1D_0.dat");
	if (load_1D_binary(E_1D, n, file_name))
	{
		free(file_name);
		return 2; //faild to load data
	}

	//load Hx_1D
	strcpy(file_name,path);
	strcat(file_name,"/H_1D_0.dat");
	if (load_1D_binary(H_1D, n-1, file_name))
	{
		free(file_name);
		return 2; //faild to load data
	}

	free(file_name);

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Free the allocated memory
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_PML_LOSS::Free_Mem(void)
{
	if (E_1D)
	{
		free(E_1D);
		E_1D = NULL;
	}

	/*if (F_1)
	{
		free(F_1);
		F_1 = NULL;
	}*/

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

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

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

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

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

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