📄 fdtd_3d_xyzpml_threads_decomp.cpp
字号:
long Save_nr_W_DIV_xORz = Save_nr_W_DIV_x*Save_nr_W_DIV_z;
#endif
long l_x, l_y, l_z;
long n_TS_xa_MIN_1, n_TS_ya_MIN_1, n_TS_za_MIN_1;
double cos_teta, sin_teta, cos_phi, sin_phi, sin_gamma, cos_gamma;
long l_y_MIN_1__OR_l_z, l_y_OR__l_z_MIN_1, l_x_MIN_1__OR_l_z;
long l_x_OR__l_z_MIN_1, l_x_MIN_1__OR_l_y, l_x_OR__l_y_MIN_1;
//long n_TS_xaMIN1, n_TS_yaMIN1, n_TS_ybMIN1, n_TS_zbPL1, n_TS_ybMINn_TS_yaPL1;
//double coef_TS_Hx, coef_TS_Hy;
if (jel_plane_wave == 1)
{
//the size of the total field zone
l_x = n_TS_xb - n_TS_xa + 1;
l_y = n_TS_yb - n_TS_ya + 1;
l_z = n_TS_zb - n_TS_za + 1;
cos_teta = cos(teta);
sin_teta = sin(teta);
cos_phi = cos(phi);
sin_phi = sin(phi);
cos_gamma = cos(gamma);
sin_gamma = sin(gamma);
n_TS_xa_MIN_1 = n_TS_xa - 1;
n_TS_ya_MIN_1 = n_TS_ya - 1;
n_TS_za_MIN_1 = n_TS_za - 1;
l_y_MIN_1__OR_l_z = (l_y-1)*l_z;
l_y_OR__l_z_MIN_1 = l_y*(l_z-1);
l_x_MIN_1__OR_l_z = (l_x-1)*l_z;
l_x_OR__l_z_MIN_1 = l_x*(l_z-1);
l_x_MIN_1__OR_l_y = (l_x-1)*l_y;
l_x_OR__l_y_MIN_1 = l_x*(l_y-1);
//n_TS_xaMIN1 = n_TS_xa - 1;
//n_TS_yaMIN1 = n_TS_ya - 1;
//n_TS_ybMIN1 = n_TS_yb - 1;
//n_TS_zbPL1 = n_TS_zb + 1;
//n_TS_ybMINn_TS_yaPL1 = n_TS_yb - n_TS_ya + 1;
//coef_TS_Hx = dt/(mu_0*dy);
//coef_TS_Hy = dt/(mu_0*dx);
}
/////////////////////////////////////////////////////////////////////////////
//Initialization of the FDTD algorithm
/////////////////////////////////////////////////////////////////////////////
double ***Ex = NULL;
double ***Fx_1 = NULL, ***Fx_2 = NULL, ***Fx_3 = NULL;
double ***Fx_4 = NULL, ***Fx_5 = NULL, ***Fx_6 = NULL;
double ***Fx_7 = NULL, ***Fx_8 = NULL, ***Fx_9 = NULL;
double ***Fx_10 = NULL, ***Fx_11 = NULL, ***Fx_12 = NULL;
double ***Fx_13 = NULL, ***Fx_15 = NULL;
double ***Fx_16 = NULL, ***Fx_17 = NULL, ***Fx_18 = NULL;
double ***Fx_19 = NULL, ***Fx_20 = NULL, ***Fx_21 = NULL;
double ***Fx_22 = NULL, ***Fx_23 = NULL, ***Fx_24 = NULL;
double ***Fx_25 = NULL, ***Fx_26 = NULL, ***Fx_27 = NULL;
double ***Gx_1 = NULL, ***Gx_2 = NULL, ***Gx_3 = NULL;
double ***Gx_7 = NULL, ***Gx_8 = NULL, ***Gx_9 = NULL;
double ***Gx_10 = NULL, ***Gx_12 = NULL;
double ***Gx_16 = NULL, ***Gx_18 = NULL;
double ***Gx_19 = NULL, ***Gx_20 = NULL, ***Gx_21 = NULL;
double ***Gx_25 = NULL, ***Gx_26 = NULL, ***Gx_27 = NULL;
double ***Ey = NULL;
double ***Fy_1 = NULL, ***Fy_2 = NULL, ***Fy_3 = NULL;
double ***Fy_4 = NULL, ***Fy_5 = NULL, ***Fy_6 = NULL;
double ***Fy_7 = NULL, ***Fy_8 = NULL, ***Fy_9 = NULL;
double ***Fy_10 = NULL, ***Fy_11 = NULL, ***Fy_12 = NULL;
double ***Fy_13 = NULL, ***Fy_15 = NULL;
double ***Fy_16 = NULL, ***Fy_17 = NULL, ***Fy_18 = NULL;
double ***Fy_19 = NULL, ***Fy_20 = NULL, ***Fy_21 = NULL;
double ***Fy_22 = NULL, ***Fy_23 = NULL, ***Fy_24 = NULL;
double ***Fy_25 = NULL, ***Fy_26 = NULL, ***Fy_27 = NULL;
double ***Gy_1 = NULL, ***Gy_2 = NULL, ***Gy_3 = NULL;
double ***Gy_4 = NULL, ***Gy_6 = NULL;
double ***Gy_7 = NULL, ***Gy_8 = NULL, ***Gy_9 = NULL;
double ***Gy_19 = NULL, ***Gy_20 = NULL, ***Gy_21 = NULL;
double ***Gy_22 = NULL, ***Gy_24 = NULL;
double ***Gy_25 = NULL, ***Gy_26 = NULL, ***Gy_27 = NULL;
double ***Ez = NULL;
double ***Fz_1 = NULL, ***Fz_2 = NULL, ***Fz_3 = NULL;
double ***Fz_4 = NULL, ***Fz_5 = NULL, ***Fz_6 = NULL;
double ***Fz_7 = NULL, ***Fz_8 = NULL, ***Fz_9 = NULL;
double ***Fz_10 = NULL, ***Fz_11 = NULL, ***Fz_12 = NULL;
double ***Fz_13 = NULL, ***Fz_15 = NULL;
double ***Fz_16 = NULL, ***Fz_17 = NULL, ***Fz_18 = NULL;
double ***Fz_19 = NULL, ***Fz_20 = NULL, ***Fz_21 = NULL;
double ***Fz_22 = NULL, ***Fz_23 = NULL, ***Fz_24 = NULL;
double ***Fz_25 = NULL, ***Fz_26 = NULL, ***Fz_27 = NULL;
double ***Gz_1 = NULL, ***Gz_3 = NULL;
double ***Gz_4 = NULL, ***Gz_6 = NULL;
double ***Gz_7 = NULL, ***Gz_9 = NULL;
double ***Gz_10 = NULL, ***Gz_12 = NULL;
double ***Gz_16 = NULL, ***Gz_18 = NULL;
double ***Gz_19 = NULL, ***Gz_21 = NULL;
double ***Gz_22 = NULL, ***Gz_24 = NULL;
double ***Gz_25 = NULL, ***Gz_27 = NULL;
double ***Hx = NULL;
double ***Bx_1 = NULL, ***Bx_2 = NULL, ***Bx_3 = NULL;
double ***Bx_4 = NULL, ***Bx_6 = NULL;
double ***Bx_7 = NULL, ***Bx_8 = NULL, ***Bx_9 = NULL;
double ***Bx_10 = NULL, ***Bx_12 = NULL;
double ***Bx_13 = NULL, ***Bx_15 = NULL;
double ***Bx_16 = NULL, ***Bx_18 = NULL;
double ***Bx_19 = NULL, ***Bx_20 = NULL, ***Bx_21 = NULL;
double ***Bx_22 = NULL, ***Bx_24 = NULL;
double ***Bx_25 = NULL, ***Bx_26 = NULL, ***Bx_27 = NULL;
double ***Hy = NULL;
double ***By_1 = NULL, ***By_2 = NULL, ***By_3 = NULL;
double ***By_4 = NULL, ***By_6 = NULL;
double ***By_7 = NULL, ***By_8 = NULL, ***By_9 = NULL;
double ***By_10 = NULL, ***By_11 = NULL, ***By_12 = NULL;
double ***By_16 = NULL, ***By_17 = NULL, ***By_18 = NULL;
double ***By_19 = NULL, ***By_20 = NULL, ***By_21 = NULL;
double ***By_22 = NULL, ***By_24 = NULL;
double ***By_25 = NULL, ***By_26 = NULL, ***By_27 = NULL;
double ***Hz = NULL;
double ***Bz_1 = NULL, ***Bz_2 = NULL, ***Bz_3 = NULL;
double ***Bz_4 = NULL, ***Bz_5 = NULL, ***Bz_6 = NULL;
double ***Bz_7 = NULL, ***Bz_8 = NULL, ***Bz_9 = NULL;
double ***Bz_10 = NULL, ***Bz_12 = NULL;
double ***Bz_16 = NULL, ***Bz_18 = NULL;
double ***Bz_19 = NULL, ***Bz_20 = NULL, ***Bz_21 = NULL;
double ***Bz_22 = NULL, ***Bz_23 = NULL, ***Bz_24 = NULL;
double ***Bz_25 = NULL, ***Bz_26 = NULL, ***Bz_27 = NULL;
double *K_a = NULL, *K_b = NULL;
//Coefficients containing the PML boundary parameters
//Electric field
double *K_Gx_a_1 = NULL, *K_Gx_b_1 = NULL;
double *K_Ex_a_1 = NULL, *K_Ex_b_1 = NULL, *K_Ex_c_1 = NULL, *K_Ex_d_1 = NULL;
double *K_Gx_a_2 = NULL, *K_Gx_b_2 = NULL;
double *K_Ex_a_2 = NULL, *K_Ex_b_2 = NULL, *K_Ex_c_2 = NULL, *K_Ex_d_2 = NULL;
double *K_Gy_a_1 = NULL, *K_Gy_b_1 = NULL;
double *K_Ey_a_1 = NULL, *K_Ey_b_1 = NULL, *K_Ey_c_1 = NULL, *K_Ey_d_1 = NULL;
double *K_Gy_a_2 = NULL, *K_Gy_b_2 = NULL;
double *K_Ey_a_2 = NULL, *K_Ey_b_2 = NULL, *K_Ey_c_2 = NULL, *K_Ey_d_2 = NULL;
double *K_Gz_a_1 = NULL, *K_Gz_b_1 = NULL;
double *K_Ez_a_1 = NULL, *K_Ez_b_1 = NULL, *K_Ez_c_1 = NULL, *K_Ez_d_1 = NULL;
double *K_Gz_a_2 = NULL, *K_Gz_b_2 = NULL;
double *K_Ez_a_2 = NULL, *K_Ez_b_2 = NULL, *K_Ez_c_2 = NULL, *K_Ez_d_2 = NULL;
//Magnetic field
double *K_Bx_a_1 = NULL, *K_Bx_b_1 = NULL;
double *K_Hx_a_1 = NULL, *K_Hx_b_1 = NULL, *K_Hx_c_1 = NULL, *K_Hx_d_1 = NULL;
double *K_Bx_a_2 = NULL, *K_Bx_b_2 = NULL;
double *K_Hx_a_2 = NULL, *K_Hx_b_2 = NULL, *K_Hx_c_2 = NULL, *K_Hx_d_2 = NULL;
double *K_By_a_1 = NULL, *K_By_b_1 = NULL;
double *K_Hy_a_1 = NULL, *K_Hy_b_1 = NULL, *K_Hy_c_1 = NULL, *K_Hy_d_1 = NULL;
double *K_By_a_2 = NULL, *K_By_b_2 = NULL;
double *K_Hy_a_2 = NULL, *K_Hy_b_2 = NULL, *K_Hy_c_2 = NULL, *K_Hy_d_2 = NULL;
double *K_Bz_a_1 = NULL, *K_Bz_b_1 = NULL;
double *K_Hz_a_1 = NULL, *K_Hz_b_1 = NULL, *K_Hz_c_1 = NULL, *K_Hz_d_1 = NULL;
double *K_Bz_a_2 = NULL, *K_Bz_b_2 = NULL;
double *K_Hz_a_2 = NULL, *K_Hz_b_2 = NULL, *K_Hz_c_2 = NULL, *K_Hz_d_2 = NULL;
/////////////////////////
//plane wave excitation with total field - scatterd field formulation
/////////////////////////
CFDTD_1D_EH_PML_LOSS fdtd_1D;
long n_1D, n_1D_MIN_1, n_PML_1D = 20;
double d_1D;
//the interface between the total field-scattered field zone
double ct_Ex_1, ct_Ex_2, ct_Ey_1, ct_Ey_2, ct_Ez_1, ct_Ez_2;
double ct_Hx_1, ct_Hx_2, ct_Hy_1, ct_Hy_2, ct_Hz_1, ct_Hz_2;
//the incident field for total field scattered field formulation
double *E_1D = NULL, *H_1D = NULL;
double *ll_1D_E = NULL, *ll_1D_H = NULL;
//the recuired incident magnetic fields to compute the electric field components
double **Hz_i0 = NULL, **Hy_i0 = NULL, **Hz_i1 = NULL, **Hy_i1 = NULL;
double **Hz_j0 = NULL, **Hx_j0 = NULL, **Hz_j1 = NULL, **Hx_j1 = NULL;
double **Hy_k0 = NULL, **Hx_k0 = NULL, **Hy_k1 = NULL, **Hx_k1 = NULL;
double **face_Hz_i0 = NULL, **face_Hy_i0 = NULL, **face_Hz_i1 = NULL, **face_Hy_i1 = NULL;
double **face_Hz_j0 = NULL, **face_Hx_j0 = NULL, **face_Hz_j1 = NULL, **face_Hx_j1 = NULL;
double **face_Hy_k0 = NULL, **face_Hx_k0 = NULL, **face_Hy_k1 = NULL, **face_Hx_k1 = NULL;
//the recuired incident electric fields to compute the magnetic field components
double **Ey_i0 = NULL, **Ez_i0 = NULL, **Ey_i1 = NULL, **Ez_i1 = NULL;
double **Ex_j0 = NULL, **Ez_j0 = NULL, **Ex_j1 = NULL, **Ez_j1 = NULL;
double **Ex_k0 = NULL, **Ey_k0 = NULL, **Ex_k1 = NULL, **Ey_k1 = NULL;
double **face_Ey_i0 = NULL, **face_Ez_i0 = NULL, **face_Ey_i1 = NULL, **face_Ez_i1 = NULL;
double **face_Ex_j0 = NULL, **face_Ez_j0 = NULL, **face_Ex_j1 = NULL, **face_Ez_j1 = NULL;
double **face_Ex_k0 = NULL, **face_Ey_k0 = NULL, **face_Ex_k1 = NULL, **face_Ey_k1 = NULL;
/////////////////////////
//Contains the followed field components
/////////////////////////
double **Ex_Foll = NULL, **Ey_Foll = NULL, **Ez_Foll = NULL;
double **Hx_Foll = NULL, **Hy_Foll = NULL, **Hz_Foll = NULL;
double *Ex_Foll_av = NULL, *Ey_Foll_av = NULL, *Ez_Foll_av = NULL;
double *Hx_Foll_av = NULL, *Hy_Foll_av = NULL, *Hz_Foll_av = NULL;
double **Ex_Br_av_re, **Ex_Br_av_im, **Ey_Br_av_re, **Ey_Br_av_im, **Ez_Br_av_re, **Ez_Br_av_im;
double **Hx_Br_av_re, **Hx_Br_av_im, **Hy_Br_av_re, **Hy_Br_av_im, **Hz_Br_av_re, **Hz_Br_av_im;
double *Ex_Foll_av_WignerC = NULL, *Ey_Foll_av_WignerC = NULL, *Ez_Foll_av_WignerC = NULL;
double *Hx_Foll_av_WignerC = NULL, *Hy_Foll_av_WignerC = NULL, *Hz_Foll_av_WignerC = NULL;
/////////////////////////
//Fourier transform in a subvolume
/////////////////////////
double **Ex_fourier_TR_re = NULL, **Ex_fourier_TR_im = NULL;
double **Ey_fourier_TR_re = NULL, **Ey_fourier_TR_im = NULL;
double **Ez_fourier_TR_re = NULL, **Ez_fourier_TR_im = NULL;
double **Hx_fourier_TR_re = NULL, **Hx_fourier_TR_im = NULL;
double **Hy_fourier_TR_re = NULL, **Hy_fourier_TR_im = NULL;
double **Hz_fourier_TR_re = NULL, **Hz_fourier_TR_im = NULL;
double *fourier_omega = NULL;
//to save the field components
char *Path = NULL;
char *Path_Data = NULL;
char *path_name_E_1D = NULL;
char *path_name_Ex = NULL, *path_name_Ey = NULL, *path_name_Ez = NULL;
char *path_name_Hx = NULL, *path_name_Hy = NULL, *path_name_Hz = NULL;
bool er_Ex, er_Ey, er_Ez, er_Hx, er_Hy, er_Hz;
bool er_PML_Ex, er_PML_Ey, er_PML_Ez, er_PML_Hx, er_PML_Hy, er_PML_Hz;
bool er_TS; //plane wave excitation with total field - scattered field formulation
//contains the id of the threads
//computing threads
pthread_t *field_threads = NULL;
//data save threads
//xy - slice
pthread_t *data_save_threads_Ex_xy = NULL, *data_save_threads_Ey_xy = NULL,
*data_save_threads_Ez_xy = NULL, *data_save_threads_Hx_xy = NULL,
*data_save_threads_Hy_xy = NULL, *data_save_threads_Hz_xy = NULL;
//yz - slice
pthread_t *data_save_threads_Ex_yz = NULL, *data_save_threads_Ey_yz = NULL,
*data_save_threads_Ez_yz = NULL, *data_save_threads_Hx_yz = NULL,
*data_save_threads_Hy_yz = NULL, *data_save_threads_Hz_yz = NULL;
//xz - slice
pthread_t *data_save_threads_Ex_xz = NULL, *data_save_threads_Ey_xz = NULL,
*data_save_threads_Ez_xz = NULL, *data_save_threads_Hx_xz = NULL,
*data_save_threads_Hy_xz = NULL, *data_save_threads_Hz_xz = NULL;
#ifdef eval_W
double ***W = NULL;
double *W_aver = NULL;
double *W_aver_xz = NULL; //avereged in x and z direction
double *W_aver_xzt = NULL;//averaged in x and z direction and time
pthread_t *eval_W_thread = NULL;
Data_En *Data_W = NULL;
double *local_W_aver = NULL;
char *path_name_W = NULL;
pthread_t *data_save_threads_W_xy = NULL;
pthread_t *data_save_threads_W_yz = NULL;
pthread_t *data_save_threads_W_xz = NULL;
Data_Save_Slice *Save_W_xy = NULL;
Data_Save_Slice *Save_W_yz = NULL;
Data_Save_Slice *Save_W_xz = NULL;
bool jel_W_data_save_join = false;
#endif
//used to pass the data to threads
//E components
Data_E_comp_PML *Data_Ex_PML = NULL; //1-13 15-27 - inside of the PML
Data_E_comp *Data_Ex = NULL; //outside of the PML -- 14 zone
Data_E_comp_PML *Data_Ey_PML = NULL; //1-13 15-27 - inside of the PML
Data_E_comp *Data_Ey = NULL; //outside of the PML -- 14 zone
Data_E_comp_PML *Data_Ez_PML = NULL; //1-13 15-27 - inside of the PML
Data_E_comp *Data_Ez = NULL; //outside of the PML -- 14 zone
//H components
Data_H_comp_PML *Data_Hx_PML = NULL; //1-13 15-27 - inside of the PML
Data_H_comp *Data_Hx = NULL; //outside of the PML -- 14 zone
Data_H_comp_PML *Data_Hy_PML = NULL; //1-13 15-27 - inside of the PML
Data_H_comp *Data_Hy = NULL; //outside of the PML -- 14 zone
Data_H_comp_PML *Data_Hz_PML = NULL; //1-13 15-27 - inside of the PML
Data_H_comp *Data_Hz = NULL; //outside of the PML -- 14 zone
//used to pass data to save threads
//xy - slice
Data_Save_Slice *Save_Ex_xy = NULL;
Data_Save_Slice *Save_Ey_xy = NULL;
Data_Save_Slice *Save_Ez_xy = NULL;
Data_Save_Slice *Save_Hx_xy = NULL;
Data_Save_Slice *Save_Hy_xy = NULL;
Data_Save_Slice *Save_Hz_xy = NULL;
//yz - slice
Data_Save_Slice *Save_Ex_yz = NULL;
Data_Save_Slice *Save_Ey_yz = NULL;
Data_Save_Slice *Save_Ez_yz = NULL;
Data_Save_Slice *Save_Hx_yz = NULL;
Data_Save_Slice *Save_Hy_yz = NULL;
Data_Save_Slice *Save_Hz_yz = NULL;
//xz - slice
Data_Save_Slice *Save_Ex_xz = NULL;
Data_Save_Slice *Save_Ey_xz = NULL;
Data_Save_Slice *Save_Ez_xz = NULL;
Data_Save_Slice *Save_Hx_xz = NULL;
Data_Save_Slice *Save_Hy_xz = NULL;
Data_Save_Slice *Save_Hz_xz = NULL;
double **K_a_thread = NULL; double **K_b_thread = NULL;
double **mu_thread = NULL;
///////////////////////////////////////////////////////////////////////////////////////
//Allocate memory
///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////
//start parallel memory allocation region
///////////////////////////////////////////////////////////////////////////////////////
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -