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.