Gauss Siedel iterative solver using OpenMP does not converge

I want to write a Gauss Siedel iterative solver using OpenMP. My solver does not converge to the correct result, and I could not figure it out why.

The governing equation is a steady-state heat equation: del^2(T)=S, where S is the heat source function. S=-35*pi/2*cos(pi*x/2)*cos(3*pi*y/2)*cos(5*pi*z/2)

I implement the OpenMP inside the dowhile loop since it does not let me start parallel form a do while loop. Is there any way to change the do while loop to a do loop?

The result converges without parallel computing, but after adding the openmp, then it does not converge anymore.

Here is my code:

    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

Solution

Gauss-Seidel is a sequential algorithm that cannot be parallelized easily. If you take a look at the update of the T array, you'll see that you're reading values from other threads that may or may not have been updated when the current thread tries to process them. That a typical race condition.

You have basically two options:

use loop skewing to "turn" the loop nest by 45 degrees and use a wave front to go through the grid. That way updated cells will be available when the current threads wants to read the updated value.
use the OpenMP 4.5 feature 'ordered depend` to express the data dependency in your code and let the OpenMP compiler add the proper synchronization to avoid the race condition.