Segmentation fault with openmp in Fortran 90 program

Question

I am writing a small program to practice Fortran 90 and openmp. The following example code returns a segmentation fault when I compile it with gfortran -fopenmp heat.f90 -o heat but runs just fine when compiled without the -fopenmp option. Can someone help me understand why this is happening?

I am setting export OMP_NUM_THREADS=4 in my shell (Bash on Ubuntu on Windows 10)

program heat

    use omp_lib
    implicit none

    integer, parameter :: nx = 1000
    integer, parameter :: ny = 800
    real, parameter :: H = 1.
    real, parameter :: alpha = 1.e-4
    real, parameter :: dt = 0.001
    real, parameter :: physical_time = 50.

    real aspect, L, dx, dy
    real, dimension(nx,ny) :: T

    aspect = real(nx)/real(ny)
    L = H*aspect
    dx = L/real(nx)
    dy = H/real(ny)

    T = initialize_T()
    call evolve_field()

    contains

        function initialize_T() result(T)
            implicit none
            real, parameter :: sigma = 0.2
            integer i, j
            real x, y
            real, dimension(nx,ny) :: T
            do i=1, nx
                do j=1, ny
                    x = real(i)/real(nx)*L
                    y = real(j)/real(ny)*H
                    T(i,j) = 10. + &
                             2.* ( &
                                   1./(2.*3.14*sigma**2.) * &
                                   exp(-1. * ( (x-L/2.)**2. + (y-H/2.)**2. ) / (2.*sigma**2.)) &
                                 )
                enddo
            enddo
        end function initialize_T

        subroutine heat_eqn()
            implicit none
            real, dimension(nx,ny) :: Tn
            real d2Tdx2, d2Tdy2
            integer i, j

            Tn(:,:) = T(:,:)

            !$omp parallel shared(T) private(i,j,d2Tdx2,d2Tdy2)
            !$omp do
            do i=2, nx-1
                do j=2, ny-1            
                    d2Tdx2 = ( Tn(i+1,j) - 2.*Tn(i,j) + Tn(i-1,j) ) / dx**2.
                    d2Tdy2 = ( Tn(i,j+1) - 2.*Tn(i,j) + Tn(i,j-1) ) / dy**2.
                    T(i,j) =  Tn(i,j) + dt*(alpha*(d2Tdx2 + d2Tdy2))
                end do
            end do
            !$omp end do                     
            !$omp end parallel

            T(:, 1)  = T(:, 2)
            T(:, ny) = T(:, ny-1)
            T(1, :)  = T(2, :)
            T(nx, :) = T(nx-1, :)

        end subroutine heat_eqn

        subroutine evolve_field()
            implicit none
            integer ts, total_ts, frames_total, output_period_ts, pic_counter
            real progress_pct
            character(len=16) :: pic_counter_str

            total_ts = ceiling(physical_time/dt)
            frames_total = 30*20
            output_period_ts = int(total_ts/frames_total)

            pic_counter = 0

            do ts = 0, total_ts

                if (mod(ts,output_period_ts) .eq. 0) then

                    progress_pct = 100.*real(ts)/real(total_ts)
                    print '(I8, F7.2, A)', ts, progress_pct, "%"

                    !! -- for plotting

                    !open(3, file="T.dat", access="stream")
                    !write(3) T(:,:)
                    !close(3)

                    !write (pic_counter_str,'(I5.5)') pic_counter
                    !call system('gnuplot -c plot3d.gnu '//trim(pic_counter_str))
                    !pic_counter = pic_counter + 1

                end if

                ! -----

                call heat_eqn()

                ! -----

            end do

        end subroutine evolve_field

end program heat

You have likely exhausted available stack memory. I don't use Ubuntu, so cannot help with how to increase the stack size. Running Ubuntu under Windows 10 may exacerbate the problem. — evets, Nov 23 '19 at 20:06
@IanBush good catch! I've updated the question to not have that error. Unfortunately that doesn't fix the segfault though. — HotDogCannon, Nov 23 '19 at 21:05

score 4 · Answer 1 · answered Nov 23 '19 at 23:37

As suggested by evets in the comments, the program crashes due to the stack memory shortage. That is, you created too many too large automatic (stack-allocated) arrays. As an indication you can run the program in Valgrind

> gfortran -fopenmp heat.f90
> valgrind ./a.out

Before printing a lot of "invalid write" errors, it will output something like

==31723== Warning: client switching stacks?  SP change: 0x1ffefffbc0 --> 0x1ffecf2740
==31723==          to suppress, use: --max-stackframe=3200128 or greater

The warning about "client switching stacks" is often present when running out of stack memory.

In a real-world program, you would use allocatable arrays, which are not artificially limited. However, in this simple case I had success with the flag -fno-automatic, which avoids the use of automatic (stack-allocated) arrays, placing them into a different kind of memory instead:

> gfortran -fno-automatic -fopenmp heat.f90
> ./a.out
       0   0.00%
     500   0.17%
    1000   0.33%
    1500   0.50%
...

You can read details about that option in the man page.

Note: I tested the program in a different Linux distribution (not in Windows).

And a working example of stack overflow should gets two points! It is pretty common when one gets a program that is relatively large... And clearly super large compared to the days of old. — Holmz, Nov 24 '19 at 09:18
Upvote for the allocatable arrays comment - this is the proper way to fix this — Ian Bush, Nov 24 '19 at 11:31

Segmentation fault with openmp in Fortran 90 program

1 Answers1