OpenMP Block gives false results

Jack_S_ · ‎12-02-2014

Hi all,

I would appreciate your point of view where I might did wrong using OpenMP. I parallelized this code pretty straight forward - yet even with single thread (i.e., call omp_set_num_threads(1)) I get wrong results.

I have checked with Intel Inspector, and I do not have a race condition, yet the Inspector tool indicated (as a warning) that a thread might approach other thread stack (I have this warning in other code I have, and it runs well with OpenMP). I'm pretty sure this is not relate to the problem.

Thanks, Jack.

SUBROUTINE GR(NUMBER_D, RAD_D, RAD_CC, SPECT)

  use TERM,only: DENSITY, TEMPERATURE, VISCOSITY, WATER_DENSITY, & 
                 PRESSURE, D_HOR, D_VER, D_TEMP, QQQ, UMU
  use SATUR,only: FF, A1, A2, AAA, BBB, SAT
  use DELTA,only: DDM, DT
  use CONST,only: PI, G

  IMPLICIT NONE

  INTEGER,INTENT(IN) :: NUMBER_D
  DOUBLE PRECISION,INTENT(IN) :: RAD_CC(NUMBER_D), SPECT(NUMBER_D)
  DOUBLE PRECISION,INTENT(INOUT) :: RAD_D(NUMBER_D)

  DOUBLE PRECISION :: R3, DR3, C2, C0, P, Q, RAD_CR, SAT_CR, C4, A, & 
                      C, D, CC, DD, CC2, DD2, RAD_ST, DRAA, DRA, DM, X1                    
  INTEGER :: I

  DDM = 0.0D0

  !$OMP PARALLEL DO DEFAULT(SHARED) & 
  !$OMP PRIVATE(I,R3,DR3,C2,C0,P,Q,SAT,SAT_CR,C4,A) &
  !$OMP PRIVATE (C,D,CC,DD,CC2,DD2,RAD_ST,DRAA,DRA,DM,RAD_CR,X1) &
  !$OMP REDUCTION (+:DDM)
  DO I=1,NUMBER_D
     R3   = RAD_CC(I)**3
     DR3  = RAD_D(I)**3-R3
     IF(DR3.LT.1.0D-100) DR3 = 1.0D-100

     C2 = -DSQRT(3.0D0*BBB*R3/AAA)
     C0 = -R3
     P  = -0.3333333333D0*C2**2
     Q  = C0+0.074074074D0*C2**3
     CALL CUBIC(P, Q, RAD_CR)
     RAD_CR = RAD_CR - 0.3333333333D0*C2
     SAT_CR = DEXP(AAA/RAD_CR-BBB*R3/(RAD_CR**3-R3))-1.0D0

     DRA  = DT*(SAT+1.0D0-DEXP(AAA/RAD_DROP(I)-BBB*R3/DR3))/ &
           (FF*RAD_D(I))

     IF(SAT.LT.SAT_CR) THEN

        IF(DABS(SAT).LT.1.0D-10) THEN
           P = -BBB*R3/AAA
           Q = -R3
           CALL CUBIC(P, Q, RAD_ST)
           GO TO 22
        END IF


        C4  = DLOG(SAT+1.0D0)
        A   = -AAA/C4
        C   = (BBB-C4)*R3/C4
        D   = -A*R3
        P   = A*C-4.0D0*D
        Q   = -(A**2*D+C**2)

        CALL CUBIC(P, Q, X1)

        CC  = DSQRT(A**2+4.D0*X1)
        DD  = DSQRT(X1**2-4.D0*D)
        CC2 = 0.5D0*(A-CC)
        IF(SAT.LT.0.0D0) THEN
           DD2    = 0.5D0*(X1-DD)
           RAD_ST = 0.5D0*(-CC2+DSQRT(CC2**2-4.0D0*DD2))
        ELSE
           DD2    = 0.5D0*(X1+DD)
           RAD_ST = 0.5D0*(-CC2-DSQRT(CC2**2-4.0D0*DD2))
        END IF

22       CONTINUE

        DRAA = RAD_ST-RAD_D(I)
        IF(ABS(DRAA).LT.ABS(DRA)) THEN
           DRA = DRAA
           DM  = 1.3333333333333333D0*PI*WATER_DENSITY* &
                (RAD_ST**3-RAD_D(I)**3)
        ELSE
           DM  = 4.0D0*PI*WATER_DENSITY*RAD_D(I)**2*DRA
        END IF
        DDM = DDM+SPECT(I)*DM
        RAD_D(I) = RAD_D(I) + DRA
     ELSE
        DM  = 4.0D0*PI*WATER_DENSITY*RAD_D(I)**2*DRA
        DDM = DDM+SPECT(I)*DM
        RAD_D(I) = RAD_D(I) + DRA
     END IF

  END DO
  !$OMP END PARALLEL DO

  RETURN
  END SUBROUTINE GR

Jack_S_ · ‎12-02-2014

The problem is solved.

I'm keeping it as a reference to other people. The argument 'SAT' is set automatically (by default) to zero if declared 'PRIVATE'. It has a value before the parallel block and its 'read only' argument, so it should declared as ''SHARED'. That was the problem here.

Jack.