Итерационный решатель Fortran Gauss Siedel, использующий OpenMP, не сходится - PullRequest
Новая
       13

Итерационный решатель Fortran Gauss Siedel, использующий OpenMP, не сходится

0 голосов
/ 07 апреля 2019

Я хочу написать итерационный решатель Gauss Siedel с использованием OpenMP.Мой решатель не сходится к правильному результату, и я не мог понять, почему.

Основным уравнением является стационарное уравнение теплопроводности: del ^ 2 (T) = S, где S - функция источника тепла.S = -35 * pi / 2 * cos (pi * x / 2) * cos (3 * pi * y / 2) * cos (5 * pi * z / 2)

Я реализую OpenMP внутрицикл Dowhile, поскольку он не позволяет начать параллельный цикл do while.Есть ли способ изменить цикл do while на цикл do?

Результат сходится без параллельных вычислений, но после добавления openmp он больше не сходится.

Вот мой код: 

    PROGRAM GAUSS_MP
        USE omp_lib
        IMPLICIT NONE

        REAL*8, DIMENSION(:,:,:), ALLOCATABLE :: S, T
        REAL*8 :: X, Y, Z
        REAL*8, PARAMETER :: PI=2*ASIN(1.0)
        REAL*8 :: DX                     ! STEP SIZE DX=DY=DZ
        REAL*8, PARAMETER :: TOL=1.0E-5  ! TOLERANCE 
        REAL*8 :: DUMAX
        REAL*8 :: T_OLD
        REAL*8 :: T1,T2

        INTEGER, PARAMETER :: N=100        ! GRID NUMBER, START FROM 1
        INTEGER, PARAMETER :: ITERMAX=1E5   ! MAXIMUM ITERATION
        INTEGER :: I, J, ITER, K
        INTEGER :: POINT_NUM
        INTEGER, PARAMETER :: NUM_THREADS=32
    !    INTEGER :: A

        ! INITIALIZE OPENMP THREADS
        CALL OMP_SET_NUM_THREADS(NUM_THREADS)


        ! COMPUTE STEP SIZE
        DX=2.0/REAL(N-1, KIND=8)  
    !    PRINT *, "DX=", DX

        ! INITIALIZE THE T ARRAYS AND S(I)
        ALLOCATE(S(N,N,N), T(N,N,N))

            ! INITIAL GUESS
            T=1.0

        ! BOUNDARY CONDITION
        T(1,:,:)=0.0; T(N,:,:)=0.0; T(:,:,1)=0.0; T(:,:,N)=0.0;
        T(:,1,:)=0.0; T(:,N,:)=0.0;

        X=0.0D0 ! COORDINATES
        Y=0.0D0

        S=0.0D0 ! SOURCE

            ! SOURCE MATRIX
            !$OMP PARALLEL DO PRIVATE(I,J,K)
             DO K=2,N-1
                Z=-1.0+(K-1)*DX
                DO I=2,N-1
                    Y=-1.0+(I-1)*DX
                    DO J=2,N-1
                        X=-1.0+(J-1)*DX
                        S(I,J,K)=-35.0*PI/2.0*COS(PI*X/2.0)*COS(3.0*PI*Y/2.0)&
                                *COS(5.0*PI*Z/2.0)
                    END DO
                END DO
            END DO
            !$OMP END PARALLEL DO

        ! GAUSS-SEIDEL ITERATION
        PRINT *, 'PARALLEL START'
        T1=OMP_GET_WTIME()

        ITER=0
        DUMAX=1.0D0 ! UPDATE DIFFERENCE
        DO WHILE(ITER <= ITERMAX .AND. DUMAX >= TOL)
            ITER=ITER+1
            DUMAX=0.0D0
            !$OMP PARALLEL PRIVATE(T_OLD, K, I, J, DUMAX)
            !$OMP DO REDUCTION(MAX:DUMAX)
            DO K=2,N-1
                DO I=2, N-1
                    DO J=2, N-1
                        T_OLD=T(I,J,K)
                        T(I,J,K)=1.0/6.0*(T(I+1,J,K)+T(I-1,J,K)+T(I,J-1,K)+T(I,J+1,K) &
                                    +T(I,J,K+1)+T(I,J,K-1) &
                                    -DX**2*S(I,J,K))
                        DUMAX=MAX(DUMAX, ABS(T_OLD-T(I,J,K)))
                    END DO
                END DO
            END DO
            !$OMP END DO
           !$OMP END PARALLEL
        END DO

        T2=OMP_GET_WTIME()

    END PROGRAM GAUSS_MP

1 Ответ

0 голосов
/ 07 апреля 2019

Gauss-Seidel - это последовательный алгоритм, который нельзя легко распараллелить.Если вы посмотрите на обновление массива T, вы увидите, что вы читаете значения из других потоков, которые могли или не могли быть обновлены, когда текущий поток пытается их обработать.Это типичное условие гонки.

У вас есть в основном два варианта:

  • используйте наклон петли, чтобы «повернуть» гнездо петли на 45 градусов, и использовать волновой фронт, чтобы идтичерез сетку.Таким образом, обновленные ячейки будут доступны, когда текущие потоки захотят прочитать обновленное значение.

  • использовать функцию OpenMP 4.5 'упорядоченная зависимость', чтобы выразить зависимость данных в вашем коде и позволитьКомпилятор OpenMP добавляет правильную синхронизацию, чтобы избежать состояния гонки.

...