Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

[Fortran90] fdtd in polar coordinates, got infinity output

  1. Jul 23, 2013 #1
    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

    Code (Text):

    !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

     
     
  2. jcsd
  3. Jul 24, 2013 #2
    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.
     
  4. Jul 24, 2013 #3
    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 dont need initialization here? but it change nothing even i disable it.

    Code (Text):


     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
     
  5. Jul 25, 2013 #4
    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?
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: [Fortran90] fdtd in polar coordinates, got infinity output
  1. Fortran90 help (Replies: 4)

  2. Fortran90 help (Replies: 5)

Loading...