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

In summary, the conversation is about writing programming code in Fortran 90 for a simulation of a TMz mode incident plane in 2D. The code includes parameters, initializations, and time loops for the simulation, as well as coefficients and properties for the 2D wave. The purpose of the code is to model electromagnetic fields.
  • #1
s_hy
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
 
Technology news on Phys.org
  • #2
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
I tried to use <> to indicate the position of the error but posting pushed it back to the left margin.
 
  • #4
I'm stumped as to why your Fortran compiler didn't catch this error.
 
  • #5
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 [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
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...
 

1. What is Fortran 90 and how is it used in scientific research?

Fortran 90 is a high-level programming language commonly used in scientific and engineering applications. It is known for its efficient handling of mathematical and scientific operations, making it a popular choice for numerical simulations and data analysis.

2. What does it mean to "call back a value" in a Fortran 90 loop?

Calling back a value in a Fortran 90 loop refers to using a variable or data from a previous iteration of the loop in the current iteration. This is often done in order to update or use the value in some way in the current iteration.

3. Why might someone encounter problems when calling back a value from a previous loop in Fortran 90?

There are a few potential reasons for encountering problems when calling back a value in Fortran 90. It could be due to incorrect usage of variables or arrays, improper indexing, or a logical error in the code.

4. How can I troubleshoot issues with calling back values in a Fortran 90 loop?

If you are experiencing problems with calling back values in a Fortran 90 loop, it is important to carefully review your code and check for any errors or inconsistencies. You may also want to consult online resources or seek assistance from a more experienced Fortran programmer.

5. Are there any best practices for calling back values in Fortran 90 loops?

Yes, there are some best practices to keep in mind when calling back values in Fortran 90 loops. It is important to properly initialize variables and arrays, use appropriate data types, and avoid overwriting values unintentionally. It is also helpful to use comments and clear variable names to make the code more readable and easier to troubleshoot.

Similar threads

Fortran [Fortran90] how to write several output into one file
    • Programming and Computer Science
Replies
4
Views
3K
Fortran [fortran90] have problem in real array index
    • Programming and Computer Science
Replies
4
Views
12K
Fortran [Fortran90] how to write increment in dimension
    • Programming and Computer Science
Replies
4
Views
3K
Fortran [Fortran90] fdtd in polar coordinates, got infinity output
    • Programming and Computer Science
Replies
4
Views
3K
Fortran Spring Pendulum system in FORTRAN
    • Programming and Computer Science
Replies
8
Views
1K
Fortran Getting my method of bisection to output the correct value
    • Programming and Computer Science
Replies
11
Views
1K
[FDTD/Fotran] detected reflection near boundary but don't know why
    • Programming and Computer Science
Replies
4
Views
2K
How to put mur's absorbing boundary condition in 1d fdtd maxwell equation
    • Programming and Computer Science
Replies
1
Views
4K
Fortran [Fortran90]Problem with Gaussian distribution
    • Programming and Computer Science
Replies
8
Views
4K
Fortran Having problem to understand kind notation in fortran 90
    • Programming and Computer Science
Replies
8
Views
6K

Hot Threads

Recent Insights

Back
Top