fdtd_3d.cpp

fdtd 3D xyzPML MPI OpenMP

								  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[i][j][k]];
				}
			}
		}

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

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

	////////////////////////////////////////////////////////////////////
	//Total field scattered field formulation
	////////////////////////////////////////////////////////////////////
	if (jel_plane_wave == 1 && jel_TS == 1) 
	{
		//the parallelization can be done this way because in case of 
		//magnetic field there aro no overlapping elements on the interface
		#pragma omp parallel sections default(shared)
		{
			#pragma omp section
			{
				//face j0
				if (jel_TS_planes[2] == 1)
					TS_Hx_j0();
			}
			#pragma omp section
			{
				//face j1
				if (jel_TS_planes[3] == 1)
					TS_Hx_j1();
			}
			#pragma omp section
			{
				//face k0
				if (jel_TS_planes[4] == 1)
					TS_Hx_k0();
			}
			#pragma omp section
			{
				//face k1
				if (jel_TS_planes[5] == 1)
					TS_Hx_k1();
			}
		}
	}
}

///////////////////////////////////////////////////////////////////////////////
//Update the Hx field components which will be send to other processes
///////////////////////////////////////////////////////////////////////////////
void CFDTD_3D::Update_Hx_send()
{
	long i, j, k;

	//update Hx_send_j
	if ( nly_Hx > 0 )
	{
		for (i = 0; i < nlx_Hx; i++)
		{
			for (k = 0; k < nlz_Hx; k++)
			{
				Hx_send_j[i][k] = Hx[i][nly_HxMIN1][k];
			}
		}
	}
	//update Hx_send_k
	if ( nlz_Hx > 0 )
	{
		for (i = 0; i < nlx_Hx; i++)
		{
			for (j = 0; j < nly_Hx; j++)
			{
				Hx_send_k[i][j] = Hx[i][j][nlz_HxMIN1];
			}
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate the Hy field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D::Calc_Hy(long  nlx, long  nly, long  nlz)
{
	long  i, j, k;
	double By_r;

	long nlxMIN1 = nlx - 1;
	long nlzMIN1 = nlz - 1;

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

				By[i][j][k] = K_By_a[k]*By[i][j][k] +
							  K_By_b[k]*( (Ez_recv_i[j][k] - Ez[i][j][k])*inv_dx -
										  (Ex_recv_k[i][j] - Ex[i][j][k])*inv_dz );
				
				Hy[i][j][k] = K_Hy_a[i]*Hy[i][j][k] + 
							  K_Hy_b[i]*( K_Hy_c[j]*By[i][j][k] + 
										  K_Hy_d[j]*By_r )/mu_r[Ind[i][j][k]];
			}
		}

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

					By[i][j][k] = K_By_a[k]*By[i][j][k] +
								  K_By_b[k]*( (Ez_recv_i[j][k] - Ez[i][j][k])*inv_dx -
											  (Ex[i][j][k+1] - Ex[i][j][k])*inv_dz );
					
					Hy[i][j][k] = K_Hy_a[i]*Hy[i][j][k] + 
								  K_Hy_b[i]*( K_Hy_c[j]*By[i][j][k] + 
											  K_Hy_d[j]*By_r )/mu_r[Ind[i][j][k]];
				}
			}
		}

		k = nlzMIN1;
		if ( nlz > 0 )
		{
			#pragma omp for schedule(dynamic,nr_threads) nowait
			for (i = 0; i < nlxMIN1; i++)
			{	
				for (j = 0; j < nly; j++)
				{	
					By_r = By[i][j][k];

					By[i][j][k] = K_By_a[k]*By[i][j][k] +
								  K_By_b[k]*( (Ez[i+1][j][k] - Ez[i][j][k])*inv_dx -
											  (Ex_recv_k[i][j] - Ex[i][j][k])*inv_dz );
					
					Hy[i][j][k] = K_Hy_a[i]*Hy[i][j][k] + 
								  K_Hy_b[i]*( K_Hy_c[j]*By[i][j][k] + 
											  K_Hy_d[j]*By_r )/mu_r[Ind[i][j][k]];
				}
			}
		}

		#pragma omp for schedule(dynamic,nr_threads)
		for (i = 0; i < nlxMIN1; i++)
		{	
			for (j = 0; j < nly; j++)
			{	
				for (k = 0; k < nlzMIN1; k++)
				{
					By_r = By[i][j][k];

					By[i][j][k] = K_By_a[k]*By[i][j][k] +
								  K_By_b[k]*( (Ez[i+1][j][k] - Ez[i][j][k])*inv_dx -
											  (Ex[i][j][k+1] - Ex[i][j][k])*inv_dz );
					
					Hy[i][j][k] = K_Hy_a[i]*Hy[i][j][k] + 
								  K_Hy_b[i]*( K_Hy_c[j]*By[i][j][k] + 
											  K_Hy_d[j]*By_r )/mu_r[Ind[i][j][k]];
				}
			}
		}
	}

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

	////////////////////////////////////////////////////////////////////
	//Total field scattered field formulation
	////////////////////////////////////////////////////////////////////
	if (jel_plane_wave == 1 && jel_TS == 1) 
	{
		//the parallelization can be done this way because in case of 
		//magnetic field there aro no overlapping elements on the interface
		#pragma omp parallel sections default(shared)
		{
			#pragma omp section
			{
				//i0 face
				if(jel_TS_planes[0] == 1)
					TS_Hy_i0();
			}
			#pragma omp section
			{
				//i1 face
				if(jel_TS_planes[1] == 1)
					TS_Hy_i1();
			}
			#pragma omp section
			{
				//k0 face
				if(jel_TS_planes[4] == 1)
 					TS_Hy_k0();
			}
			#pragma omp section
			{
				//k1 face
				if(jel_TS_planes[5] == 1)
					TS_Hy_k1();
			}
		}
	}
}

