Я столкнулся с этой ошибкой при решении уравнения диффузии на фортране с использованием метода неявного переменного направления (ADI).
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,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
Что меня интересует, так это то, что сначала это хорошо.
Сначала я установил imax=200 (максимальный индекс координаты x).
Вот исходный код.
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
Это хорошо работает.
Что не так в моем коде?
