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

Large/Infinite results

  1. Jun 18, 2015 #1
    Hello to all!
    I have written this 1D FDTD code but the results I take are very large and after some timestep they get to infinite/NaN... I really don't know what's wrong. The expressions I have used are the exact given solution for the Maxwell equations.

    Thank you for your help.
    Code (Fortran):
    program yee
        implicit none


        real :: Delta, tau,omega,n,mu,time,Eprev,Hprev
        real, dimension(:), allocatable :: sigma,eps,A,B,C,D,J
        real, dimension(:), allocatable :: E,H
        real, parameter :: pi=3.1415926
        integer :: i,m=10000, Length=5000,ls,k


        Allocate(E(Length),H(Length),sigma(Length),eps(Length))
        Allocate(A(Length),B(Length),C(Length),D(Length))
        Allocate(J(Length))

    !Initialization
        open(unit=1,file="Elektrfeld.txt")
        Delta=0.02d0
        tau=0.9d0*Delta
        omega=2.0d0*pi
        n=1.46d0
        mu=1
        ls=1000
        E(:)=0.0d0
        H(:)=0.0d0
        A(:)=0.0d0
        B(:)=0.0d0
        C(:)=0.0d0
        D(:)=0.0d0
        J(:)=0.0d0
        Eprev=0.0d0
        Hprev=0.0d0


        do i=1,Length
            if (i.le.6/Delta) then
                sigma(i)=1.0d0
            else if (i.gt.6/Delta .and. i.lt.((Length*Delta-6.0d0)/Delta))then
                sigma(i)=0.0d0
            else if (i.ge.(Length*Delta-6.0d0)/Delta .and. i.le.Length) then
                sigma(i)=1.0d0
            end if
        end do


        do i=1,Length
            if (i.lt.Length*0.5d0) then
                eps(i)=1.
            else if (i.ge.1100 .and. i.lt.(1100+2)) then
                eps(i)=n**2.
            else if (i.ge.(1100+2)) then
                eps(i)=1.
            end if
        end do


        do i=1,Length
            A(i)=(1.-sigma(i)*tau/(2.*mu))/(1.+sigma(i)*tau/(2.*mu))
            B(i)=(tau/mu)/(1.+(sigma(i)*tau)/(2.*mu))
            C(i)=(1.-sigma(i)*tau/(2.*eps(i)))/(1.+sigma(i)*tau/(2.*eps(i)))
            D(i)=(tau/eps(i))/(1.+sigma(i)*tau/(2.*eps(i)))
        end do


        write(*,*) "time"
        read(*,*) time

    !----------------Main calculation----------------
        do k=1,m        !time
            E(1)=0.0d0
            E(Length)=0.0d0
            J(ls)=sin(2.*pi*k*tau)*exp(-((k*tau-30.)/10.)**2.)
            do i=2,Length
                Eprev=E(i)
                E(i)=D(i)*(H(i)-H(i-1))/Delta+C(i)*E(i)-D(i)*J(i)
            end do
            do i=1,Length
                Hprev=H(i)
                H(i)=B(i)*(E(i)-E(i-1))/Delta+A(i)*Hprev
            end do
            if (k.eq.time ) then
                write(1,*) "-----",k
                do i=800,1200
                    write(1,*) E(i)
                end do
            end if
        end do

        close(1)
    end program yee
     
    Last edited: Jun 18, 2015
  2. jcsd
  3. Jun 18, 2015 #2

    phinds

    User Avatar
    Gold Member
    2016 Award

    we have "code" tags that get your code properly formatted. It would be a good idea if you modified your post to use them.
     
  4. Jun 18, 2015 #3
    Well, at quick glance, the second i-loop goes from 1 to Length and it contains the term E(i-1), this means that for i=1, you are calling on E(0) which is not defined...it is outside your program altogether given the way you allocated the arrays, their indexes start at 1.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Large/Infinite results
  1. Infinite Loop (Replies: 0)

  2. Infinite Databases? (Replies: 10)

Loading...