[FDTD/FORTRAN] problem with tfsf boundary and berenger's pml

[s_hy]

hi all...

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.

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

 

[jedishrfu]

Not understanding what you're trying to do but trying to help.

Often when a program mostly works but has some problem on the fringes of data computation means that perhaps your for loop indexes are off by 1 so say you want to loop ten times because you have ten values to average but an index is off by one then you actually loop by 9 or by 11 times meaning either one number was not averaged in or some undefined number was averaged in causing strange artifacts in your data as you used the average to compute over things.

So the bottom line is check your loops to be sure they are iterating over your data correctly.

 

[SteamKing]

You have some weird code:

 imp = sqrt(miu0/eps0)
 smax = -3.0*eps0*c*log(1e-5)/(2.0*dx*npml)
 rmax = smax*(imp**2)

Why not the following:

 imp = sqrt(miu0/eps0)
 smax = -3.0*eps0*c*log(1e-5)/(2.0*dx*npml)
 rmax = smax*(miu0/eps0)

Taking square roots and later squaring the result for re-use kills some of the precision of the original value.

 

[s_hy]

Steamking,

I have changed the weird part, but the output is still same.

 

[DrDu]

Do you really expect that everybody who might help you with your programming problem knows what abbreviations like "2d fdtd", "tfsf", "berenger's pml" means?
Have a look here:
https://www.physicsforums.com/showthread.php?t=617567

 

[s_hy]

sorry...I expected anyone who know FDTD are familiar with term tfsf and berenger's pml because you can find it in fdtd textbook everywhere. Sorry, I will more concern about this. Thank you for your advise.

FDTD : finite difference time domain
TFSf : total field scattered field
PML : perfectly matched layer

 

[DrDu]

I see lots of fishy stuff, which may make problems. In defining double precision constants, why do you use "e" instead of "d"? Pi is a double precision number, but you only assign it to 3.14159. Hence, cosphi and sinphi should both be equal to sqrt(2.0d0)/2.0d0. Are they?
Always write e.g. "1.0d0" instead of "1.0", only.

 

[s_hy]

I am newbie in Fortran. I may mislook about these things. I noted that pi should be assigned as pi = 3.14159265358979. for double precision constant, as example c = 3e8 means that c = light velocity = 3x10**8. about cosphi and sinphi, I don't understand why did you suggested that?

 

[DrDu]

. for double precision constant, as example c = 3e8 means that c = light velocity = 3x10**8. about cosphi and sinphi, I don't understand why did you suggested that?

I mean you should write c=3.0d8 instead of 3e8. d is for double precision numbers, e for single precision.
This cosphi and sinphi are used very often and imprecisions in pi may matter.