[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?

 

[LudusRex]

end subroutine test3


I would first suggest checking the boundary conditions in the code. It is possible that the infinite output is due to incorrect boundary conditions, which can cause waves to reflect infinitely within the simulation.

Additionally, I would recommend checking the values of the physical parameters used in the code (such as f0, miu0, eps0, etc.) to ensure they are appropriate for the problem being solved.

Another approach would be to compare the code with the reference mentioned (Numerical Electromagnetic: The FDTD method) and see if there are any discrepancies.

Lastly, I would suggest reaching out to the author of the reference or other experts in the field for further guidance and advice on troubleshooting the code.