///////////////////////////////////////////////////////////////////////////////
//Update the Hy field components which will be send to other processes
///////////////////////////////////////////////////////////////////////////////
void CFDTD_3D::Update_Hy_send()
{
	long i, j, k;

	//update Hy_send_i
	if ( nlx_Hy > 0 )
	{
		for (j = 0; j < nly_Hy; j++)
		{
			for (k = 0; k < nlz_Hy; k++)
			{
				Hy_send_i[j][k] = Hy[nlx_HyMIN1][j][k];
			}
		}
	}
	//update Hy_send_k
	if ( nlz_Hy > 0 )
	{
		for (i = 0; i < nlx_Hy; i++)
		{
			for (j = 0; j < nly_Hy; j++)
			{
				Hy_send_k[i][j] = Hy[i][j][nlz_HyMIN1];
			}
		}
	}	
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate the Hz field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D::Calc_Hz(long  nlx, long  nly, long  nlz)
{
	long  i, j, k;
	double Bz_r;
	
	long nlxMIN1 = nlx - 1;
	long nlyMIN1 = nly - 1;

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

				Bz[i][j][k] = K_Bz_a[i]*Bz[i][j][k] +  
							  K_Bz_b[i]*( (Ex_recv_j[i][k] - Ex[i][j][k])*inv_dy - 
										  (Ey_recv_i[j][k] - Ey[i][j][k])*inv_dx );

				Hz[i][j][k] = K_Hz_a[j]*Hz[i][j][k] + 
							  K_Hz_b[j]*( K_Hz_c[k]*Bz[i][j][k] + 
										  K_Hz_d[k]*Bz_r )/mu_r[Ind[i][j][k]];
			}
		}

		j = nlyMIN1;
		if ( nly > 0 )
		{
			#pragma omp for schedule(dynamic,nr_threads) nowait
			for (i = 0; i < nlxMIN1; i++)
			{	
				for (k = 0; k < nlz; k++)
				{
					Bz_r = Bz[i][j][k];

					Bz[i][j][k] = K_Bz_a[i]*Bz[i][j][k] +  
								  K_Bz_b[i]*( (Ex_recv_j[i][k] - Ex[i][j][k])*inv_dy - 
											(  Ey[i+1][j][k] - Ey[i][j][k])*inv_dx );

					Hz[i][j][k] = K_Hz_a[j]*Hz[i][j][k] + 
								  K_Hz_b[j]*( K_Hz_c[k]*Bz[i][j][k] + 
											  K_Hz_d[k]*Bz_r )/mu_r[Ind[i][j][k]];
				}
			}
		}

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

					Bz[i][j][k] = K_Bz_a[i]*Bz[i][j][k] +  
								  K_Bz_b[i]*( (Ex[i][j+1][k] - Ex[i][j][k])*inv_dy - 
											  (Ey_recv_i[j][k] - Ey[i][j][k])*inv_dx );

					Hz[i][j][k] = K_Hz_a[j]*Hz[i][j][k] + 
								  K_Hz_b[j]*( K_Hz_c[k]*Bz[i][j][k] + 
											  K_Hz_d[k]*Bz_r )/mu_r[Ind[i][j][k]];
				}
			}
		}

		#pragma omp for schedule(dynamic,nr_threads)
		for (i = 0; i < nlxMIN1; i++)
		{	
			for (j = 0; j < nlyMIN1; j++)
			{	
				for (k = 0; k < nlz; k++)
				{
					Bz_r = Bz[i][j][k];

					Bz[i][j][k] = K_Bz_a[i]*Bz[i][j][k] +  
								  K_Bz_b[i]*( (Ex[i][j+1][k] - Ex[i][j][k])*inv_dy - 
											  (Ey[i+1][j][k] - Ey[i][j][k])*inv_dx );

					Hz[i][j][k] = K_Hz_a[j]*Hz[i][j][k] + 
								  K_Hz_b[j]*( K_Hz_c[k]*Bz[i][j][k] + 
											  K_Hz_d[k]*Bz_r )/mu_r[Ind[i][j][k]];
				}
			}
		}
	}

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

	////////////////////////////////////////////////////////////////////
	//Total field scattered field formulation
	////////////////////////////////////////////////////////////////////
	if (jel_plane_wave == 1 && jel_TS == 1) 
	{
		//the parallelization can be done this way because in case of 
		//magnetic field there aro no overlapping elements on the interface
		#pragma omp parallel sections default(shared)
		{
			#pragma omp section
			{
				//face i0
				if(jel_TS_planes[0] == 1)
					TS_Hz_i0();
			}
			#pragma omp section
			{
				//face i1
				if(jel_TS_planes[1] == 1)
					TS_Hz_i1();
			}
			#pragma omp section
			{
				//face j0
				if(jel_TS_planes[2] == 1)
					TS_Hz_j0();
			}
			#pragma omp section
			{
				//face j1
				if(jel_TS_planes[3] == 1)
					TS_Hz_j1();
			}
		}
	}
}

