📄 fdtd_3d.cpp
字号:
#include "fdtd_3d.h"
#include "Matrix.h"
#include "Load_Save_File_Data.h"
#include <math.h>
///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices in x directions
///////////////////////////////////////////////////////////////////////////////////////
void Init_PML_Par_x(double *K_Ey_a, double *K_Ey_b, double *K_Gz_a, double *K_Gz_b,
double *K_Hx_c, double *K_Hx_d, double *K_Ex_c, double *K_Ex_d,
double *K_Hy_a, double *K_Hy_b, double *K_Bz_a, double *K_Bz_b,
double eps_r_x, double mu_r_x, long nPML_x, double dx, double dt)
{
long i;
//permittivity of free space
double eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]
//permeability of free space
double mu_0 = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6; // [H/m]
//PML_x parameters
double ka_max = 1.0;
int exponent = 4;
double R_err = 1e-16;
double eta = sqrt(mu_0*mu_r_x/eps_0/eps_r_x);
double sigma_x, ka_x;
double sigma_max= -(exponent+1)*log(R_err)/(2*eta*nPML_x*dx);
for (i = 0; i < nPML_x; i++)
{
//i
sigma_x = sigma_max*pow( (nPML_x - i)/((double) nPML_x) ,exponent);
ka_x = 1.0 + (ka_max - 1.0)*pow( (nPML_x - i)/((double) nPML_x) ,exponent);
K_Ey_a[i] = (2.0*eps_0*ka_x - sigma_x*dt)/(2.0*eps_0*ka_x + sigma_x*dt);
K_Ey_b[i] = 1.0/(2.0*eps_0*ka_x + sigma_x*dt);
K_Gz_a[i] = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
K_Gz_b[i] = 2.0*eps_0/(2.0*eps_0*ka_x + sigma_x*dt);
K_Hx_c[i] = (2.0*eps_0*ka_x + sigma_x*dt)/mu_0;
K_Hx_d[i] = -(2.0*eps_0*ka_x - sigma_x*dt)/mu_0;
//i+0.5
sigma_x = sigma_max*pow( (nPML_x - i - 0.5)/nPML_x ,exponent);
ka_x = 1.0 + (ka_max - 1.0)*pow( (nPML_x - i - 0.5)/nPML_x ,exponent);
K_Ex_c[i] = 2.0*eps_0*ka_x + sigma_x*dt;
K_Ex_d[i] = -(2.0*eps_0*ka_x - sigma_x*dt);
K_Hy_a[i] = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
K_Hy_b[i] = 1.0/(2.0*eps_0*ka_x+sigma_x*dt);
K_Bz_a[i] = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
K_Bz_b[i] = (2.0*eps_0*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
}
}
///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices in y directions
///////////////////////////////////////////////////////////////////////////////////////
void Init_PML_Par_y(double *K_Gx_a, double *K_Gx_b, double *K_Ez_a, double *K_Ez_b,
double *K_Hy_c, double *K_Hy_d, double *K_Ey_c, double *K_Ey_d,
double *K_Bx_a, double *K_Bx_b, double *K_Hz_a, double *K_Hz_b,
double eps_r_y, double mu_r_y, long nPML_y, double dy, double dt)
{
long j;
//permittivity of free space
double eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]
//permeability of free space
double mu_0 = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6; // [H/m]
//PML_y parameters
double ka_max = 1.0;
int exponent = 4;
double R_err = 1e-16;
//PML_y parameters
double eta = sqrt(mu_0*mu_r_y/eps_0/eps_r_y);
double sigma_y, ka_y;
double sigma_max= -(exponent+1)*log(R_err)/(2*eta*nPML_y*dy);
for (j = 0; j < nPML_y; j++)
{
//j
sigma_y = sigma_max*pow( (nPML_y - j)/((double) nPML_y) ,exponent);
ka_y = 1.0 + (ka_max - 1.0)*pow( (nPML_y - j)/((double) nPML_y) ,exponent);
K_Gx_a[j] = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Gx_b[j] = 2.0*eps_0/(2.0*eps_0*ka_y+sigma_y*dt);
K_Ez_a[j] = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Ez_b[j] = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);
K_Hy_c[j] = (2.0*eps_0*ka_y + sigma_y*dt)/mu_0;
K_Hy_d[j] = -(2.0*eps_0*ka_y - sigma_y*dt)/mu_0;
//j+0.5
sigma_y = sigma_max*pow( (nPML_y - j - 0.5)/nPML_y ,exponent);
ka_y = 1.0 + (ka_max - 1.0)*pow( (nPML_y - j - 0.5)/nPML_y ,exponent);
K_Ey_c[j] = 2.0*eps_0*ka_y + sigma_y*dt;
K_Ey_d[j] = -(2.0*eps_0*ka_y - sigma_y*dt);
K_Bx_a[j] = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Bx_b[j] = (2.0*eps_0*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Hz_a[j] = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Hz_b[j] = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);
}
}
///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices in z directions - 1 layer
///////////////////////////////////////////////////////////////////////////////////////
void Init_PML_Par_z(double *K_Ex_a, double *K_Ex_b, double *K_Gy_a, double *K_Gy_b,
double *K_Hz_c, double *K_Hz_d, double *K_Ez_c, double *K_Ez_d,
double *K_Hx_a, double *K_Hx_b, double *K_By_a, double *K_By_b,
double eps_r_z, double mu_r_z, long nPML_z, double dz, double dt)
{
long k;
//permittivity of free space
double eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]
//permeability of free space
double mu_0 = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6; // [H/m]
//PML_y parameters
double ka_max = 1.0;
int exponent = 4;
double R_err = 1e-16;
//PML_z parameters
double eta = sqrt(mu_0*mu_r_z/eps_0/eps_r_z);
double sigma_z, ka_z;
double sigma_max= -(exponent+1)*log(R_err)/(2*eta*nPML_z*dz);
for (k = 0; k < nPML_z; k++)
{
//k
sigma_z = sigma_max*pow( (nPML_z - k)/((double) nPML_z) ,exponent);
ka_z = 1.0 + (ka_max - 1.0)*pow( (nPML_z - k)/((double) nPML_z) ,exponent);
K_Ex_a[k] = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
K_Ex_b[k] = 1.0/(2.0*eps_0*ka_z + sigma_z*dt);
K_Gy_a[k] = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
K_Gy_b[k] = 2.0*eps_0/(2.0*eps_0*ka_z + sigma_z*dt);
K_Hz_c[k] = (2.0*eps_0*ka_z + sigma_z*dt)/mu_0;
K_Hz_d[k] = -(2.0*eps_0*ka_z - sigma_z*dt)/mu_0;
//k+0.5
sigma_z = sigma_max*pow( (nPML_z - k - 0.5)/nPML_z ,exponent);
ka_z = 1.0 + (ka_max - 1.0)*pow( (nPML_z - k - 0.5)/nPML_z ,exponent);
K_Ez_c[k] = 2.0*eps_0*ka_z + sigma_z*dt;
K_Ez_d[k] = -(2.0*eps_0*ka_z - sigma_z*dt);
K_Hx_a[k] = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
K_Hx_b[k] = 1.0/(2.0*eps_0*ka_z + sigma_z*dt);
K_By_a[k] = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
K_By_b[k] = (2.0*eps_0*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
}
}
///////////////////////////////////////////////////////////////////////////////////////
//Curent source J
///////////////////////////////////////////////////////////////////////////////////////
void ptSource_J(double ***&xx, double **&Par_ptS, long **Coord_ptSource, double time,
double dt, long n_ptSource, int source_type, double const_alfa)
{
long i, i1, j1, k1;
double J0, omega, phi, t0, tw, alfa;
//double eps_0 = 8.8541878176203898505365630317107e-12; // [F/m]
//double dtDIVeps0 = dt/eps_0;
//time = time - dt/2;
for (i = 0; i < n_ptSource; i++ )
{
i1 = Coord_ptSource[i][0];
j1 = Coord_ptSource[i][1];
k1 = Coord_ptSource[i][2];
J0 = Par_ptS[i][0];//*dtDIVeps0/Mat[Ind[i1][j1][k1]][0];
switch (source_type)
{
case 1: //Gaussian pulse
tw = Par_ptS[i][1];
t0 = Par_ptS[i][2];
xx[i1][j1][k1] += J0*exp( -pow( (time - t0)/tw ,2) );
break;
case 2: //Sinusoidal plane wave
omega = Par_ptS[i][1];
phi = Par_ptS[i][2];
alfa = omega*const_alfa;
xx[i1][j1][k1] += J0*(1.0-exp(-alfa*time))*cos(omega*time + phi);
break;
case 3: //Wave packet(sinusoidal modulated Gaussian pulse)
tw = Par_ptS[i][1];
t0 = Par_ptS[i][2];
omega = Par_ptS[i][3];
phi = Par_ptS[i][4];
xx[i1][j1][k1] += J0*cos( omega*(time-t0) + phi)*
exp( -pow( (time - t0)/tw ,2) );
}
}
}
///////////////////////////////////////////////////////////////////////////////////////
//Add the incident field to Ey on total-field scattered field boundary - y direction
///////////////////////////////////////////////////////////////////////////////////////
void TS_Ey_ydir(double ***Ey, long ***Ind, double *K_b, double *Hx_1D, long n_TS_xa,
long n_TS_xb, long n_TS_ya, long n_TS_ybMIN1, long n_TS_za, long n_TS_zbPL1,
long n_TS_yaMIN1, double inv_dz)
{
long i, j, jj;
for (i = n_TS_xa; i <= n_TS_xb; i++)
{
for (j = n_TS_ya; j <= n_TS_ybMIN1; j++)
{
jj = j - n_TS_yaMIN1;
Ey[i][j][n_TS_za] -= K_b[Ind[i][j][n_TS_za]]*inv_dz*Hx_1D[jj];
Ey[i][j][n_TS_zbPL1] += K_b[Ind[i][j][n_TS_zbPL1]]*inv_dz*Hx_1D[jj];
}
}
}
///////////////////////////////////////////////////////////////////////////////////////
//Add the incident field to Ez on total-field scattered field boundary - y direction
///////////////////////////////////////////////////////////////////////////////////////
void TS_Ez_ydir(double ***Ez, long ***Ind, double *K_b, double *Hx_1D, long n_TS_xa,
long n_TS_xb, long n_TS_ya, long n_TS_yb, long n_TS_za, long n_TS_zb,
long n_TS_ybMINn_TS_yaPL1, double inv_dy)
{
long i, k;
for (i = n_TS_xa; i <= n_TS_xb; i++)
{
for (k = n_TS_za; k <= n_TS_zb; k++)
{
Ez[i][n_TS_ya][k] += K_b[Ind[i][n_TS_ya][k]]*inv_dy*Hx_1D[0];
Ez[i][n_TS_yb][k] -= K_b[Ind[i][n_TS_yb][k]]*inv_dy*Hx_1D[n_TS_ybMINn_TS_yaPL1];
}
}
}
///////////////////////////////////////////////////////////////////////////////////////
//Add the incident field to Hx on total-field scattered field boundary - y direction
///////////////////////////////////////////////////////////////////////////////////////
void TS_Hx_ydir(double ***Hx, long ***Ind, double **Mat, double *Ez_1D, long n_TS_xa,
long n_TS_xb, long n_TS_yaMIN1, long n_TS_yb, long n_TS_za, long n_TS_zb,
long n_TS_ybMINn_TS_yaPL1, double coef_TS_Hx)
{
long i, k;
double mu_r;
for ( i = n_TS_xa; i <= n_TS_xb; i++ )
{
for ( k = n_TS_za; k <= n_TS_zb; k++ )
{
mu_r = Mat[Ind[i][n_TS_yaMIN1][k]][1];
Hx[i][n_TS_yaMIN1][k] += coef_TS_Hx*Ez_1D[1]/mu_r;
mu_r = Mat[Ind[i][n_TS_yb][k]][1];
Hx[i][n_TS_yb][k] -= coef_TS_Hx*Ez_1D[n_TS_ybMINn_TS_yaPL1]/mu_r;
}
}
}
///////////////////////////////////////////////////////////////////////////////////////
//Add the incident field to Hy on total-field scattered field boundary - y direction
///////////////////////////////////////////////////////////////////////////////////////
void TS_Hy_ydir(double ***Hy, long ***Ind, double **Mat, double *Ez_1D, long n_TS_xaMIN1,
long n_TS_xb, long n_TS_ya, long n_TS_yb, long n_TS_za, long n_TS_zb,
long n_TS_yaMIN1, double coef_TS_Hy)
{
long j, k;
double mu_r;
for ( j = n_TS_ya; j <= n_TS_yb; j++ )
{
for ( k = n_TS_za; k <= n_TS_zb; k++ )
{
mu_r = Mat[Ind[n_TS_xaMIN1][j][k]][1];
Hy[n_TS_xaMIN1][j][k] -= coef_TS_Hy*Ez_1D[j-n_TS_yaMIN1]/mu_r;
mu_r = Mat[Ind[n_TS_xb][j][k]][1];
Hy[n_TS_xb][j][k] += coef_TS_Hy*Ez_1D[j-n_TS_yaMIN1]/mu_r;
}
}
}
int Init_TS_Param(double **ll_1D_E, double **ll_1D_H,
double ***Hz_i0, double ***Hy_i0, double ***Hz_i1, double ***Hy_i1,
double ***Hz_j0, double ***Hx_j0, double ***Hz_j1, double ***Hx_j1,
double ***Hy_k0, double ***Hx_k0, double ***Hy_k1, double ***Hx_k1,
double ***face_Hz_i0, double ***face_Hy_i0, double ***face_Hz_i1,
double ***face_Hy_i1, double ***face_Hz_j0, double ***face_Hx_j0,
double ***face_Hz_j1, double ***face_Hx_j1, double ***face_Hy_k0,
double ***face_Hx_k0, double ***face_Hy_k1, double ***face_Hx_k1,
double ***Ey_i0, double ***Ez_i0, double ***Ey_i1, double ***Ez_i1,
double ***Ex_j0, double ***Ez_j0, double ***Ex_j1, double ***Ez_j1,
double ***Ex_k0, double ***Ey_k0, double ***Ex_k1, double ***Ey_k1,
double ***face_Ey_i0, double ***face_Ez_i0, double ***face_Ey_i1, double ***face_Ez_i1,
double ***face_Ex_j0, double ***face_Ez_j0, double ***face_Ex_j1, double ***face_Ez_j1,
double ***face_Ex_k0, double ***face_Ey_k0, double ***face_Ex_k1, double ***face_Ey_k1,
long l_1D, long l_x, long l_y, long l_z, long n_PML, double teta, double phi,
double gamma, double dx, double dy, double dz, double d, double dt,
double &ct_Ex_1, double &ct_Ex_2, double &ct_Ey_1, double &ct_Ey_2,
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -