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

  1. 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.

    Code (Text):


    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

     
    [​IMG]
     
  2. jcsd
  3. jedishrfu

    Staff: Mentor

    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.
     
  4. SteamKing

    SteamKing 8,285
    Staff Emeritus
    Science Advisor
    Homework Helper

    You have some weird code:
    Code (Text):

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

    Code (Text):

     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.
     
  5. Steamking,

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

    DrDu 4,088
    Science Advisor

  7. 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
     
  8. DrDu

    DrDu 4,088
    Science Advisor

    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.
     
  9. 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?
     
  10. DrDu

    DrDu 4,088
    Science Advisor

    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.
     
Know someone interested in this topic? Share a link to this question via email, Google+, Twitter, or Facebook