Weird issues with private variables

Accelerated Computing HPC Compilers Legacy PGI Compilers
#1

I’ll list out the part of the code that I am having trouble with.

   subroutine acc_getVisDeri(Nelements, Npts, Nvar, &
                                     u, ux, uy, uz,       &
                                     gradU,               &
                                     sub1,                &
                                     elem, pt,            &
                                     Cv,                  &
                                     S1, S2, S3, gradT, div)
      !Get all viscous derivatives
      !Strain, T_x, divergence

      !$acc routine seq

      integer(c_int), intent(in)  :: Nelements, Npts, Nvar
      integer(c_int), intent(in)  :: sub1, elem, pt
      real(c_double), intent(in)  :: u(Npts, Nvar, Nelements)  
      real(c_double), intent(in)  :: ux(Npts, Nvar, Nelements)
      real(c_double), intent(in)  :: uy(Npts, Nvar, Nelements)
      real(c_double), intent(in)  :: uz(Npts, Nvar, Nelements)
      real(c_double), intent(in)  :: gradU(Npts, Nvar, Nelements)
      real(c_double), intent(in)  :: Cv

      real(c_double), intent(out) :: S1, S2, S3
      real(c_double), intent(out)  :: gradT
      real(c_double), intent(out)  :: div

      real(c_double)             :: rho, vel, velx, vely, velz, E
      real(c_double)             :: gradM1, gradM2, gradM3, gradE
      real(c_double)             :: gradVel1, gradVel2, gradVel3
      real(c_double)             :: gradRho1, gradRho2, gradRho3

      !Splitting up calculation of strain rates

      rho    = u(pt, 1, elem)
      velx   = u(pt, 2, elem)/rho
      vely   = u(pt, 3, elem)/rho
      velz   = u(pt, 4, elem)/rho
      E      = u(pt, 5, elem)
      gradRho1 = gradU(pt, 1, elem)
      gradM1   = gradU(pt, 2, elem)
      gradM2   = gradU(pt, 3, elem)
      gradM3   = gradU(pt, 4, elem)
      gradE    = gradU(pt, 5, elem)

      gradVel1 = ONE/rho*(gradM1 - velx*gradRho1)
      gradVel2 = ONE/rho*(gradM2 - vely*gradRho1)
      gradVel3 = ONE/rho*(gradM3 - velz*gradRho1)

      S1  = HALF*(gradVel1)
      S2  = HALF*(gradVel2)
      S3  = HALF*(gradVel3)

      gradT = ONE/(rho*Cv)*(gradE - E/rho*gradRho1 - &
          rho*(velx*gradVel1 + vely*gradVel2 + velz*gradVel3))

      vel       = u(pt, sub1, elem)/u(pt, 1, elem)

      gradRho1  = ux(pt,    1, elem)
      gradM1    = ux(pt, sub1, elem)
      gradRho2  = uy(pt,    1, elem)
      gradM2    = uy(pt, sub1, elem)
      gradRho3  = uz(pt,    1, elem)
      gradM3    = uz(pt, sub1, elem)

      gradVel1 = ONE/rho*(gradM1 - vel*gradRho1)
      gradVel2 = ONE/rho*(gradM2 - vel*gradRho2)
      gradVel3 = ONE/rho*(gradM3 - vel*gradRho3)

      S1  = S1 + HALF*(gradVel1)
      S2  = S2 + HALF*(gradVel2)
      S3  = S3 + HALF*(gradVel3)

      gradM1    = ux(pt, 2_c_int, elem)
      gradM2    = uy(pt, 3_c_int, elem)
      gradM3    = uz(pt, 4_c_int, elem)

      gradVel1 = ONE/rho*(gradM1 - velx*gradRho1)
      gradVel2 = ONE/rho*(gradM2 - vely*gradRho2)
      gradVel3 = ONE/rho*(gradM3 - velz*gradRho3)

      div = gradVel1 + gradVel2 + gradVel3

    end subroutine acc_getVisDeri


    subroutine acc_evaluate_vis_flux(Nelements, Npts, Nvar, &
                                      u, ux, uy, uz, &
                                      Fv, Gv, Hv, &
                                      mu, Cv, Cp, Pr)

      integer(c_int), intent(in) :: Nelements, Npts, Nvar
      real(c_double), intent(in) :: u(Npts, Nvar, Nelements)
      real(c_double), intent(in) :: ux(Npts, Nvar, Nelements)
      real(c_double), intent(in) :: uy(Npts, Nvar, Nelements)
      real(c_double), intent(in) :: uz(Npts, Nvar, Nelements)

      real(c_double), intent(out) :: Fv(Npts, Nvar, Nelements)
      real(c_double), intent(out) :: Gv(Npts, Nvar, Nelements)
      real(c_double), intent(out) :: Hv(Npts, Nvar, Nelements)

      real(c_double), intent(in) :: mu, Cv, Cp, Pr

      real(c_double) :: rho, Mx, My, Mz, E
      real(c_double) :: velx, vely, velz

      real(c_double) :: T_x, T_y, T_z

      real(c_double) :: div
      real(c_double) :: tauxx, tauxy, tauxz
      real(c_double) :: tauyy, tauyz, tauzz

      real(c_double) :: Sdot

      real(c_double) :: S1, S2, S3

      integer(c_int) :: elem_index, p_indx, sub1

      !$acc enter data present_or_copyin(Cv, Nelements, Nvar, Npts)

      !$acc parallel loop independent gang  &
      !$acc present(u)                              &
      !$acc present(ux, uy, uz)                     &
      !$acc present(Cv, Nelements, Nvar, Npts)      &
      !$acc present(Fv, Gv, Hv)                     
      do elem_index = 1, Nelements
          !$acc loop independent vector private(S1, S2, S3) &
          !$acc private(T_x, T_y, T_z, div)
          do p_indx = 1, Npts
              rho = u(p_indx, 1, elem_index)
              Mx  = u(p_indx, 2, elem_index)
              My  = u(p_indx, 3, elem_index)
              E   = u(p_indx, 4, elem_index)
              Mz  = ZERO
#ifdef DIM3
              Mz  = u(p_indx, 4, elem_index)
              E   = u(p_indx, 5, elem_index)
#endif
              ! velocity 
              velx = Mx/rho
              vely = My/rho
              velz = Mz/rho

              sub1 = 2_c_int
              call acc_getVisDeri(Nelements, Npts, Nvar,  &
                                     u, ux, uy, uz,       &
                                     ux,                  &
                                     sub1,                &
                                     elem_index, p_indx,  &
                                     Cv,                  &
                                     S1, S2, S3, T_x, div)
              S1 = S1 + HALF*( - TWO3RD*div)

              tauxx = 2*(mu)*S1
              tauxy = 2*(mu)*S2
              tauxz = 2*(mu)*S3

              Fv(p_indx, 1, elem_index) = ZERO
              Fv(p_indx, 2, elem_index) = -tauxx
              Fv(p_indx, 3, elem_index) = -tauxy
              Fv(p_indx, 4, elem_index) = -(velx*tauxx + vely*tauxy)  - mu*Cp/Pr*T_x
#ifdef DIM3
              Fv(p_indx, 4, elem_index) = -tauxz
              Fv(p_indx, 5, elem_index) = -(velx*tauxx + vely*tauxy + velz*tauxz) - mu*Cp/Pr*T_x
#endif
              sub1 = 3_c_int
              call acc_getVisDeri(Nelements, Npts, Nvar, &
                                u, ux, uy, uz, uy,       &
                                sub1,                    &
                                elem_index, p_indx,      &
                                Cv,                      &
                                S1, S2, S3, T_y, div)

              S2 = S2 + HALF*( - TWO3RD*div)

              tauyy = 2*(mu)*S2
              tauyz = 2*(mu)*S3

              !Y-viscous fluxes 
              Gv(p_indx, 1, elem_index) = ZERO
              Gv(p_indx, 2, elem_index) = -tauxy
              Gv(p_indx, 3, elem_index) = -tauyy
              Gv(p_indx, 4, elem_index) = -(velx*tauxy + vely*tauyy) - mu*Cp/Pr*T_y

#ifdef DIM3
              Gv(p_indx, 4, elem_index) = -tauyz
              Gv(p_indx, 5, elem_index) = -(velx*tauxy + vely*tauyy + velz*tauyz) - mu*Cp/Pr*T_y
#endif

#ifdef DIM3
              call acc_getVisDeri(Nelements, Npts, Nvar, &
                                u, ux, uy, uz, uz,       &
                                4_c_int,                 &
                                elem_index, p_indx,      &
                                Cv,                      &
                                S1, S2, S3, T_z, div)

              S3 = S3 + HALF*( - TWO3RD*div)

              tauzz = 2*(mu)*S3
              !Z-viscous fluxes 
              Hv(p_indx, 1, elem_index) = ZERO
              Hv(p_indx, 2, elem_index) = -tauxz
              Hv(p_indx, 3, elem_index) = -tauyz
              Hv(p_indx, 4, elem_index) = -tauzz
              Hv(p_indx, 5, elem_index) = -(velx*tauxz + vely*tauyz + velz*tauzz) - mu*Cp/Pr*T_z
#else
              do var = 1, nvar
                  Hv(var) = ZERO
              end do
#endif

            end do
        end do
        !$acc end parallel

        !$acc exit data delete(Cv, Nelements, Nvar, Npts)

    end subroutine acc_evaluate_vis_flux

Sorry for the long copy paste.

Now, what is weird is that I have to specify S1, S2, S3 as private even though I thought by default openacc makes all scalars private. In addition even though, I have made T_x, T_y, T_z, div private, I see that if I don’t make them private the code works fine.

Another problem is that the compiler prompts me to copyin Nelements, Npts, Nvar even though, again, they are scalars.

#2

Hi vsingh96824

Now, what is weird is that I have to specify S1, S2, S3 as private even though I thought by default openacc makes all scalars private. In addition even though, I have made T_x, T_y, T_z, div private, I see that if I don’t make them private the code works fine.

Another problem is that the compiler prompts me to copyin Nelements, Npts, Nvar even though, again, they are scalars

Yes, this is one of the few cases where you need to manually privatize scalars.

In Fortran, arguments are passed by address. This causes these variables to “escape”, meaning that the compiler can no longer analyze how they are being used and hence must assume the worst case scenario that other global references can now point to them. The same would be true in C if the scalars were passed by reference.

The easiest solution to this issue is to use the F2003 “value” attribute when declaring the variables in “acc_getVisDeri”. This tells the compiler to pass the variables by value rather than address.

For example:

integer(c_int), intent(in), value  :: Nelements, Npts, Nvar

Hope this helps,
Mat

#3

Thanks for the reply Mat.

When you say this being one of the few cases, what specifically defines the case. Is it the calling of a subroutine from within an openacc loop?

Vikram

#4

Hi Viram,

When you say this being one of the few cases, what specifically defines the case. Is it the calling of a subroutine from within an openacc loop?

When passing a scalar variable by address/reference to an OpenACC routine from within an OpenACC compute region, the compiler can no longer make the scalar implicitly private.

  • Mat