fdtd_3d_lorentz.cpp

三维FDTD

				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;
				}
			}
		}

		i = 0;
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (j = 0; j < nly; j++)
		{	
			for (k = 1; k < nlz; k++)
			{
				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[i][j][k-1])*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;
				}
			}
		}

		#pragma omp for schedule(dynamic,nr_threads)
		for (i = 1; i < nlx; i++)
		{	
			for (j = 0; j < nly; j++)
			{	
				for (k = 1; k < nlz; k++)
				{
					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[i][j][k-1])*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);

					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;
					}
				}
			}
		}
	}

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

///////////////////////////////////////////////////////////////////////////////
//Update the Ey field components which will be send to other processes
///////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_LORENTZ::Update_Ey_send()
{
	long i, j, k;
	//update Ey_send_i
	for (j = 0; j < nly_Ey; j++)
	{
		for (k = 0; k < nlz_Ey; k++)
		{
			Ey_send_i[j][k] = Ey[0][j][k];
		}
	}
	//update Ey_send_k
	for (i = 0; i < nlx_Ey; i++)
	{
		for (j = 0; j < nly_Ey; j++)
		{
			Ey_send_k[i][j] = Ey[i][j][0];
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate the Ez field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_LORENTZ::Calc_Ez(long  nlx, long  nly, long  nlz)
{
	long  i, j, k, kk, Nr_Lorentz;
	double Gz_r, Sz_r, Sum_S;
	double eps_0; 

	if (myrank_i == iprocs-1)
	{
		nlx--;
	}
	if (myrank_j == jprocs-1)
	{
		nly--;
	}
	
	#pragma omp parallel default(shared) private(i,j,k,kk,Gz_r,Sz_r,Nr_Lorentz,Sum_S,eps_0)
	{
	    eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]

		i = 0;
		j = 0;
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (k = 0; k < nlz; k++)
		{
			Gz_r = Gz[i][j][k];

			Gz[i][j][k] = K_Gz_a[i]*Gz[i][j][k] + 
						  K_Gz_b[i]*( (Hy[i][j][k] - Hy_recv_i[j][k])*inv_dx -
									  (Hx[i][j][k] - Hx_recv_j[i][k])*inv_dy );
		
			//Total field scattered field formulation
			if (jel_plane_wave == 1 && jel_TS == 1) 
			{
				//i0 face
				if(jel_TS_planes[0] == 1)
					TS_Gz_i0(i,j,k);
				//i1 face
				if(jel_TS_planes[1] == 1)
					TS_Gz_i1(i,j,k);
				//j0 face
				if(jel_TS_planes[2] == 1)
					TS_Gz_j0(i,j,k);
				//j1 face
				if(jel_TS_planes[3] == 1)
					TS_Gz_j1(i,j,k);
			}

			Dz[i][j][k] = K_Dz_a[j]*(K_Dz_b[j]*Dz[i][j][k] + 
						             K_Dz_c[k]*Gz[i][j][k] + K_Dz_d[k]*Gz_r);

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

			for (kk = 0; kk < Nr_Lorentz; kk++)
			{
				Sz_r = Sz[i][j][k][kk]; 
				Sz[i][j][k][kk]   = K_a[Ind[i][j][k]][kk] * Sz_r +
				 					K_b[Ind[i][j][k]][kk] * Sz_2[i][j][k][kk] + 
									K_c[Ind[i][j][k]][kk] * Ez[i][j][k];
				Sz_2[i][j][k][kk] = Sz_r;
			}
		}
	

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

				Gz[i][j][k] = K_Gz_a[i]*Gz[i][j][k] + 
						      K_Gz_b[i]*( (Hy[i][j][k] - Hy_recv_i[j][k])*inv_dx -
										  (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );
			
				//Total field scattered field formulation
				if (jel_plane_wave == 1 && jel_TS == 1) 
				{
					//i0 face
					if(jel_TS_planes[0] == 1)
						TS_Gz_i0(i,j,k);
					//i1 face
					if(jel_TS_planes[1] == 1)
						TS_Gz_i1(i,j,k);
					//j0 face
					if(jel_TS_planes[2] == 1)
						TS_Gz_j0(i,j,k);
					//j1 face
					if(jel_TS_planes[3] == 1)
						TS_Gz_j1(i,j,k);
				}

				Dz[i][j][k] = K_Dz_a[j]*(K_Dz_b[j]*Dz[i][j][k] + 
						                 K_Dz_c[k]*Gz[i][j][k] + K_Dz_d[k]*Gz_r);

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

				for (kk = 0; kk < Nr_Lorentz; kk++)
				{
					Sz_r = Sz[i][j][k][kk]; 
					Sz[i][j][k][kk]   = K_a[Ind[i][j][k]][kk] * Sz_r +
				 					    K_b[Ind[i][j][k]][kk] * Sz_2[i][j][k][kk] + 
									    K_c[Ind[i][j][k]][kk] * Ez[i][j][k];
					Sz_2[i][j][k][kk] = Sz_r;
				}
			}
		}

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

				Gz[i][j][k] = K_Gz_a[i]*Gz[i][j][k] + 
						      K_Gz_b[i]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
										  (Hx[i][j][k] - Hx_recv_j[i][k])*inv_dy );
			
				//Total field scattered field formulation
				if (jel_plane_wave == 1 && jel_TS == 1) 
				{
					//i0 face
					if(jel_TS_planes[0] == 1)
						TS_Gz_i0(i,j,k);
					//i1 face
					if(jel_TS_planes[1] == 1)
						TS_Gz_i1(i,j,k);
					//j0 face
					if(jel_TS_planes[2] == 1)
						TS_Gz_j0(i,j,k);
					//j1 face
					if(jel_TS_planes[3] == 1)
						TS_Gz_j1(i,j,k);
				}

				Dz[i][j][k] = K_Dz_a[j]*(K_Dz_b[j]*Dz[i][j][k] + 
						                 K_Dz_c[k]*Gz[i][j][k] + K_Dz_d[k]*Gz_r);

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

				for (kk = 0; kk < Nr_Lorentz; kk++)
				{
					Sz_r = Sz[i][j][k][kk]; 
					Sz[i][j][k][kk]   = K_a[Ind[i][j][k]][kk] * Sz_r +
				 					    K_b[Ind[i][j][k]][kk] * Sz_2[i][j][k][kk] + 
									    K_c[Ind[i][j][k]][kk] * Ez[i][j][k];
					Sz_2[i][j][k][kk] = Sz_r;
				}
			}
		}

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

					Gz[i][j][k] = K_Gz_a[i]*Gz[i][j][k] + 
						          K_Gz_b[i]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
											  (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );
				
					//Total field scattered field formulation
					if (jel_plane_wave == 1 && jel_TS == 1) 
					{
						//i0 face
						if(jel_TS_planes[0] == 1)
							TS_Gz_i0(i,j,k);
						//i1 face
						if(jel_TS_planes[1] == 1)
							TS_Gz_i1(i,j,k);
						//j0 face
						if(jel_TS_planes[2] == 1)
							TS_Gz_j0(i,j,k);
						//j1 face
						if(jel_TS_planes[3] == 1)
							TS_Gz_j1(i,j,k);
					}

					Dz[i][j][k] = K_Dz_a[j]*(K_Dz_b[j]*Dz[i][j][k] + 
						                     K_Dz_c[k]*Gz[i][j][k] + K_Dz_d[k]*Gz_r);

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

					for (kk = 0; kk < Nr_Lorentz; kk++)
					{
						Sz_r = Sz[i][j][k][kk]; 
						Sz[i][j][k][kk]   = K_a[Ind[i][j][k]][kk] * Sz_r +
				 					        K_b[Ind[i][j][k]][kk] * Sz_2[i][j][k][kk] + 
									        K_c[Ind[i][j][k]][kk] * Ez[i][j][k];
					    Sz_2[i][j][k][kk] = Sz_r;
					}
				}
			}
		}
	}

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

