使用 openmp 并行化 N 体模拟

问题描述 投票：0回答：1

我正在尝试并行化 Fortran 代码以进行 N 体模拟，但与顺序版本相比，它总是需要更长的执行时间，并且执行时间随着用于并行化的线程数量的增加而增加。

我尝试一次并行化一个子例程，但每个子例程都会出现问题，我尝试使用所有归约子句、Flush、Barrier 来避免数组中的竞争条件，但没有得到好的结果。即使在计算标量势能时也会出现同样的问题。即使在主文件中并行化时间步骤。有一点不明白希望你能帮助我。

module functions
    implicit none
    contains
  subroutine initial_speed(nb_particles,v,position)
     implicit none
     integer, parameter :: dp = kind(1.0d0)
     integer :: i, nb_particles
     real(kind=dp), ALLOCATABLE:: position(:,:) , omega(:,:) , v(:,:)
     !$OMP parallel shared(position,v) private(i)
    ! $OMP do schedule (dynamic)

     do i = 1, nb_particles 
        v(1, i) =  - position(2,i)
        v(2, i) =  position(1, i) 
        v(3, i) =  0
     end do
     !$OMP END DO
     !$OMP END PARALLEL
  end subroutine initial_speed

 subroutine acceleration(nb_particles,position,ai1,m,G )
     implicit none
     integer :: i,j,nb_particles
     integer, parameter :: dp = kind(1.0d0)
     real(kind=dp) :: G,m ,x ,y, z , r, dist(3)
     real(kind=dp), ALLOCATABLE:: position(:,:), ai1(:,:)
     ai1(:,:) =0.
     !$OMP parallel shared (position ,G,m,ai1) private(dist,i,j) 
      !$OMP do  schedule (DYNAMIC) reduction(+:ai1)
     do i= 1, nb_particles
        do j= i+1, nb_particles
             !  dist(1) = position(1,j) - position(1,i)
              ! dist(2) = position(2,j) - position(2,i)
              ! dist(3) = position(3,j) - position(3,i)
                dist(:) = position(:,j) - position(:,i)

                !ai1(1,i) =ai1(1,i) +(G*m*(x))/((sqrt(abs(r)**2+0.02**2))**3 )
                !ai1(2,i) = ai1(2,i)+(G*m*(y))/((sqrt(abs(r)**2+0.02**2))**3 )
                !ai1(3,i) = ai1(3,i)+(G*m*(z))/((sqrt(abs(r)**2+0.02**2))**3 )
               
                ai1(:,i) =ai1(:,i) +(G*m*(dist(:)))/(dot_product(dist,dist)+0.04**2)**1.5 
                ai1(:,j) =ai1(:,j) -(G*m*(dist(:)))/(dot_product(dist,dist)+0.04**2)**1.5
                              
        enddo
     enddo
    !$OMP END DO
      !$OMP END PARALLEL
    end subroutine acceleration


    subroutine leapfrog(nb_particles,position,v,dt,a,ai1,m,G)
     implicit none
     integer, parameter :: dp = kind(1.0d0)
     integer :: i ,nb_particles
     real(kind=dp) :: dt, G,m
     real(kind=dp), ALLOCATABLE:: position(:,:), v(:,:), a(:,:),ai1(:,:)
     call acceleration(nb_particles,position,a,m,G)
     !$OMP PARALLEL SHARED(position,a, dt) Private(i)
     !$OMP DO SCHEDULE(DYNAMIC)
     do i = 1, nb_particles
         !position(1,i)=position(1,i)+v(1,i)*dt + 0.5*a(1,i)*(dt**2)
         !position(2,i)=position(2,i)+v(2,i)*dt + 0.5*a(2,i)*(dt**2)
         !position(3,i)=position(3,i)+v(3,i)*dt + 0.5*a(3,i)*(dt**2)
         position(:,i)=position(:,i)+v(:,i)*dt + 0.5*a(:,i)*(dt**2)
     end do
    !$OMP END DO 
    !$OMP END PARALLEL
     
     call  acceleration(nb_particles,position,ai1,m,G)
     
     !$OMP PARALLEL SHARED(v,a,ai1, dt) Private(i)
     !$OMP do schedule (Dynamic) 
          do  i = 1, nb_particles
         !v(1,i) = v(1,i) + 0.5 * (a(1,i)+ai1(1,i))*dt
         !v(2,i) = v(2,i) + 0.5 * (a(2,i)+ai1(2,i))*dt
         !v(3,i) = v(3,i) + 0.5 * (a(3,i)+ai1(3,i))*dt
          v(:,i) = v(:,i) + 0.5 * (a(:,i)+ai1(:,i))*dt
     end do
     !$OMP END DO 
     !$OMP END PARALLEL
     a = ai1

    end subroutine

    subroutine KE(kin_e,nb_particles,v,m)
        implicit none
        integer, parameter :: dp = kind(1.0d0)
        integer:: nb_particles,i
        real(dp) :: m, kin_e ,vv 
        real(dp), allocatable :: v(:,:)
         kin_e = kin_e + 0.5 * m * sum(v**2)
    end subroutine

    subroutine PE(pot,nb_particles,G,m,position)
      implicit none
      integer, parameter :: dp = kind(1.0d0)
      integer :: nb_particles, i, j 
      real(dp) :: G, m, pot,x,y,z ,r
      real(dp), allocatable :: position(:,:)
      pot = 0.0
      !$OMP PARALLEL SHARED (position) private(x,y,z,i,j)
      !$OMP DO SCHEDULE (dynamic) reduction(- : pot)
      do i = 1, nb_particles
         do j = i,nb_particles
            if (j/=i) then
              x = position(1,j)- position(1,i)
               y = position(2,j)- position(2,i)
               z = position(3,j)- position(3,i)

               pot = pot - (m**2)*G/ dsqrt(x**2+y**2+z**2 +0.04**2)
            end if
         enddo
      enddo
      !$OMP END DO
      !$OMP END PARALLEL
   end subroutine

  program nbody
    use functions
    implicit none
    integer, parameter :: dp = kind(1.0d0)
    integer :: nb_particles=2000, radius, nb_steps=1000 ,t1,t2,rate
    real(dp):: x,y,z,r, G = 1 ,m, dt=0.01,kin_e=0.0,pot=0.0, time, dist(3)
    integer :: i,t = 0,j
    real(dp), ALLOCATABLE:: position(:,:)  , v(:,:), a(:,:), ai1(:,:), ai1_private(:,:)

    radius = 1
    G = 1.
    dt=0.01
    m = 1./real(nb_particles,kind=dp)
    open(unit=1, file='position_ini.txt', status='old', action='read')
    
    open(unit=100,file="Energy.txt")
    allocate(position(3,nb_particles)) 
    do i =1,nb_particles
        read(1,*) position(1,i), position(2,i), position(3,i)    
    enddo
    
    

    allocate(v(3,nb_particles))
    call initial_speed(nb_particles,v,position)
 
    allocate(a(3,nb_particles))    
    allocate(ai1(3,nb_particles))
 
   call KE(kin_e,nb_particles,v,m)
    call PE(pot,nb_particles,G,m,position)
    write(100,*) kin_e , pot , kin_e+pot, t*dt
  !  !$OMP parallel DO SCHEDULE(runtime)  private(t)
    do t = 1 , 1000
     call leapfrog(nb_particles,position,v,dt,a,ai1,m,G)
     call KE(kin_e,nb_particles,v,m)
     call PE(pot,nb_particles,G,m,position)
     write(100,*) kin_e , pot , kin_e+pot , t
    end do
    !!$OMP END PARALLEL DO
     
    close(100)
    close(1)
  

   


