fdtd_3d_lorentz.cpp

三维FDTD

		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 = nlz-1;
		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_Dx_a[kk]     = 1/(2.0*eps_0*ka_z + sigma_z*dt);
		
			K_Dx_b[kk]     = 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*dt/(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_Dx_a[kk]     = 1/(2.0*eps_0*ka_z + sigma_z*dt);
		
			K_Dx_b[kk]     = 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*dt/(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_Dz_c[kk]     = 2.0*eps_0*ka_z + sigma_z*dt;
		
			K_Dz_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_Dz_c[kk]     = 2.0*eps_0*ka_z + sigma_z*dt;
				
				K_Dz_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_Dz_c[nlz-1]     = 2.0*eps_0*ka_z + sigma_z*dt;
		
		K_Dz_d[nlz-1]     = -2.0*eps_0*ka_z + sigma_z*dt;
		
		K_Hx_a[nlz-1]     = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
		
		K_Hx_b[nlz-1]     = 1.0/(2.0*eps_0*ka_z + sigma_z*dt);
		
		K_By_a[nlz-1]     = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
		
		K_By_b[nlz-1]     = (2.0*eps_0*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
	}

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate the Ex field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_LORENTZ::Calc_Ex(long  nlx, long  nly, long  nlz)
{
	long  i, j, k, kk, Nr_Lorentz;
	double Gx_r, Sx_r, Sum_S;
	double eps_0; 

	if (myrank_j == jprocs-1)
	{
		nly--;
	}
	if (myrank_k == kprocs-1)
	{
		nlz--;
	}
	
	#pragma omp parallel default(shared) private(i,j,k,kk,Gx_r,Sx_r,Nr_Lorentz,Sum_S,eps_0)
	{
	    eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]

		j = 0;
		k = 0;
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (i = 0; i < nlx; i++)
		{	
			Gx_r = Gx[i][j][k];

			Gx[i][j][k] = K_Gx_a[j]*Gx[i][j][k] + 
						  K_Gx_b[j]*( (Hz[i][j][k] - Hz_recv_j[i][k])*inv_dy - 
									  (Hy[i][j][k] - Hy_recv_k[i][j])*inv_dz );
		
			//Total field scattered field formulation
			if (jel_plane_wave == 1 && jel_TS == 1) 
			{
				//j0 face
				if(jel_TS_planes[2] == 1)
					TS_Gx_j0(i,j,k);
				//j1 face
				if(jel_TS_planes[3] == 1)
					TS_Gx_j1(i,j,k);
				//k0 face
				if(jel_TS_planes[4] == 1)
					TS_Gx_k0(i,j,k);
				//k1 face
				if(jel_TS_planes[5] == 1)
					TS_Gx_k1(i,j,k);
			}

			Dx[i][j][k] = K_Dx_a[k]*(K_Dx_b[k]*Dx[i][j][k] + 
						             K_Dx_c[i]*Gx[i][j][k] + K_Dx_d[i]*Gx_r);

			Nr_Lorentz = (long) Mat[Ind[i][j][k]][0];
			Sum_S = 0;
			for (kk = 0; kk < Nr_Lorentz; kk++)
			{
				Sum_S = Sum_S + Sx[i][j][k][kk];
			}
			Ex[i][j][k] = (Dx[i][j][k] - eps_0*Sum_S)/(eps_0*Mat[Ind[i][j][k]][1]);

			for (kk = 0; kk < Nr_Lorentz; kk++)
			{
				Sx_r = Sx[i][j][k][kk]; 
				Sx[i][j][k][kk]   = K_a[Ind[i][j][k]][kk] * Sx_r +
				 					K_b[Ind[i][j][k]][kk] * Sx_2[i][j][k][kk] + 
									K_c[Ind[i][j][k]][kk] * Ex[i][j][k];
				Sx_2[i][j][k][kk] = Sx_r;
			}
		}

		j = 0;
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (i = 0; i < nlx; i++)
		{	
			for (k = 1; k < nlz; k++)
			{		
				Gx_r = Gx[i][j][k];

				Gx[i][j][k] = K_Gx_a[j]*Gx[i][j][k] + 
						      K_Gx_b[j]*( (Hz[i][j][k] - Hz_recv_j[i][k])*inv_dy - 
										  (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );
			
				//Total field scattered field formulation
				if (jel_plane_wave == 1 && jel_TS == 1) 
				{
					//j0 face
					if(jel_TS_planes[2] == 1)
						TS_Gx_j0(i,j,k);
					//j1 face
					if(jel_TS_planes[3] == 1)
						TS_Gx_j1(i,j,k);
					//k0 face
					if(jel_TS_planes[4] == 1)
						TS_Gx_k0(i,j,k);
					//k1 face
					if(jel_TS_planes[5] == 1)
						TS_Gx_k1(i,j,k);
				}

				Dx[i][j][k] = K_Dx_a[k]*(K_Dx_b[k]*Dx[i][j][k] + 
						                 K_Dx_c[i]*Gx[i][j][k] + K_Dx_d[i]*Gx_r);

				Nr_Lorentz = (long) Mat[Ind[i][j][k]][0];
				Sum_S = 0;
				for (kk = 0; kk < Nr_Lorentz; kk++)
				{
					Sum_S = Sum_S + Sx[i][j][k][kk];
				}
				Ex[i][j][k] = (Dx[i][j][k] - eps_0*Sum_S)/(eps_0*Mat[Ind[i][j][k]][1]);

				for (kk = 0; kk < Nr_Lorentz; kk++)
				{
					Sx_r = Sx[i][j][k][kk]; 
					Sx[i][j][k][kk]   = K_a[Ind[i][j][k]][kk] * Sx_r +
				 					    K_b[Ind[i][j][k]][kk] * Sx_2[i][j][k][kk] + 
									    K_c[Ind[i][j][k]][kk] * Ex[i][j][k];
					Sx_2[i][j][k][kk] = Sx_r;
				}
			}
		}
		
		k = 0;
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (i = 0; i < nlx; i++)
		{	
			for (j = 1; j < nly; j++)
			{	
				Gx_r = Gx[i][j][k];

				Gx[i][j][k] = K_Gx_a[j]*Gx[i][j][k] + 
						        K_Gx_b[j]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
										    (Hy[i][j][k] - Hy_recv_k[i][j])*inv_dz );
			
				//Total field scattered field formulation
				if (jel_plane_wave == 1 && jel_TS == 1) 
				{
					//j0 face
					if(jel_TS_planes[2] == 1)
						TS_Gx_j0(i,j,k);
					//j1 face
					if(jel_TS_planes[3] == 1)
						TS_Gx_j1(i,j,k);
					//k0 face
					if(jel_TS_planes[4] == 1)
						TS_Gx_k0(i,j,k);
					//k1 face
					if(jel_TS_planes[5] == 1)
						TS_Gx_k1(i,j,k);
				}

				Dx[i][j][k] = K_Dx_a[k]*(K_Dx_b[k]*Dx[i][j][k] +
						                 K_Dx_c[i]*Gx[i][j][k] + K_Dx_d[i]*Gx_r);

				Nr_Lorentz = (long) Mat[Ind[i][j][k]][0];
				Sum_S = 0;
				for (kk = 0; kk<Nr_Lorentz; kk++)
				{
					Sum_S = Sum_S + Sx[i][j][k][kk];
				}
				Ex[i][j][k] = (Dx[i][j][k] - eps_0*Sum_S)/(eps_0*Mat[Ind[i][j][k]][1]);

				for (kk = 0; kk < Nr_Lorentz; kk++)
				{
					Sx_r = Sx[i][j][k][kk]; 
					Sx[i][j][k][kk]   = K_a[Ind[i][j][k]][kk] * Sx_r +
				 					    K_b[Ind[i][j][k]][kk] * Sx_2[i][j][k][kk] + 
									    K_c[Ind[i][j][k]][kk] * Ex[i][j][k];
					Sx_2[i][j][k][kk] = Sx_r;
				}
			}
		}
	
		#pragma omp for schedule(dynamic,nr_threads)
		for (i = 0; i < nlx; i++)
		{	
			for (j = 1; j < nly; j++)
			{	
				for (k = 1; k < nlz; k++)
				{		
					Gx_r = Gx[i][j][k];

					Gx[i][j][k] = K_Gx_a[j]*Gx[i][j][k] + 
						          K_Gx_b[j]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
											  (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );
				
					//Total field scattered field formulation
					if (jel_plane_wave == 1 && jel_TS == 1) 
					{
						//j0 face
						if(jel_TS_planes[2] == 1)
							TS_Gx_j0(i,j,k);
						//j1 face
						if(jel_TS_planes[3] == 1)
							TS_Gx_j1(i,j,k);
						//k0 face
						if(jel_TS_planes[4] == 1)
							TS_Gx_k0(i,j,k);
						//k1 face
						if(jel_TS_planes[5] == 1)
							TS_Gx_k1(i,j,k);
					}

					Dx[i][j][k] = K_Dx_a[k]*(K_Dx_b[k]*Dx[i][j][k] + 
						                     K_Dx_c[i]*Gx[i][j][k] + K_Dx_d[i]*Gx_r);

					Nr_Lorentz = (long) Mat[Ind[i][j][k]][0];
					Sum_S = 0;
					for (kk = 0; kk < Nr_Lorentz; kk++)
					{
						Sum_S = Sum_S + Sx[i][j][k][kk];
					}
					Ex[i][j][k] = (Dx[i][j][k] - eps_0*Sum_S)/(eps_0*Mat[Ind[i][j][k]][1]);

					for (kk = 0; kk < Nr_Lorentz; kk++)
					{
						Sx_r = Sx[i][j][k][kk]; 
						Sx[i][j][k][kk]   = K_a[Ind[i][j][k]][kk] * Sx_r +
				 					        K_b[Ind[i][j][k]][kk] * Sx_2[i][j][k][kk] + 
									        K_c[Ind[i][j][k]][kk] * Ex[i][j][k];
					    Sx_2[i][j][k][kk] = Sx_r;
					}
				}
			}
		}
	}

	////////////////////////////////////////////////////////////////////
	//Point source
	////////////////////////////////////////////////////////////////////
	if (jel_plane_wave == 0 && pt_source_Ex == 1 && n_Coord_ptSource>0 && iter <= switch_off_time)
	{
		PtSource_J(Ex, time);
	}
}

///////////////////////////////////////////////////////////////////////////////
//Update the Ex field components which will be send to other processes
///////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_LORENTZ::Update_Ex_send()
{
	long i, j, k;
	//update Ex_send_j
	for (i = 0; i < nlx_Ex; i++)
	{
		for (k = 0; k < nlz_Ex; k++)
		{
			Ex_send_j[i][k] = Ex[i][0][k];
		}
	}
	//update  Ex_send_k
	for (i = 0; i < nlx_Ex; i++)
	{
		for (j = 0; j < nly_Ex; j++)
		{
			Ex_send_k[i][j] = Ex[i][j][0];
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate the Ey field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_LORENTZ::Calc_Ey(long  nlx, long  nly, long  nlz)
{
	long  i, j, k, kk, Nr_Lorentz;
	double Gy_r, Sy_r, Sum_S;
	double eps_0; 
	
	if (myrank_i == iprocs-1)
	{
		nlx--;
	}
	if (myrank_k == kprocs-1)
	{
		nlz--;
	}

	#pragma omp parallel default(shared) private(i,j,k,kk,Gy_r,Sy_r,Nr_Lorentz,Sum_S,eps_0)
	{
	    eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]

		i = 0;
		k = 0;
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (j = 0; j < nly; j++)
		{	
			Gy_r = Gy[i][j][k];

			Gy[i][j][k] = K_Gy_a[k]*Gy[i][j][k] + 
						  K_Gy_b[k]*( (Hx[i][j][k] - Hx_recv_k[i][j])*inv_dz - 
							  		  (Hz[i][j][k] - Hz_recv_i[j][k])*inv_dx );
		
			//Total field scattered field formulation
			if (jel_plane_wave == 1 && jel_TS == 1) 
			{
				//i0 face
				if(jel_TS_planes[0] == 1)
					TS_Gy_i0(i,j,k);
				//i1 face
				if(jel_TS_planes[1] == 1)
					TS_Gy_i1(i,j,k);
				//k0 face
				if(jel_TS_planes[4] == 1)
					TS_Gy_k0(i,j,k);
				//k1 face
				if(jel_TS_planes[5] == 1)
					TS_Gy_k1(i,j,k);
			}

			Dy[i][j][k] = K_Dy_a[i]*(K_Dy_b[i]*Dy[i][j][k] + 
						             K_Dy_c[j]*Gy[i][j][k] + K_Dy_d[j]*Gy_r);

			Nr_Lorentz = (long) Mat[Ind[i][j][k]][0];
			Sum_S = 0;
			for (kk = 0; kk < Nr_Lorentz; kk++)
			{
				Sum_S = Sum_S + Sy[i][j][k][kk];
			}
			Ey[i][j][k] = (Dy[i][j][k] - eps_0*Sum_S)/(eps_0*Mat[Ind[i][j][k]][1]);

			for (kk = 0; kk < Nr_Lorentz; kk++)
			{
				Sy_r = Sy[i][j][k][kk]; 
				Sy[i][j][k][kk]   = K_a[Ind[i][j][k]][kk] * Sy_r +
				 					K_b[Ind[i][j][k]][kk] * Sy_2[i][j][k][kk] + 
									K_c[Ind[i][j][k]][kk] * Ey[i][j][k];
				Sy_2[i][j][k][kk] = Sy_r;
			}
		}


		k = 0;
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (i = 1; i < nlx; i++)
		{	
			for (j = 0; j < nly; j++)
			{	
				Gy_r = Gy[i][j][k];

				Gy[i][j][k] = K_Gy_a[k]*Gy[i][j][k] + 
						      K_Gy_b[k]*( (Hx[i][j][k] - Hx_recv_k[i][j])*inv_dz - 
								  	      (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );
			
				//Total field scattered field formulation
				if (jel_plane_wave == 1 && jel_TS == 1) 
				{
					//i0 face
					if(jel_TS_planes[0] == 1)
						TS_Gy_i0(i,j,k);
					//i1 face
					if(jel_TS_planes[1] == 1)
						TS_Gy_i1(i,j,k);
					//k0 face
					if(jel_TS_planes[4] == 1)
						TS_Gy_k0(i,j,k);
					//k1 face
					if(jel_TS_planes[5] == 1)
						TS_Gy_k1(i,j,k);
				}

				Dy[i][j][k] = K_Dy_a[i]*(K_Dy_b[i]*Dy[i][j][k] + 
						                 K_Dy_c[j]*Gy[i][j][k] + K_Dy_d[j]*Gy_r);