- #1

s_hy

- 61

- 0

I have written codes for 2d fdtd tfsf with berenger's pml absorbing boundary. But I have serious leakage at front and right boundary. At the first place I think the problem is because the pml, but pml is working perfectly in the left and bottom boundary. I need an advice whether the problem is coming from tfsf boundary or absorbing boundary condition. Attached here is my code as well as output produced.

Code:

```
module data_module
implicit none
integer :: i,j,k,ie,je,ia,ib,ja,jb,hinit,hlast, n,nt,m0,iinit,jinit
double precision,parameter :: f0 = 3.0e7,c = 3.e8,eps0 = 8.8e-12,pi = 3.14159,miu0 = 4*pi*1e-7
double precision,parameter :: s = -0.477369,p = -1.0/3.0,q = -miu0*c/6,r = -miu0*c/2 ! mur absorbing boundary condition
double precision,parameter :: lamda = c/f0,dx = lamda/10,dt = dx/(6e8)
double precision,parameter :: A0 = 1.0
double precision,parameter :: Cb = 1/(2*eps0*c),Db = 1/(2*miu0*c)
double precision :: epsr,miur,phi,cosphi,sinphi
double precision,dimension(-10:1000) :: Einc0,Einc1,Hinc0
double precision,dimension(1000,1000) :: Ez,Hy,Hx
double precision,dimension(1000,1000) :: Ez1,Ez2,Ez3,Ez4,Ez0
double precision,dimension(-10:1000,-10:1000) :: Ezinc,Hxinc,Hyinc
double precision,dimension(1000,1000) :: Cax,Cbx,Dax,Dbx
double precision,dimension(1000,1000) :: miu,eps,sigma, sigmat
character(len=20) :: filename
end module data_module
! #################################################################################################
subroutine tfsfpml
use data_module
implicit none
integer :: npml,m,ip,jp
double precision :: smax,rmax,imp ! sigma max, rho max, impedance
double precision :: t,d,dp
double precision :: re,rm
double precision,dimension (500) :: ga,gb,ha,hb,sig,rho
phi = pi/4
cosphi = cos(phi)
sinphi = sin(phi)
m0 = 1
! model size
ie = 400
je = 400
hinit = -10
hlast = 500
iinit = 1
jinit = 1
! tfsf boundary
!ia = ie/4
!ja = je/4
!ib = ie*3/4
!jb = ib
ia = 20
ja = 20
ib = 280
jb = 280
! initialization for incident wave
Hinc0 = 0
Einc0 = 0
! initialization for field evolution
Hx = 0
Hy = 0
Ez = 0
Ez1 = 0
Ez2 = 0
Ez3 = complex(0.0,0.0)
Ez4 = 0
! :::::::::::::::::BERENGER'S ABC PML:::::::::::::::::::::::::::::::::::::::::::
! set up for berenger pml abc material
npml = 10
ip = ie - npml
jp = je - npml
imp = sqrt(miu0/eps0)
smax = -3.0*eps0*c*log(1e-5)/(2.0*dx*npml)
rmax = smax*(imp**2)
do m = 1,npml
sig(m)=smax*((m-0.5)/(npml+0.5))**2
rho(m)=rmax*(m/(npml+0.5))**2
end do
! set up constant needed in fdtd equations
do m = 1,npml+5
re = sig(m)*dt/eps0
rm = rho(m)*dt/miu0
ga(m) = exp(-re)
gb(m) = -(exp(-re)-1.0)/sig(m)/dx
ha(m) = exp(-rm)
hb(m) = -(exp(-rm)-1.0)/rho(m)/dx
end do
! medium's properties
miur = 1.0
epsr = 1.0
! initialize permittivity and permeability
do i = 1,ie
do j = 1,je
eps(i,j) = epsr*eps0
miu(i,j) = miur*miu0
end do
end do
! initialize electric conductivity & magnetic conductivity
do i = 1,ie
do j = 1,je
sigma(i,j) = 0.0
sigmat(i,j) = 0.0
Cax(i,j) = (1.- sigma(i,j)*dt/(2.*eps(i,j)))/(1. + (sigma(i, j)*dt)/(2.*eps(i,j)))
Cbx(i,j) = (dt/(eps(i,j)*dx))/(1.+sigma(i,j)*dt/(2.*eps(i,j)))
Dax(i,j) = (1.-sigmat(i,j)*dt/(2.*miu(i,j)))/(1. +(sigmat(i,j)*dt)/(2.*miu(i,j)))
Dbx(i,j) = (dt/(miu(i,j)*dx))/(1. + sigmat(i,j)*dt/(2.*miu(i,j)))
!print*, 'Cax(i,j)=',Cax(i,j),'Cbx(i,j)=',Cbx(i,j),'Dax(i,j)=',Dax(i,j),'Dbx(i,j)=',Dbx(i,j)
end do
end do
!::::: initialize all of the matrices for the Berenger pml absorbing boundaries ::::::::::::::
! Ez field - left and right pml regions
do i = 2,ie+50
do j = 2,npml
m = npml + 2 - j
Cax(i,j) = ga(m)
Cbx(i,j) = gb(m)
end do
do j = jp-1,je+5
m = j-jp
Cax(i,j) = ga(m)
Cbx(i,j) = gb(m)
end do
end do
! Ez field - back and front pml regions
do j = 2,je+50
do i = 2,npml
m = npml + 2 - i
Cax(i,j) = ga(m)
Cbx(i,j) = gb(m)
end do
do i = ip-1,ie+5
m = i-ip
Cax(i,j) = ga(m)
Cbx(i,j) = gb(m)
end do
end do
! Hx field - left and right pml regions
do i = 2,ie+50
do j = 1,npml
m = npml + 1-j
Dax(i,j) = ha(m)
Dbx(i,j) = hb(m)
end do
do j = jp-1,je+5
m = j-jp
Dax(i,j) = ha(m)
Dbx(i,j) = hb(m)
end do
end do
! Hy field - back and front pml regions
do j = 2,je+50
do i =1,npml
m = npml+1-i
Dax(i,j) = ha(m)
Dbx(i,j) = hb(m)
end do
do i = ip-1,ie+5
m = i-ip
Dax(i,j) = ha(m)
Dbx(i,j) = hb(m)
end do
end do
!-------------------------------------------------------------------------------------------------------------
! time steps start
do n = 1,1000
write (filename, "('data',I5.4,'.dat')") n
open (unit=31,file=filename)
! incident wave
t = n*dt
Einc1 = Einc0
do k = hinit,hlast
Hinc0(k) = Hinc0(k) - Db*(Einc1(k+1)-Einc1(k))
end do
do k = hinit+1,hlast-1
Einc0(k) = Einc1(k) - Cb*(Hinc0(k)-Hinc0(k-1))
end do
do i = ia,ib
do j = ja,jb
d = cosphi*(i - ia) + sinphi*(j - ja)
dp = d - int(d)
Ezinc(i,j) = (1-dp)*Einc0(m0+int(d))+dp*Einc0(m0+int(d+1))
end do
end do
Einc0(2) = A0*sin (2*pi*f0*n*dt)
k = hlast
Einc0(k) = Einc1(k-1)-1.0/3.0*(Einc0(k-1)-Einc1(k))
k = hinit
Einc0(k) = Einc1(k+1)-1.0/3.0*(Einc0(k+1)-Einc1(k))
! Hxincident field
do i = ia,ib
j = ja
d = cosphi*(i-ia)+sinphi*(j-0.5-ja)+0.5
dp = d - int(d)
Hxinc(i,j-1) = (dp*Hinc0(m0-1+int(d)+1)+(1-dp)*Hinc0(m0-1+int(d)))*sinphi
end do
do i = ia,ib
j = jb
d = cosphi*(i-ia)+sinphi*(j+0.5-ja)+0.5
dp = d - int(d)
Hxinc(i,j) = (dp*Hinc0(m0-1+int(d)+1)+(1-dp)*Hinc0(m0-1+int(d)))*sinphi
end do
! Hyincident field
do j = ja,jb
i = ia
d = cosphi*(i-0.5-ia) + sinphi*(j-ja)+0.5
dp = d - int(d)
Hyinc(i-1,j) = ((dp*Hinc0(m0-1+int(d)+1)) + (1-dp)*Hinc0(m0-1+int(d)))*(-cosphi)
end do
do j = ja,jb
i = ib
d = cosphi*(i+0.5-ia) + sinphi*(j-ja)+0.5
dp = d - int(d)
Hyinc(i,j) = ((dp*Hinc0(m0-1+int(d)+1))+(1-dp)*Hinc0(m0-1+int(d)))*(-cosphi)
end do
! update Hx field-------------------------------------------------------------------------------------------------------
do i = 1,ie
do j = 1,je-1
Hx(i,j) = Hx(i,j) - Dbx(i,jb)*(Ez(i,j+1)-Ez(i,j))
end do
end do
! tfsf for hx
do i = ia,ib
Hx(i,ja-1) = Hx(i,ja-1) + Dbx(i,ja-1)*Ezinc(i,ja)
end do
do i = ia,ib
Hx(i,jb+1) = Hx(i,jb+1) - Dbx(i,jb+1)*Ezinc(i,jb)
end do
! update Hy field-----------------------------------------------------------------------------------------------------
do i = 1,ie-1
do j = 1,je
Hy(i,j) = Hy(i,j) + Dbx(i,j)*(Ez(i+1,j)-Ez(i,j))
end do
end do
! tfsf for hy
do j = ja,jb
Hy(ia-1,j) = Hy(ia-1,j) - Dbx(ia-1,j)*Ezinc(ia,j)
end do
do j = ja,jb
Hy(ib+1,j) = Hy(ib+1,j) + Dbx(ib+1,j)*Ezinc(ib,j)
end do
! update for Ez-----------------------------------------------------------------------------------------------------------
do i = 2, ie-1
do j = 2,je-1
Ez(i,j) = Cax(i,j)*Ez(i,j) + Cbx(i,j)*(Hy(i,j) - Hy(i-1,j) - Hx(i,j) + Hx(i,j-1))
write (31,*) i,j,Ez(i,j)
if (j == 399) write (31,*) ' '
end do
end do
! update tfsf for Ez
do j = ja+1,jb-1
Ez(ia,j) = Ez(ia,j) - Cbx(ia,j)*Hyinc(ia-1,j)
end do
do j = ja+1,jb-1
Ez(ib,j) = Ez(ib,j) + Cbx(ib,j)*Hyinc(ib,j)
end do
do i = ia+1,ib-1
Ez(i,ja) = Ez(i,ja) + Cbx(i,ja)*Hxinc(i,ja-1)
end do
do i = ia+1,ib-1
Ez(i,jb) = Ez(i,jb) - Cbx(i,jb)*Hxinc(i,jb)
end do
close (unit=31)
end do
end subroutine tfsfpml
```