///////////////////////////////////////////////////////////////////////////////
//Update the Ez field components which will be send to other processes
///////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_LORENTZ::Update_Ez_send()
{
	long i, j, k;

	//update Ez_send_i
	for (j = 0; j < nly_Ez; j++)
	{
		for (k = 0; k < nlz_Ez; k++)
		{
			Ez_send_i[j][k] = Ez[0][j][k];
		}
	}
	//update Ez_send_j
	for (i = 0; i < nlx_Ez; i++)
	{
		for (k = 0; k < nlz_Ez; k++)
		{
			Ez_send_j[i][k] = Ez[i][0][k];
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate the Hx field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_LORENTZ::Calc_Hx(long  nlx, long  nly, long  nlz)
{
	long  i, j, k;
	double Bx_r;

	long nlyMIN1 = nly - 1;
	long nlzMIN1 = nlz - 1;

	#pragma omp parallel default(shared) private(i,j,k,Bx_r)
	{
		j = nlyMIN1;
		k = nlzMIN1;
		if ( nly > 0 && nlz >0 )
		{
			#pragma omp for schedule(dynamic,nr_threads) nowait
			for (i = 0; i < nlx; i++)
			{	
				Bx_r = Bx[i][j][k];

				Bx[i][j][k] = K_Bx_a[j]*Bx[i][j][k] + 
							  K_Bx_b[j]*( (Ey_recv_k[i][j] - Ey[i][j][k])*inv_dz -
										  (Ez_recv_j[i][k] - Ez[i][j][k])*inv_dy);

				Hx[i][j][k] = K_Hx_a[k]*Hx[i][j][k] + 
							  K_Hx_b[k]*( K_Hx_c[i]*Bx[i][j][k] + 
										  K_Hx_d[i]*Bx_r )/mu_r[Ind[i][j][k]];
			}
		}

		k = nlzMIN1;
		if ( nlz > 0 )
		{
			#pragma omp for schedule(dynamic,nr_threads) nowait
			for (i = 0; i < nlx; i++)
			{	
				for (j = 0; j < nlyMIN1; j++)
				{	
					Bx_r = Bx[i][j][k];

					Bx[i][j][k] = K_Bx_a[j]*Bx[i][j][k] + 
								  K_Bx_b[j]*( (Ey_recv_k[i][j] - Ey[i][j][k])*inv_dz -
											  (Ez[i][j+1][k] - Ez[i][j][k])*inv_dy);

					Hx[i][j][k] = K_Hx_a[k]*Hx[i][j][k] + 
								  K_Hx_b[k]*( K_Hx_c[i]*Bx[i][j][k] + 
											  K_Hx_d[i]*Bx_r )/mu_r[Ind[i][j][k]];
				}
			}
		}

		j = nlyMIN1;
		if ( nly > 0 )
		{
			#pragma omp for schedule(dynamic,nr_threads) nowait
			for (i = 0; i < nlx; i++)
			{	
				for (k = 0; k < nlzMIN1; k++)
				{
					Bx_r = Bx[i][j][k];

					Bx[i][j][k] = K_Bx_a[j]*Bx[i][j][k] + 
								  K_Bx_b[j]*( (Ey[i][j][k+1] - Ey[i][j][k])*inv_dz -
											  (Ez_recv_j[i][k] - Ez[i][j][k])*inv_dy);

					Hx[i][j][k] = K_Hx_a[k]*Hx[i][j][k] + 
								  K_Hx_b[k]*( K_Hx_c[i]*Bx[i][j][k] + 
											  K_Hx_d[i]*Bx_r )/mu_r[Ind