OpenMP 不一致的段错误

计算科学 并行计算 内存管理
2021-12-06 10:56:05

我用 openMP 并行化了我的代码,现在我的代码中出现了这个非常奇怪的错误。如果这是我第一次在计算机上运行,​​它会出现段错误,但如果我再次运行该程序,它运行良好。我从来没有见过这样的事情,并且很好奇是否有人以前见过这个或者有很好的资源或关于如何解决这个问题的提示。这是一个非常简单的循环,它会出现段错误,我将每个循环并行化,默认为 none。下面是循环后面的变量声明。为了便于阅读,我之前发布了循环的精简版本,但为了更清楚起见,要求我填写它:

 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)
   !DIR$ IVDEP:LOOP
   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
           else
               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
   !$OMP END PARALLEL DO

我可以将完整的代码放在代码块中,但为了便于阅读,我认为最好不要这样做。如果有人有任何建议,我将不胜感激。

谢谢!

0个回答
没有发现任何回复~