?? fdtd_3d.cpp
字號:
K_Hx_d[ii] = -(2.0*eps_0*ka_x - sigma_x*dt)/mu_0;
}
//i+0.5
if (jel1 == 1)
{
sigma_x = sigma_max_1*pow( (nPML_x_1 - i - 0.5)/nPML_x_1 ,exponent);
ka_x = 1.0 + (ka_max - 1.0)*pow( (nPML_x_1 - i - 0.5)/nPML_x_1 ,exponent);
ii = cik;
K_Ex_c[ii] = 2.0*eps_0*ka_x + sigma_x*dt;
K_Ex_d[ii] = -(2.0*eps_0*ka_x - sigma_x*dt);
K_Hy_a[ii] = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
K_Hy_b[ii] = 1.0/(2.0*eps_0*ka_x+sigma_x*dt);
K_Bz_a[ii] = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
K_Bz_b[ii] = (2.0*eps_0*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
}
if (jel2 == 1)
{
sigma_x = sigma_max_2*pow( (nPML_x_2 - i - 0.5)/nPML_x_2 ,exponent);
ka_x = 1.0 + (ka_max - 1.0)*pow( (nPML_x_2 - i - 0.5)/nPML_x_2 ,exponent);
ii = nlx-cik-2;
if (ii >=0 )
{
K_Ex_c[ii] = 2.0*eps_0*ka_x + sigma_x*dt;
K_Ex_d[ii] = -(2.0*eps_0*ka_x - sigma_x*dt);
K_Hy_a[ii] = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
K_Hy_b[ii] = 1.0/(2.0*eps_0*ka_x+sigma_x*dt);
K_Bz_a[ii] = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
K_Bz_b[ii] = (2.0*eps_0*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
}
}
cik++;
}
if ( (iend < nx-1) && (iend >= nx - nPML_x_2 - 1) )
{
i = n1;
if (ii < 0)
{
i--;
}
sigma_x = sigma_max_2*pow( (nPML_x_2 - i - 0.5)/nPML_x_2 ,exponent);
ka_x = 1.0 + (ka_max - 1.0)*pow( (nPML_x_2 - i - 0.5)/nPML_x_2 ,exponent);
K_Ex_c[nlxMIN1] = 2.0*eps_0*ka_x + sigma_x*dt;
K_Ex_d[nlxMIN1] = -(2.0*eps_0*ka_x - sigma_x*dt);
K_Hy_a[nlxMIN1] = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
K_Hy_b[nlxMIN1] = 1.0/(2.0*eps_0*ka_x+sigma_x*dt);
K_Bz_a[nlxMIN1] = (2.0*eps_0*ka_x-sigma_x*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
K_Bz_b[nlxMIN1] = (2.0*eps_0*dt)/(2.0*eps_0*ka_x+sigma_x*dt);
}
return 0;
}
///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices in x directions
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_3D::Init_PML_Par_y(double eps_r_y_1, double mu_r_y_1, double eps_r_y_2,
double mu_r_y_2)
{
long j, jj;
long nly = jend - jsta + 1;
long nlyMIN1 = nly -1;
//Outside of the PML region
//j
for (j = 0; j < nlyMIN1; j++)
{
K_Gx_a[j] = 1.0;
K_Gx_b[j] = 1.0;
K_Ey_c[j] = 2.0*eps_0;
K_Ey_d[j] = -2.0*eps_0;
K_Ez_a[j] = 1.0;
K_Ez_b[j] = 1.0/(2.0*eps_0);
K_Bx_a[j] = 1.0;
K_Bx_b[j] = dt;
K_Hy_c[j] = 2*eps_0/mu_0;
K_Hy_d[j] = -2*eps_0/mu_0;
K_Hz_a[j] = 1.0;
K_Hz_b[j] = 1.0/(2.0*eps_0);
}
K_Gx_a[nlyMIN1] = 1.0;
K_Gx_b[nlyMIN1] = 1.0;
K_Ez_a[nlyMIN1] = 1.0;
K_Ez_b[nlyMIN1] = 1.0/(2.0*eps_0);
K_Hy_c[nlyMIN1] = 2*eps_0/mu_0;
K_Hy_d[nlyMIN1] = -2*eps_0/mu_0;
if ( jend < ny-1 )
{
K_Ey_c[nlyMIN1] = 2.0*eps_0;
K_Ey_d[nlyMIN1] = -2.0*eps_0;
K_Bx_a[nlyMIN1] = 1.0;
K_Bx_b[nlyMIN1] = dt;
K_Hz_a[nlyMIN1] = 1.0;
K_Hz_b[nlyMIN1] = 1.0/(2.0*eps_0);
}
//PML_y parameters
double ka_max = 1;
int exponent = 4;
double R_err = 1e-16;
double eta_1 = sqrt(mu_0*mu_r_y_1/eps_0/eps_r_y_1);
double eta_2 = sqrt(mu_0*mu_r_y_2/eps_0/eps_r_y_2);
double sigma_y, ka_y;
double sigma_max_1= -(exponent+1.0)*log(R_err)/(2.0*eta_1*nPML_y_1*dy);
double sigma_max_2= -(exponent+1.0)*log(R_err)/(2.0*eta_2*nPML_y_2*dy);
long n1 = 1, n2 = 0, jel1 = 0, jel2 = 0;
if (jend <= nPML_y_1-1)
{
n1 = jsta;
n2 = jend;
jel1 = 1;
}
if ( (jsta <= nPML_y_1-1) && (jend >= nPML_y_1-1) )
{
n1 = jsta;
n2 = nPML_y_1-1;
jel1 = 1;
}
if ( (jsta <= ny - nPML_y_2) && (jend > ny - nPML_y_2) )
{
n1 = ( ny - jend-1);
n2 = nPML_y_2-1;
jel2 = 1;
}
if ( jsta > ny - nPML_y_2)
{
n1 = 0;
n2 = nlyMIN1;
jel2 = 1;
}
long cik = 0;
for (j = n1; j <= n2; j++)
{
//j
if (jel1 == 1)
{
sigma_y = sigma_max_1*pow( (nPML_y_1 - j)/((double) nPML_y_1) ,exponent);
ka_y = 1.0 + (ka_max - 1.0)*pow( (nPML_y_1 - j)/((double) nPML_y_1) ,exponent);
jj = cik;
K_Gx_a[jj] = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Gx_b[jj] = 2.0*eps_0/(2.0*eps_0*ka_y+sigma_y*dt);
K_Ez_a[jj] = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Ez_b[jj] = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);
K_Hy_c[jj] = (2.0*eps_0*ka_y + sigma_y*dt)/mu_0;
K_Hy_d[jj] = -(2.0*eps_0*ka_y - sigma_y*dt)/mu_0;
}
if (jel2 == 1)
{
sigma_y = sigma_max_2*pow( (nPML_y_2 - j)/((double) nPML_y_2) ,exponent);
ka_y = 1.0 + (ka_max - 1.0)*pow( (nPML_y_2 - j)/((double) nPML_y_2) ,exponent);
jj = nly - cik - 1;
K_Gx_a[jj] = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Gx_b[jj] = 2.0*eps_0/(2.0*eps_0*ka_y+sigma_y*dt);
K_Ez_a[jj] = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Ez_b[jj] = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);
K_Hy_c[jj] = (2.0*eps_0*ka_y + sigma_y*dt)/mu_0;
K_Hy_d[jj] = -(2.0*eps_0*ka_y - sigma_y*dt)/mu_0;
}
//j+0.5
if (jel1 == 1)
{
sigma_y = sigma_max_1*pow( (nPML_y_1 - j - 0.5)/nPML_y_1 ,exponent);
ka_y = 1.0 + (ka_max - 1.0)*pow( (nPML_y_1 - j - 0.5)/nPML_y_1 ,exponent);
jj = cik;
K_Ey_c[jj] = 2.0*eps_0*ka_y + sigma_y*dt;
K_Ey_d[jj] = -(2.0*eps_0*ka_y - sigma_y*dt);
K_Bx_a[jj] = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Bx_b[jj] = (2.0*eps_0*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Hz_a[jj] = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Hz_b[jj] = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);
}
if (jel2 == 1)
{
sigma_y = sigma_max_2*pow( (nPML_y_2 - j - 0.5)/nPML_y_2 ,exponent);
ka_y = 1.0 + (ka_max - 1.0)*pow( (nPML_y_2 - j - 0.5)/nPML_y_2 ,exponent);
jj = nly-cik-2;
if (jj >=0 )
{
K_Ey_c[jj] = 2.0*eps_0*ka_y + sigma_y*dt;
K_Ey_d[jj] = -(2.0*eps_0*ka_y - sigma_y*dt);
K_Bx_a[jj] = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Bx_b[jj] = (2.0*eps_0*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Hz_a[jj] = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Hz_b[jj] = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);
}
}
cik++;
}
if ( (jend < ny-1) && (jend >= ny - nPML_y_2 - 1) )
{
j = n1;
if (jj < 0)
{
j--;
}
sigma_y = sigma_max_2*pow( (nPML_y_2 - j - 0.5)/nPML_y_2 ,exponent);
ka_y = 1.0 + (ka_max - 1.0)*pow( (nPML_y_2 - j - 0.5)/nPML_y_2 ,exponent);
K_Ey_c[nlyMIN1] = 2.0*eps_0*ka_y + sigma_y*dt;
K_Ey_d[nlyMIN1] = -(2.0*eps_0*ka_y - sigma_y*dt);
K_Bx_a[nlyMIN1] = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Bx_b[nlyMIN1] = (2.0*eps_0*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Hz_a[nlyMIN1] = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
K_Hz_b[nlyMIN1] = 1.0/(2.0*eps_0*ka_y + sigma_y*dt);
}
return 0;
}
///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices in z directions
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_3D::Init_PML_Par_z(double eps_r_z_1, double mu_r_z_1, double eps_r_z_2,
double mu_r_z_2)
{
long k, kk;
long nlz = kend - ksta + 1;
long nlzMIN1 = nlz -1;
//Outside of the PML region
//k
for (k = 0; k < nlzMIN1; k++)
{
K_Ex_a[k] = 1.0;
K_Ex_b[k] = 1.0/(2.0*eps_0);
K_Gy_a[k] = 1.0;
K_Gy_b[k] = 1.0;
K_Ez_c[k] = 2.0*eps_0;
K_Ez_d[k] = -2.0*eps_0;
K_Hx_a[k] = 1.0;
K_Hx_b[k] = 1.0/(2.0*eps_0);
K_By_a[k] = 1.0;
K_By_b[k] = dt;
K_Hz_c[k] = 2.0*eps_0/mu_0;
K_Hz_d[k] = -2.0*eps_0/mu_0;
}
K_Ex_a[nlzMIN1] = 1.0;
K_Ex_b[nlzMIN1] = 1.0/(2.0*eps_0);
K_Gy_a[nlzMIN1] = 1.0;
K_Gy_b[nlzMIN1] = 1.0;
K_Hz_c[nlzMIN1] = 2.0*eps_0/mu_0;
K_Hz_d[nlzMIN1] = -2.0*eps_0/mu_0;
if ( kend < nz-1 )
{
K_Ez_c[nlzMIN1] = 2.0*eps_0;
K_Ez_d[nlzMIN1] = -2.0*eps_0;
K_Hx_a[nlzMIN1] = 1.0;
K_Hx_b[nlzMIN1] = 1.0/(2.0*eps_0);
K_By_a[nlzMIN1] = 1.0;
K_By_b[nlzMIN1] = dt;
}
//PML_z parameters
double ka_max = 1;
int exponent = 4;
double R_err = 1e-16;
double eta_1 = sqrt(mu_0*mu_r_z_1/eps_0/eps_r_z_1);
double eta_2 = sqrt(mu_0*mu_r_z_2/eps_0/eps_r_z_2);
double sigma_z, ka_z;
double sigma_max_1= -(exponent+1.0)*log(R_err)/(2.0*eta_1*nPML_z_1*dz);
double sigma_max_2= -(exponent+1.0)*log(R_err)/(2.0*eta_2*nPML_z_2*dz);
long n1 = 1, n2 = 0, jel1 = 0, jel2 = 0;
if (kend <= nPML_z_1-1)
{
n1 = ksta;
n2 = kend;
jel1 = 1;
}
if ( (ksta <= nPML_z_1-1) && (kend >= nPML_z_1-1) )
{
n1 = ksta;
n2 = nPML_z_1-1;
jel1 = 1;
}
if ( (ksta <= nz - nPML_z_2) && (kend > nz - nPML_z_2) )
{
n1 = ( nz - kend - 1);
n2 = nPML_z_2-1;
jel2 = 1;
}
if ( ksta > nz - nPML_z_2)
{
n1 = 0;
n2 = nlzMIN1;
jel2 = 1;
}
long cik = 0;
for (k = n1; k <= n2; k++)
{
//k
if (jel1 == 1)
{
sigma_z = sigma_max_1*pow( (nPML_z_1 - k)/((double) nPML_z_1) ,exponent);
ka_z = 1.0 + (ka_max - 1.0)*pow( (nPML_z_1 - k)/((double) nPML_z_1) ,exponent);
kk = cik;
K_Ex_a[kk] = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
K_Ex_b[kk] = 1.0/(2.0*eps_0*ka_z + sigma_z*dt);
K_Gy_a[kk] = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
K_Gy_b[kk] = 2.0*eps_0/(2.0*eps_0*ka_z + sigma_z*dt);
K_Hz_c[kk] = (2.0*eps_0*ka_z + sigma_z*dt)/mu_0;
K_Hz_d[kk] = -(2.0*eps_0*ka_z - sigma_z*dt)/mu_0;
}
if (jel2 == 1)
{
sigma_z = sigma_max_2*pow( (nPML_z_2 - k)/((double) nPML_z_2) ,exponent);
ka_z = 1.0 + (ka_max - 1.0)*pow( (nPML_z_2 - k)/((double) nPML_z_2) ,exponent);
kk = nlz - cik - 1;
K_Ex_a[kk] = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
K_Ex_b[kk] = 1.0/(2.0*eps_0*ka_z + sigma_z*dt);
K_Gy_a[kk] = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
K_Gy_b[kk] = 2.0*eps_0/(2.0*eps_0*ka_z + sigma_z*dt);
K_Hz_c[kk] = (2.0*eps_0*ka_z + sigma_z*dt)/mu_0;
K_Hz_d[kk] = -(2.0*eps_0*ka_z - sigma_z*dt)/mu_0;
}
//k+0.5
if (jel1 == 1)
{
sigma_z = sigma_max_1*pow( (nPML_z_1 - k - 0.5)/nPML_z_1 ,exponent);
ka_z = 1.0 + (ka_max - 1.0)*pow( (nPML_z_1 - k - 0.5)/nPML_z_1 ,exponent);
kk = cik;
K_Ez_c[kk] = 2.0*eps_0*ka_z + sigma_z*dt;
K_Ez_d[kk] = -(2.0*eps_0*ka_z - sigma_z*dt);
K_Hx_a[kk] = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
K_Hx_b[kk] = 1.0/(2.0*eps_0*ka_z + sigma_z*dt);
K_By_a[kk] = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
K_By_b[kk] = (2.0*eps_0*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
}
if (jel2 == 1)
{
sigma_z = sigma_max_2*pow( (nPML_z_2 - k - 0.5)/nPML_z_2 ,exponent);
ka_z = 1.0 + (ka_max - 1.0)*pow( (nPML_z_2 - k - 0.5)/nPML_z_2 ,exponent);
kk = nlz - cik - 2;
if (kk >= 0 )
{
K_Ez_c[kk] = 2.0*eps_0*ka_z + sigma_z*dt;
K_Ez_d[kk] = -(2.0*eps_0*ka_z - sigma_z*dt);
K_Hx_a[kk] = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
K_Hx_b[kk] = 1.0/(2.0*eps_0*ka_z + sigma_z*dt);
K_By_a[kk] = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
K_By_b[kk] = (2.0*eps_0*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
}
}
cik++;
}
if ( (kend < nz-1) && (kend >= nz - nPML_z_2 - 1) )
{
k = n1;
if (kk < 0)
{
k--;
}
sigma_z = sigma_max_2*pow( (nPML_z_2 - k - 0.5)/nPML_z_2 ,exponent);
ka_z = 1.0 + (ka_max - 1.0)*pow( (nPML_z_2 - k - 0.5)/nPML_z_2 ,exponent);
K_Ez_c[nlzMIN1] = 2.0*eps_0*ka_z + sigma_z*dt;
K_Ez_d[nlzMIN1] = -(2.0*eps_0*ka_z - sigma_z*dt);
K_Hx_a[nlzMIN1] = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
K_Hx_b[nlzMIN1] = 1.0/(2.0*eps_0*ka_z + sigma_z*dt);
K_By_a[nlzMIN1] = (2.0*eps_0*ka_z - sigma_z*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
K_By_b[nlzMIN1] = (2.0*eps_0*dt)/(2.0*eps_0*ka_z + sigma_z*dt);
}
return 0;
}
///////////////////////////////////////////////////////////////////////////////////////
//Calculate the Ex field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D::Calc_Ex(long nlx, long nly, long nlz)
{
double Fx_r, Gx_r;
long i, j, k;
if (myrank_j == jprocsMIN1)
{
nly--;
}
if (myrank_k == kprocsMIN1)
{
nlz--;
}
#pragma omp parallel default(shared) private(i,j,k,Fx_r,Gx_r)
{
j = 0;
k = 0;
#pragma omp for schedule(dynamic,nr_threads) nowait
for (i = 0; i < nlx; i++)
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -