📄 fdtd_3d_xyzpml_threads_decomp.cpp
字号:
phase = Inp_D->phase;
break;
}
pt_source_Ez = Inp_D->pt_source_Ez;
pt_source_Hz = Inp_D->pt_source_Hz;
switch (source_type)
{
case 1:
for (i = 0; i<n_Coord; i++)
{
//Gaussian source
Param_Source[i][0] = X0; //J0
Param_Source[i][1] = tw; //tw
Param_Source[i][2] = t0; //t0
switch_off_time = num_iter;
}
break;
case 2:
for (i = 0; i<n_Coord; i++)
{
//Sinusoidal source
Param_Source[i][0] = X0; //J0
Param_Source[i][1] = omega; //omega
Param_Source[i][2] = phase; //phase
switch_off_time = num_iter;//(int) ( (n_Per*pi/2.0 - Param_Source[i][2])/(omega*dt) );
}
break;
case 3:
for (i = 0; i<n_Coord; i++)
{
//Sinusoidal modulated Gaussian source
Param_Source[i][0] = X0; //J0
Param_Source[i][1] = tw; //tw
Param_Source[i][2] = t0; //t0
Param_Source[i][3] = omega; //omega
Param_Source[i][4] = phase; //phase
switch_off_time = num_iter;
}
break;
}
break;
}
case 1: //plane wave excitation with total field - scattered field formulation
//the direction and polarization of the incident field
teta = Inp_D->teta*pi/180; //angle relative to +z axis 0 < teta < 180
phi = Inp_D->phi*pi/180; //angle relative to +x axis 0 <= phi < 360
gamma = Inp_D->gamma*pi/180; //polarization
cout << "teta = " << teta << endl;
cout << "phi = " << phi << endl;
cout << "gamma = " << gamma << endl;
X0 = Inp_D->X0; //amplituide
switch ( Inp_D->source_type)
{
case 1: //Gaussian source
tw = Inp_D->tw;
t0 = Inp_D->t0;
break;
case 2: //Sin source
omega = Inp_D->omega;
phase = Inp_D->phase;
const_alfa = Inp_D->const_alfa;
break;
case 3://Gaussian-Sin source
tw = Inp_D->tw;
t0 = Inp_D->t0;
omega = Inp_D->omega;
phase = Inp_D->phase;
break;
}
n_TS = Inp_D->n_TS;
//the position of the total field-scattered field region
n_TS_xa = nPML_x_1 + n_TS;
n_TS_xb = nx - nPML_x_2 - n_TS;
n_TS_ya = nPML_y_1 + n_TS;
n_TS_yb = ny - nPML_y_2 - n_TS;
n_TS_za = nPML_z_1 + n_TS;
n_TS_zb = nz - nPML_z_2 - n_TS;
break;
}
///////////////////////////////////////////////////////////////////////////////////////
//average over an arbitrary volume to compute the specular and/or the Bragg reflectance
///////////////////////////////////////////////////////////////////////////////////////
long **Mask_Wigner = NULL;
long n_Mask_Wigner;
double **Gxyz = NULL;
long n_Gxyz;
double kx, ky, kz;
if ( (Inp_D->aver_field_volume_WignerUC == 1 || Inp_D->aver_Bragg_WignerUC == 1)
&& jel_plane_wave == 1)
{
n_Mask_Wigner = Inp_D->n_Wigner_Cell;
Mask_Wigner = Init_Matrix_2D<long>(n_Mask_Wigner,3);
if (!Mask_Wigner)
{
//cout << "Memory allocation problem - **Mat" << endl;
printf("Memory allocation problem - **Mask_Wigner\n");
return -1;
}
switch (load_2D_long(Mask_Wigner, n_Mask_Wigner, 3, Inp_D->path_name_Wigner_Cell))
{
case 0:
printf("Mask_Wigner file loaded successful\n");
break;
case 1:
printf("faild to open the Mask_Wigner file \n");
return 2;
case 2:
printf("wrong Mask_Wigner file content\n");
return 3;
case 3:
printf("wrong Mask_Wigner file dimension\n");
return 3;
}
//average over a Wigner cell to compute the Bragg reflectance
if (Inp_D->aver_Bragg_WignerUC == 1)
{
n_Gxyz = Inp_D->n_Gxyz;
Gxyz = Init_Matrix_2D<double>(n_Gxyz,3);
if (!Gxyz)
{
cout << "Memory allocation problem - **Gxyz" << endl;
return 1;
}
switch (load_2D(Gxyz, n_Gxyz, 3, Inp_D->path_name_Gxyz))
{
case 0:
printf("Gxyz file loaded successful\n");
break;
case 1:
printf("faild to open the Gxyz file \n");
return 2;
case 2:
printf("wrong Gxyz file content\n");
return 3;
case 3:
printf("wrong Gxyz file dimension\n");
return 3;
}
for (i = 0; i< n_Gxyz; i++)
{
cout << "Gx = " << Gxyz[i][0] << " Gy = " << Gxyz[i][1]
<< " Gz = " << Gxyz[i][2] << endl;
}
kx = sin(teta)*cos(phi);
ky = sin(teta)*sin(phi);
kz = cos(teta);
}
}
/////////////////////////////////////////////////////////////////////////////
//Initialization of the geometry
/////////////////////////////////////////////////////////////////////////////
long ***Index = NULL;
Index = Init_Matrix_3D<long>(nx,ny,nz);
if (!Index)
{
printf("Memory allocation problem - ***Index\n");
return -1;
}
//Reads from file the elements of the Geometry matrix
switch ( load_3D_Geom_long(Index, nx, ny, nz, Inp_D->path_name_Index) )
{
case 0:
printf("Geom file loaded successful\n");
break;
case 1:
printf("faild to open the Geom file \n");
return 2;
case 2:
printf("wrong Geom file content\n");
return 3;
}
/////////////////////////////////////////////////////////////////////////////
double **Mat = NULL;
long n_Mat = Inp_D->n_Mat;
Mat = Init_Matrix_2D<double>(n_Mat,3);
if (!Mat)
{
//cout << "Memory allocation problem - **Mat" << endl;
printf("Memory allocation problem - **Mat\n");
return -1;
}
switch (load_2D(Mat, n_Mat, 3, Inp_D->path_name_MatParam))
{
case 0:
printf("Mat file loaded successful\n");
break;
case 1:
printf("faild to open the Mat file \n");
return 2;
case 2:
printf("wrong Mat file content\n");
return 3;
case 3:
printf("wrong Mat file dimension\n");
return 3;
}
//Set the points for the followed field components
long n_Ind_F;
long **Ind_F = NULL;
if (Inp_D->load_foll_field_points == 1)
{
n_Ind_F = Inp_D->n_Ind_F;
Ind_F = Init_Matrix_2D<long >(n_Ind_F,3);
if(!Ind_F)
{
cout << "Memory allocation problem - Ind_F" << endl;
return -1;
}
switch (load_2D_long(Ind_F, n_Ind_F, 3, Inp_D->path_name_foll_field_points))
{
case 0:
printf("Ind_F loaded successful\n");
break;
case 1:
printf("faild to open the Ind_F file \n");
return 2;
case 2:
printf("wrong Ind_F file content\n");
return 3;
default:
printf("error reading Ind_F from file\n");
return 4;
}
}
else //if there are no external Followed field points specified
{
n_Ind_F = 60;
Ind_F = Init_Matrix_2D<long >(n_Ind_F,3);
if(!Ind_F)
{
cout << "Memory allocation problem - Ind_F" << endl;
return -1;
}
long step_x = (nx-nPML_x_1-nPML_x_2-20)/19;
for (i = 0; i<20; i++)
{
Ind_F[i][0] = nPML_x_1 + n_TS + 1 + i*step_x; Ind_F[i][1] = ny/2; Ind_F[i][2] = nz/2;
}
long step_y = (ny-nPML_y_1-nPML_y_2-20)/19;
for (i = 20; i<40; i++)
{
Ind_F[i][0] = nx/2; Ind_F[i][1] = nPML_y_1 + n_TS + 1 + (i-20)*step_y; Ind_F[i][2] = nz/2;
}
long step_z = (nz-nPML_z_1-nPML_z_2-20)/19;
for (i = 40; i<60; i++)
{
Ind_F[i][0] = nx/2; Ind_F[i][1] = ny/2; Ind_F[i][2] = nPML_z_1 + n_TS + 1 + (i-40)*step_z;
}
}
/////////////////////////////////////////////////////////////////////////////
//auxiliary parameters to reduce the computational cost
/////////////////////////////////////////////////////////////////////////////
double TwoORpi = 2.0*pi;
double TwoOREp0 = 2.0*eps_0;
double inv_TwoOREp0 = 1.0/TwoOREp0;
double TwoOREp0DIVMu0 = TwoOREp0/mu_0;
double inv_dx = 1.0/dx;
double inv_dy = 1.0/dy;
double inv_dz = 1.0/dz;
double dtDIVMu0 = dt/mu_0;
double TwoOREp0DIVMu0ORdt = TwoOREp0DIVMu0*dt;
double inv_dxORmu_0 = 1.0/(dx*mu_0);
double inv_dyORmu_0 = 1.0/(dy*mu_0);
double inv_dzORmu_0 = 1.0/(dz*mu_0);
double dtDIVdx = dt/dx;
double dtDIVdy = dt/dy;
double dtDIVdz = dt/dz;
double dtDIVMu0DIVdx = dtDIVMu0/dx;
double dtDIVMu0DIVdy = dtDIVMu0/dy;
double dtDIVMu0DIVdz = dtDIVMu0/dz;
double TwoOREp0DIVMu0ORdtDIVdx = TwoOREp0DIVMu0ORdt/dx;
double TwoOREp0DIVMu0ORdtDIVdy = TwoOREp0DIVMu0ORdt/dy;
double TwoOREp0DIVMu0ORdtDIVdz = TwoOREp0DIVMu0ORdt/dz;
long nx_MIN_1 = nx - 1;
long ny_MIN_1 = ny - 1;
long nz_MIN_1 = nz - 1;
long nx_MIN_nPML_x_2 = nx - nPML_x_2;
long ny_MIN_nPML_y_2 = ny - nPML_y_2;
long nz_MIN_nPML_z_2 = nz - nPML_z_2;
long nx_MIN_1_MIN_nPML_x_2 = nx - 1 - nPML_x_2;
long ny_MIN_1_MIN_nPML_y_2 = ny - 1 - nPML_y_2;
long nz_MIN_1_MIN_nPML_z_2 = nz - 1 - nPML_z_2;
long nPML_x_2_MIN_1 = nPML_x_2 - 1;
long nPML_y_2_MIN_1 = nPML_y_2 - 1;
long nPML_z_2_MIN_1 = nPML_z_2 - 1;
long nxDIV2 = nx/2;
long nyDIV2 = ny/2;
long nzDIV2 = nz/2;
long Save_nr_DIV_xORy = Save_nr_DIV_x*Save_nr_DIV_y;
long Save_nr_DIV_yORz = Save_nr_DIV_y*Save_nr_DIV_z;
long Save_nr_DIV_xORz = Save_nr_DIV_x*Save_nr_DIV_z;
#ifdef eval_W
long Save_nr_W_DIV_xORy = Save_nr_W_DIV_x*Save_nr_W_DIV_y;
long Save_nr_W_DIV_yORz = Save_nr_W_DIV_y*Save_nr_W_DIV_z;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -