Fortran: 2D array taking value not assigned

Question

I am creating a 2D array phi and initializing it to 0. I am then going through a pair of do loops to modify some of the values of phi. However, some of the values are being assigned when they shouldn't be.

Specifically when I run the example below (compiling with either ifort 14 or 15) about 1/3rd of the bottom row are -20, when I compile with gfortran about 1/3rd of the bottom row are -20, and another 1/3rd are 10.

When I uncomment the phi(i,j) = 0.0 before the if-elseif statement, the value of phi(i,j) stays 0.0. So it's not that the if-elseif is changing its value.

Does anyone know what is causing this?

As a side issue when I uncomment close (42), nothing I place after it will execute. So I'm not sure if that's related.

Thank you,

John

program test
implicit none
real, parameter, dimension(2) :: delta = (/ 1E-7, 1E-7 /)
real,parameter,dimension(2) ::  Length = (/ 6E-6, 5E-6 /)
integer, PARAMETER, dimension(2) :: cells = Length / delta
real, dimension(cells(1), cells(2)) :: phi

REAL,parameter :: V0 = -20
REAL,parameter :: VL = 10
integer :: i,j

phi = 0.0

open(unit=42,file='ChangeValuePhi.csv',form='formatted')
do i = 0, size(phi,1)
    do j = 0, size(phi,2)
            if (j .ne. 0) then
                write(42,'(A)',Advance="NO") ","
            endif

!           phi(i,j) = 0.0
            if (j <= nint(lowerBoundary(i * delta(1)) / delta(2))) then     !If at minimum z feature
                phi(i,j) = V0
            elseif (j >= nint(upperBoundary(i * delta(1)) / delta(2))) then     !If at maximum z feature
                phi(i,j) = VL
            end if

            write(42,'(F6.2)',Advance="NO") phi(i,j)
    enddo
    write(42,*) 
enddo
!close(42)

contains

function lowerBoundary(r) result(z)
    real :: r,z
    real, parameter :: RL = 3E-6
    real, parameter :: HL = 3E-6

    if (r < RL) then
        z = HL - HL * (r / RL)
    elseif (r >= RL) then
        z = 0.0
    endif
end function lowerBoundary

function upperBoundary(r) result(z)
    real :: r,z
    real, parameter :: RU = 0.0
    real, parameter :: HU = 0.0

    if (r <= RU) then
        z = Length(2) - (HU - HU * (r / RU))
    elseif (r > RU) then
        z = Length(2) - 0.0
    end if
end function upperBoundary


end program test

are you not getting any runtime error with division by zero in upperBoundary? you are dividing by RU that is zero. — innoSPG
– innoSPG, Commented Aug 23, 2015 at 22:55
No because r>= 0, but you're right that I should change where the "or equal" is to avoid a division by zero error. — user1543042
– user1543042, Commented Aug 23, 2015 at 23:09

casey · Accepted Answer · 2015-08-23 23:14:14Z

1

You are accessing phi out of bounds. In your loop,

do i = 0, size(phi,1)
    do j = 0, size(phi,2)
...
!           phi(i,j) = 0.0
            if (j <= nint(lowerBoundary(i * delta(1)) / delta(2))) then     !If at minimum z feature
                phi(i,j) = V0
            elseif (j >= nint(upperBoundary(i * delta(1)) / delta(2))) then     !If at maximum z feature
                phi(i,j) = VL
            end if
...
enddo

you are accessing phi with either i or j equal to 0. Fortran arrays are 1-indexed unless you specifically declare the lower bound to be something else (e.g. 0). The incorrect output you are seeing is due to the out of bounds access.

You can change the beginning of your do loop to

  do i = 1, size(phi,1)
     do j = 1, size(phi,2)
        if (j .ne. 1) then
           write(42,'(A)',Advance="NO") ","
        endif

to avoid the out of bounds access.

answered Aug 23, 2015 at 23:14

casey

6,9251 gold badge26 silver badges38 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

user1543042 Over a year ago

Okay, I've apparently been staring at my computer for to many hours. Thank you.

casey Over a year ago

@user1543042 I know that feeling all too well.

Collectives™ on Stack Overflow

Fortran: 2D array taking value not assigned

1 Answer 1

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related