end program

已编辑

这是完整的代码，我使用了 

 integer, parameter :: dp = kind(1.0d0)
而不是 
real(8)
，时间步进无与伦比。，也使用了 
iso_fortran_env
但没有任何积极的结果。

我使用 gfortran 编译器（gcc 11.4）和 openmp 4.5。 CPU是i7 13600 H。我使用time命令测量执行时间 这是我编译和运行程序的方式 

 gfortran main.f90 functions.f90 -fopenmp -O3  -o Nbody

time ./Nbody

fortran openmp
1个回答
0
投票

OpenMP 并不是一根可以加速一切的魔杖，有些算法本身就不能很好地并行化。当他们这样做时，可能有很多原因导致 OpenMP 代码没有预期那么快，甚至比串行代码慢。

如评论中所述，嵌套 OpenMP 区域通常效率不高，如果不是真的有必要，请尝试避免这种情况。

总工作负载必须足够大才能获得并行化的好处并隐藏开销。也许您的代码可以使用更多粒子来更好地扩展。

在您的代码中，您到处都使用 

schedule(dynamic)
，这是一个坏主意，特别是当并行循环每次迭代的工作负载非常低时，如 
initial_speed()
和 
leapfrog
所示：每个线程一次仅执行迭代，大部分时间都花在将迭代分配给线程上。此外，在您的代码中，它会导致错误共享，这是性能杀手。除非您有充分的理由这样做，否则请勿更改默认计划，即 
static
。

由于迭代不平衡，

schedule(dynamic)
唯一有意义的地方是在
acceleration()
例程中。然而，即使在那里，也不能保证真正有效。也许比默认值（1）更大的块会有所帮助：
SCHEDULE(DYNAMIC,10)

我在类似案例中的经验表明，从工作量最少的迭代开始并使用

guided
时间表实际上更高效：

     !$OMP PARALLEL PRIVATE(dist,i,j,tmp) 
     !$OMP DO SCHEDULE(GUIDED) REDUCTION(+:ai1)
     do i= nb_particles, 1, -1
        do j= i+1, nb_particles
           dist(:) = position(:,j) - position(:,i)
           tmp = G*m/(dot_product(dist,dist)+0.04**2)**1.5 
           ai1(:,i) = ai1(:,i) + tmp*dist(:) 
           ai1(:,j) = ai1(:,j) - tmp*dist(:)        
        enddo
     enddo
     !$OMP END DO
     !$OMP END PARALLEL

减少也可能会损害性能，以及内循环为外循环的每次迭代更新整个 

ai1(:,:)
数组的事实，以及一些缓存未命中。您可以通过基本上仅更新 
ai1(:,i)
来避免这种情况，但是您必须扫描内部循环中的整个索引。有些计算会进行两次，但也许这是更好的并行化方法。另一个优点是可以使用默认计划（并且不需要减少，因为现在可以共享
ai1
：

     !$OMP PARALLEL PRIVATE(dist,i,j,tmp) 
     !$OMP DO
     do i= 1, nb_particles
        do j= 1, nb_particles
           dist(:) = position(:,j) - position(:,i)
           tmp = G*m/(dot_product(dist,dist)+0.04**2)**1.5 
           ai1(:,i) = ai1(:,i) + tmp*dist(:) 
        enddo
     enddo
     !$OMP END DO
     !$OMP END PARALLEL
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.