[Fortran90] fdtd in polar coordinates, got infinity output

[s_hy]

hi all,

attached here is my code for 2d fdtd in polar coordinates, from 'numerical electromagnetic: the fdtd method (umran s inan, pg 94-96) written in fortran90. I have try a few approach I could think about to troubleshoot this code but the output is still infinity. Anybody here can give me advice. Thank you

!2d polar coordinates
!reference: Numerical electromagnetic: fdtd method (Umran S Inan, R Marshall)

subroutine test3
implicit none

double precision                   :: f0,miu0,eps0
double precision                   :: delta, deltat,delr,delphi
double precision                   :: c0
integer                                 :: i,j
integer                                 :: n,M
double precision,dimension(600,600):: Ez,Hp,Hr    !Ez,Hy,Hx
double precision, parameter             :: pi = 3.14159265
character(len=20)                  :: filename


 f0   = 1.e6
 c0   = 3.e8
 miu0 = 4.0*pi*1.0e-7
 eps0 = 8.854e-12
! deltat = 1/sqrt(2.0)*(delta/c0)

 delphi = 0.0314    !radian
 delr   = 15  
 !delt   = 1/sqrt(2.0)*delr*delphi/c0
 deltat  = delr*delphi/(c0*2)              
 !print*,'delt=',delt

 M = 200
!print*, 'see me?'

!initialize Ez, Hx,Hy to zero at t=0
do i = 1,M
 do j = 1,M
  Hp(i,j) = 0
  Hr(i,j) = 0
  Ez(i,j) = 0
 end do
end do

do n = 1,500
   
  write (filename, "('data',I3.3,'.dat')") n
  open (unit=130,file=filename)
 
!initiate sinusoidal wavepulse at center
  !Ez(1,1) = sin((n+1)*omega*deltat)
  Ez(100,100) = Ez(100,100) + sin(2*pi*f0*n*deltat)
  !print *, 'Ez(100,100)=',Ez(100,100)


!update electric field equation
  do i = 2,M
   do j = 2,M    
    Ez(i,j) = Ez(i,j)+(deltat/eps0)*((delr*(i)*Hp(i,j)-delr*(i-1)*Hp(i-1,j))/(delr*(i)*delr))-(Hr(i,j)-Hr(i,j-1)/(delr*(i)*delphi))
    write (130,*) i,j,Ez(i,j)
    if (j == 199) write (130,*) ' '
    !print *, 'i=',i,'j=',j,'Ez(i,j)=',Ez(i,j)
   end do
  end do


!update magnetic field equation
  do i = 1,M
   do j = 1,M
    Hr(i,j) = Hr(i,j) - (deltat/miu0)*((Ez(i,j+1) - Ez(i,j))/delr*(i)*delphi)
    print *,'i=',i,'j=',j,'Hr(i,j)=',Hr(i,j)
   end do
  end do

  do i = 1,M
   do j = 1,M
    Hp(i,j) = Hp(i,j) + (deltat/miu0)*((-Ez(i+1,j)+Ez(i,j))/delr)
    !print *,'i=',i,'j=',j,'Hp(i,j)=',Hp(i,j)
   end do
  end do


close (unit=130) 

end do !n

 

[gsal]

baby steps...baby steps!

...before you go on and open 500 files, make sure your program is working; try n=1,5, first.

One of the lines in your code is tooooo long, break it down with continuation marks...with fortran 90, I think you just need to stay within 132 columns wide.

Also, you double nested do-loop where you calculate Ez, it shows to be dependent on Hp and Hr, but they both are zero; so, I don't see how if ever you will get out of zero.

 

[s_hy]

eh...such a silly mistake...

changed...continuation marks and so on, the output is a little bit differ...as example -378 into -062...:D

is it i don't need initialization here? but it change nothing even i disable it.

 M = 200
!print*, 'see me?'

!initialize Ez, Hx,Hy to zero at t=0
do i = 1,M
 do j = 1,M
  Hp(i,j) = 0
  Hr(i,j) = 0
  Ez(i,j) = 0
 end do
end do

do n = 1,100
   
  write (filename, "('data',I3.3,'.dat')") n
  open (unit=130,file=filename)
 
!initiate sinusoidal wavepulse at center
  Ez(100,100) = sin(2*pi*f0*n*deltat)

!update electric field equation
  do i = 3,M
   do j = 3,M    
    Ez(i,j) = Ez(i,j)+(deltat/eps0)*((delr*(i)*Hp(i,j)-delr*(i-1)*Hp(i-1,j))&
/(delr*(i)*delr))-(Hr(i,j)-Hr(i,j-1)/(delr*(i)*delphi))
    write (130,*) i,j,Ez(i,j)
    if (j == 200) write (130,*) ' '
   end do
  end do

!update magnetic field equation
  do i = 2,M
   do j = 2,M
    Hr(i,j) = Hr(i,j) - (deltat/miu0)*((Ez(i,j+1) - Ez(i,j))/delr*(i)*delphi)
   end do
  end do

  do i = 2,M
   do j = 2,M
    Hp(i,j) = Hp(i,j) + (deltat/miu0)*((-Ez(i+1,j)+Ez(i,j))/delr)
   end do
  end do

close (unit=130) 

end do !n

end SUBROUTINE test3

thanks for your attention

 

[gsal]

You have not look into the zero problem.

Just because you are using a computer, it does not mean you are not allowed to use your own cpu (brain)...go ahead and mentally evaluate your Ez(i,j) for the very first iteration, i=3,j=3...what do you get?