fortran : invalid memory reference
I encountered this error while solving the diffusion equation on fortran using alternating directional implicit(ADI) method.
Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
#0 0xffffffff
#1 0xffffffff
#2 0xffffffff
#3 0xffffffff
#4 0xffffffff
#5 0xffffffff
#6 0xffffffff
#7 0xffffffff
#8 0xffffffff
#9 0xffffffff
#10 0xffffffff
#11 0xffffffff
#12 0xffffffff
#13 0xffffffff
#14 0xffffffff
And here is my code
program HW1_1_ADI
implicit none
real :: delx, dely, delt, dx, dy
real, parameter :: a=0.7, m=0.1, eta=0.001, pi=3.141592
real, dimension(401,101,50) :: psi
real, dimension(201,101,50) :: psi2, psi_half
real , dimension(401) :: x
real, dimension(400) :: a_D, b_D, C_D, v_D, o_D
real, dimension(100) :: a_D2, b_D2, C_D2, v_D2, o_D2
integer :: i, j, n, imax, jmax, nmax, ihalf
open(1,file='HW1.txt')
2 format (101((400(e10.3,','), e10.3), /))
imax = 401
jmax = 101
nmax = 50
ihalf = (imax+1)/2
delx = 2.0/(imax-1)
dely = 1.0/(jmax-1)
delt = 0.2 * 10**3
dx = eta * delt / delx**2
dy = eta * delt / dely**2
! x-coordinate
do i = 1, imax
x(i) = -1.0 + (i-1)*delx
end do
! initial condition
do i = 1, imax
do j = 1, jmax
psi(i,j,1) = cos(m*pi*x(i))-a*cos(3.0*m*pi*x(i))
end do
end do
do i = 1, ihalf
do j = 1, jmax
psi2(i,j,1) = cos(m*pi*x(i))-a*cos(3.0*m*pi*x(i))
end do
end do
! ADI
do n = 1, nmax-1
do j=2, jmax
do i = 2, ihalf
b_D(i-1) = 1.0 + dx
a_D(i-1) = -0.5*dx
c_D(i-1) = -0.5*dx
if (i == ihalf) then
a_D(i-1) = -dx
end if
if (j==jmax) then
if (i==2) then
v_D(i-1) = dy*psi2(i,jmax-1,n) + (1.0-dy)*psi2(i,jmax,n) + 0.5*dx*psi2(1,jmax,n)
else
v_D(i-1) = dy*psi2(i,jmax-1,n) + (1.0-dy)*psi2(i,jmax,n)
end if
else
if (i==2) then
v_D(i-1) = 0.5*dy*psi2(2,j-1,n) + (1.0-dy)*psi2(2,j,n) + 0.5*dy*psi2(2,j+1,n)&
+0.5*dx*psi2(1,j,n)
else
v_D(i-1) = 0.5*dy*psi2(i,j-1,n) + (1.0-dy)*psi2(i,j,n) + 0.5*dy*psi2(i,j+1,n)
end if
end if
end do
call tridag(a_D,b_D,c_D,v_D,o_D,ihalf-1)
do i = 1, ihalf-1
psi_half(i+1,j,n) = o_D(i)
end do
end do
! boundary condition
! y=0
do i = 1, ihalf
psi_half(i,1,:) = cos(m*pi*x(i))-a*cos(3.0*m*pi*x(i))
end do
! x=-1
do j = 1, jmax
psi_half(1,j,:) = cos(m*pi) - a*cos(3.0*m*pi)
end do
do i = 2, ihalf
do j = 2, jmax
b_D2(j-1) = 1.0 + dy
a_D2(j-1) = -0.5*dy
c_D2(j-1) = -0.5*dy
if (j == jmax) then
a_D2(j-1) = -dy
end if
if (i==ihalf) then
if (j==2) then
v_D2(j-1) = dx*psi_half(ihalf-1,j,n) + (1.0-dx)*psi_half(ihalf,j,n) + 0.5*dy*psi_half(ihalf,1,n)
else
v_D2(j-1) = dx*psi_half(ihalf-1,j,n) + (1.0-dx)*psi_half(ihalf,j,n)
end if
else
if (j==2) then
v_D2(j-1) = 0.5*dx*psi_half(i-1,j,n) + (1.0-dx)*psi_half(i,j,n) + 0.5*dx*psi_half(i+1,j,n)&
+0.5*dy*psi_half(i,1,n)
else
v_D2(j-1) = 0.5*dx*psi_half(i-1,j,n) + (1.0-dx)*psi_half(i,j,n) + 0.5*dx*psi_half(i+1,j,n)
end if
end if
end do
call tridag(a_D2,b_D2,c_D2,v_D2,o_D2,jmax-1)
do j = 1, jmax-1
psi2(i,j+1,n+1) = o_D2(j)
end do
end do
! boundary condition
! y=0
do i = 1, ihalf
psi2(i,1,n+1) = cos(m*pi*x(i))-a*cos(3.0*m*pi*x(i))
end do
! x=-1
do j = 1, jmax
psi2(1,j,n+1) = cos(m*pi) - a*cos(3.0*m*pi)
end do
! reflection condition
do i = 1, ihalf
do j = 1, jmax
psi(i,j,n+1) = psi2(i,j,n+1)
end do
end do
end do
do i = ihalf+1, imax
do j = 1, jmax
psi(i,j,:) = psi(imax+1-i,j,:)
end do
end do
write(1,2) psi(:,:,2)
contains
subroutine tridag(a,b,c,r,u,n)
implicit none
integer n, nMAX
real a(n), b(n), c(n), r(n), u(n)
parameter (nMAX = 100)
integer j
real bet, gam(nMAX)
if (b(1) == 0.0) print *, 'tridag: rewrite equations'
bet = b(1)
u(1) = r(1)/bet
do j = 2, n
gam(j) = c(j-1)/bet
bet = b(j) - a(j)*gam(j)
if (bet == 0.0) print *, 'tridag failed'
u(j) = (r(j)-a(j)*u(j-1)) / bet
end do
do j = n-1, 1, -1
u(j) = u(j) - gam(j+1)*u(j+1)
end do
return
end subroutine
end program
What I'm curious about is that it's good at first.
At first, I set it up imax=200(maximum index of x-coordinate).
Here is original code.
program HW1_1_ADI
implicit none
real :: delx, dely, delt, dx, dy
real, parameter :: a=0.7, m=0.1, eta=0.001, pi=3.141592
real, dimension(201,101,50) :: psi
real, dimension(101,101,50) :: psi2, psi_half
real , dimension(201) :: x
real, dimension(100) :: a_D, b_D, C_D, v_D, o_D
real, dimension(100) :: a_D2, b_D2, C_D2, v_D2, o_D2
integer :: i, j, n, imax, jmax, nmax, ihalf
open(1,file='HW1.txt')
2 format (101((200(e10.3,','), e10.3), /))
imax = 201
jmax = 101
nmax = 50
ihalf = (imax+1)/2
delx = 2.0/(imax-1)
dely = 1.0/(jmax-1)
delt = 0.2 * 10**3
dx = eta * delt / delx**2
dy = eta * delt / dely**2
! x-coordinate
do i = 1, imax
x(i) = -1.0 + (i-1)*delx
end do
! initial condition
do i = 1, imax
do j = 1, jmax
psi(i,j,1) = cos(m*pi*x(i))-a*cos(3.0*m*pi*x(i))
end do
end do
do i = 1, ihalf
do j = 1, jmax
psi2(i,j,1) = cos(m*pi*x(i))-a*cos(3.0*m*pi*x(i))
end do
end do
! ADI
do n = 1, nmax-1
do j=2, jmax
do i = 2, ihalf
b_D(i-1) = 1.0 + dx
a_D(i-1) = -0.5*dx
c_D(i-1) = -0.5*dx
if (i == ihalf) then
a_D(i-1) = -dx
end if
if (j==jmax) then
if (i==2) then
v_D(i-1) = dy*psi2(i,jmax-1,n) + (1.0-dy)*psi2(i,jmax,n) + 0.5*dx*psi2(1,jmax,n)
else
v_D(i-1) = dy*psi2(i,jmax-1,n) + (1.0-dy)*psi2(i,jmax,n)
end if
else
if (i==2) then
v_D(i-1) = 0.5*dy*psi2(2,j-1,n) + (1.0-dy)*psi2(2,j,n) + 0.5*dy*psi2(2,j+1,n)&
+0.5*dx*psi2(1,j,n)
else
v_D(i-1) = 0.5*dy*psi2(i,j-1,n) + (1.0-dy)*psi2(i,j,n) + 0.5*dy*psi2(i,j+1,n)
end if
end if
end do
call tridag(a_D,b_D,c_D,v_D,o_D,ihalf-1)
do i = 1, ihalf-1
psi_half(i+1,j,n) = o_D(i)
end do
end do
! boundary condition
! y=0
do i = 1, ihalf
psi_half(i,1,:) = cos(m*pi*x(i))-a*cos(3.0*m*pi*x(i))
end do
! x=-1
do j = 1, jmax
psi_half(1,j,:) = cos(m*pi) - a*cos(3.0*m*pi)
end do
do i = 2, ihalf
do j = 2, jmax
b_D2(j-1) = 1.0 + dy
a_D2(j-1) = -0.5*dy
c_D2(j-1) = -0.5*dy
if (j == jmax) then
a_D2(j-1) = -dy
end if
if (i==ihalf) then
if (j==2) then
v_D2(j-1) = dx*psi_half(ihalf-1,j,n) + (1.0-dx)*psi_half(ihalf,j,n) + 0.5*dy*psi_half(ihalf,1,n)
else
v_D2(j-1) = dx*psi_half(ihalf-1,j,n) + (1.0-dx)*psi_half(ihalf,j,n)
end if
else
if (j==2) then
v_D2(j-1) = 0.5*dx*psi_half(i-1,j,n) + (1.0-dx)*psi_half(i,j,n) + 0.5*dx*psi_half(i+1,j,n)&
+0.5*dy*psi_half(i,1,n)
else
v_D2(j-1) = 0.5*dx*psi_half(i-1,j,n) + (1.0-dx)*psi_half(i,j,n) + 0.5*dx*psi_half(i+1,j,n)
end if
end if
end do
call tridag(a_D2,b_D2,c_D2,v_D2,o_D2,jmax-1)
do j = 1, jmax-1
psi2(i,j+1,n+1) = o_D2(j)
end do
end do
! boundary condition
! y=0
do i = 1, ihalf
psi2(i,1,n+1) = cos(m*pi*x(i))-a*cos(3.0*m*pi*x(i))
end do
! x=-1
do j = 1, jmax
psi2(1,j,n+1) = cos(m*pi) - a*cos(3.0*m*pi)
end do
! reflection condition
do i = 1, ihalf
do j = 1, jmax
psi(i,j,n+1) = psi2(i,j,n+1)
end do
end do
end do
do i = ihalf+1, imax
do j = 1, jmax
psi(i,j,:) = psi(imax+1-i,j,:)
end do
end do
write(1,2) psi(:,:,2)
contains
subroutine tridag(a,b,c,r,u,n)
implicit none
integer n, nMAX
real a(n), b(n), c(n), r(n), u(n)
parameter (nMAX = 100)
integer j
real bet, gam(nMAX)
if (b(1) == 0.0) print *, 'tridag: rewrite equations'
bet = b(1)
u(1) = r(1)/bet
do j = 2, n
gam(j) = c(j-1)/bet
bet = b(j) - a(j)*gam(j)
if (bet == 0.0) print *, 'tridag failed'
u(j) = (r(j)-a(j)*u(j-1)) / bet
end do
do j = n-1, 1, -1
u(j) = u(j) - gam(j+1)*u(j+1)
end do
return
end subroutine
end program
This works well.
What's wrong in my code?
asked on Stack Overflow Apr 5, 2021 by 
이종훈
이종훈1 Answer
I suggest you learn how to use your compiler to help you diagnose these problems - when developing always turn on run time error checks. Here is what gfortran can tell you, note the -fcheck=all and -g flags, though I would recommend all the ones I use:
ijb@ijb-Latitude-5410:~/work/stack$ gfortran --version
GNU Fortran (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
ijb@ijb-Latitude-5410:~/work/stack$ gfortran -Wall -Wextra -fcheck=all -O -g -std=f2008 imax.f90
imax.f90:160:12:
160 | if (b(1) == 0.0) print *, 'tridag: rewrite equations'
| 1
Warning: Equality comparison for REAL(4) at (1) [-Wcompare-reals]
imax.f90:168:16:
168 | if (bet == 0.0) print *, 'tridag failed'
| 1
Warning: Equality comparison for REAL(4) at (1) [-Wcompare-reals]
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
At line 166 of file imax.f90
Fortran runtime error: Index '101' of dimension 1 of array 'gam' above upper bound of 100
Error termination. Backtrace:
#0 0x7f998e9abd01 in ???
#1 0x7f998e9ac849 in ???
#2 0x7f998e9acec6 in ???
#3 0x5562d1681337 in tridag
at /home/ijb/work/stack/imax.f90:166
#4 0x5562d1681971 in hw1_1_adi
at /home/ijb/work/stack/imax.f90:72
#5 0x5562d1681f56 in main
at /home/ijb/work/stack/imax.f90:149
From this it looks like you haven't adjusted the nMax parameter in tridag correctly - I suggest you also learn about allocatable arrays and how you can avoid this problem in a way which means you don't have to change your code for every different size you try.
answered on Stack Overflow Apr 5, 2021 by 
Ian Bush
Ian BushUser contributions licensed under CC BY-SA 3.0