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

  • Context: Fortran 
  • Thread starter Thread starter s_hy
  • Start date Start date
  • Tags Tags
    Fortran90 Loop Value
Click For Summary

Discussion Overview

The discussion revolves around a Fortran 90 code implementation for a 2D Finite-Difference Time-Domain (FDTD) simulation, specifically focusing on the Total Field/Scattered Field (TFSF) boundary condition. Participants are examining issues related to the interaction of the field evolution with the boundary conditions in the simulation.

Discussion Character

  • Technical explanation
  • Debate/contested

Main Points Raised

  • The original poster describes a problem where the field evolution does not seem to interact with the TFSF boundary condition, resulting in the incident wave being zero.
  • One participant points out a potential error in the code, suggesting that there is a missing operator in the expression for Hy in the field evolution loop.
  • Another participant expresses confusion about why the Fortran compiler did not catch the mentioned error.
  • A later reply attempts to clarify the location of the missing operator in the code snippet provided by the original poster.

Areas of Agreement / Disagreement

Participants have identified a specific issue in the code, but there is no consensus on whether this is the sole reason for the problem described. The discussion remains unresolved regarding the overall impact of the identified error on the simulation's behavior.

Contextual Notes

There are indications of potential misunderstandings regarding the compiler's error detection capabilities and the specific coding practices in Fortran that may lead to such issues.

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...
 

Similar threads

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