Identical function call leading to segmentation fault on 1 of 6 processors -- MPI

Help
1

I am writing an MPI program for a project in a course on parallel computing. The program is supposed to solve Poisson’s equation using FEM with a mesh of square 9-node Lagrange elements. One part of the code is of course to build the force/right hand vector. In this part of the code I get a segmentation fault on one of 6 processors (when running on such a configuration).

This is the code:

    ! Construct right hand vector (source contributions)
    ! This is only for myElements (excluding overlap)
    allocate(gb(intNodes))
    gb=0.0_dp
    do k=1,sd*dimension
        k2 = myElemNoOverLap(k)
        do j=1,9
            M  = elementNodes(k2,j)
            if (nodeTypes(mod(M,nodesOnEdge), ceiling(real(M)/nodesOnEdge)) == 1) then
                M2 = GlobalToLocal(M, nodesOnEdge)
                ! if (myrank==5) print *, 
                rho = GaussQuad(sourceTerm, 1.0_dp, 0.0_dp, 0.0_dp, 1)
                ! gb(M2) = gb(M2) + GaussQuad(sourceTerm, elemArea, xCenters(k2), yCenters(k2), j)
                ! gb(M2) = gb(M2) + GaussQuad(sourceTerm, elemArea, 0.0_dp, 0.0_dp, j)
            end if
        end do
    end do
    print *, rho

The row I actually want to execute is the one with xCenters etc., but I have reduced the problem so that the same thing happens on every processor. The assignment itself seems to work, as commenting out the print statement makes the segmentation fault go away, so the segmentation fault doesn’t occur until I try to do anything with the assigned value.

I am at a complete loss as to why the segmentation fault only happens on one processor, since the function call is now identical, so I was hoping someone here could have some insight. The system I am running on is not my own but a HPC cluster owned by my university.

The function being called, in case it is necessary:

function GaussQuad(func, area, centerX, centerY, nodeIndex) result(b)
    ! This function takes a function "func" to integrate over a square
    ! domain of area "area" centered at (centerX,centerY).
    ! This is accomplished through Gauss-quadrature of order 8
    ! If the optional argument "nodeIndex" is present, 
    ! the function computes the integral of the function times its shape-function.
    ! This is useful for computing contributions from source functions.
    interface 
        pure function func(xin,yin)
            double precision             :: func
            double precision, intent(in) :: xin, yin
        end function func
    end interface
    real(dp), intent(in)               :: area, centerX, centerY
    real(dp)                           :: b
    real(dp), parameter, dimension(8)  :: w  = [0.3626837833783620, &
                                                0.3626837833783620, &
                                                0.3137066458778873, &
                                                0.3137066458778873, &
                                                0.2223810344533745, &
                                                0.2223810344533745, &
                                                0.1012285362903763, &
                                                0.1012285362903763]
    real(dp), parameter, dimension(8)  :: xi = [-0.1834346424956498, &
                                                 0.1834346424956498, &
                                                -0.5255324099163290, &
                                                 0.5255324099163290, &
                                                -0.7966664774136267, &
                                                 0.7966664774136267, &
                                                -0.9602898564975363, &
                                                 0.9602898564975363]
    real(dp)                           :: scale, x, y
    integer(1)                         :: i,j
    integer, optional    :: nodeIndex
    scale=sqrt(area)*0.5_dp
    b=0.0_dp
    if (present(nodeIndex)) then
        do i=1,8
            x=scale*xi(i)+centerX
            do j=1,8
                y=scale*xi(j)+centerY
                
                b=b+w(i)*w(j)*func(x,y)*N(xi(i), xi(j), mod(nodeIndex,3), ceiling(real(nodeIndex)/3))
            end do
        end do
    else
        do i=1,8
            x=scale*xi(i)+centerX
            do j=1,8
                y=scale*xi(j)+centerY
                
                b=b+w(i)*w(j)*func(x,y)
            end do
        end do
    end if
    b=b*area*0.25_dp
end function GaussQuad
3

Welcome to Fortran discourse.

If the segmentation fault can be triggered by adding or removing the print statement, then it may be that the problem is elsewhere in the code – bugs involving undefined memory accesses can do this.

Perhaps try running the code with bounds checking and similar compiler options.

BTW, unrelated to this, your code GaussQuad will have less precision than expected. This line

real(dp), parameter, dimension(8)  :: w  = [0.3626837833783620, &
    0.3626837833783620, &
    ! ... elided
    0.1012285362903763]

will cause w to be populated with single precision values (converted to dp) because the right hand side array is single precision, notwithstanding the many significant figures. In Fortran one should specify the precision of constants to avoid truncation to single precision, like

real(dp), parameter, dimension(8)  :: w  = [0.3626837833783620_dp, &
    0.3626837833783620_dp, &
    ! ... elided
    0.1012285362903763_dp]

and similarly for xi.

4

Thanks for the welcoming, and pointing out the precision issue!

Okay, in that case I will have to keep searching for the error someplace earlier in the code, thank you!