///////////////////////////////////////////////////////////////////////////////
//Update the Hz field components which will be send to other processes
///////////////////////////////////////////////////////////////////////////////
void CFDTD_3D::Update_Hz_send()
{
	long i, j, k;

	long nlx_HzMIN1 = nlx_Hz - 1;
	long nly_HzMIN1 = nly_Hz - 1;

	//update Hz_send_i
	if ( nlx_Hz > 0 )
	{
		for (j = 0; j < nly_Hz; j++)
		{
			for (k = 0; k < nlz_Hz; k++)
			{
				Hz_send_i[j][k] = Hz[nlx_HzMIN1][j][k];
			}
		}
	}
	//update Hz_send_j
	if ( nly_Hz > 0 )
	{
		for (i = 0; i < nlx_Hz; i++)
		{
			for (k = 0; k < nlz_Hz; k++)
			{
				Hz_send_j[i][k] = Hz[i][nly_HzMIN1][k];
			}
		}
	}
}

void CFDTD_3D::Get_SendRecv(double **&hz_recv_j, double **&hy_recv_k, 
							double **&ex_send_j, double **&ex_send_k,
							double **&hx_recv_k, double **&hz_recv_i, 
							double **&ey_send_i, double **&ey_send_k,
							double **&hy_recv_i, double **&hx_recv_j,
							double **&ez_send_i, double **&ez_send_j,
							double **&ey_recv_k, double **&ez_recv_j,
							double **&hx_send_j, double **&hx_send_k,
							double **&ez_recv_i, double **&ex_recv_k, 
							double **&hy_send_i, double **&hy_send_k,
							double **&ex_recv_j, double **&ey_recv_i, 
							double **&hz_send_i, double **&hz_send_j)
{
	//Ex
	hz_recv_j = Hz_recv_j;   hy_recv_k = Hy_recv_k; 
	ex_send_j = Ex_send_j;   ex_send_k = Ex_send_k;
	//Ey
	hx_recv_k = Hx_recv_k;   hz_recv_i = Hz_recv_i; 
	ey_send_i = Ey_send_i;   ey_send_k = Ey_send_k;
	//Ez
	hy_recv_i = Hy_recv_i;   hx_recv_j = Hx_recv_j;
	ez_send_i = Ez_send_i;   ez_send_j = Ez_s