fdtd_3d.cpp

fdtd 3D xyzPML MPI OpenMP

			
			K_Hx_d[ii]     = -(2.0*eps_0*ka_x - sigma_x*dt)/mu_0;
		}
		
		//i+0.5
		if (jel1 == 1)
		{
			sigma_x         = sigma_max_1*pow( (nPML_x_1 - i - 0.5)/nPML_x_1 ,exponent);
			ka_x            = 1.0 + (ka_max - 1.0)*pow( (nPML_x_1 - i - 0.5)/nPML_x_1 ,exponent);

			ii = cik;
			K_Ex_c[ii]     = 2.0*eps_0*ka_x + sigma_x*dt;
				
			K_Ex_d[ii]     = -(2.0*eps_0*ka_x - sigma_x*dt);
			
			K_Hy_a[ii]     = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
			
			K_Hy_b[ii]     = 1.0/(2.0*eps_0*ka_x+sigma_x*dt);
			
			K_Bz_a[ii]     = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
			
			K_Bz_b[ii]     = (2.0*eps_0*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
		}

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

			ii = nlx-cik-2;
			if (ii >=0 )
			{
				K_Ex_c[ii]     = 2.0*eps_0*ka_x + sigma_x*dt;
				
				K_Ex_d[ii]     = -(2.0*eps_0*ka_x - sigma_x*dt);
				
				K_Hy_a[ii]     = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
				
				K_Hy_b[ii]     = 1.0/(2.0*eps_0*ka_x+sigma_x*dt);
				
				K_Bz_a[ii]     = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
				
				K_Bz_b[ii]     = (2.0*eps_0*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
			}
		}
	
		cik++;
	}
	
	if ( (iend < nx-1) && (iend >= nx - nPML_x_2 - 1) )
	{
		i = n1;
		if (ii < 0)
		{
			i--;
		}

		sigma_x         = sigma_max_2*pow( (nPML_x_2 - i - 0.5)/nPML_x_2 ,exponent);
		ka_x            = 1.0 + (ka_max - 1.0)*pow( (nPML_x_2 - i - 0.5)/nPML_x_2 ,exponent);

		K_Ex_c[nlxMIN1]     = 2.0*eps_0*ka_x + sigma_x*dt;
		
		K_Ex_d[nlxMIN1]     = -(2.0*eps_0*ka_x - sigma_x*dt);
		
		K_Hy_a[nlxMIN1]     = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
		
		K_Hy_b[nlxMIN1]     = 1.0/(2.0*eps_0*ka_x+sigma_x*dt);
		
		K_Bz_a[nlxMIN1]     = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
		
		K_Bz_b[nlxMIN1]     = (2.0*eps_0*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
	}
	
	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices in x directions
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_3D::Init_PML_Par_y(double eps_r_y_1, double mu_r_y_1, double eps_r_y_2, 
					         double mu_r_y_2)
{
	long  j, jj; 
	long nly = jend - jsta + 1;
	long nlyMIN1 = nly -1;

	//Outside of the PML region
	//j
	for (j = 0; j < nlyMIN1; j++)
	{
		K_Gx_a[j] = 1.0;
		K_Gx_b[j] = 1.0;

		K_Ey_c[j] =  2.0*eps_0;
		K_Ey_d[j] = -2.0*eps_0;

		K_Ez_a[j] = 1.0;
		K_Ez_b[j] = 1.0/(2.0*eps_0);

		K_Bx_a[j] = 1.0;
		K_Bx_b[j] = dt;

		K_Hy_c[j] = 2*eps_0/mu_0;
		K_Hy_d[j] = -2*eps_0/mu_0;

		K_Hz_a[j] = 1.0;
		K_Hz_b[j] = 1.0/(2.0*eps_0);
	}
	K_Gx_a[nlyMIN1] =  1.0;
	K_Gx_b[nlyMIN1] =  1.0;
	K_Ez_a[nlyMIN1] =  1.0;
	K_Ez_b[nlyMIN1] =  1.0/(2.0*eps_0);
	K_Hy_c[nlyMIN1] =  2*eps_0/mu_0;
	K_Hy_d[nlyMIN1] = -2*eps_0/mu_0;
	if ( jend < ny-1 )
	{
		K_Ey_c[nlyMIN1] =  2.0*eps_0;
		K_Ey_d[nlyMIN1] = -2.0*eps_0;
		K_Bx_a[nlyMIN1] = 1.0;
		K_Bx_b[nlyMIN1] = dt;
		K_Hz_a[nlyMIN1] = 1.0;
		K_Hz_b[nlyMIN1] = 1.0/(2.0*eps_0);
	}

    //PML_y parameters
	double ka_max = 1;
	int  exponent = 4;
	double R_err = 1e-16;

	double eta_1 = sqrt(mu_0*mu_r_y_1/eps_0/eps_r_y_1);
	double eta_2 = sqrt(mu_0*mu_r_y_2/eps_0/eps_r_y_2);
	
	double sigma_y, ka_y;
	double sigma_max_1= -(exponent+1.0)*log(R_err)/(2.0*eta_1*nPML_y_1*dy);
	double sigma_max_2= -(exponent+1.0)*log(R_err)/(2.0*eta_2*nPML_y_2*dy);

	long n1 = 1, n2 = 0, jel1 = 0, jel2 = 0;
	if (jend <= nPML_y_1-1)
	{
		n1 = jsta;
		n2 = jend;
		jel1 = 1;
	}
	if ( (jsta <= nPML_y_1-1) && (jend >= nPML_y_1-1) )
	{
		n1 = jsta;
		n2 = nPML_y_1-1;
		jel1 = 1;
	}
	if ( (jsta <= ny - nPML_y_2) && (jend > ny - nPML_y_2) )
	{
		n1 = ( ny - jend-1);
		n2 = nPML_y_2-1;
		jel2 = 1;
	}
	if ( jsta > ny - nPML_y_2)
	{
		n1 = 0;
		n2 = nlyMIN1;
		jel2 = 1;
	}
    
	long cik = 0;
	for (j = n1; j <= n2; j++)
	{
		//j
		if (jel1 == 1)
		{
			sigma_y         = sigma_max_1*pow( (nPML_y_1 - j)/((double) nPML_y_1) ,exponent);
			ka_y            = 1.0 + (ka_max - 1.0)*pow( (nPML_y_1 - j)/((double) nPML_y_1) ,exponent);

			jj = cik;
			K_Gx_a[jj]     = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
			
			K_Gx_b[jj]     = 2.0*eps_0/(2.0*eps_0*ka_y+sigma_y*dt);
			
			K_Ez_a[jj]     = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
			
			K_Ez_b[jj]     = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);
			
			K_Hy_c[jj]     = (2.0*eps_0*ka_y + sigma_y*dt)/mu_0;
			
			K_Hy_d[jj]     = -(2.0*eps_0*ka_y - sigma_y*dt)/mu_0;
		}

		if (jel2 == 1)
		{
			sigma_y         = sigma_max_2*pow( (nPML_y_2 - j)/((double) nPML_y_2) ,exponent);
			ka_y            = 1.0 + (ka_max - 1.0)*pow( (nPML_y_2 - j)/((double) nPML_y_2) ,exponent);

			jj = nly - cik - 1;
			K_Gx_a[jj]     = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
			
			K_Gx_b[jj]     = 2.0*eps_0/(2.0*eps_0*ka_y+sigma_y*dt);
			
			K_Ez_a[jj]     = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
			
			K_Ez_b[jj]     = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);
			
			K_Hy_c[jj]     = (2.0*eps_0*ka_y + sigma_y*dt)/mu_0;
			
			K_Hy_d[jj]     = -(2.0*eps_0*ka_y - sigma_y*dt)/mu_0;
		}
		
		//j+0.5
		if (jel1 == 1)
		{
			sigma_y         = sigma_max_1*pow( (nPML_y_1 - j - 0.5)/nPML_y_1 ,exponent);
			ka_y            = 1.0 + (ka_max - 1.0)*pow( (nPML_y_1 - j - 0.5)/nPML_y_1 ,exponent);
		
			jj = cik;
			K_Ey_c[jj]     = 2.0*eps_0*ka_y + sigma_y*dt;
				
			K_Ey_d[jj]     = -(2.0*eps_0*ka_y - sigma_y*dt);
			
			K_Bx_a[jj]     = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
			
			K_Bx_b[jj]     = (2.0*eps_0*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
			
			K_Hz_a[jj]     = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
			
			K_Hz_b[jj]     = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);
		}

		if (jel2 == 1)
		{
			sigma_y         = sigma_max_2*pow( (nPML_y_2 - j - 0.5)/nPML_y_2 ,exponent);
			ka_y            = 1.0 + (ka_max - 1.0)*pow( (nPML_y_2 - j - 0.5)/nPML_y_2 ,exponent);

			jj = nly-cik-2;
			if (jj >=0 )
			{
				K_Ey_c[jj]     = 2.0*eps_0*ka_y + sigma_y*dt;
				
				K_Ey_d[jj]     = -(2.0*eps_0*ka_y - sigma_y*dt);
				
				K_Bx_a[jj]     = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
				
				K_Bx_b[jj]     = (2.0*eps_0*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
				
				K_Hz_a[jj]     = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
				
				K_Hz_b[jj]     = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);
			}
		}
		cik++;
	}

	if ( (jend < ny-1) && (jend >= ny - nPML_y_2 - 1) )
	{
		j = n1;
		if (jj < 0)
		{
		  j--;
		}

		sigma_y         = sigma_max_2*pow( (nPML_y_2 - j - 0.5)/nPML_y_2 ,exponent);
		ka_y            = 1.0 + (ka_max - 1.0)*pow( (nPML_y_2 - j - 0.5)/nPML_y_2 ,exponent);

		K_Ey_c[nlyMIN1]     = 2.0*eps_0*ka_y + sigma_y*dt;
			
		K_Ey_d[nlyMIN1]     = -(2.0*eps_0*ka_y - sigma_y*dt);
		
		K_Bx_a[nlyMIN1]     = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
		
		K_Bx_b[nlyMIN1]     = (2.0*eps_0*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
		
		K_Hz_a[nlyMIN1]     = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
		
		K_Hz_b[nlyMIN1]     = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);
	}
	
	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices in z directions
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_3D::Init_PML_Par_z(double eps_r_z_1, double mu_r_z_1, double eps_r_z_2, 
					         double mu_r_z_2)
{
	long  k, kk; 
	long nlz = kend - ksta + 1;
	long nlzMIN1 = nlz -1;

	//Outside of the PML region
	//k
	for (k = 0; k < nlzMIN1; k++)
	{
		K_Ex_a[k] = 1.0;
		K_Ex_b[k] = 1.0/(2.0*eps_0);

		K_Gy_a[k] = 1.0;
		K_Gy_b[k] = 1.0;

		K_Ez_c[k] =  2.0*eps_0;
		K_Ez_d[k] = -2.0*eps_0;

		K_Hx_a[k] = 1.0;
		K_Hx_b[k] = 1.0/(2.0*eps_0);

		K_By_a[k] = 1.0;
		K_By_b[k] = dt;

		K_Hz_c[k] =  2.0*eps_0/mu_0;
		K_Hz_d[k] = -2.0*eps_0/mu_0;
	}
	K_Ex_a[nlzMIN1] =  1.0;
	K_Ex_b[nlzMIN1] =  1.0/(2.0*eps_0);
	K_Gy_a[nlzMIN1] =  1.0;
	K_Gy_b[nlzMIN1] =  1.0;
	K_Hz_c[nlzMIN1] =  2.0*eps_0/mu_0;
	K_Hz_d[nlzMIN1] = -2.0*eps_0/mu_0;
	if ( kend < nz-1 )
	{
		K_Ez_c[nlzMIN1] =  2.0*eps_0;
		K_Ez_d[nlzMIN1] = -2.0*eps_0;
		K_Hx_a[nlzMIN1] = 1.0;
		K_Hx_b[nlzMIN1] = 1.0/(2.0*eps_0);
		K_By_a[nlzMIN1] = 1.0;
		K_By_b[nlzMIN1] = dt;
	}
		
	//PML_z parameters
	double ka_max = 1;
	int exponent = 4;
	double R_err = 1e-16;

	double eta_1 = sqrt(mu_0*mu_r_z_1/eps_0/eps_r_z_1);
	double eta_2 = sqrt(mu_0*mu_r_z_2/eps_0/eps_r_z_2);
	
	double sigma_z, ka_z;
	double sigma_max_1= -(exponent+1.0)*log(R_err)/(2.0*eta_1*nPML_z_1*dz);
	double sigma_max_2= -(exponent+1.0)*log(R_err)/(2.0*eta_2*nPML_z_2*dz);
	
	long n1 = 1, n2 = 0, jel1 = 0, jel2 = 0;
	if (kend <= nPML_z_1-1)
	{
		n1 = ksta;
		n2 = kend;
		jel1 = 1;
	}
	if ( (ksta <= nPML_z_1-1) && (kend >= nPML_z_1-1) )
	{
		n1 = ksta;
		n2 = nPML_z_1-1;
		jel1 = 1;
	}
	if ( (ksta <= nz - nPML_z_2) && (kend > nz - nPML_z_2) )
	{
		n1 = ( nz - kend - 1);
		n2 = nPML_z_2-1;
		jel2 = 1;
	}
	if ( ksta > nz - nPML_z_2)
	{
		n1 = 0;
		n2 = nlzMIN1;
		jel2 = 1;
	}

	long cik = 0;
	for (k = n1; k <= n2; k++)
	{
		//k
		if (jel1 == 1)
		{
			sigma_z         = sigma_max_1*pow( (nPML_z_1 - k)/((double) nPML_z_1) ,exponent);
			ka_z            = 1.0 + (ka_max - 1.0)*pow( (nPML_z_1 - k)/((double) nPML_z_1) ,exponent);

			kk = cik;

			K_Ex_a[kk]     = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
		
			K_Ex_b[kk]     = 1.0/(2.0*eps_0*ka_z + sigma_z*dt);
			
			K_Gy_a[kk]     = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
			
			K_Gy_b[kk]     = 2.0*eps_0/(2.0*eps_0*ka_z + sigma_z*dt);
			
			K_Hz_c[kk]     = (2.0*eps_0*ka_z + sigma_z*dt)/mu_0;
			
			K_Hz_d[kk]     = -(2.0*eps_0*ka_z - sigma_z*dt)/mu_0;
		}

		if (jel2 == 1)
		{
			sigma_z         = sigma_max_2*pow( (nPML_z_2 - k)/((double) nPML_z_2) ,exponent);
			ka_z            = 1.0 + (ka_max - 1.0)*pow( (nPML_z_2 - k)/((double) nPML_z_2) ,exponent);

			kk = nlz - cik - 1;

			K_Ex_a[kk]     = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
		
			K_Ex_b[kk]     = 1.0/(2.0*eps_0*ka_z + sigma_z*dt);
			
			K_Gy_a[kk]     = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
			
			K_Gy_b[kk]     = 2.0*eps_0/(2.0*eps_0*ka_z + sigma_z*dt);
			
			K_Hz_c[kk]     = (2.0*eps_0*ka_z + sigma_z*dt)/mu_0;
			
			K_Hz_d[kk]     = -(2.0*eps_0*ka_z - sigma_z*dt)/mu_0;
		}
		
		//k+0.5
		if (jel1 == 1)
		{
			sigma_z         = sigma_max_1*pow( (nPML_z_1 - k - 0.5)/nPML_z_1 ,exponent);
			ka_z            = 1.0 + (ka_max - 1.0)*pow( (nPML_z_1 - k - 0.5)/nPML_z_1 ,exponent);
	
			kk = cik;
			K_Ez_c[kk]     = 2.0*eps_0*ka_z + sigma_z*dt;
		
			K_Ez_d[kk]     = -(2.0*eps_0*ka_z - sigma_z*dt);
			
			K_Hx_a[kk]     = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
			
			K_Hx_b[kk]     = 1.0/(2.0*eps_0*ka_z + sigma_z*dt);
	        
			K_By_a[kk]     = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
			
			K_By_b[kk]     = (2.0*eps_0*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
		}

		if (jel2 == 1)
		{
			sigma_z         = sigma_max_2*pow( (nPML_z_2 - k - 0.5)/nPML_z_2 ,exponent);
			ka_z            = 1.0 + (ka_max - 1.0)*pow( (nPML_z_2 - k - 0.5)/nPML_z_2 ,exponent);

			kk = nlz - cik - 2;
			if (kk >= 0 )
			{
				K_Ez_c[kk]     = 2.0*eps_0*ka_z + sigma_z*dt;
				
				K_Ez_d[kk]     = -(2.0*eps_0*ka_z - sigma_z*dt);
				
				K_Hx_a[kk]     = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
				
				K_Hx_b[kk]     = 1.0/(2.0*eps_0*ka_z + sigma_z*dt);
		        
				K_By_a[kk]     = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
				
				K_By_b[kk]     = (2.0*eps_0*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
			}
		}
		cik++;
	}

	if ( (kend < nz-1) && (kend >= nz - nPML_z_2 - 1) )
	{
		k = n1;
		if (kk < 0)
		{
		  k--;
		}

		sigma_z         = sigma_max_2*pow( (nPML_z_2 - k - 0.5)/nPML_z_2 ,exponent);
		ka_z            = 1.0 + (ka_max - 1.0)*pow( (nPML_z_2 - k - 0.5)/nPML_z_2 ,exponent);

		K_Ez_c[nlzMIN1]     = 2.0*eps_0*ka_z + sigma_z*dt;
		
		K_Ez_d[nlzMIN1]     = -(2.0*eps_0*ka_z - sigma_z*dt);
		
		K_Hx_a[nlzMIN1]     = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
		
		K_Hx_b[nlzMIN1]     = 1.0/(2.0*eps_0*ka_z + sigma_z*dt);
		
		K_By_a[nlzMIN1]     = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
		
		K_By_b[nlzMIN1]     = (2.0*eps_0*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
	}

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate the Ex field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D::Calc_Ex(long  nlx, long  nly, long  nlz)
{
	double Fx_r, Gx_r;
	long  i, j, k;

	if (myrank_j == jprocsMIN1)
	{
		nly--;
	}
	if (myrank_k == kprocsMIN1)
	{
		nlz--;
	}
	
	#pragma omp parallel default(shared) private(i,j,k,Fx_r,Gx_r)
	{
		j = 0;
		k = 0;
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (i = 0; i < nlx; i++)