📄 fdtd_3d_xyzpml_threads_decomp.cpp
字号:
Save_Hz_xz = (Data_Save_Slice *)calloc(Save_nr_DIV_xORz, sizeof(Data_Save_Slice) );
}
if (jel_plane_wave == 1)
{
n_1D = (int) (l_x*abs(sin_teta)*abs(cos_phi) +
l_y*abs(sin_teta)*abs(sin_phi) + l_z*abs(cos_teta) ) +
10 + 2*n_PML_1D;
n_1D_MIN_1 = n_1D - 1;
d_1D = dx*sin_teta*sin_teta*cos_phi*cos_phi + dy*sin_teta*sin_teta*sin_phi*sin_phi +
dz*cos_teta*cos_teta;
er_TS = Init_TS_Param(&ll_1D_E, &ll_1D_H,
&Hz_i0, &Hy_i0, &Hz_i1, &Hy_i1, &Hz_j0, &Hx_j0, &Hz_j1,
&Hx_j1, &Hy_k0, &Hx_k0, &Hy_k1, &Hx_k1,
&face_Hz_i0, &face_Hy_i0, &face_Hz_i1, &face_Hy_i1,
&face_Hz_j0, &face_Hx_j0, &face_Hz_j1, &face_Hx_j1,
&face_Hy_k0, &face_Hx_k0, &face_Hy_k1, &face_Hx_k1,
&Ey_i0, &Ez_i0, &Ey_i1, &Ez_i1, &Ex_j0, &Ez_j0, &Ex_j1,
&Ez_j1, &Ex_k0, &Ey_k0, &Ex_k1, &Ey_k1,
&face_Ey_i0, &face_Ez_i0, &face_Ey_i1, &face_Ez_i1,
&face_Ex_j0, &face_Ez_j0, &face_Ex_j1, &face_Ez_j1,
&face_Ex_k0, &face_Ey_k0, &face_Ex_k1, &face_Ey_k1,
n_1D, l_x, l_y, l_z, n_PML_1D, teta, phi, gamma,
dx, dy, dz, d_1D, dt,
ct_Ex_1, ct_Ex_2, ct_Ey_1, ct_Ey_2, ct_Ez_1, ct_Ez_2,
ct_Hx_1, ct_Hx_2, ct_Hy_1, ct_Hy_2, ct_Hz_1, ct_Hz_2);
}
//Join the memory allocation threads
pthread_join(thread_FGE_x,NULL);
pthread_join(thread_FGE_y,NULL);
pthread_join(thread_FGE_z,NULL);
pthread_join(thread_BH_x,NULL);
pthread_join(thread_BH_y,NULL);
pthread_join(thread_BH_z,NULL);
//clear up the memory and exit if the memory allocation is unsuccesfull
if ( !(Allocate_FGE_x->er && Allocate_FGE_y->er && Allocate_FGE_z->er &&
Allocate_BH_x->er && Allocate_BH_y->er && Allocate_BH_z->er && er_PML_Ex
&& er_PML_Ey && er_PML_Ez && er_PML_Hx && er_PML_Hy && er_PML_Hz)
&& !K_a && !K_b && !Path && !Path_Data && !path_name_E_1D &&
!path_name_Ex && !path_name_Ey && !path_name_Ez && !path_name_Hx &&
!path_name_Hy && !path_name_Hz && !Ex_Foll && !Ey_Foll && !Ez_Foll &&
!Hx_Foll && !Hy_Foll && !Hz_Foll && !field_threads && !Data_Ex_PML &&
!Data_Ex && !Data_Ey_PML && !Data_Ey && !Data_Ez_PML && !Data_Ez &&
!K_a_thread && !K_b_thread && !mu_thread && (!data_save_threads_Ex_xy &&
!data_save_threads_Ey_xy && !data_save_threads_Ez_xy &&
!data_save_threads_Hx_xy && !data_save_threads_Hy_xy &&
!data_save_threads_Hz_xy && !data_save_threads_Ex_yz &&
!data_save_threads_Ey_yz && !data_save_threads_Ez_yz &&
!data_save_threads_Hx_yz && !data_save_threads_Hy_yz &&
!data_save_threads_Hz_yz && !data_save_threads_Ex_xz &&
!data_save_threads_Ey_xz && !data_save_threads_Ez_xz &&
!data_save_threads_Hx_xz && !data_save_threads_Hy_xz &&
!data_save_threads_Hz_xz && !Save_Ex_xy && !Save_Ey_xy &&
!Save_Ez_xy && !Save_Hx_xy && !Save_Hy_xy && !Save_Hz_xy &&
!Save_Ex_yz && !Save_Ey_yz && !Save_Ez_yz && !Save_Hx_yz &&
!Save_Hy_yz && !Save_Hz_yz && !Save_Ex_xz && !Save_Ey_xz &&
!Save_Ez_xz && !Save_Hx_xz && !Save_Hy_xz && !Save_Hz_xz ||
!save_field) && jel_alloc_W || er_TS)
{
Free_FxGxEx(&Fx_1, &Gx_1, &Fx_2, &Gx_2, &Fx_3, &Gx_3, &Fx_4, &Fx_5,
&Fx_6, &Fx_7, &Gx_7, &Fx_8, &Gx_8, &Fx_9, &Gx_9, &Fx_10,
&Gx_10, &Fx_11, &Fx_12, &Gx_12, &Fx_13, &Fx_15, &Fx_16, &Gx_16,
&Fx_17, &Fx_18, &Gx_18, &Fx_19, &Gx_19, &Fx_20, &Gx_20, &Fx_21,
&Gx_21, &Fx_22, &Fx_23, &Fx_24, &Fx_25, &Gx_25, &Fx_26, &Gx_26,
&Fx_27, &Gx_27, &Ex, nPML_x_1, nPML_x_2, nx);
Free_FyGyEy(&Fy_1, &Gy_1, &Fy_2, &Gy_2, &Fy_3, &Gy_3, &Fy_4, &Gy_4,
&Fy_5, &Fy_6, &Gy_6, &Fy_7, &Gy_7, &Fy_8, &Gy_8, &Fy_9,
&Gy_9, &Fy_10, &Fy_11, &Fy_12, &Fy_13, &Fy_15, &Fy_16, &Fy_17,
&Fy_18, &Fy_19, &Gy_19, &Fy_20, &Gy_20, &Fy_21, &Gy_21, &Fy_22,
&Gy_22, &Fy_23, &Fy_24, &Gy_24, &Fy_25, &Gy_25, &Fy_26, &Gy_26,
&Fy_27, &Gy_27, &Ey, nPML_x_1, nPML_x_2, nx);
Free_FzGzEz(&Fz_1, &Gz_1, &Fz_2, &Fz_3, &Gz_3, &Fz_4, &Gz_4, &Fz_5,
&Fz_6, &Gz_6, &Fz_7, &Gz_7, &Fz_8, &Fz_9, &Gz_9, &Fz_10,
&Gz_10, &Fz_11, &Fz_12, &Gz_12, &Fz_13, &Fz_15, &Fz_16, &Gz_16,
&Fz_17, &Fz_18, &Gz_18, &Fz_19, &Gz_19, &Fz_20, &Fz_21, &Gz_21,
&Fz_22, &Gz_22, &Fz_23, &Fz_24, &Gz_24, &Fz_25, &Gz_25, &Fz_26,
&Fz_27, &Gz_27, &Ez, nPML_x_1, nPML_x_2, nx);
Free_BxHx(&Bx_1, &Bx_2, &Bx_3, &Bx_4, &Bx_6, &Bx_7, &Bx_8, &Bx_9,
&Bx_10, &Bx_12, &Bx_13, &Bx_15,&Bx_16, &Bx_18, &Bx_19, &Bx_20,
&Bx_21, &Bx_22, &Bx_24, &Bx_25,&Bx_26, &Bx_27, &Hx, nPML_x_1,
nPML_x_2, nx);
Free_ByHy(&By_1, &By_2, &By_3, &By_4, &By_6, &By_7, &By_8, &By_9,
&By_10, &By_11, &By_12, &By_16, &By_17, &By_18, &By_19, &By_20,
&By_21, &By_22, &By_24, &By_25, &By_26, &By_27, &Hy, nPML_x_1,
nPML_x_2, nx);
Free_BzHz(&Bz_1, &Bz_2, &Bz_3, &Bz_4, &Bz_5, &Bz_6, &Bz_7, &Bz_8,
&Bz_9, &Bz_10, &Bz_12, &Bz_16, &Bz_18, &Bz_19, &Bz_20, &Bz_21,
&Bz_22, &Bz_23, &Bz_24, &Bz_25, &Bz_26, &Bz_27, &Hz, nPML_x_1,
nPML_x_2, nx);
Free_PML_Ex(&K_Gx_a_1, &K_Gx_b_1, &K_Ex_a_1, &K_Ex_b_1, &K_Ex_c_1,
&K_Ex_d_1, &K_Gx_a_2, &K_Gx_b_2, &K_Ex_a_2, &K_Ex_b_2,
&K_Ex_c_2, &K_Ex_d_2, nPML_x_1, nPML_x_2);
Free_PML_Ey(&K_Gy_a_1, &K_Gy_b_1, &K_Ey_a_1, &K_Ey_b_1, &K_Ey_c_1,
&K_Ey_d_1, &K_Gy_a_2, &K_Gy_b_2, &K_Ey_a_2, &K_Ey_b_2,
&K_Ey_c_2, &K_Ey_d_2, nPML_x_1, nPML_x_2);
Free_PML_Ez(&K_Gz_a_1, &K_Gz_b_1, &K_Ez_a_1, &K_Ez_b_1, &K_Ez_c_1,
&K_Ez_d_1, &K_Gz_a_2, &K_Gz_b_2, &K_Ez_a_2, &K_Ez_b_2,
&K_Ez_c_2, &K_Ez_d_2, nPML_x_1, nPML_x_2);
Free_PML_Hx(&K_Bx_a_1, &K_Bx_b_1, &K_Hx_a_1, &K_Hx_b_1, &K_Hx_c_1,
&K_Hx_d_1, &K_Bx_a_2, &K_Bx_b_2, &K_Hx_a_2, &K_Hx_b_2,
&K_Hx_c_2, &K_Hx_d_2, nPML_x_1, nPML_x_2);
Free_PML_Hy(&K_By_a_1, &K_By_b_1, &K_Hy_a_1, &K_Hy_b_1, &K_Hy_c_1,
&K_Hy_d_1, &K_By_a_2, &K_By_b_2, &K_Hy_a_2, &K_Hy_b_2,
&K_Hy_c_2, &K_Hy_d_2, nPML_x_1, nPML_x_2);
Free_PML_Hz(&K_Bz_a_1, &K_Bz_b_1, &K_Hz_a_1, &K_Hz_b_1, &K_Hz_c_1,
&K_Hz_d_1, &K_Bz_a_2, &K_Bz_b_2, &K_Hz_a_2, &K_Hz_b_2,
&K_Hz_c_2, &K_Hz_d_2, nPML_x_1, nPML_x_2);
Free_K_a_K_b(&K_a, &K_b);
Free_path_name(&Path, &Path_Data, &path_name_E_1D, &path_name_Ex,
&path_name_Ey, &path_name_Ez, &path_name_Hx, &path_name_Hy,
&path_name_Hz);
Free_Foll(&Ex_Foll, &Ey_Foll, &Ez_Foll, &Hx_Foll, &Hy_Foll, &Hz_Foll);
Free_Threads_Param(&field_threads, &Data_Ex_PML, &Data_Ex, &Data_Ey_PML,
&Data_Ey, &Data_Ez_PML, &Data_Ez, &Data_Hx_PML, &Data_Hx,
&Data_Hy_PML, &Data_Hy, &Data_Hz_PML, &Data_Hz, &K_a_thread,
&K_b_thread, &mu_thread);
Free_Data_SaveThreads_Param(&data_save_threads_Ex_xy, &data_save_threads_Ey_xy,
&data_save_threads_Ez_xy, &data_save_threads_Hx_xy,
&data_save_threads_Hy_xy, &data_save_threads_Hz_xy,
&data_save_threads_Ex_yz, &data_save_threads_Ey_yz,
&data_save_threads_Ez_yz, &data_save_threads_Hx_yz,
&data_save_threads_Hy_yz, &data_save_threads_Hz_yz,
&data_save_threads_Ex_xz, &data_save_threads_Ey_xz,
&data_save_threads_Ez_xz, &data_save_threads_Hx_xz,
&data_save_threads_Hy_xz, &data_save_threads_Hz_xz,
&Save_Ex_xy, &Save_Ey_xy, &Save_Ez_xy,
&Save_Hx_xy, &Save_Hy_xy, &Save_Hz_xy,
&Save_Ex_yz, &Save_Ey_yz, &Save_Ez_yz,
&Save_Hx_yz, &Save_Hy_yz, &Save_Hz_yz,
&Save_Ex_xz, &Save_Ey_xz, &Save_Ez_xz,
&Save_Hx_xz, &Save_Hy_xz, &Save_Hz_xz);
#ifdef eval_W
Free_W_Param(&W, &W_aver, &W_aver_xz, &W_aver_xzt, &eval_W_thread, &Data_W,
&local_W_aver, &path_name_W,nx_W);
#endif
Free_TS_Param(&ll_1D_E, &ll_1D_H,
&Hz_i0, &Hy_i0, &Hz_i1, &Hy_i1, &Hz_j0, &Hx_j0, &Hz_j1,
&Hx_j1, &Hy_k0, &Hx_k0, &Hy_k1, &Hx_k1,
&face_Hz_i0, &face_Hy_i0, &face_Hz_i1, &face_Hy_i1,
&face_Hz_j0, &face_Hx_j0, &face_Hz_j1, &face_Hx_j1,
&face_Hy_k0, &face_Hx_k0, &face_Hy_k1, &face_Hx_k1,
&Ey_i0, &Ez_i0, &Ey_i1, &Ez_i1, &Ex_j0, &Ez_j0, &Ex_j1,
&Ez_j1, &Ex_k0, &Ey_k0, &Ex_k1, &Ey_k1,
&face_Ey_i0, &face_Ez_i0, &face_Ey_i1, &face_Ez_i1,
&face_Ex_j0, &face_Ez_j0, &face_Ex_j1, &face_Ez_j1,
&face_Ex_k0, &face_Ey_k0, &face_Ex_k1, &face_Ey_k1);
cout << "Memory allocation problem FDTD parameters." << endl;
return -1;
}
///////////////////////////////////////////////////////////////////////////////////////
//End main memory allocatin region
///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////
//Initialize K_a and K_b - the contribution of materials
///////////////////////////////////////////////////////////////////////////////////////
double eps_r, mu_r, sigma;
for (i =0; i<n_Mat; i++)
{
eps_r = Mat[i][0];
sigma = Mat[i][2];
K_a[i] = ( 2.0*eps_0*eps_r - sigma*dt )/( 2.0*eps_0*eps_r + sigma*dt );
K_b[i] = 2.0*dt/( 2.0*eps_0*eps_r + sigma*dt );
}
long j;
//reproduce the K_a, K_b and mu_r to avoid concurent memory reads by threads
for (i = 0; i<nr_Threads; i++)
{
for (j = 0; j<n_Mat; j++)
{
K_a_thread[i][j] = K_a[j];
K_b_thread[i][j] = K_b[j];
mu_thread[i][j] = Mat[j][1]; //mu_r
}
}
///////////////////////////////////////////////////////////////////////////////////////
//Initialize the Fourier transform
///////////////////////////////////////////////////////////////////////////////////////
if(Inp_D->fourier_transf_vol == 1)
{
double delta_omega = TwoORpi*(frec_2 - frec_1)/(n_frec-1);
fourier_omega[0] = TwoORpi*frec_1;
for (i = 1; i< n_frec; i++)
{
fourier_omega[i] = fourier_omega[i-1] + delta_omega;
}
}
///////////////////////////////////////////////////////////////////////////////////////
//Initialize the PML parameters
///////////////////////////////////////////////////////////////////////////////////////
//in x direction
double PML_eps_r_x_1 = Inp_D->PML_eps_r_x_1;
double PML_mu_r_x_1 = Inp_D->PML_mu_r_x_1;
double PML_eps_r_x_2 = Inp_D->PML_eps_r_x_2;
double PML_mu_r_x_2 = Inp_D->PML_mu_r_x_2;
//in y direction
double PML_eps_r_y_1 = Inp_D->PML_eps_r_y_1;
double PML_mu_r_y_1 = Inp_D->PML_mu_r_y_1;
double PML_eps_r_y_2 = Inp_D->PML_eps_r_y_2;
double PML_mu_r_y_2 = Inp_D->PML_mu_r_y_2;
//in z direction
double PML_eps_r_z_1 = Inp_D->PML_eps_r_z_1;
double PML_mu_r_z_1 = Inp_D->PML_mu_r_z_1;
double PML_eps_r_z_2 = Inp_D->PML_eps_r_z_2;
double PML_mu_r_z_2 = Inp_D->PML_mu_r_z_2;
//in x direction
//1
Init_PML_Par_x(K_Ey_a_1, K_Ey_b_1, K_Gz_a_1, K_Gz_b_1, K_Hx_c_1, K_Hx_d_1,
K_Ex_c_1, K_Ex_d_1, K_Hy_a_1, K_Hy_b_1, K_Bz_a_1, K_Bz_b_1,
PML_eps_r_x_1, PML_mu_r_x_1, nPML_x_1, dx, dt);
//2
Init_PML_Par_x(K_Ey_a_2, K_Ey_b_2, K_Gz_a_2, K_Gz_b_2, K_Hx_c_2, K_Hx_d_2,
K_Ex_c_2, K_Ex_d_2, K_Hy_a_2, K_Hy_b_2, K_Bz_a_2, K_Bz_b_2,
PML_eps_r_x_2, PML_mu_r_x_2, nPML_x_2, dx, dt);
//in y direction
//1
Init_PML_Par_y(K_Gx_a_1, K_Gx_b_1, K_Ez_a_1, K_Ez_b_1, K_Hy_c_1, K_Hy_d_1,
K_Ey_c_1, K_Ey_d_1, K_Bx_a_1, K_Bx_b_1, K_Hz_a_1, K_Hz_b_1,
PML_eps_r_y_1, PML_mu_r_y_1, nPML_y_1, dy, dt);
//2
Init_PML_Par_y(K_Gx_a_2, K_Gx_b_2, K_Ez_a_2, K_Ez_b_2, K_Hy_c_2, K_Hy_d_2,
K_Ey_c_2, K_Ey_d_2, K_Bx_a_2, K_Bx_b_2, K_Hz_a_2, K_Hz_b_2,
PML_eps_r_y_2, PML_mu_r_y_2, nPML_y_2, dy, dt);
//in z direction
//1
Init_PML_Par_z(K_Ex_a_1, K_Ex_b_1, K_Gy_a_1, K_Gy_b_1, K_Hz_c_1, K_Hz_d_1,
K_Ez_c_1, K_Ez_d_1, K_Hx_a_1, K_Hx_b_1, K_By_a_1, K_By_b_1,
PML_eps_r_z_1, PML_mu_r_z_1, nPML_z_1, dz, dt);
//2
Init_PML_Par_z(K_Ex_a_2, K_Ex_b_2, K_Gy_a_2, K_Gy_b_2, K_Hz_c_2, K_Hz_d_2,
K_Ez_c_2, K_Ez_d_2, K_Hx_a_2, K_Hx_b_2, K_By_a_2, K_By_b_2,
PML_eps_r_z_2, PML_mu_r_z_2, nPML_z_2, dz, dt);
///////////////////////////////////////////////////////////////////////////////////////
//Initializations in case of plane wave excitation
///////////////////////////////////////////////////////////////////////////////////////
double phase_vel, phase_vel_ref, corr_term;
if (jel_plane_wave == 1)
{
eps_r = Mat[0][0];
mu_r = Mat[0][1];
sigma = Mat[0][2];
if (source_type > 1) //Sin or SinGauss
{
//correction term to shyncronize the 1D and 2D wave propagations
phase_velocity(dx, dy, dz, dt, Mat[0][0], Mat[0][1], omega, teta, phi, &phase_vel);
if (abs(sin_teta*cos_phi) >= abs(sin_teta*sin_phi) &&
abs(sin_teta*cos_phi) >= abs(cos_teta) )
{
phase_velocity(dx, dy, dz, dt, Mat[0][0], Mat[0][1], omega, pi/2.0, 0.0, &phase_vel_ref);
}
if (abs(sin_teta*sin_phi) >= abs(sin_teta*cos_phi) &&
abs(sin_teta*sin_phi) >= abs(cos_teta) )
{
phase_velocity(dx, dy, dz, dt, Mat[0][0], Mat[0][1], omega, pi/2.0, pi/2.0, &phase_vel_ref);
}
if (abs(cos_teta) >= abs(sin_teta*cos_phi) &&
abs(cos_teta) >= abs(sin_teta*sin_phi) )
{
phase_velocity(dx, dy, dz, dt, Mat[0][0], Mat[0][1], omega, 0, pi/2.0, &phase_vel_ref);
}
corr_term = phase_vel_ref/phase_vel;
}
switch (source_type)
{
case 1: //Gauss
fdtd_1D.Init(n_1D, n_PML_1D, dt, d_1D, eps_r, mu_r, sigma);
fdtd_1D.Init_Gauss_Source(X0, t0, tw);
break;
case 2://Sin
fdtd_1D.Init(n_1D, n_PML_1D, dt, d_1D, corr_term*eps_r, corr_term*mu_r, sigma);
fdtd_1D.Init_Sin_Source(X0, omega, phase);
break;
case 3://SinGauss
fdtd_1D.Init(n_1D, n_PML_1D, dt, d_1D, corr_term*eps_r, corr_term*mu_r, sigma);
fdtd_1D.Init_GaussSin_Source(X0, omega, phase, t0, tw);
break;
}
fdtd_1D.Init_PML_Param();
fdtd_1D.Get_Data(E_1D,H_1D);
}
/////////////////////////////////////////////////////////////////////////////
//Initializations for data save
/⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -