fdtd_3d_lorentz.cpp

三维FDTD

	//memory allocations to compute Hx
	///////////////////////////////////////////////////////////////////////////////////////
	ista_Hx = ista;	   iend_Hx = iend;    jsta_Hx = jsta;
	jend_Hx = jend;	   ksta_Hx = ksta;	  kend_Hx = kend;
	nlx_Hx = nlx;	   nly_Hx = nly;	  nlz_Hx = nlz;
	if (myrank_j == jprocs-1)
	{
		jend_Hx--;
		nly_Hx--;
	}
	if (myrank_k == kprocs-1)
	{
		kend_Hx--;
		nlz_Hx--;
	}

	Hx = Init_Matrix_3D<double>(nlx_Hx, nly_Hx, nlz_Hx);
	if (!Hx)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hx - ");
		return 1;
	}

	Bx = Init_Matrix_3D<double>(nlx_Hx, nly_Hx, nlz_Hx);
	if (!Bx)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Bx - ");
		return 1;
	}

	Ey_recv_k = Init_Matrix_2D<double>(nlx_Hx, nly_Hx);
	if (!Ey_recv_k)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ey_recv_k - ");
		return 1;
	}

	Ez_recv_j = Init_Matrix_2D<double>(nlx_Hx, nlz_Hx);
	if (!Ez_recv_j)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ez_recv_j - ");
		return 1;
	}

	Hx_send_j = Init_Matrix_2D<double>(nlx_Hx, nlz_Hx);
	if (!Hx_send_j)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hx_send_j - ");
		return 1;
	}

	Hx_send_k = Init_Matrix_2D<double>(nlx_Hx, nly_Hx);
	if (!Hx_send_k)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hx_send_k - ");
		return 1;
	}
	
	///////////////////////////////////////////////////////////////////////////////////////
	//memory allocations to compute Hy
	///////////////////////////////////////////////////////////////////////////////////////
	ista_Hy = ista;	   iend_Hy = iend;    jsta_Hy = jsta;
	jend_Hy = jend;	   ksta_Hy = ksta;	  kend_Hy = kend;
	nlx_Hy = nlx;	   nly_Hy = nly;	  nlz_Hy = nlz;
	if (myrank_i == iprocs-1)
	{
		iend_Hy--;
		nlx_Hy--;
	}
	if (myrank_k == kprocs-1)
	{
		kend_Hy--;
		nlz_Hy--;
	}

	Hy = Init_Matrix_3D<double>(nlx_Hy, nly_Hy, nlz_Hy);
	if (!Hy)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hy - ");
		return 1;
	}

	By = Init_Matrix_3D<double>(nlx_Hy, nly_Hy, nlz_Hy);
	if (!By)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - By - ");
		return 1;
	}

	Ez_recv_i = Init_Matrix_2D<double>(nly_Hy, nlz_Hy);
	if (!Ez_recv_i)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ez_recv_i - ");
		return 1;
	}

	Ex_recv_k = Init_Matrix_2D<double>(nlx_Hy, nly_Hy);
	if (!Ex_recv_k)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ex_recv_k - ");
		return 1;
	}

	Hy_send_i = Init_Matrix_2D<double>(nly_Hy, nlz_Hy);
	if (!Hy_send_i)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hy_send_i - ");
		return 1;
	}

	Hy_send_k = Init_Matrix_2D<double>(nlx_Hy, nly_Hy);
	if (!Hy_send_k)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hy_send_k - ");
		return 1;
	}

	///////////////////////////////////////////////////////////////////////////////////////
	//memory allocations to compute Hz
	///////////////////////////////////////////////////////////////////////////////////////
	ista_Hz = ista;	   iend_Hz = iend;    jsta_Hz = jsta;
	jend_Hz = jend;	   ksta_Hz = ksta;	  kend_Hz = kend;
	nlx_Hz = nlx;	   nly_Hz = nly;	  nlz_Hz = nlz;
	if (myrank_i == iprocs-1)
	{
		iend_Hz--;
		nlx_Hz--;
	}
	if (myrank_j == jprocs-1)
	{
		jend_Hz--;
		nly_Hz--;
	}
	
	Hz = Init_Matrix_3D<double>(nlx_Hz, nly_Hz, nlz_Hz);
	if (!Hz)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hz - ");
		return 1;
	}

	Bz = Init_Matrix_3D<double>(nlx_Hz, nly_Hz, nlz_Hz);
	if (!Bz)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Bz - ");
		return 1;
	}

	Ex_recv_j = Init_Matrix_2D<double>(nlx_Hz, nlz_Hz);
	if (!Ex_recv_j)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ex_recv_j - ");
		return 1;
	}

	Ey_recv_i = Init_Matrix_2D<double>(nly_Hz, nlz_Hz);
	if (!Ey_recv_i)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ey_recv_i - ");
		return 1;
	}

	Hz_send_i = Init_Matrix_2D<double>(nly_Hz, nlz_Hz);
	if (!Hz_send_i)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hz_send_i - ");
		return 1;
	}

	Hz_send_j = Init_Matrix_2D<double>(nlx_Hz, nlz_Hz);
	if (!Hz_send_j)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hz_send_j - ");
		return 1;
	}
	
	///////////////////////////////////////////////////////////////////////////////////////
	//Coefficients containing the PML boundary parameters
	///////////////////////////////////////////////////////////////////////////////////////
	//PML-Ex field
	K_Gx_a = (double *)calloc(nly_Ex,sizeof(double)); 
	if (!K_Gx_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Gx_a - ");
		return 1;
	}

	K_Gx_b = (double *)calloc(nly_Ex,sizeof(double)); 
	if (!K_Gx_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Gx_b - ");
		return 1;
	}

	K_Dx_a = (double *)calloc(nlz_Ex,sizeof(double)); 
	if (!K_Dx_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Dx_a - ");
		return 1;
	}

	K_Dx_b = (double *)calloc(nlz_Ex,sizeof(double)); 
	if (!K_Dx_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Dx_b - ");
		return 1;
	}

	K_Dx_c = (double *)calloc(nlx_Ex,sizeof(double)); 
	if (!K_Dx_c)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Dx_c - ");
		return 1;
	}

	K_Dx_d = (double *)calloc(nlx_Ex,sizeof(double)); 
	if (!K_Dx_d)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Dx_d - ");
		return 1;
	}

	//PML-Ey field
	K_Gy_a = (double *)calloc(nlz_Ey,sizeof(double)); 
	if (!K_Gy_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Gy_a - ");
		return 1;
	}

	K_Gy_b = (double *)calloc(nlz_Ey,sizeof(double)); 
	if (!K_Gy_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Gy_b - ");
		return 1;
	}

	K_Dy_a = (double *)calloc(nlx_Ey,sizeof(double)); 
	if (!K_Dy_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Dy_a - ");
		return 1;
	}

	K_Dy_b = (double *)calloc(nlx_Ey,sizeof(double)); 
	if (!K_Dy_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Dy_b - ");
		return 1;
	}

	K_Dy_c = (double *)calloc(nly_Ey,sizeof(double)); 
	if (!K_Dy_c)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Dy_c - ");
		return 1;
	}

	K_Dy_d = (double *)calloc(nly_Ey,sizeof(double)); 
	if (!K_Dy_d)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Dy_d - ");
		return 1;
	}

	//PML-Ez field
	K_Gz_a = (double *)calloc(nlx_Ez,sizeof(double)); 
	if (!K_Gz_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Gz_a - ");
		return 1;
	}

	K_Gz_b = (double *)calloc(nlx_Ez,sizeof(double)); 
	if (!K_Gz_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Gz_b - ");
		return 1;
	}

	K_Dz_a = (double *)calloc(nly_Ez,sizeof(double)); 
	if (!K_Dz_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Dz_a - ");
		return 1;
	}

	K_Dz_b = (double *)calloc(nly_Ez,sizeof(double)); 
	if (!K_Dz_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Dz_b - ");
		return 1;
	}

	K_Dz_c = (double *)calloc(nlz_Ez,sizeof(double)); 
	if (!K_Dz_c)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Dz_c - ");
		return 1;
	}

	K_Dz_d = (double *)calloc(nlz_Ez,sizeof(double)); 
	if (!K_Dz_d)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Dz_d - ");
		return 1;
	}

	//PML-Hx field
	K_Bx_a = (double *)calloc(nly_Hx,sizeof(double)); 
	if (!K_Bx_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Bx_a - ");
		return 1;
	}

	K_Bx_b = (double *)calloc(nly_Hx,sizeof(double)); 
	if (!K_Bx_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Bx_b - ");
		return 1;
	}

	K_Hx_a = (double *)calloc(nlz_Hx,sizeof(double)); 
	if (!K_Hx_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hx_a - ");
		return 1;
	}

	K_Hx_b = (double *)calloc(nlz_Hx,sizeof(double)); 
	if (!K_Hx_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hx_b - ");
		return 1;
	}

	K_Hx_c = (double *)calloc(nlx_Hx,sizeof(double)); 
	if (!K_Hx_c)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hx_c - ");
		return 1;
	}

	K_Hx_d = (double *)calloc(nlx_Hx,sizeof(double)); 
	if (!K_Hx_d)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hx_d - ");
		return 1;
	}

	//PML-Hy field
	K_By_a = (double *)calloc(nlz_Hy,sizeof(double)); 
	if (!K_By_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_By_a - ");
		return 1;
	}

	K_By_b = (double *)calloc(nlz_Hy,sizeof(double)); 
	if (!K_By_b )
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_By_b  - ");
		return 1;
	}

	K_Hy_a = (double *)calloc(nlx_Hy,sizeof(double)); 
	if (!K_Hy_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hy_a - ");
		return 1;
	}

	K_Hy_b = (double *)calloc(nlx_Hy,sizeof(double)); 
	if (!K_Hy_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hy_b - ");
		return 1;
	}

	K_Hy_c = (double *)calloc(nly_Hy,sizeof(double)); 
	if (!K_Hy_c)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hy_c - ");
		return 1;
	}

	K_Hy_d = (double *)calloc(nly_Hy,sizeof(double)); 
	if (!K_Hy_d)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hy_d - ");
		return 1;
	}
	
	//PML-Hz field
	K_Bz_a = (double *)calloc(nlx_Hz,sizeof(double)); 
	if (!K_Bz_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Bz_a - ");
		return 1;
	}

	K_Bz_b = (double *)calloc(nlx_Hz,sizeof(double)); 
	if (!K_Bz_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Bz_b - ");
		return 1;
	}

	K_Hz_a = (double *)calloc(nly_Hz,sizeof(double)); 
	if (!K_Hz_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hz_a - ");
		return 1;
	}

	K_Hz_b = (double *)calloc(nly_Hz,sizeof(double)); 
	if (!K_Hz_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hz_b - ");
		return 1;
	}

	K_Hz_c = (double *)calloc(nlz_Hz,sizeof(double)); 
	if (!K_Hz_c)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hz_c - ");
		return 1;
	}

	K_Hz_d = (double *)calloc(nlz_Hz,sizeof(double)); 
	if (!K_Hz_d)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_Hz_d - ");
		return 1;
	}
	
	///////////////////////////////////////////////////////////////////////////////////////
	//Materials matrices
	///////////////////////////////////////////////////////////////////////////////////////
	long *n_lorentz = NULL;
	n_lorentz = (long *)calloc(n_Mat,sizeof(int));
	if (!n_lorentz)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - n_lorentz - ");
		return 1;
	}

	for (i =0; i<n_Mat; i++)
	{
		n_lorentz[i] = (long) Mat[i][0];
	}

	K_a = Init_Matrix_2D<double>(n_Mat,n_lorentz);
	if (!K_a)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_a - ");
		return 1;
	}

	K_b = Init_Matrix_2D<double>(n_Mat,n_lorentz);
	if (!K_b)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_b - ");
		return 1;
	}

	K_c = Init_Matrix_2D<double>(n_Mat,n_lorentz);
	if (!K_c)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - K_c - ");
		return 1;
	}
	
	if (n_lorentz)
	{
		free(n_lorentz);
		n_lorentz = NULL;
	}

	///////////////////////////////////////////////////////////////////////////////////////
	//Initialize K_a and K_b K_c - the contribution of Lorentz dispersive materials
	///////////////////////////////////////////////////////////////////////////////////////
	double am, delta_0, omega_0;
	for (i = 0; i<n_Mat; i++)
	{
		for (j = 0; j<Mat[i][0]; j++)
		{
			am = Mat[i][2+3*j];
			delta_0 = Mat[i][2+3*j+1];
			omega_0 = Mat[i][2+3*j+2];
			
			K_a[i][j] = 2*(2 - omega_0*omega_0*dt*dt)/