问题描述
我在使用交替方向隐式 (ADI) 方法求解 fortran 上的扩散方程时遇到了这个错误。
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
这是我的代码
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,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,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,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,n) + 0.5*dx*psi2(1,n)
else
v_D(i-1) = dy*psi2(i,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,n) + 0.5*dy*psi2(2,j+1,n)&
+0.5*dx*psi2(1,n)
else
v_D(i-1) = 0.5*dy*psi2(i,n) + 0.5*dy*psi2(i,n)
end if
end if
end do
call tridag(a_D,c_D,o_D,ihalf-1)
do i = 1,ihalf-1
psi_half(i+1,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,:) = 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,n) + (1.0-dx)*psi_half(ihalf,n) + 0.5*dy*psi_half(ihalf,n)
else
v_D2(j-1) = dx*psi_half(ihalf-1,n)
end if
else
if (j==2) then
v_D2(j-1) = 0.5*dx*psi_half(i-1,n) + (1.0-dx)*psi_half(i,n) + 0.5*dx*psi_half(i+1,n)&
+0.5*dy*psi_half(i,n)
else
v_D2(j-1) = 0.5*dx*psi_half(i-1,n)
end if
end if
end do
call tridag(a_D2,c_D2,o_D2,jmax-1)
do j = 1,jmax-1
psi2(i,n+1) = o_D2(j)
end do
end do
! boundary condition
! y=0
do i = 1,ihalf
psi2(i,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,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,n+1) = psi2(i,n+1)
end do
end do
end do
do i = ihalf+1,:) = psi(imax+1-i,:)
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
u(j) = u(j) - gam(j+1)*u(j+1)
end do
return
end subroutine
end program
我很好奇的是一开始就很好。
一开始，我设置了imax=200（x坐标的最大索引）。
这是原始代码。
program HW1_1_ADI
implicit none
real :: delx,dimension(101,dimension(201) :: x
real,dimension(100) :: a_D,file='HW1.txt')
2 format (101((200(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,-1
u(j) = u(j) - gam(j+1)*u(j+1)
end do
return
end subroutine
end program
效果很好。
我的代码有什么问题？
解决方法
我建议你学习如何使用你的编译器来帮助你诊断这些问题——在开发时总是打开运行时错误检查。以下是 gfortran 可以告诉您的信息，请注意 -fcheck=all 和 -g 标志，但我会推荐我使用的所有标志：
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
由此看来，您似乎没有正确调整 nMax 中的 tridag 参数 - 我建议您还了解可分配数组以及如何以某种方式避免此问题，这意味着您不会不必为您尝试的每种不同尺寸更改代码。