[Fortran90] fdtd in polar coordinates, got infinity output

In summary, the conversation is about a person seeking advice on a code they wrote for 2D FDTD in polar coordinates using Fortran90. They have tried different approaches to troubleshoot the code, but it still outputs infinity. The person shares their code and mentions that they are getting an error for a line that is too long. They also mention that they are unsure if they need to initialize some variables or not. Lastly, they acknowledge a mistake in their code and provide an updated version.
  • #1
s_hy
61
0
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:
!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
 
Technology news on Phys.org
  • #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.
 
  • #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 don't need initialization here? but it change nothing even i disable it.

Code:
 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
 
  • #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?
 
  • #5


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.
 

Related to [Fortran90] fdtd in polar coordinates, got infinity output

1. What is FDTD in polar coordinates?

FDTD (Finite-Difference Time-Domain) is a numerical method used to solve electromagnetic field equations. In polar coordinates, the equations are solved in a spherical coordinate system, rather than the more commonly used Cartesian coordinate system.

2. Why am I getting an infinity output when running my FDTD simulation?

This could be due to a number of reasons, such as incorrect boundary conditions, improper initialization of the simulation, or numerical instability. It is important to carefully check your code and make sure all inputs and parameters are correct.

3. How can I troubleshoot my FDTD simulation in polar coordinates?

First, check your code for any errors or potential sources of numerical instability. You can also try adjusting your boundary conditions or initial conditions to see if that affects the output. Additionally, using smaller time steps or grid sizes may help resolve the issue.

4. Is FDTD in polar coordinates better than FDTD in Cartesian coordinates?

It depends on the specific problem being solved. FDTD in polar coordinates may be more suitable for certain types of simulations, such as those involving spherical objects or radiation patterns. However, it may also be more computationally intensive and require more complex code compared to FDTD in Cartesian coordinates.

5. Are there any resources available to help me with FDTD in polar coordinates?

Yes, there are many online resources, tutorials, and textbooks available that cover FDTD in polar coordinates. It may also be helpful to consult with other scientists or researchers who have experience with FDTD simulations.

Similar threads

  • Programming and Computer Science
Replies
4
Views
3K
  • Programming and Computer Science
Replies
4
Views
2K
Replies
1
Views
4K
  • Programming and Computer Science
Replies
4
Views
3K
  • Programming and Computer Science
Replies
5
Views
2K
  • Programming and Computer Science
Replies
2
Views
4K
  • Programming and Computer Science
Replies
8
Views
6K
  • Programming and Computer Science
Replies
11
Views
2K
  • Atomic and Condensed Matter
Replies
8
Views
4K
  • Programming and Computer Science
Replies
2
Views
1K
Back
Top