fdtd_3d.cpp

fdtd 3D xyzPML MPI OpenMP

		ErrorMessage(1, myrank, " -- Memory allocation problem - Hy_send_k - ");
		return 1;
	}

	///////////////////////////////////////////////////////////////////////////////////////
	//memory allocations to compute Hz
	///////////////////////////////////////////////////////////////////////////////////////
	ista_Hz = ista;	   iend_Hz = iend;    jsta_Hz = jsta;
	jend_Hz = jend;	   ksta_Hz = ksta;	  kend_Hz = kend;
	nlx_Hz = nlx;	   nly_Hz = nly;	  nlz_Hz = nlz;
	if (myrank_i == iprocs-1)
	{
		iend_Hz--;
		nlx_Hz--;
	}
	if (myrank_j == jprocs-1)
	{
		jend_Hz--;
		nly_Hz--;
	}

	nlx_HzMIN1 = nlx_Hz - 1;
	nly_HzMIN1 = nly_Hz - 1;

	Hz = Init_Matrix_3D<double>(nlx_Hz, nly_Hz, nlz_Hz);
	if (!Hz)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hz - ");
		return 1;
	}

	Bz = Init_Matrix_3D<double>(nlx_Hz, nly_Hz, nlz_Hz);
	if (!Bz)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Bz - ");
		return 1;
	}

	Ex_recv_j = Init_Matrix_2D<double>(nlx_Hz, nlz_Hz);
	if (!Ex_recv_j)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ex_recv_j - ");
		return 1;
	}

	Ey_recv_i = Init_Matrix_2D<double>(nly_Hz, nlz_Hz);
	if (!Ey_recv_i)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ey_recv_i - ");
		return 1;
	}

	Hz_send_i = Init_Matrix_2D<double>(nly_Hz, nlz_Hz);
	if (!Hz_send_i)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hz_send_i - ");
		return 1;
	}

	Hz_send_j = Init_Matrix_2D<double>(nlx_Hz, nlz_Hz);
	if (!Hz_send_j)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hz_send_j - ");
		return 1;
	}
	
	///////////////////////////////////////////////////////////////////////////////////////
	//Coefficients containing the PML boundary parameters
	///////////////////////////////////////////////////////////////////////////////////////
	//PML-Ex field
	K_Gx_a = (double *)calloc(nly_Ex,sizeof(double)); 
	if (!K_Gx_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Gx_a - ");
		return 1;
	}

	K_Gx_b = (double *)calloc(nly_Ex,sizeof(double)); 
	if (!K_Gx_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Gx_b - ");
		return 1;
	}

	K_Ex_a = (double *)calloc(nlz_Ex,sizeof(double)); 
	if (!K_Ex_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Ex_a - ");
		return 1;
	}

	K_Ex_b = (double *)calloc(nlz_Ex,sizeof(double)); 
	if (!K_Ex_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Ex_b - ");
		return 1;
	}

	K_Ex_c = (double *)calloc(nlx_Ex,sizeof(double)); 
	if (!K_Ex_c)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Ex_c - ");
		return 1;
	}

	K_Ex_d = (double *)calloc(nlx_Ex,sizeof(double)); 
	if (!K_Ex_d)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Ex_d - ");
		return 1;
	}

	//PML-Ey field
	K_Gy_a = (double *)calloc(nlz_Ey,sizeof(double)); 
	if (!K_Gy_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Gy_a - ");
		return 1;
	}

	K_Gy_b = (double *)calloc(nlz_Ey,sizeof(double)); 
	if (!K_Gy_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Gy_b - ");
		return 1;
	}

	K_Ey_a = (double *)calloc(nlx_Ey,sizeof(double)); 
	if (!K_Ey_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Ey_a - ");
		return 1;
	}

	K_Ey_b = (double *)calloc(nlx_Ey,sizeof(double)); 
	if (!K_Ey_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Ey_b - ");
		return 1;
	}

	K_Ey_c = (double *)calloc(nly_Ey,sizeof(double)); 
	if (!K_Ey_c)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Ey_c - ");
		return 1;
	}

	K_Ey_d = (double *)calloc(nly_Ey,sizeof(double)); 
	if (!K_Ey_d)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Ey_d - ");
		return 1;
	}

	//PML-Ez field
	K_Gz_a = (double *)calloc(nlx_Ez,sizeof(double)); 
	if (!K_Gz_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Gz_a - ");
		return 1;
	}

	K_Gz_b = (double *)calloc(nlx_Ez,sizeof(double)); 
	if (!K_Gz_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Gz_b - ");
		return 1;
	}

	K_Ez_a = (double *)calloc(nly_Ez,sizeof(double)); 
	if (!K_Ez_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Ez_a - ");
		return 1;
	}

	K_Ez_b = (double *)calloc(nly_Ez,sizeof(double)); 
	if (!K_Ez_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Ez_b - ");
		return 1;
	}

	K_Ez_c = (double *)calloc(nlz_Ez,sizeof(double)); 
	if (!K_Ez_c)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Ez_c - ");
		return 1;
	}

	K_Ez_d = (double *)calloc(nlz_Ez,sizeof(double)); 
	if (!K_Ez_d)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Ez_d - ");
		return 1;
	}

	//PML-Hx field
	K_Bx_a = (double *)calloc(nly_Hx,sizeof(double)); 
	if (!K_Bx_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Bx_a - ");
		return 1;
	}

	K_Bx_b = (double *)calloc(nly_Hx,sizeof(double)); 
	if (!K_Bx_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Bx_b - ");
		return 1;
	}

	K_Hx_a = (double *)calloc(nlz_Hx,sizeof(double)); 
	if (!K_Hx_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hx_a - ");
		return 1;
	}

	K_Hx_b = (double *)calloc(nlz_Hx,sizeof(double)); 
	if (!K_Hx_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hx_b - ");
		return 1;
	}

	K_Hx_c = (double *)calloc(nlx_Hx,sizeof(double)); 
	if (!K_Hx_c)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hx_c - ");
		return 1;
	}

	K_Hx_d = (double *)calloc(nlx_Hx,sizeof(double)); 
	if (!K_Hx_d)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hx_d - ");
		return 1;
	}

	//PML-Hy field
	K_By_a = (double *)calloc(nlz_Hy,sizeof(double)); 
	if (!K_By_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_By_a - ");
		return 1;
	}

	K_By_b = (double *)calloc(nlz_Hy,sizeof(double)); 
	if (!K_By_b )
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_By_b  - ");
		return 1;
	}

	K_Hy_a = (double *)calloc(nlx_Hy,sizeof(double)); 
	if (!K_Hy_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hy_a - ");
		return 1;
	}

	K_Hy_b = (double *)calloc(nlx_Hy,sizeof(double)); 
	if (!K_Hy_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hy_b - ");
		return 1;
	}

	K_Hy_c = (double *)calloc(nly_Hy,sizeof(double)); 
	if (!K_Hy_c)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hy_c - ");
		return 1;
	}

	K_Hy_d = (double *)calloc(nly_Hy,sizeof(double)); 
	if (!K_Hy_d)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hy_d - ");
		return 1;
	}
	
	//PML-Hz field
	K_Bz_a = (double *)calloc(nlx_Hz,sizeof(double)); 
	if (!K_Bz_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Bz_a - ");
		return 1;
	}

	K_Bz_b = (double *)calloc(nlx_Hz,sizeof(double)); 
	if (!K_Bz_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Bz_b - ");
		return 1;
	}

	K_Hz_a = (double *)calloc(nly_Hz,sizeof(double)); 
	if (!K_Hz_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hz_a - ");
		return 1;
	}

	K_Hz_b = (double *)calloc(nly_Hz,sizeof(double)); 
	if (!K_Hz_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hz_b - ");
		return 1;
	}

	K_Hz_c = (double *)calloc(nlz_Hz,sizeof(double)); 
	if (!K_Hz_c)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hz_c - ");
		return 1;
	}

	K_Hz_d = (double *)calloc(nlz_Hz,sizeof(double)); 
	if (!K_Hz_d)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hz_d - ");
		return 1;
	}
	
	///////////////////////////////////////////////////////////////////////////////////////
	//Materials matrices
	///////////////////////////////////////////////////////////////////////////////////////
	//Materials matrices
	K_a = (double *)calloc(n_Mat,sizeof(double)); 
	if (!K_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ka - ");
		return 1;
	}

	K_b = (double *)calloc(n_Mat,sizeof(double)); 
	if (!K_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Kb - ");
		return 1;
	}

	///////////////////////////////////////////////////////////////////////////////////////
	//Initialize K_a and K_b - the contribution of materials
	///////////////////////////////////////////////////////////////////////////////////////
	double eps_r, sigma;
	long i;
	for (i =0; i<n_Mat; i++)
	{
		eps_r = Mat[i][0];
		sigma = Mat[i][2];
		K_a[i] = ( 2.0*eps_0*eps_r - sigma*dt )/( 2.0*eps_0*eps_r + sigma*dt );
		K_b[i] = 2.0*dt/( 2.0*eps_0*eps_r + sigma*dt );
	}

	//Init the magnetic permittivity
	mu_r = (double *) calloc(n_Mat,sizeof(double));
	if (!mu_r)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Mu_r - ");
		return 1;
	}

	for (i = 0; i<n_Mat; i++)
	{
		mu_r[i] = Mat[i][1];
	}

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices in x directions
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_3D::Init_PML_Par_x(double eps_r_x_1, double mu_r_x_1, double eps_r_x_2, 
							 double mu_r_x_2)
{
	long  i, ii; 
	long nlx = iend - ista + 1;
	long nlxMIN1 = nlx -1;

	//Outside of the PML region
	//i
	for (i = 0; i < nlxMIN1; i++)
	{
		K_Ex_c[i] =  2.0*eps_0;
		K_Ex_d[i] = -2.0*eps_0;

		K_Ey_a[i] = 1.0;
		K_Ey_b[i] = 1.0/(2.0*eps_0);

		K_Gz_a[i] = 1.0;
		K_Gz_b[i] = 1.0;

		K_Hx_c[i] =  2.0*eps_0/mu_0;
		K_Hx_d[i] = -2.0*eps_0/mu_0;

		K_Hy_a[i] = 1.0;
		K_Hy_b[i] = 1.0/(2.0*eps_0);

		K_Bz_a[i] = 1.0;
		K_Bz_b[i] = dt;
	}
	K_Ey_a[nlxMIN1] =  1.0;
	K_Ey_b[nlxMIN1] =  1.0/(2.0*eps_0);
	K_Gz_a[nlxMIN1] =  1.0;
	K_Gz_b[nlxMIN1] =  1.0;
	K_Hx_c[nlxMIN1] =  2.0*eps_0/mu_0;
	K_Hx_d[nlxMIN1] = -2.0*eps_0/mu_0;
	if ( iend < nx-1 )
	{
		K_Ex_c[nlxMIN1] =  2.0*eps_0;
		K_Ex_d[nlxMIN1] = -2.0*eps_0;
		K_Hy_a[nlxMIN1] = 1.0;
		K_Hy_b[nlxMIN1] = 1.0/(2.0*eps_0);
		K_Bz_a[nlxMIN1] = 1.0;
		K_Bz_b[nlxMIN1] = dt;
	}

	//PML_x parameters
	double ka_max = 1.0;
	int  exponent = 4;
	double R_err = 1e-16;
	
	double eta_1 = sqrt(mu_0*mu_r_x_1/eps_0/eps_r_x_1);
	double eta_2 = sqrt(mu_0*mu_r_x_2/eps_0/eps_r_x_2);
	
	double sigma_x, ka_x;
	double sigma_max_1= -(exponent+1.0)*log(R_err)/(2.0*eta_1*nPML_x_1*dx);
	double sigma_max_2= -(exponent+1.0)*log(R_err)/(2.0*eta_2*nPML_x_2*dx);

	long n1 = 1, n2 = 0, jel1 = 0, jel2 = 0;
	if (iend <= nPML_x_1-1)
	{
		n1 = ista;
		n2 = iend;
		jel1 = 1;
	}
	if ( (ista <= nPML_x_1-1) && (iend >= nPML_x_1-1) )
	{
		n1 = ista;
		n2 = nPML_x_1-1;
		jel1 = 1;
	}
	if ( (ista <= nx - nPML_x_2) && (iend > nx - nPML_x_2) )
	{
		n1 = ( nx - iend-1);
		n2 = nPML_x_2-1;
		jel2 = 1;
	}
	if ( ista > nx - nPML_x_2)
	{
		n1 = 0;
		n2 = nlxMIN1;
		jel2 = 1;
	}
    
	long cik = 0;
	for (i = n1; i <= n2; i++)
	{
		//i
		if (jel1 == 1)
		{
			sigma_x         = sigma_max_1*pow( (nPML_x_1 - i)/((double) nPML_x_1) ,exponent);
			ka_x            = 1.0 + (ka_max - 1.0)*pow( (nPML_x_1 - i)/((double) nPML_x_1) ,exponent);

			ii = cik;

			K_Ey_a[ii]     = (2.0*eps_0*ka_x - sigma_x*dt)/(2.0*eps_0*ka_x + sigma_x*dt);
			
			K_Ey_b[ii]     = 1.0/(2.0*eps_0*ka_x + sigma_x*dt);
			
			K_Gz_a[ii]     = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
			
			K_Gz_b[ii]     = 2.0*eps_0/(2.0*eps_0*ka_x + sigma_x*dt);
			
			K_Hx_c[ii]     = (2.0*eps_0*ka_x + sigma_x*dt)/mu_0;
			
			K_Hx_d[ii]     = -(2.0*eps_0*ka_x - sigma_x*dt)/mu_0;
		}

		if (jel2 == 1)
		{
			sigma_x         = sigma_max_2*pow( (nPML_x_2 - i)/((double) nPML_x_2) ,exponent);
			ka_x            = 1.0 + (ka_max - 1.0)*pow( (nPML_x_2 - i)/((double) nPML_x_2) ,exponent);

			ii = nlx - cik - 1;
			K_Ey_a[ii]     = (2.0*eps_0*ka_x - sigma_x*dt)/(2.0*eps_0*ka_x + sigma_x*dt);
			
			K_Ey_b[ii]     = 1.0/(2.0*eps_0*ka_x + sigma_x*dt);
			
			K_Gz_a[ii]     = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
			
			K_Gz_b[ii]     = 2.0*eps_0/(2.0*eps_0*ka_x + sigma_x*dt);
			
			K_Hx_c[ii]     = (2.0*eps_0*ka_x + sigma_x*dt)/mu_0;