fdtd_3d.cpp

fdtd 3D xyzPML MPI OpenMP

		{	
			Fx_r = Fx[i][j][k];
			
			Fx[i][j][k] = K_a[Ind[i][j][k]]*Fx[i][j][k] +
						  K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz_recv_j[i][k])*inv_dy - 
											  (Hy[i][j][k] - Hy_recv_k[i][j])*inv_dz );

			Gx_r = Gx[i][j][k];

			Gx[i][j][k] = K_Gx_a[j]*Gx[i][j][k] + K_Gx_b[j]*(Fx[i][j][k] - Fx_r);

			Ex[i][j][k] = K_Ex_a[k]*Ex[i][j][k] + 
						  K_Ex_b[k]*( K_Ex_c[i]*Gx[i][j][k] + K_Ex_d[i]*Gx_r);
		}

		j = 0;
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (i = 0; i < nlx; i++)
		{	
			for (k = 1; k < nlz; k++)
			{
				Fx_r = Fx[i][j][k];
				
				Fx[i][j][k] = K_a[Ind[i][j][k]]*Fx[i][j][k] +
							  K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz_recv_j[i][k])*inv_dy - 
												  (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx[i][j][k];

				Gx[i][j][k] = K_Gx_a[j]*Gx[i][j][k] + K_Gx_b[j]*(Fx[i][j][k] - Fx_r);

				Ex[i][j][k] = K_Ex_a[k]*Ex[i][j][k] + 
							  K_Ex_b[k]*( K_Ex_c[i]*Gx[i][j][k] + K_Ex_d[i]*Gx_r);
			}
		}
		
		k = 0;
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (i = 0; i < nlx; i++)
		{	
			for (j = 1; j < nly; j++)
			{	
				Fx_r = Fx[i][j][k];
				
				Fx[i][j][k] = K_a[Ind[i][j][k]]*Fx[i][j][k] +
							  K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
												  (Hy[i][j][k] - Hy_recv_k[i][j])*inv_dz );

				Gx_r = Gx[i][j][k];

				Gx[i][j][k] = K_Gx_a[j]*Gx[i][j][k] + K_Gx_b[j]*(Fx[i][j][k] - Fx_r);

				Ex[i][j][k] = K_Ex_a[k]*Ex[i][j][k] + 
							  K_Ex_b[k]*( K_Ex_c[i]*Gx[i][j][k] + K_Ex_d[i]*Gx_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++)
				{
					Fx_r = Fx[i][j][k];
					
					Fx[i][j][k] = K_a[Ind[i][j][k]]*Fx[i][j][k] +
								  K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
													  (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

					Gx_r = Gx[i][j][k];

					Gx[i][j][k] = K_Gx_a[j]*Gx[i][j][k] + 
								  K_Gx_b[j]*(Fx[i][j][k] - Fx_r);

					Ex[i][j][k] = K_Ex_a[k]*Ex[i][j][k] + 
								  K_Ex_b[k]*( K_Ex_c[i]*Gx[i][j][k] + K_Ex_d[i]*Gx_r);
				}
			}
		}
	}

	////////////////////////////////////////////////////////////////////
	//Total field scattered field formulation
	////////////////////////////////////////////////////////////////////
	if (jel_plane_wave == 1 && jel_TS == 1) 
	{
		//j0 face
		if(jel_TS_planes[2] == 1)
			TS_Ex_j0();
		//j1 face
		if(jel_TS_planes[3] == 1)
			TS_Ex_j1();
		//k0 face
		if(jel_TS_planes[4] == 1)
			TS_Ex_k0();
		//k1 face
		if(jel_TS_planes[5] == 1)
			TS_Ex_k1();
	}

	////////////////////////////////////////////////////////////////////
	//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::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::Calc_Ey(long  nlx, long  nly, long  nlz)
{
	long  i, j, k;
	double Fy_r, Gy_r;
	
	if (myrank_i == iprocsMIN1)
	{
		nlx--;
	}
	if (myrank_k == kprocsMIN1)
	{
		nlz--;
	}

	#pragma omp parallel default(shared) private(i,j,k,Fy_r,Gy_r)
	{
		i = 0;
		k = 0;
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (j = 0; j < nly; j++)
		{	
			Fy_r = Fy[i][j][k];
			
			Fy[i][j][k] = K_a[Ind[i][j][k]]*Fy[i][j][k] + 
						  K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx_recv_k[i][j])*inv_dz - 
											  (Hz[i][j][k] - Hz_recv_i[j][k])*inv_dx );

			Gy_r = Gy[i][j][k];

			Gy[i][j][k] = K_Gy_a[k]*Gy[i][j][k] + K_Gy_b[k]*(Fy[i][j][k] - Fy_r);

			Ey[i][j][k] = K_Ey_a[i]*Ey[i][j][k] + 
						  K_Ey_b[i]*( K_Ey_c[j]*Gy[i][j][k] + K_Ey_d[j]*Gy_r);
		}


		k = 0;
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (i = 1; i < nlx; i++)
		{	
			for (j = 0; j < nly; j++)
			{	
				Fy_r = Fy[i][j][k];
				
				Fy[i][j][k] = K_a[Ind[i][j][k]]*Fy[i][j][k] + 
							  K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx_recv_k[i][j])*inv_dz - 
												  (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

				Gy_r = Gy[i][j][k];

				Gy[i][j][k] = K_Gy_a[k]*Gy[i][j][k] + 
							  K_Gy_b[k]*(Fy[i][j][k] - Fy_r);

				Ey[i][j][k] = K_Ey_a[i]*Ey[i][j][k] + 
							  K_Ey_b[i]*( K_Ey_c[j]*Gy[i][j][k] + K_Ey_d[j]*Gy_r);
			}
		}

		i = 0;
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (j = 0; j < nly; j++)
		{	
			for (k = 1; k < nlz; k++)
			{
				Fy_r = Fy[i][j][k];
				
				Fy[i][j][k] = K_a[Ind[i][j][k]]*Fy[i][j][k] + 
							  K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
							  					  (Hz[i][j][k] - Hz_recv_i[j][k])*inv_dx );

				Gy_r = Gy[i][j][k];

				Gy[i][j][k] = K_Gy_a[k]*Gy[i][j][k] + 
							  K_Gy_b[k]*(Fy[i][j][k] - Fy_r);

				Ey[i][j][k] = K_Ey_a[i]*Ey[i][j][k] + 
							  K_Ey_b[i]*( K_Ey_c[j]*Gy[i][j][k] + K_Ey_d[j]*Gy_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++)
				{
					Fy_r = Fy[i][j][k];
					
					Fy[i][j][k] = K_a[Ind[i][j][k]]*Fy[i][j][k] + 
								  K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
								  					  (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

					Gy_r = Gy[i][j][k];

					Gy[i][j][k] = K_Gy_a[k]*Gy[i][j][k] + 
								  K_Gy_b[k]*(Fy[i][j][k] - Fy_r);

					Ey[i][j][k] = K_Ey_a[i]*Ey[i][j][k] + 
								  K_Ey_b[i]*( K_Ey_c[j]*Gy[i][j][k] + K_Ey_d[j]*Gy_r);
				}
			}
		}
	}

	////////////////////////////////////////////////////////////////////
	//Total field scattered field formulation
	////////////////////////////////////////////////////////////////////
	if (jel_plane_wave == 1 && jel_TS == 1) 
	{
		//i0 face
		if(jel_TS_planes[0] == 1)
			TS_Ey_i0();
		//i1 face
		if(jel_TS_planes[1] == 1)
			TS_Ey_i1();
		//k0 face
		if(jel_TS_planes[4] == 1)
			TS_Ey_k0();
		//k1 face
		if(jel_TS_planes[5] == 1)
			TS_Ey_k1();
	}

	////////////////////////////////////////////////////////////////////
	//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::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::Calc_Ez(long  nlx, long  nly, long  nlz)
{
	long  i, j, k;
	double Fz_r, Gz_r;

	if (myrank_i == iprocsMIN1)
	{
		nlx--;
	}
	if (myrank_j == jprocsMIN1)
	{
		nly--;
	}
	
	#pragma omp parallel default(shared) private(i,j,k,Fz_r,Gz_r)
	{
		i = 0;
		j = 0;
		#pragma omp for schedule(dynamic,nr_threads) nowait
		for (k = 0; k < nlz; k++)
		{
			Fz_r = Fz[i][j][k];

			Fz[i][j][k] = K_a[Ind[i][j][k]]*Fz[i][j][k] +
						  K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy_recv_i[j][k])*inv_dx -
											  (Hx[i][j][k] - Hx_recv_j[i][k])*inv_dy );

			Gz_r = Gz[i][j][k];

			Gz[i][j][k] = K_Gz_a[i]*Gz[i][j][k] + K_Gz_b[i]*(Fz[i][j][k] - Fz_r);

			Ez[i][j][k] = K_Ez_a[j]*Ez[i][j][k] + 
						  K_Ez_b[j]*( K_Ez_c[k]*Gz[i][j][k] + K_Ez_d[k]*Gz_r);
		}
	

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

				Fz[i][j][k] = K_a[Ind[i][j][k]]*Fz[i][j][k] +
							  K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy_recv_i[j][k])*inv_dx -
												  (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

				Gz_r = Gz[i][j][k];

				Gz[i][j][k] = K_Gz_a[i]*Gz[i][j][k] + K_Gz_b[i]*(Fz[i][j][k] - Fz_r);

				Ez[i][j][k] = K_Ez_a[j]*Ez[i][j][k] + 
							  K_Ez_b[j]*( K_Ez_c[k]*Gz[i][j][k] + K_Ez_d[k]*Gz_r);
			}
		}

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

				Fz[i][j][k] = K_a[Ind[i][j][k]]*Fz[i][j][k] +
							K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
												(Hx[i][j][k] - Hx_recv_j[i][k])*inv_dy );

				Gz_r = Gz[i][j][k];

				Gz[i][j][k] = K_Gz_a[i]*Gz[i][j][k] + K_Gz_b[i]*(Fz[i][j][k] - Fz_r);

				Ez[i][j][k] = K_Ez_a[j]*Ez[i][j][k] + 
							K_Ez_b[j]*( K_Ez_c[k]*Gz[i][j][k] + K_Ez_d[k]*Gz_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++)
				{
					Fz_r = Fz[i][j][k];

					Fz[i][j][k] = K_a[Ind[i][j][k]]*Fz[i][j][k] +
								  K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
													  (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

					Gz_r = Gz[i][j][k];

					Gz[i][j][k] = K_Gz_a[i]*Gz[i][j][k] + 
								  K_Gz_b[i]*(Fz[i][j][k] - Fz_r);

					Ez[i][j][k] = K_Ez_a[j]*Ez[i][j][k] + 
								  K_Ez_b[j]*( K_Ez_c[k]*Gz[i][j][k] + K_Ez_d[k]*Gz_r);
				}
			}
		}
	}

	////////////////////////////////////////////////////////////////////
	//Total field scattered field formulation
	////////////////////////////////////////////////////////////////////
	if (jel_plane_wave == 1 && jel_TS == 1) 
	{
		//i0 face
		if(jel_TS_planes[0] == 1)
			TS_Ez_i0();
		//i1 face
		if(jel_TS_planes[1] == 1)
			TS_Ez_i1();
		//j0 face
		if(jel_TS_planes[2] == 1)
			TS_Ez_j0();
		//j1 face
		if(jel_TS_planes[3] == 1)
			TS_Ez_j1();
	}

	////////////////////////////////////////////////////////////////////
	//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::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::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] +