fdtd_3d_lorentz.cpp

三维FDTD

			              (2 + delta_0*omega_0*dt);
			K_b[i][j] = (-2 + delta_0*omega_0*dt)/
			            ( 2 + delta_0*omega_0*dt);
			K_c[i][j] = 2*am*omega_0*omega_0*dt*dt/
						( 2 + delta_0*omega_0*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][3*((long)Mat[i][0])+2];
	}

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices in x directions
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_3D_LORENTZ::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;

	//Outside of the PML region
	//i
	for (i = 0; i < nlx-1; i++)
	{
		K_Dx_c[i] =  2.0*eps_0;
		K_Dx_d[i] = -2.0*eps_0;

		K_Dy_a[i] = 1.0/(2.0*eps_0);
		K_Dy_b[i] = 2.0*eps_0;
		
		K_Gz_a[i] = 1.0;
		K_Gz_b[i] = dt;

		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_Dy_a[nlx-1] =  1.0/(2.0*eps_0);
	K_Dy_b[nlx-1] =  2.0*eps_0;
	K_Gz_a[nlx-1] =  1.0;
	K_Gz_b[nlx-1] =  dt;
	K_Hx_c[nlx-1] =  2.0*eps_0/mu_0;
	K_Hx_d[nlx-1] = -2.0*eps_0/mu_0;
	if ( iend < nx-1 )
	{
		K_Dx_c[nlx-1] =  2.0*eps_0;
		K_Dx_d[nlx-1] = -2.0*eps_0;
		K_Hy_a[nlx-1] = 1.0;
		K_Hy_b[nlx-1] = 1.0/(2.0*eps_0);
		K_Bz_a[nlx-1] = 1.0;
		K_Bz_b[nlx-1] = 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 = nlx-1;
		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_Dy_a[ii]     = 1.0/(2.0*eps_0*ka_x + sigma_x*dt);

			K_Dy_b[ii]     = 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*dt/(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_Dy_a[ii]     = 1.0/(2.0*eps_0*ka_x + sigma_x*dt);

			K_Dy_b[ii]     = 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*dt/(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;
		}
		
		//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_Dx_c[ii]     = 2.0*eps_0*ka_x + sigma_x*dt;
				
			K_Dx_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_Dx_c[ii]     = 2.0*eps_0*ka_x + sigma_x*dt;
				
				K_Dx_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) )
	{
		cout << n1 << endl;
		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_Dx_c[nlx-1]     = 2.0*eps_0*ka_x + sigma_x*dt;
		
		K_Dx_d[nlx-1]     = -2.0*eps_0*ka_x + sigma_x*dt;
		
		K_Hy_a[nlx-1]     = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
		
		K_Hy_b[nlx-1]     = 1.0/(2.0*eps_0*ka_x+sigma_x*dt);
		
		K_Bz_a[nlx-1]     = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
		
		K_Bz_b[nlx-1]     = (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_LORENTZ::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;

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

		K_Dy_c[j] =  2.0*eps_0;
		K_Dy_d[j] = -2.0*eps_0;

		K_Dz_a[j] = 1.0/(2.0*eps_0);
		K_Dz_b[j] = 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[nly-1] =  1.0;
	K_Gx_b[nly-1] =  dt;
	K_Dz_a[nly-1] =  1.0/(2.0*eps_0);
	K_Dz_b[nly-1] =  2.0*eps_0;
	K_Hy_c[nly-1] =  2*eps_0/mu_0;
	K_Hy_d[nly-1] = -2*eps_0/mu_0;
	if ( jend < ny-1 )
	{
		K_Dy_c[nly-1] =  2.0*eps_0;
		K_Dy_d[nly-1] = -2.0*eps_0;
		K_Bx_a[nly-1] = 1.0;
		K_Bx_b[nly-1] = dt;
		K_Hz_a[nly-1] = 1.0;
		K_Hz_b[nly-1] = 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 = nly-1;
		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*dt/(2.0*eps_0*ka_y+sigma_y*dt);
			
			K_Dz_a[jj]     = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);

			K_Dz_b[jj]     = 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*dt/(2.0*eps_0*ka_y+sigma_y*dt);
			
			K_Dz_a[jj]     = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);

			K_Dz_b[jj]     = 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_Dy_c[jj]     = 2.0*eps_0*ka_y + sigma_y*dt;
				
			K_Dy_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_Dy_c[jj]     = 2.0*eps_0*ka_y + sigma_y*dt;
				
				K_Dy_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_Dy_c[nly-1]     = 2.0*eps_0*ka_y + sigma_y*dt;
			
		K_Dy_d[nly-1]     = -2.0*eps_0*ka_y + sigma_y*dt;
		
		K_Bx_a[nly-1]     = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
		
		K_Bx_b[nly-1]     = (2.0*eps_0*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
		
		K_Hz_a[nly-1]     = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
		
		K_Hz_b[nly-1]     = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);
	}
	
	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices in z directions
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_3D_LORENTZ::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;

	//Outside of the PML region
	//k
	for (k = 0; k < nlz-1; k++)
	{
		K_Dx_a[k] = 1.0/(2.0*eps_0);
		K_Dx_b[k] = 2.0*eps_0;

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

		K_Dz_c[k] =  2.0*eps_0;
		K_Dz_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_Dx_a[nlz-1] =  1.0/(2.0*eps_0);
	K_Dx_b[nlz-1] =  2.0*eps_0;
	K_Gy_a[nlz-1] =  1.0;
	K_Gy_b[nlz-1] =  dt;
	K_Hz_c[nlz-1] =  2.0*eps_0/mu_0;
	K_Hz_d[nlz-1] = -2.0*eps_0/mu_0;
	if ( kend < nz-1 )
	{
		K_Dz_c[nlz-1] =  2.0*eps_0;
		K_Dz_d[nlz-1] = -2.0*eps_0;
		K_Hx_a[nlz-1] = 1.0;
		K_Hx_b[nlz-1] = 1.0/(2.0*eps_0);
		K_By_a[nlz-1] = 1.0;
		K_By_b[nlz-1] = 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;