📄 fdtd_3d_bloch_pml.cpp
字号:
/////////////////////////////////////////////////////////////////
//main program
/////////////////////////////////////////////////////////////////
#include "run_enviroment.h"
#ifdef USE_MPI
//included due to a bug in MPI-2 standard
#undef SEEK_SET
#undef SEEK_END
#undef SEEK_CUR
#include "mpi.h"
#else
#include <ctime>
#define MPI_PROC_NULL 0
#endif
#include "read_input_file.h"
#include "FDTD_3D_COMPLEX.h"
int main(int argc, char* argv[])
{
char *file_name = NULL;
file_name = (char *) calloc(512,sizeof(char));
if (!file_name)
{
cout << "Memory allocation problem - *file_name" << endl;
return 1;
}
if (argc > 1)
{
strcpy(file_name,argv[1]);
}
else
{
strcpy(file_name,"Input_Data.txt");
}
///////////////////////////////////////////////////////
//initialization to read from file the input data
///////////////////////////////////////////////////////
Input_Data *Inp_D = NULL;
Inp_D = (Input_Data *) calloc(1,sizeof(Input_Data));
Inp_D->Path_Kx_Ky_Kz = NULL;
Inp_D->Path_Load_Workspace_Data = NULL;
Inp_D->path_name_CoordPointSource = NULL;
Inp_D->path_name_Index = NULL;
Inp_D->path_name_MatParam = NULL;
Inp_D->path_name_ParamPointSource = NULL;
if (readFile_fdtd_data(Inp_D,file_name))
{
cout << "Error reading Input_Data.dat" << endl;
return 1;
}
else
{
cout << "Input file: " << file_name << " loaded succesfull" << endl;
}
///////////////////////////////////////////////////////
///////////////////////////////////////////////////////
//initialization to start the MPI
///////////////////////////////////////////////////////
#ifdef USE_MPI
int myrank, nprocs;
MPI_Init(&argc, &argv); /* start MPI */
MPI_Comm_rank(MPI_COMM_WORLD, &myrank); /* get my proc id */
MPI_Comm_size(MPI_COMM_WORLD, &nprocs); /* get no.r of procs */
if (Inp_D->n_kxkykz != nprocs)
{
cout << "Wrong nr of processes\n" << endl;
MPI_Finalize(); //finish MPI
return 1;
}
cout << "MPI Process number " << myrank << " of " << nprocs << " is alive" << endl;
#endif
/////////////////////////////////////////////////////////////////////////////
//to measure the ellapsed time
/////////////////////////////////////////////////////////////////////////////
#ifdef run_enviroment //WIN
clock_t start, finish;
double duration;
start = clock();
#else //IBM AIX
double time_start, time_end, el_time;
struct timeval tv;
struct timezone tz;
gettimeofday(&tv, &tz);
time_start = (double) tv.tv_sec + (double) tv.tv_usec / 1000000.0;
#endif
double pi = 180.0*atan(1.0)/45.0;
//permittivity of free space
double eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]
//permeability of free space
double mu_0 = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6; // [H/m]
//the speed of the light in the vacuum
double speed_c = 1.0/sqrt(eps_0*mu_0); // [m/s]
long nx = Inp_D->nx; //total number of cells in x direction
long ny = Inp_D->ny; //total number of cells in y direction
long nz = Inp_D->nz; //total number of cells in z direction
//setting n_PML == 1 - PEC
long n_PML_x_1 = Inp_D->n_PML_x_1; //[cells] the size of PML-1 zone in x
long n_PML_x_2 = Inp_D->n_PML_x_2; //[cells] the size of PML-2 zone in x
long n_PML_y_1 = Inp_D->n_PML_y_1; //[cells] the size of PML-2 zone in y
long n_PML_y_2 = Inp_D->n_PML_y_2; //[cells] the size of PML-2 zone in y
long n_PML_z_1 = Inp_D->n_PML_z_1; //[cells] the size of PML-2 zone in z
long n_PML_z_2 = Inp_D->n_PML_z_2; //[cells] the size of PML-2 zone in z
double dx = Inp_D->dx;
double dy = Inp_D->dy;
double dz = Inp_D->dz;
double dt = ( 1.0/speed_c/sqrt(1.0/(dx*dx)+1.0/(dy*dy)+1.0/(dz*dz)) )*0.9;
cout << "dt = " << dt << endl;
long num_iter = Inp_D->num_iter; //the number of FDTD iterations
//the excitation
long source_type = Inp_D->source_type; // 1- Gauss; 2- Sin; 3- Gauss-Sin
long jel_plane_wave = Inp_D->jel_plane_wave; //the type of the excitation 0 - point source
//Reads from file the inverse lattice points
long n_kxkykz;
double **Kx_Ky_Kz = NULL;
long n_Coord;
long **Coord_ptSource = NULL;
double **Param_Source = NULL;
double X0 = 0.0, tw = 1.0, t0 = 0.0, omega = 0.0, phase = 0.0;
double teta = 90.0, phi = 0.0, gamma = 90.0;
double const_alfa = 0.0;
long i, ii = 0;
long n_TS = 0;
switch (jel_plane_wave)
{
case 0:
n_kxkykz = Inp_D->n_kxkykz;
Kx_Ky_Kz = Init_Matrix_2D<double>(n_kxkykz,3);
if (!Kx_Ky_Kz)
{
#ifdef USE_MPI
cout << "Error in the process : " << myrank << " - Memory allocation problem - **Kx_Ky_Kz" << endl;
MPI_Finalize();
#else
cout << "Memory allocation problem - **Kx_Ky_Kz" << endl;
#endif
return 1;
}
switch ( load_2D(Kx_Ky_Kz, n_kxkykz, 3, Inp_D->Path_Kx_Ky_Kz) )
{
case 0:
#ifdef USE_MPI
cout << "Process " << myrank << " : Inverse lattice points file loaded successfull" << endl;
#else
cout << "Inverse lattice points file loaded successfull" << endl;
#endif
break;
case 1:
#ifdef USE_MPI
cout << "Error in the process : " << myrank << " - Faild to open the file Kx_Ky_Kz" << endl;
MPI_Finalize();
#else
cout << "Faild to open the file Kx_Ky_Kz" << endl;
#endif
return 2;
case 2:
#ifdef USE_MPI
cout << "Error in the process : " << myrank << " - Wrong file content - Kx_Ky_Kz" << endl;
MPI_Finalize();
#else
cout << "Wrong file content - Kx_Ky_Kz" << endl;
#endif
return 3;
case 3:
#ifdef USE_MPI
cout << "Error in the process : " << myrank << " - Wrong dimension - Kx_Ky_Kz" << endl;
MPI_Finalize();
#else //IBM AIX
cout << "Wrong dimension - Kx_Ky_Kz" << endl;
#endif
return 4;
}
n_Coord = Inp_D->n_Coord;
Coord_ptSource = Init_Matrix_2D<long>(n_Coord,3);
if (!Coord_ptSource)
{
#ifdef USE_MPI
cout << "Error in the process : " << myrank << " - Memory allocation problem - **Coord_ptSource" << endl;
MPI_Finalize();
#else //IBM AIX
cout << "Memory allocation problem - **Coord_ptSource" << endl;
#endif
return 1;
}
switch ( load_2D_long(Coord_ptSource, n_Coord, 3, Inp_D->path_name_CoordPointSource ) )
{
case 0:
#ifdef USE_MPI
cout << "Process " << myrank << " : Coord_ptSource file loaded successfull" << endl;
#else //IBM AIX
cout << "Coord_ptSource file loaded successfull" << endl;
#endif
break;
case 1:
#ifdef USE_MPI
cout << "Error in the process : " << myrank << " - Failed to open the file Coord_ptSource" << endl;
MPI_Finalize();
#else //IBM AIX
cout << "Failed to open the file Coord_ptSource" << endl;
#endif
return 2;
case 2:
#ifdef USE_MPI
cout << "Error in the process : " << myrank << " - Wrong Coord_ptSource file content" << endl;
MPI_Finalize();
#else //IBM AIX
cout << "Wrong Coord_ptSource file content" << endl;
#endif
return 3;
case 3:
#ifdef USE_MPI
cout << "Error in the process : " << myrank << " - Wrong Coord_ptSource file dimension" << endl;
MPI_Finalize();
#else //IBM AIX
cout << "Wrong Coord_ptSource file dimension" << endl;
#endif
return 4;
}
switch (source_type)
{
case 1:
Param_Source = Init_Matrix_2D<double>(n_Coord,3);
break;
case 2:
Param_Source = Init_Matrix_2D<double>(n_Coord,4);
break;
case 3:
Param_Source = Init_Matrix_2D<double>(n_Coord,5);
break;
}
if (!Param_Source)
{
#ifdef USE_MPI
cout << "Error in the process : " << myrank << " - Memory allocation problem - **Param_Source" << endl;
MPI_Finalize();
#else //IBM AIX
cout << "Memory allocation problem - **Param_Source" << endl;
#endif
return 1;
}
X0 = Inp_D->X0; //amplitude
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;
break;
case 3://Gaussian-Sin source
tw = Inp_D->tw;
t0 = Inp_D->t0;
omega = Inp_D->omega;
break;
}
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
}
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] = 0; //phase
Param_Source[i][3] = omega*3.0/(2.0*pi); //const_alfa
}
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] = 0; //phase
}
break;
}
break;
case 1: //plane wave excitation with total - scattered field formulation
cout << "The plane wave excitation with the total Field - Scattered Field formulation is not implemented yet - Use point sources instead" << endl;
#ifdef USE_MPI
MPI_Finalize();
#endif
return 2;
//the direction and polarization of the incident field
teta = Inp_D->teta*pi/180.0; //angle relative to +z axis 0 < teta < 180
phi = Inp_D->phi*pi/180.0; //angle relative to +x axis 0 <= phi < 360
gamma = Inp_D->gamma*pi/180.0; //polarization
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;
const_alfa = Inp_D->const_alfa;
break;
case 3://Gaussian-Sin source
tw = Inp_D->tw;
t0 = Inp_D->t0;
omega = Inp_D->omega;
break;
}
n_TS = Inp_D->n_TS;
n_kxkykz = 1;
Kx_Ky_Kz = Init_Matrix_2D<double>(n_kxkykz,3);
if (!Kx_Ky_Kz)
{
#ifdef USE_MPI
cout << "Error in the process : " << myrank << " - Memory allocation problem - **Kx_Ky_Kz" << endl;
MPI_Finalize();
#else
cout << "Memory allocation problem - **Kx_Ky_Kz" << endl;
#endif
return 1;
}
Kx_Ky_Kz[0][0] = omega/speed_c*sin(teta)*cos(phi);
Kx_Ky_Kz[0][1] = omega/speed_c*sin(teta)*sin(phi);
Kx_Ky_Kz[0][2] = omega/speed_c*cos(teta);
break;
}
/////////////////////////////////////////////////////////////////////////////
//Initialization of the geometry
/////////////////////////////////////////////////////////////////////////////
long ***Index = NULL;
Index = Init_Matrix_3D<long>(nx,ny,nz);
if (!Index)
{
#ifdef USE_MPI
cout << "Error in the process : " << myrank << " - Memory allocation problem - ***Index" << endl;
MPI_Finalize();
#else //IBM AIX
cout << "Memory allocation problem - ***Index" << endl;
#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -