Fortran [fortran90] having problem to call back the value from previous loop

AI Thread Summary
The discussion revolves around a Fortran 90 code implementation for a 2D Total Field/Scattered Field (TFSF) simulation in the context of electromagnetic wave propagation. The user reports an issue where the field evolution does not respond correctly to the TFSF boundary conditions, resulting in an incident wave value of zero. A key point raised is the identification of a missing operator in the field evolution loop, specifically in the equation for updating the electric field (Ez). Suggestions include improving code readability by adding spaces for clarity and using more descriptive variable names to facilitate debugging and collaboration. The user indicates that they have addressed the missing part, suggesting progress in resolving the issue.
s_hy
Messages
57
Reaction score
0
hi all,

I wrote programming code for fdtd 2d tfsf in fortran 90 as below:

Code: 
!simulation for tmz mode incident plane in 2d

subroutine tfsf2d
implicit none

double precision             :: f0,miu,eps,delta,S,c,lambda,E0,N0  
double precision             :: deltat
integer                            :: minit,pinit
integer                            :: iinit,ilast
integer                            :: jinit,jlast
double precision,parameter         :: pi = 3.14159265
double precision             :: phi,cosphi,sinphi
integer                            :: h,m
double precision,dimension(300)    :: Hm,Em
integer                            :: y
double precision             :: f,G,A,B
double precision,dimension(100)    :: k 
double precision,dimension(100)    :: Vp
double precision              :: Ca,Da,Cb,Db
integer                            :: i,j
integer                            :: n,p
double precision             :: d,dp,dpp
double precision             :: Einc,Hinc
double precision,dimension(300,300):: Ezinc,Hxinc,Hyinc
character(len=20)          :: filename

double precision             :: miu0,eps0
double precision             :: deltax, deltay
double precision             :: sigma,sigmax,sigmay,sigmat,sigmatx,sigmaty
integer                            :: iinitp,ilastp
integer                            :: jinitp,jlastp
double precision,dimension(300,300) :: Cax,Cay,Dax,Day,Cbx,Cby,Dbx,Dby
double precision,dimension(300,300) :: Ez
double precision,dimension(300,300) :: Hy,Hx

!## parameters ##
 f0 = 100.0
 miu = 1.0
 eps = 1.0
 delta = 1.e-3
 S = 1.0001
 c = 1.0
 lambda = c/f0
 E0 = 1.0D0
 N0 = lambda/delta
 deltat = (delta/c)/sqrt(2.0) 

 minit = 1
 pinit = 1
 
 phi = pi/4.0
 cosphi = cos(phi)
 sinphi = sin(phi)
 !print *, 'cosphi=',cosphi,'sinphi=',sinphi

 miu0 = 1.0
 eps0 = 1.0
 sigma = 1.0
 !print *, 'lambda=',lambda,'deltat=',deltat

 deltax = delta
 deltay = delta

 iinit = 10
 ilast = 100
 jinit = 10
 jlast = 100
 iinitp = 1
 ilastp = 110
 jinitp = 1
 jlastp = 110


 sigmax = sigma
 sigmay = sigma
 sigmat = sigma
 sigmatx = sigma
 sigmaty = sigma

!initialization for 1d grid
 do h = 1, 300
    Hm(h) = 0
    Em(h) = 0
 end do

!initialization for incident wave
do p = 1,300
 do m = 1,300
      Hxinc(m-1,p+2) = 0
      Hyinc(m,p+1) = 0

      Ezinc(m-1,p+1) = 0
 end do
end do

!initialization for field evolution 
 do p = 1,300
  do m = 1,300
     Ez(m-1,p+1)=0

     Hx(m-1,p+2) = 0
     Hy(m,p+1) = 0
  end do
 end do 
 
!properties coefficient for 2d wave
 do y = 1,45
  f = 0.5
  k(1) = 2*pi
  G = (1/f**2)*sin(f*pi/N0)**2             
  A = (delta*cos(y*pi/180))/2
  B = (delta*sin(y*pi/180))/2
  k(y+1) = k(y)-(sin(A*k(y))**2+sin(B*k(y))**2)-G/(A*sin(2*A*k(y))+B*sin(2*B*k(y)))
  Vp(y) = 2*pi/k(y+1)

  Ca = 1.0
  Cb = deltat/((Vp(1)/Vp(y))*eps*delta)
 
  Da = 1.0
  Db = deltat/((Vp(1)/Vp(y))*miu*delta)
  !Print *,'Ca=',Ca,'Cb=',Cb,'Da=',Da,'Db=',Db
 end do


!property coefficient 
 do m = 1,300
  do p = 1,300
   Cax(m,p) = exp(-sigma*deltat/eps)
   Cbx(m,p) = (1-Cax(m,p))/(sigma*deltax)
   Cay(m,p) = exp(-sigma*deltat/eps)
   Cby(m,p) = (1-Cay(m,p))/(sigma*deltay)
   Dax(m,p) = exp(-sigmatx*deltat/miu)
   Dbx(m,p) = (1-Dax(m,p))/(sigmatx*deltax)
   Day(m,p) = exp(-sigmaty*deltat/miu)
   Dby(m,p) = (1-Day(m,p))/(sigmaty*deltay)
  end do
 end do

!time loop
do n = 1,200
  
  write (filename, "('data',I3.3,'.dat')") n
  open (unit=31,file=filename)

!1d background wave
 do h = 1, 300
  if (Mod(h,2) /= 0) then  !j is odd; 
   Em(h) = Ca*Em(h) + Cb*(Hm(h+1)-Hm(h-1))
   !print *,'h=',h,'Em(h)=',Em(h)
  end if
 end do !h

 do h = 1, 300
  if (Mod(h,2) == 0) then  !j is even;
   Hm(h) = Da*Hm(h) + Db*(Em(h+1)-Em(h-1)) 
   !print *,'h=',h,'Hm(h)=',Hm(h)
  end if 
 end do !h

!plane wave source
  Em(1) = E0*sin (2*pi*f0*n*deltat)

!incident field 
  do m = 2,2*ilast+2           !odd
   do p = 2,2*jlast+2          !odd
    if (Mod(m,2)/=0) then
     if (Mod(p,2)/=0) then
      d = cosphi*(m - minit) + sinphi*(p - pinit)
      dp = d - int(d)     
      Einc = (1-dp)*Em(int(d))+dp*Em(int(d+1))
      Ezinc(m,p+2) = Einc
     end if
    end if
   end do
  end do

  do m = 2,2*ilast+2
   do p = 2,2*jlast+2
    if (Mod(m,2) == 0) then     !even
     if (Mod(p,2) == 0) then    !even
      d = cosphi*(m - minit) + sinphi*(p - pinit)
      dpp = int(d)+0.5
      dp = dpp-int(dpp)           
      Hinc = (1-dp)*Hm(int(dpp-0.5))+dp*Hm(int(dpp+0.5))
      Hxinc(m,p+1) = sinphi*Hinc
      Hyinc(m-1,p+2) = -cosphi*Hinc
      !print *,'m=',m,'p=',p,'Hyinc(m-1,p+2)=',Hyinc(m-1,p+2)
      !end if
     end if
    end if
   end do
 end do

!##TFSF boundary## 
 do p = jinit-1,jlast-1
  if (Mod(p,2) == 0) then 
   Hy(10,p)=Hy(10,p)+Dby(10,p)*Ezinc(11,p+1)
   !print *,'p=',p,'Hy(10,p)=',Hy(10,p)
  end if
 end do

 do p = jinit-1,jlast-1
  if (Mod(p,2) == 0) then
   Hy(200,p)=Hy(200,p)+Dby(200,p)*Ezinc(199,p+1)
   !print *,'p=',p,'Hy(200),=',Hy(200,p)
  end if
 end do

 do m = iinit-1,ilast-1 !iinit+1,2*ilast-1
  if (Mod(m,2) == 0) then
   Hx(m,10)=Hx(m,10)+Dbx(m,10)*Ezinc(m+1,11)
    !print *,'Hx(m,8)=',Hx(m,8)
  end if
 end do

 do m = iinit-1,ilast-1
  if (Mod(m,2) == 0) then
   Hx(m,200)=Hx(m,200)-Dbx(m,200)*Ezinc(m+1,199)
   !print *,'Hx(m,200)=',Hx(m,200)
  end if
 end do
  
 do p = jinit-1,jlast-1
  if (Mod(p,2)/=0) then
   Ez(11,p+1) = Ez(11,p+1)-Cby(11,p)*Hyinc(10,p+1)
   !print *,'p=',p,'Ez(11,p+1)=',Ez(11,p+1)
  end if
 end do

 do p = jinit-1,jlast-1 
  if (Mod(p,2)/=0) then
   Ez(199,p+1) = Ez(199,p+1)+Cby(199,p)*Hyinc(201,p+1)
   !print *,'p=',p,'Ez(199,p+1)=',Ez(199,p+1)
  end if
 end do

 do m = iinit-1,ilast-1
  if (Mod(m,2)/=0) then
   Ez(m+1,11) = Ez(m+1,11)+Cbx(m,11)*Hxinc(m+1,9)
   !print *,'m=',m,'Ez(m+1,11)=',Ez(m+1,11)
  end if
 end do

 do m = iinit-1,ilast-1
  if (Mod(m,2)/=0) then
   Ez(m+1,199) = Ez(m+1,199)-Cby(m,199)*Hxinc(m+1,199)
   !print *,'m=',m,'Ez(m+1,199)=',Ez(m+1,199)
  end if
 end do
!## end tfsf boundary##


!## field evolution ## 
 do m = iinitp,ilastp
  do p = jinitp,jlastp
   if (Mod(m,2)/=0) then
    if (Mod(p,2)/=0) then
     Ez(m,p+2) = Cax(m,p)*Ez(m,p+2)+Cbx(m,p)*(Hy(m+1,p+2)-Hy(m  1,p+2)+Hx(m,p+1)-Hx(m,p+3))
     !print *,'n=',n,'m=',m,'p=',p,'Ez(m,p+2)=',Ez(m,p+2)
     !write (31,*) m,p,Ez(m,p+2)
     !if (p == 219) write (31,*) ' ' 
    end if
   end if
  end do
 end do

 do m = iinitp,ilastp
  do p = iinitp,jlastp
   if (Mod(m,2) == 0) then
    if (Mod(p,2) == 0) then
     Hx(m-1,p+2) = Dax(m,p)*Hx(m-1,p+2)+Dbx(m,p)*(Ez(m-2,p+1)-Ez(m-1,p+3))
     Hy(m,p+1) = Day(m,p)*Hy(m,p+1)+Dby(m,p)*(Ez(m+1,p+1)-Ez(m-1,p+1))
     !print *, 'n=',n,'m=',m,'p=',p,'Hy(m,p+1)=',Hy(m,p+1)
    end if
   end if
  end do
 end do

close (unit=31)
end do !n

return
end subroutine

i have a problem to call the tfsf boundary condition in field evolution. it's seems that the field evolution didn't react with the boundary (the incident wave = 0)...anyone have any idea which part did i make mistake or other idea how to approach this problem.

thank you
 
Technology news on Phys.org
You have a missing operator in Hy(m ?? 1, p+2) in your Field Evolution Loop:

<>
Ez(m,p+2) = Cax(m,p)*Ez(m,p+2)+Cbx(m,p)*(Hy(m+1,p+2)-Hy(m 1,p+2)+Hx(m,p+1)-Hx(m,p+3))
 
I tried to use <> to indicate the position of the error but posting pushed it back to the left margin.
 
I'm stumped as to why your Fortran compiler didn't catch this error.
 
This should be clearer as to where the missing operator is.
SteamKing said:
You have a missing operator in Hy(m ?? 1, p+2) in your Field Evolution Loop:

Code: 
Ez(m,p+2) = Cax(m,p)*Ez(m,p+2)+Cbx(m,p)*(Hy(m+1,p+2)-Hy(m <missing operator>[/color]  1,p+2)
   +Hx(m,p+1)-Hx(m,p+3))

s_hy, your code would be easier to read by humans if you added many more spaces. In the line above, the only spaces are around the '='. When everything is jammed together, it's harder for us humans to find things.

Ideally, the line above would look something like this:
Code: 
Ez(m, p+2) = Cax(m, p) * Ez(m, p+2) + Cbx(m, p) * (Hy(m+1, p+2) - Hy(m <?>[/color]1, p+2)
   + Hx(m, p+1)-Hx(m, p+3))

Also, and this is some very common amongst Fortran programmers, there are a lot of names used that seem arcane to the point of gobbledygook. For example:
tmz
tfsf2d
TFSF
Einc
Hinc
Ezinc
Hxinc
Hyinc
Cax
Cay
Dax
Day
Cbx
Cby
Dbx
Dby
If you're writing a program that only you will debug and maintain, I suppose that's some justification for doing something quick and dirty. But when your code doesn't work, and you need to get help from others, it makes it much more difficult for others than it needs to be if your variable names have little or no connection with what they represent.
 
Last edited:
Code: 
Ez(m,p+2) = Cax(m,p)*Ez(m,p+2)+Cbx(m,p)*(Hy(m+1,p+2)-Hy(m - 1,p+2)
   +Hx(m,p+1)-Hx(m,p+3))

i have fix the missing part...
 

Thread 'Learning Assembly and computer architecture for x86'

Dear Peeps I have posted a few questions about programing on this sectio of the PF forum. I want to ask you veterans how you folks learn program in assembly and about computer architecture for the x86 family. In addition to finish learning C, I am also reading the book From bits to Gates to C and Beyond. In the book, it uses the mini LC3 assembly language. I also have books on assembly programming and computer architecture. The few famous ones i have are Computer Organization and...
View full post »

Thread 'Learning data structures and algorithms in different programming languages'

I have a quick questions. I am going through a book on C programming on my own. Afterwards, I plan to go through something call data structures and algorithms on my own also in C. I also need to learn C++, Matlab and for personal interest Haskell. For the two topic of data structures and algorithms, I understand there are standard ones across all programming languages. After learning it through C, what would be the biggest issue when trying to implement the same data...
View full post »

Similar threads

Fortran [Fortran90] how to write several output into one file
Replies
3
Views
3K
Fortran [fortran90] have problem in real array index
Replies
3
Views
12K
Fortran [Fortran90] how to write increment in dimension
Replies
4
Views
4K
Fortran [Fortran90] fdtd in polar coordinates, got infinity output
Replies
3
Views
3K
Fortran Spring Pendulum system in FORTRAN
Replies
8
Views
2K
Fortran Having problem to understand kind notation in fortran 90
Replies
8
Views
6K
Fortran Getting my method of bisection to output the correct value
Replies
11
Views
2K
[FDTD/Fotran] detected reflection near boundary but don't know why
Replies
3
Views
2K
Fortran [Fortran90]Problem with Gaussian distribution
Replies
8
Views
4K
Fortran Simulating Contact Force in a Moving Dashpot System Using FORTRAN
Replies
26
Views
3K

Hot Threads

Recent Insights

Back
Top