Код застревает при решении уравнения Пуассона с MPI

Question

Я реализовал простой параллельный решатель уравнения Пуассона с MPI, чтобы ознакомиться с библиотекой MPI. Я разработал код для работы с неопределенным количеством процессоров (в том числе всего с 1).

Код работает и дает хорошие результаты при работе на 1 или 2 процессорах. Тем не менее, он застревает на звонках mpi_send и mpi_recv с 4 процессорами. Поэтому я ожидаю, что моя реализация обмена призракными точками неверна.

Поскольку код слишком большой, чтобы включить здесь, я только включил схему Якоби и обмен данными:

do iter=1,max_iter                      
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
    ! Initial guess, on interior points only               
    Ujacob(min_x+1:max_x-1) = 0._dp                 
    Ujacob_all(0:grid_nx-1) = 0._dp                 

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
    ! Store solution vector from last iteration              
    Uold (:) = Ujacob (:)                  
    Uold_all(:) = Ujacob_all(:)                  

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
    ! Jacobi scheme                     
    do ii=min_x+1,max_x-1                    
     !Ujacob(ii) = 0.5_dp * (Uold (ii-1) + Uold (ii+1) - grid_delta_x**2 * Urhs(ii))    
    end do                        

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
    ! Gather Ujacob vector                    
    call mpi_allgather(Ujacob(0:proc_nx-1), proc_nx, mpi_float,       &    
        & Ujacob_all,   proc_nx, mpi_float, mpi_comm_world, ierror)    

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
    ! Compute error and check if less than tolerance value            
    error = sqrt((sum(Ujacob_all - Uold_all)**2)/dble(grid_nx))          

    if(error < error_tol) return                  
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
    ! Exchange data points                    

    ! Interior processors                   
    if(Xsrc /= -1 .AND. Xdes /= -1) then 
     call mpi_send(Ujacob(  0), 1, mpi_float, Xsrc, 200, mpi_comm_world, ierror)   
     call mpi_send(Ujacob(proc_nx-1), 1, mpi_float, Xdes, 100, mpi_comm_world, ierror)   

     call mpi_recv(Ujacob( -1), 1, mpi_float, Xsrc, 100, mpi_comm_world, stat, ierror)  
     call mpi_recv(Ujacob(proc_nx), 1, mpi_float, Xdes, 200, mpi_comm_world, stat, ierror) 
    ! First processor                    
    elseif(Xsrc == -1) then                  
     call mpi_send(Ujacob(proc_nx-1), 1, mpi_float, Xdes, 100, mpi_comm_world, ierror)   
     call mpi_recv(Ujacob(proc_nx ), 1, mpi_float, Xdes, 200, mpi_comm_world, stat, ierror) 
    ! Last processor                    
    elseif(Xdes == -1) then                 
     call mpi_send(Ujacob(0), 1, mpi_float, Xsrc, 200, mpi_comm_world, ierror)     
     call mpi_recv(Ujacob(-1), 1, mpi_float, Xsrc, 100, mpi_comm_world, stat, ierror)   
    end if                           
end do 

Xsrc и Xdes устанавливаются следующим образом:

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    ! Setting the source and destination neighbors of each processor        
    if(myid == 0) then                    
     Xsrc = -1                     
     Xdes = myid + 1                   
    elseif(myid == nprocs-1) then                 
     Xsrc = myid -1                    
     Xdes = -1                     
    else                       
     Xsrc = myid - 1                   
     Xsrc = myid + 1                   
    end if 

Кроме того, я проверил, что процессор ранга 0 и nprocs-1 действительно соответствуют левым и правым ограниченным процессорам.

Я проверил, что метки хорошо установлены. Кроме того, не стесняйтесь комментировать все, что, по вашему мнению, может быть улучшено.

Answer 1

@Hristo правильно, что ваш код концептуально ошибочен в принципе. Однако почти каждая реализация MPI будет буферизовать MPI_Send для сообщения, содержащего один реальный (хотя, конечно, это не гарантируется), так что это не проблема с вашим кодом.

Я думаю, что вы несовпадающие свои метки - краевые случаи должны иметь бирки наоборот:

elseif(Xsrc == -1) then                  
     call mpi_send(Ujacob(proc_nx-1), 1, mpi_float, Xdes, 200, mpi_comm_world, ierror)   
     call mpi_recv(Ujacob(proc_nx ), 1, mpi_float, Xdes, 100, mpi_comm_world, stat, ierror) 
    ! Last processor                    
    elseif(Xdes == -1) then                 
     call mpi_send(Ujacob(0), 1, mpi_float, Xsrc, 100, mpi_comm_world, ierror)     
     call mpi_recv(Ujacob(-1), 1, mpi_float, Xsrc, 200, mpi_comm_world, stat, ierror)   
    end if  

Несколько других замечаний по коду:

• это очень неэффективно, чтобы вычислить ошибку term with allgather: вы должны суммировать только локальные элементы, а затем вычислить глобальную ошибку с помощью MPI_Allreduce;
• вы должны использовать MPI_REAL, а не MPI_FLOAT для кода Fortran;
• Я не вижу, как наш код может работать на одном процессе - здесь процесс выполнит первое предложение elseif, а затем попробует и отправит несуществующий ранг.

После того как вы подтвердили правильность ваших тегов, вы должны исправить проблемы, отмеченные @Hristo.