Fortran嵌套OpenMP循环变量作用域配置与并行失效排查

我编写了离散元法（DEM）程序中的stepper子例程，用于计算粒子i与j间的相互作用并更新力。目前尝试对核心O(N²)嵌套循环进行并行化（暂未使用更复杂的搜索算法），但始终无法成功。我知道问题源于多线程对变量的修改冲突，但不清楚如何正确处理私有/共享/归约变量，曾考虑将数组重定义为矩阵但不确定是否正确。

我通过include引入多子例程共享的变量，对此也欢迎优化建议。以下是stepper子例程代码：

subroutine stepper (tstep)

    use omp_lib

    implicit none

    include "parameter.h"
    include "CB_variables.h"
    include "CB_const.h"
    include "CB_bond.h"
    include "CB_forcings.h"

    integer :: i, j
    integer, intent(in) :: tstep

    ! reinitialize force arrays for contact and bonds
    do i = 1, n
        mc(i)  = 0d0
        mb(i)  = 0d0
        fcx(i) = 0d0
        fcy(i) = 0d0
        fbx(i) = 0d0
        fby(i) = 0d0
    end do
    
    ! put yourself in the referential of the ith particle
    ! loop through all j particles and compute interactions

    !$omp parallel do schedule(dynamic) &
    !$omp private(i,j) &
    !$omp reduction(+:tfx,tfy,fcx,fcy,fbx,fby,m,mc,mb)
    do i = 1, n
        do j = i + 1, n

            ! compute relative position and velocity
            call rel_pos_vel (i, j)

            ! bond initialization
            if ( tstep .eq. 1 ) then
                if ( -deltan(i, j) .le. 5d-1 * r(i)) then ! can be fancier
                    !bond (i, j) = 1
                end if
                call bond_properties (i, j)
            end if

            ! verify if two particles are colliding
            if ( deltan(i,j) .gt. 0 ) then

                call contact_forces (i, j)
                !call bond_creation (i, j) ! to implement

                ! update contact force on particle i by particle j
                fcx(i) = fcx(i) - fcn(i,j) * cosa(i,j)
                fcy(i) = fcy(i) - fcn(i,j) * sina(i,j)

                ! update moment on particule i by particule j due to tangent contact 
                mc(i) = mc(i) - r(i) * fct(i,j) - mcc(i,j)

                ! Newton's third law
                ! update contact force on particle j by particle i
                fcx(j) = fcx(j) + fcn(i,j) * cosa(i,j)
                fcy(j) = fcy(j) + fcn(i,j) * sina(i,j)

                ! update moment on particule j by particule i due to tangent contact 
                mc(j) = mc(j) - r(j) * fct(i,j) + mcc(i,j)

            end if

            ! compute forces from bonds between particle i and j
            if ( bond (i, j) .eq. 1 ) then

                call bond_forces (i, j)
                !call bond_breaking (i, j)

                ! update force on particle i by particle j due to bond
                fbx(i) = fbx(i) - fbn(i,j) * cosa(i,j) +    &
                                    fbt(i,j) * sina(i,j)
                fby(i) = fby(i) - fbn(i,j) * sina(i,j) -    &
                                    fbt(i,j) * cosa(i,j)

                ! update moment on particule i by particule j to to bond
                mb(i) = mb(i) - r(i) * fbt(i,j) - mbb(i, j)

                ! Newton's third law
                ! update force on particle j by particle i due to bond
                fbx(j) = fbx(j) + fbn(i,j) * cosa(i,j) -    &
                                    fbt(i,j) * sina(i,j)
                fby(j) = fby(j) + fbn(i,j) * sina(i,j) +    &
                                    fbt(i,j) * cosa(i,j)
 
                
                ! update moment on particule j by particule i to to bond
                mb(j) = mb(j) - r(i) * fbt(i,j) + mbb(j, i)

            end if

            ! compute sheltering height for particule j on particle i for air and water drag
            call sheltering(i, j)

        end do

        ! compute the total forcing from winds, currents and coriolis on particule i
        call forcing (i)
        !call coriolis(i)

        ! reinitialize total force arrays before summing everything
        m(i)   = 0d0
        tfx(i) = 0d0
        tfy(i) = 0d0

        ! sum all forces together on particule i
        tfx(i) = fcx(i) + fbx(i) + fax(i) + fwx(i)
        tfy(i) = fcy(i) + fby(i) + fay(i) + fwy(i)

        ! sum all moments on particule i together
        m(i) =  mc(i) + mb(i) + ma(i) + mw(i)
    end do
    !$omp end parallel do

    ! forces on side particles for experiments
    call experiment_forces

    ! integration in time
    call velocity
    call euler

end subroutine stepper

相关子例程说明 #

• rel_pos_vel(i, j)：计算粒子i、j间的相对位置、速度、角度等，内部使用公共块中的二维数组；
• bond_properties(i, j)：逻辑同上；
• contact_forces(i, j)：使用局部变量及公共块数组；
• bond_forces(i, j)：逻辑同上；
• sheltering和forcing：逻辑同上；

尝试的解决方法与问题 #

我尝试将tfx、tfy等修改型变量设为私有但未成功，对reduction指令的理解也不清晰；同时不清楚如何处理r(i)这类被多线程访问的常量数组。实际运行时omp_get_num_threads()输出为1，说明未开启并行。我甚至尝试并行化初始化循环：

!$omp do private(i) schedule(dynamic) num_threads(10)
print*, omp_get_num_threads()
do i = 1, n
    mc(i)  = 0d0
    mb(i)  = 0d0
    fcx(i) = 0d0
    fcy(i) = 0d0
    fbx(i) = 0d0
    fby(i) = 0d0
end do
!$omp parallel do

但仍未生效，线程数仍为1。编译采用gfortran，编译选项为-ffast-math、-O3、-fopenmp。

内容的提问来源于stack exchange，提问作者Sogapi