我用 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
我可以将完整的代码放在代码块中,但为了便于阅读,我认为最好不要这样做。如果有人有任何建议,我将不胜感激。
谢谢!