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

  • Fortran
  • Thread starter s_hy
  • Start date
  • #1
61
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
 

Answers and Replies

  • #2
SteamKing
Staff Emeritus
Science Advisor
Homework Helper
12,796
1,670
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))
 
  • #3
SteamKing
Staff Emeritus
Science Advisor
Homework Helper
12,796
1,670
I tried to use <> to indicate the position of the error but posting pushed it back to the left margin.
 
  • #4
SteamKing
Staff Emeritus
Science Advisor
Homework Helper
12,796
1,670
I'm stumped as to why your Fortran compiler didn't catch this error.
 
  • #5
35,028
6,774
This should be clearer as to where the missing operator is.
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 [color="red"]<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="red"]<???>[/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:
  • #6
61
0
Code:
Ez(m,p+2) = Cax(m,p)*Ez(m,p+2)+Cbx(m,p)*(Hy(m+1,p+2)-Hy([COLOR="Red"]m - 1[/COLOR],p+2)
   +Hx(m,p+1)-Hx(m,p+3))

i have fix the missing part...
 

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

Replies
3
Views
10K
Replies
2
Views
2K
Replies
1
Views
2K
  • Last Post
Replies
0
Views
2K
Replies
8
Views
2K
  • Last Post
Replies
8
Views
3K
Replies
5
Views
3K
Replies
4
Views
701
  • Last Post
Replies
5
Views
3K
Top