[求助] Program received signal SIGSEGV

子程序：
module mod_fdm_tool_C3
    use iso_fortran_env
    implicit none
    real(real64), parameter :: pi = 4.d0*datan(1.d0)
contains

subroutine Grid_2D_Type(nx, nt, xmin, xmax, tmin, tmax, hx, ht, x, t)
    implicit none
    integer(int32), intent(in)  :: nx, nt
    real(real64),   intent(in)  :: xmin, xmax, tmin, tmax
    real(real64),   intent(out) :: hx, ht, x(0:nx), t(0:nt)
    integer(int32) :: i, j

    hx = (xmax - xmin) / dble(nx)
    ht = (tmax - tmin) / dble(nt)

    do concurrent (i = 0 : nx)
        x(i) = xmin + dble(i) * hx
    end do
    x(nx) = xmax

    do concurrent (j = 0 : nt)
        t(j) = tmin + dble(j) * ht
    end do
    t(nt) = tmax

    return
end subroutine Grid_2D_Type



!===============================================================================
pure elemental function f_rhs(x,t) result(res)
    implicit none
    real(real64), intent(in) :: x, t
    real(real64) :: res

    res = 0

    return
end function f_rhs

!===============================================================================
pure elemental function Solu_Exact(x,t) result(res)
    implicit none
    real(real64), intent(in) :: x, t
    real(real64) :: res

    res = dexp(x+t)

    return
end function Solu_Exact





subroutine Apply_DBC(x, t,u)
    implicit none
    real(real64), intent(in)  :: x(:), t(:)
    real(real64), intent(out) :: u(:,:)
    integer(int32) :: nt

    nt = size(t)

    u(:, 1) = dexp(x)
    u(1,:) = dexp(t)
    u(nt,:) = dexp(1.0+t)



    return
end subroutine Apply_DBC

!======================================================
subroutine Gene_Alg_Sys_DBC(ht,hx,nx,r,coef1,coef2)
    ! Algebraic system, i.e., Au = f
    implicit none
    integer(int32) ,intent(in):: nx
    real(real64), intent(in)  :: hx,ht
    real(real64), intent(out) :: coef1(nx-2,3),coef2(nx-2,3),r

    r = ht/hx/hx


    ! all the following variable are arraies, not scalars
    coef1(:,1) = -r/2.d0
    coef1(:,2) = 1+r
    coef1(:,3) = -r/2.d0

    coef2(:,1) = r/2.d0
    coef2(:,2) = 1-r
    coef2(:,3) = r/2.d0



    coef1(1,1) = 0.d0
    coef1(nx-2,3) = 0.d0

    coef2(1,1) = 0.d0
    coef2(nx-2,3) = 0.d0



    return
end subroutine Gene_Alg_Sys_DBC


!======================================================
subroutine Solver_Crank_Nicolson (coef1,coef2,L,u)

    implicit none
    integer(int32) :: nx,nt,i,j
    real(real64), intent(in)  :: coef1(:,:),coef2(:,:)

    real(real64) :: L(:),r

    real(real64) :: u(:,:)

    real(real64),allocatable::g(:),w(:)




    do j = 1, nt-1
      L(1)=coef2(1,2)*u(2,j)+coef2(1,3)*u(3,j)+r*(u(1,j)+u(1,j+1))/2
      L(nx-2)=coef2(nx-2,1)*u(nx-2,j)+coef2(nx-2,2)*u(nx-1,j)+r*(u(nx,j)+u(nx,j+1))/2
      L(2:nx-3)=coef2(2:nx-3,1)*u(2:nx-3,j)+coef2(2:nx-3,2)*u(3:nx-2,j)+coef2(2:nx-3,3)*u(4:nx-1,j)

      allocate(g(nx),w(nx))

      g(1)=L(1)/coef1(1,2)
      w(1)=coef1(1,3)/coef1(1,2)

      do i=2,nx-2
        w(i)=coef1(i,2)-coef1(i,1)*w(i-1)
        g(i)=(L(i)-coef1(i,1)*g(i-1))/w(i)
        w(i)=coef1(i,3)/w(i)
      end do
      u(nx,:)=g(nx)

      do i=nx-3,2,-1
        u(i,j+1)=g(i)-w(i)*u(i+1,j+1)
      end do


    end do

     deallocate(g,w)

    return
