149 行
4.9 KiB
Fortran
149 行
4.9 KiB
Fortran
!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
|
|
!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
|
|
!Code distribution @ https://git.em3d.cn/
|
|
|
|
subroutine Get_pml_parameters
|
|
use constantparameters
|
|
USE PML_PARAMETER
|
|
implicit none
|
|
integer i,j,k,II,JJ,KK
|
|
DO i = 1,PML_X1
|
|
sig_PML_e_x1(i) = sig_x_max * ( (PML_X1 - i) / (PML_X1 - 1.0) )**ma
|
|
alpha_PML_e_x1(i) = alpha_x_max*((i-1.0)/(PML_X1-1.0))**mb
|
|
kappa_PML_e_x1(i) = 1.0+(kappa_x_max-1.0)*((PML_X1 - i) / (PML_X1 - 1.0))**ma
|
|
ENDDO
|
|
DO i = 1,PML_X1-1
|
|
sig_PML_h_x1(i) = sig_x_max * ( (PML_X1 - i - 0.5)/(PML_X1-1.0))**ma
|
|
alpha_PML_h_x1(i) = alpha_x_max*((i-0.5)/(PML_X1-1.0))**mb
|
|
kappa_PML_h_x1(i) = 1.0+(kappa_x_max-1.0)*((PML_X1 - i - 0.5) / (PML_X1 - 1.0))**ma
|
|
ENDDO
|
|
DO i = 1,PML_X2
|
|
sig_PML_e_x2(i) = sig_x_max * ( (PML_X2 - i) / (PML_X2 - 1.0) )**ma
|
|
alpha_PML_e_x2(i) = alpha_x_max*((i-1.0)/(PML_X2-1.0))**mb
|
|
kappa_PML_e_x2(i) = 1.0+(kappa_x_max-1.0)*((PML_X2 - i) / (PML_X2 - 1.0))**ma
|
|
ENDDO
|
|
DO i = 1,PML_X2-1
|
|
sig_PML_h_x2(i) = sig_x_max * ( (PML_X2 - i - 0.5)/(PML_X2-1.0))**ma
|
|
alpha_PML_h_x2(i) = alpha_x_max*((i-0.5)/(PML_X2-1.0))**mb
|
|
kappa_PML_h_x2(i) = 1.0+(kappa_x_max-1.0)*((PML_X2 - i - 0.5) / (PML_X2 - 1.0))**ma
|
|
ENDDO
|
|
|
|
!*************************************************************************************************
|
|
!y方向pml参数的求解
|
|
DO j = 1,PML_Y1
|
|
sig_PML_e_y1(j) = sig_y_max * ( (PML_Y1 - j ) / (PML_Y1 - 1.0) )**ma
|
|
alpha_PML_e_y1(j) = alpha_y_max*((j-1)/(PML_Y1-1.0))**mb
|
|
kappa_PML_e_y1(j) = 1.0+(kappa_y_max-1.0)*((PML_Y1 - j) / (PML_Y1 - 1.0))**ma
|
|
ENDDO
|
|
DO j = 1,PML_Y1-1
|
|
sig_PML_h_y1(j) = sig_y_max * ( (PML_Y1 - j - 0.5)/(PML_Y1-1.0))**ma
|
|
alpha_PML_h_y1(j) = alpha_y_max*((j-0.5)/(PML_Y1-1.0))**mb
|
|
kappa_PML_h_y1(j) = 1.0+(kappa_y_max-1.0)*((PML_Y1 - j - 0.5) / (PML_Y1 - 1.0))**ma
|
|
ENDDO
|
|
DO j = 1,PML_Y2
|
|
sig_PML_e_y2(j) = sig_y_max * ( (PML_Y2 - j ) / (PML_Y2 - 1.0) )**ma
|
|
alpha_PML_e_y2(j) = alpha_y_max*((j-1)/(PML_Y2-1.0))**mb
|
|
kappa_PML_e_y2(j) = 1.0+(kappa_y_max-1.0)*((PML_Y2 - j) / (PML_Y2 - 1.0))**ma
|
|
ENDDO
|
|
DO j = 1,PML_Y2-1
|
|
sig_PML_h_y2(j) = sig_y_max * ( (PML_Y2 - j - 0.5)/(PML_Y2-1.0))**ma
|
|
alpha_PML_h_y2(j) = alpha_y_max*((j-0.5)/(PML_Y2-1.0))**mb
|
|
kappa_PML_h_y2(j) = 1.0+(kappa_y_max-1.0)*((PML_Y2 - j - 0.5) / (PML_Y2 - 1.0))**ma
|
|
ENDDO
|
|
|
|
!*************************************************************************************************
|
|
!Z方向pml参数的求解
|
|
DO k = 1,PML_Z1
|
|
sig_PML_e_z1(k) = sig_z_max * ( (PML_Z1 - k ) / (PML_Z1 - 1.0) )**ma
|
|
alpha_PML_e_z1(k) = alpha_z_max*((k-1)/(PML_Z1-1.0))**mb
|
|
kappa_PML_e_z1(k) = 1.0+(kappa_z_max-1.0)*((PML_Z1 - k) / (PML_Z1 - 1.0))**ma
|
|
ENDDO
|
|
DO k = 1,PML_Z1-1
|
|
sig_PML_h_z1(k) = sig_z_max * ( (PML_Z1 - k - 0.5)/(PML_Z1-1.0))**ma
|
|
alpha_PML_h_z1(k) = alpha_z_max*((k-0.5)/(PML_Z1-1.0))**mb
|
|
kappa_PML_h_z1(k) = 1.0+(kappa_z_max-1.0)*((PML_Z1 - k - 0.5) / (PML_Z1 - 1.0))**ma
|
|
ENDDO
|
|
DO k = 1,PML_Z2
|
|
sig_PML_e_z2(k) = sig_z_max * ( (PML_Z2 - k ) / (PML_Z2 - 1.0) )**ma
|
|
alpha_PML_e_z2(k) = alpha_z_max*((k-1)/(PML_Z2-1.0))**mb
|
|
kappa_PML_e_z2(k) = 1.0+(kappa_z_max-1.0)*((PML_Z2 - k) / (PML_Z2 - 1.0))**ma
|
|
ENDDO
|
|
DO k = 1,PML_Z2-1
|
|
sig_PML_h_z2(k) = sig_z_max * ( (PML_Z2 - k - 0.5)/(PML_Z2-1.0))**ma
|
|
alpha_PML_h_z2(k) = alpha_z_max*((k-0.5)/(PML_Z2-1.0))**mb
|
|
kappa_PML_h_z2(k) = 1.0+(kappa_z_max-1.0)*((PML_Z2 - k - 0.5) / (PML_Z2 - 1.0))**ma
|
|
ENDDO
|
|
|
|
!求解den
|
|
!x方向
|
|
ii =PML_X2
|
|
DO i = 1,NX
|
|
if (i <= PML_X1) then
|
|
den_ex(i) = 1.0/kappa_PML_e_x1(i)
|
|
elseif (i >= NX+2-PML_X2) then
|
|
den_ex(i) = 1.0/kappa_PML_e_x2(ii)
|
|
ii = ii-1
|
|
else
|
|
den_ex(i) = 1.0
|
|
endif
|
|
ENDDO
|
|
ii =PML_X2-1
|
|
DO i = 1,NX
|
|
if (i <= PML_X1-1) then
|
|
den_hx(i) = 1.0/kappa_PML_h_x1(i)
|
|
elseif (i >= NX+2-PML_X2) then
|
|
den_hx(i) = 1.0/kappa_PML_h_x2(ii)
|
|
ii = ii-1
|
|
else
|
|
den_hx(i) = 1.0
|
|
endif
|
|
ENDDO
|
|
!y方向
|
|
jj = PML_Y2
|
|
DO j = 1,NY
|
|
if (j <= PML_Y1) then
|
|
den_ey(j) = 1.0/kappa_PML_e_y1(j)
|
|
elseif (j >= NY+2-PML_Y2) then
|
|
den_ey(j) = 1.0/kappa_PML_e_y2(jj)
|
|
jj = jj-1
|
|
else
|
|
den_ey(j) = 1.0
|
|
endif
|
|
ENDDO
|
|
jj =PML_Y2-1
|
|
DO j = 1,NY
|
|
if (j <= PML_Y1-1) then
|
|
den_hy(j) = 1.0/kappa_PML_h_y1(j)
|
|
elseif (j >= NY+2-PML_Y2) then
|
|
den_hy(j) = 1.0/kappa_PML_h_y2(jj)
|
|
jj = jj-1
|
|
else
|
|
den_hy(j) = 1.0
|
|
endif
|
|
ENDDO
|
|
!z方向
|
|
kk =PML_Z2
|
|
DO k = 1,NZ
|
|
if (k <= PML_Z1) then
|
|
den_ez(k) = 1.0/kappa_PML_e_z1(k)
|
|
elseif (k >= NZ+2-PML_Z2) then
|
|
den_ez(k) = 1.0/kappa_PML_e_z2(kk)
|
|
kk = kk - 1
|
|
else
|
|
den_ez(k) = 1.0
|
|
endif
|
|
ENDDO
|
|
kk =PML_Z2-1
|
|
DO k = 1,NZ
|
|
if (k <= PML_Z1-1) then
|
|
den_hz(k) = 1.0/kappa_PML_h_z1(k)
|
|
elseif (k >= NZ+2-PML_Z2) then
|
|
den_hz(k) = 1.0/kappa_PML_h_z2(kk)
|
|
kk = kk - 1
|
|
else
|
|
den_hz(k) = 1.0
|
|
endif
|
|
ENDDO
|
|
|
|
end subroutine Get_pml_parameters
|