I parallelized my code with openMP, and now have this bug in my code that's really odd. If it's the first time I run on a computer, it segfaults, but if I run the program again it runs fine. I've never seen anything like this and was curious if anyone has seen this before or has good resources or hints on how to fix this. It's a very simple loop that it segfaults on, and I parallelized every loop with default none. The variable declaration followed by the loop is below. I had previously posted a stripped down version of the loop for readability, but it was requested that I fill it in for the sake of more clarity :

 real(dp)    :: q, oneOverq, L, D, Pressure, u, v, En, rho, oneOverrho
 real(dp)    :: normal(Mesh%numDims), wL(Mesh%numFields), cd_dql(Mesh%numFields, Mesh%numEdges)
 real(dp)    :: norm_dx1(Mesh%numDims, Mesh%numDims), norm_dx2(Mesh%numDims, Mesh%numDims)
 real(dp)    :: L_dx(Mesh%numDims, Mesh%numNodes), D_dx(Mesh%numDims, Mesh%numNodes)
 real(dp)    :: p_dql(Mesh%numFields), L_dql(Mesh%numFields, Mesh%numEdges), Res(Mesh%numFields, Mesh%numTri)
 INTEGER(i4) :: i, j, k, m, nodecount, cellNum, edgeNum, Node1, Node2, node, cell, cells(3)

   !$OMP PARALLEL DO Reduction(+:L, D, cL_dql, L_dx, cD_dql, D_Dx) default(none) shared(Mesh, W, grad, gamma, dRdX, Jac_2nd) &
   !$OMP private(cellnum, edgenum, node1, node2, rho, oneoverrho, u, v, En, pressure, normal, norm_dx1, norm_dx2, p_dql, wL, cell, node, cells)
   do cellnum = 1, Mesh%numTri
       edgeNum = Mesh%airfoilEdges(i)
       Node1   = Mesh%edgelist(edgenum,1)
       Node2   = Mesh%edgelist(edgenum,2)
       cellNum = Mesh%edgeList(edgenum,3)
       call cell_stencil(Mesh, cellnum, cells)

       wL = W(:, cellnum)

       call linear_limited_reconstruction(cellnum, edgenum, Mesh, grad(:,:,cellNum), wL)

       rho = wL(1)
       oneOverrho = 1.0_dp/rho
       u   = wL(2)*oneOverrho
       v   = wL(3)*oneOverrho
       En  = wL(4)*oneOverrho

       call get_pressure(rho, En, u, v, pressure)

       p_dql(1) = (gamma - 1._dp)*half*( u**2 + v**2 )
       p_dql(2) = (gamma - 1._dp)*-u
       p_dql(3) = (gamma - 1._dp)*-v
       p_dql(4) = gamma-1._dp

       call norm_dx(node1, node2, Mesh, normal, norm_dx1, norm_dx2)

       L     = L + Pressure*normal(2)
       D     = D + Pressure*normal(1)

       do j = 1, 4
           if (j == 1) then 
               cell = cellnum
               cell = cells(j-1)
           end if
           if (cell /= 0) then
               cL_dql(:,cell) = cL_dql(:,cell) + matmul(P_dql(:),Jac_2nd%wL_rc_dq(:,:,j,edgenum))*normal(2)
               cD_dql(:,cell) = cD_dql(:,cell) + matmul(P_dql(:),Jac_2nd%wL_rc_dq(:,:,j,edgenum))*normal(1)
           end if 
       end do

       L_dx(:, node1)   = L_dx(:, node1) + Pressure*norm_dx1(2, :)
       L_dx(:, node2)   = L_dx(:, node2) + Pressure*norm_dx2(2, :)
       D_dx(:, node1)   = D_dx(:, node1) + Pressure*norm_dx1(1, :)
       D_dx(:, node2)   = D_dx(:, node2) + Pressure*norm_dx2(1, :)

       do j = 1, 6
           node = Mesh%NodeStencil(j, cellnum)
           if (node /= 0) then
               do k = 1, Mesh%numFields 
                   L_dx(:,node) = L_dx(:,node) + P_dql(k)*dRdX%wL_rc_dx(k,:,j,edgenum)*normal(2)              
                   D_dx(:,node) = D_dx(:,node) + P_dql(k)*dRdX%wL_rc_dx(k,:,j,edgenum)*normal(1)
               end do
           end if
       end do  

   end do

I can put in the full code up in the code block, but for the sake of readability I thought it was best not to. If anyone has any advice I'd greatly appreciate it.


  • 2
    $\begingroup$ Often these sorts of heisenbugs happens because of an out-of-bounds memory access. Depending on the operating system's whim of where arrays get allocated into heap memory, that memory access can either segfault or run over into already allocated data, the first of which is obvious and the second of which is not. You can try compiling with -fsanitize=address and see what happens; there are also thread, leak, memory sanitizers, etc. $\endgroup$ Jul 24 '20 at 15:20
  • $\begingroup$ As with all things, I posted and then found the answer (I think). It looks like default none didn't properly work. I hadn't specified two indices, i and j, that I was using to assign the data to the appropriate arrays, so I guess sometimes it corrected had them as private, and other times not so. It looks like your hunch was right, and I'll definintely use that flag in the future for my memory debugging. Do you have any other suggestions as to what sanitizers to use, etc? thanks! $\endgroup$
    – EMP
    Jul 24 '20 at 15:30
  • $\begingroup$ The address and leak sanitizers have been the most useful to me, but I code more in C/C++ than Fortran and there you have to be more careful about memory leaks. What I would definitely recommend is that you put your project is on github or gitlab, set up an automated test suite for it, and enable address sanitizer for the testing build. Once that's set up, it's easy to add more sanitizers later. $\endgroup$ Jul 24 '20 at 17:16
  • $\begingroup$ If default (none) is not working properly that is a VERY broken compiler - what is it? To be honest that would be so broken that I can't believe what you attest to be the problem. But to say more we really need to see at the very least "Computational things" and the declaration of the variables. $\endgroup$
    – Ian Bush
    Jul 24 '20 at 18:32
  • $\begingroup$ It's the intel ifort compiler. I used the i and j variables inside the computational section. From my understanding of default(none) it should have given an error since I hadn't specified either public or private. Since I made this change, I haven't had a problem. $\endgroup$
    – EMP
    Jul 24 '20 at 19:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.