Register to reply

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

Share this thread:
s_hy
#1
Mar17-14, 10:32 PM
P: 56
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
Phys.Org News Partner Physics news on Phys.org
Technique simplifies the creation of high-tech crystals
Working group explores the 'frustration' of spin glasses
New analysis of oxide glass structures could guide the forecasting of melt formation in planetary interiors
jedishrfu
#2
Mar17-14, 11:20 PM
P: 2,812
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
#3
Mar17-14, 11:55 PM
Emeritus
Sci Advisor
HW Helper
Thanks
PF Gold
P: 6,334
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
#4
Mar18-14, 12:50 AM
P: 56
[FDTD/FORTRAN] problem with tfsf boundary and berenger's pml

Steamking,

I have changed the weird part, but the output is still same.
DrDu
#5
Mar18-14, 03:11 AM
Sci Advisor
P: 3,564
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:
http://www.physicsforums.com/showthread.php?t=617567
s_hy
#6
Mar18-14, 03:18 AM
P: 56
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
#7
Mar18-14, 03:25 AM
Sci Advisor
P: 3,564
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
#8
Mar18-14, 03:51 AM
P: 56
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
#9
Mar18-14, 04:13 AM
Sci Advisor
P: 3,564
Quote Quote by s_hy View Post
. 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.


Register to reply

Related Discussions
How to Model Drude Materials using FDTD Atomic, Solid State, Comp. Physics 1
[FDTD/Fotran] detected reflection near boundary but don't know why Programming & Computer Science 3
[2D FDTD TFSF problem] got error in result but don't know why General Physics 0
How to put mur's absorbing boundary condition in 1d fdtd maxwell equation Programming & Computer Science 0
Fortran in boundary element Mechanical Engineering 0