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

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

  1. Mar 20, 2013 #1
    hi all,

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

    Code (Text):


    !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
     
  2. jcsd
  3. Mar 21, 2013 #2

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    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))
     
  4. Mar 21, 2013 #3

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    I tried to use <> to indicate the position of the error but posting pushed it back to the left margin.
     
  5. Mar 21, 2013 #4

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    I'm stumped as to why your Fortran compiler didn't catch this error.
     
  6. Mar 21, 2013 #5

    Mark44

    Staff: Mentor

    This should be clearer as to where the missing operator is.
    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 (Text):
    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: Mar 21, 2013
  7. Mar 22, 2013 #6
    Code (Text):

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




Similar Discussions: [fortran90] having problem to call back the value from previous loop
Loading...