end subroutine Solver_Crank_Nicolson

subroutine Make_Directory(dirname)
    implicit none
    character(len=:), allocatable, intent(in) :: dirname
    call system("mkdir "//trim(AdjustL(dirname)))
end subroutine Make_Directory


end module mod_fdm_tool_C3

主程序：

!=============================================================
! cLi, 2022-11-24
!
! Solving example 3.1.1 in p.67:
!
!   -u'(t) - u''(x) = 0,  x \in (0, 1),t \in(0,1]
!   u(x,0) =  exp(x),   x \in (0, 1)
!   u(0,t) = exp(t) ,u(1,t) = exp(1+t) ,t \in(0,1]
!
!
!   
!   
!   u(x,t) = exp(x+t)
!=============================================================
program Crank_Nicolson_2D_DBC

    use iso_fortran_env
    use mod_fdm_tool_C3

implicit none

integer(int32):: nx, nt, m, k, k_max, fileunit,  i,j
real(real64)   :: hx, xmin, xmax
real(real64)   :: ht, tmin, tmax

real(real64)   :: r,iter_tol, iter_fac, error = 0.d0
real(real64), allocatable :: x(:), t(:),L(:),  u(:,:), v(:,:),coef1(:,:),coef2(:,:)
character(len=:), allocatable :: dirname, filename, output_fmt
logical :: istat, debug = .false.

namelist /nml_para/ nx, nt, xmin, xmax, tmin, tmax, k_max

print '(1X,A)', 'Example_311'

filename = 'nml_para_DBC.txt'
open(newunit=fileunit, file=filename)
read(unit=fileunit, nml=nml_para)
close(fileunit)

dirname = 'Out'
inquire(file=dirname, exist=istat)
if (.not. istat) then
    call Make_Directory(dirname)
end if


filename = dirname//'/results.txt'
open(newunit=fileunit, file=trim(filename))
write(fileunit, '(a)') 'Inf error at x = 0.5, 0.5,  0.5, 0.5,  0.5,'
write(fileunit, '(a)') '         and t = 0.1, 0.2,  0.3, 0.4,  0.5 are'

do k=1,k_max
   m  = 2**(k-1)
    nx = 10*m
    nt = 10*m
    print '(1X, A, I3, A,I3, A)', 'The spatial step size of 1/', nx, ',the time step size of 1/',nt,' is running...'

    if (allocated(u)) deallocate(x, t,  u, v)
    allocate(x(0:nx), t(0:nt),  u(0:nx,0:nt), v(0:nx,0:nt))

     ! start the FDM solving procedure
    call Grid_2D_Type (nx, nt, xmin, xmax, tmin, tmax, hx, ht, x, t)

    call Apply_DBC(x, t, u)

    call Gene_Alg_Sys_DBC(ht,hx,nx,r,coef1,coef2)

    do concurrent (i = 0 : nx)

        v(i,:) = Solu_Exact(x(i),t)
    end do


    call Solver_Crank_Nicolson (coef1,coef2,L,u)

     output_fmt = '(A, I3, 2X, 6(ES9.3, 2X), F5.2)'
    write(fileunit, output_fmt) '1/', nx,'1/', nt,       &

            abs(u(nx/2*1,nt/10*1)-v(nx/2*1,nt/10*1)), &

            abs(u(nx/2*1,nt/10*2)-v(nx/2*1,nt/10*2)), &
            abs(u(nx/2*1,nt/10*3)-v(nx/2*1,nt/10*3)), &
            abs(u(nx/2*1,nt/10*4)-v(nx/2*1,nt/10*4)), &
            abs(u(nx/2*1,nt/10*5)-v(nx/2*1,nt/10*5)), &
            maxval(abs(u-v)), error/maxval(abs(u-v))

    error = maxval(abs(u-v))
end do

close(fileunit)
print *, ''
print *, 'The program is done.'

return
end program Crank_Nicolson_2D_DBC