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

 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
 P: 2,810 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.
 Emeritus Sci Advisor HW Helper Thanks PF Gold P: 6,333 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.
 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.
 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
 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
 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.
 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?