Skip to content

CCE OpenACC 19-21: an acc loop inside a routine seq plus an array-element actual argument reads garbage and drops the store #1815

Description

@sbryngelson

CCE 19.0.0 OpenACC (Frontier, -hacc -h acc_model=auto_async_none -h acc_model=no_fast_addr -O2, MFC's flags) mishandles a by-reference actual argument that is an element of a device-resident array when the callee is a !$acc routine seq that does not get inlined. Reads through such an argument return garbage; writes through one are lost. Scalars and locals are correct. Index kind (literal or device-resident derived-type member) makes no difference. CCE OpenMP offload, amdflang and nvfortran are correct on the same source.

Found on #1811, where the new EOS helpers were called from kernels with q_prim_vf%vf(i)%sf(j,k,l) inputs and a blkmod1(k,l,q) output; every test on those paths ended in NaN(s) in timestep output on the Frontier CCE gpu-acc lanes only. An in-situ check recomputed the kernel's output on the host: device 0.0 in 300/300 cells, host 1.4 from the same fields. The parent PR passed because its only such callee was a leaf that CCE inlined at device link time; the rule holds for master too, which passes field elements into device routines at several sites that work only by inlining.

Standalone reproducer (three files, one kernel per argument kind; output below is from a Frontier login node):

sanity store:      bad     0 of   300
element, type index:  bad   300 of   300  0.000000E+00  2.438493E+00
element, const index: bad   300 of   300  0.000000E+00  2.438493E+00
elements in, scalar out: bad   300 of   300           NaN  2.438493E+00
scalars in, element out: bad   300 of   300  0.000000E+00  2.438493E+00
scalar:               bad     0 of   300  2.438493E+00  2.438493E+00

Workaround, applied in #1811 and recorded in .claude/rules/common-pitfalls.md: copy the element to a scalar before the call and receive outputs into a scalar. A lint rule could enforce it (flag %sf( or a GPU_DECLAREd array element inside the argument list of a call in a kernel or device routine); it would also flag the currently-working master sites, which rely on inlining.

build.sh
#!/bin/bash
F="-hacc -h acc_model=auto_async_none -h acc_model=no_fast_addr -O2"
ftn $F -c m_eos.f90 && ftn $F -c m_rhs.f90 && ftn $F main.f90 m_eos.o m_rhs.o -o repro &&  ./repro
m_eos.f90
module m_eos
  implicit none
  integer, parameter :: wp = kind(1.0d0)
  type scalar_field
    real(wp), pointer :: sf(:, :, :) => null()
  end type
  type vector_field
    type(scalar_field), allocatable :: vf(:)
  end type
  type idx_t
    integer :: e = 1, adv = 2, cont = 3
  end type
  type(idx_t) :: eqn
  real(wp), allocatable :: gam(:), pinf(:), rho0(:), c0(:), s(:)
  integer, allocatable :: eos(:)
  logical :: state_dependent = .false.
  !$acc declare create(eqn, gam, pinf, rho0, c0, s, eos, state_dependent)
contains
  subroutine reference_curve(rho, i, p_ref, e_ref, dp, de)
    !$acc routine seq
    real(wp), intent(in)  :: rho
    integer,  intent(in)  :: i
    real(wp), intent(out) :: p_ref, e_ref, dp, de
    real(wp) :: mu, d, up, us, dus
    integer :: it
    mu = rho/rho0(i) - 1.0_wp
    select case (eos(i))
    case (1)
      if (mu < 0.0_wp) then
        p_ref = rho0(i)*c0(i)**2*mu; dp = rho0(i)*c0(i)**2
      else
        up = c0(i)*mu/(1.0_wp - (s(i) - 1.0_wp)*mu)
        !$acc loop seq
        do it = 1, 8
          us = c0(i) + s(i)*up; dus = s(i)
          up = up - (us*mu - up*(1.0_wp + mu))/(dus*mu - (1.0_wp + mu))
        end do
        p_ref = rho0(i)*(c0(i) + s(i)*up)*up; dp = rho0(i)*c0(i)**2*exp(mu)
      end if
      e_ref = p_ref*mu/(2.0_wp*rho0(i)*(1.0_wp + mu)); de = dp*mu/rho0(i)
    case default
      p_ref = 0.0_wp; e_ref = 0.0_wp; dp = 0.0_wp; de = 0.0_wp
    end select
  end subroutine
  subroutine coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi)
    !$acc routine seq
    real(wp), intent(in)  :: alpha_rho, alpha
    integer,  intent(in)  :: i
    real(wp), intent(out) :: rho, gamma, pi_inf, dpi
    real(wp) :: p_ref, e_ref, dp, de
    rho = max(alpha_rho, 1.0e-16_wp)/max(alpha, 1.0e-16_wp)
    if (state_dependent) then
      call reference_curve(rho, i, p_ref, e_ref, dp, de)
      gamma = 1.0_wp/0.4_wp; pi_inf = rho*e_ref - p_ref/0.4_wp; dpi = e_ref + rho*de - dp/0.4_wp
    else
      gamma = gam(i); pi_inf = pinf(i); dpi = 0.0_wp
    end if
  end subroutine
  subroutine bulk_modulus(pres, alpha, alpha_rho, i, blkmod)
    !$acc routine seq
    real(wp), intent(in)  :: pres, alpha, alpha_rho
    integer,  intent(in)  :: i
    real(wp), intent(out) :: blkmod
    real(wp) :: rho, gamma, pi_inf, dpi
    call coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi)
    blkmod = ((gamma + 1.0_wp)*pres + pi_inf)/gamma - rho*dpi/gamma
  end subroutine
end module
m_rhs.f90
module m_rhs
  use m_eos
  implicit none
  real(wp), allocatable :: blkmod(:, :, :)
  !$acc declare create(blkmod)
contains
  subroutine by_element_idx(q, n)          ! MFC's form: field elements indexed through a device-resident type
    type(vector_field), intent(in) :: q
    integer, intent(in) :: n
    integer :: k, l, m
    !$acc parallel loop collapse(3) private(k, l, m)
    do m = 0, 0
      do l = 0, 0
        do k = 0, n
          call bulk_modulus(q%vf(eqn%e)%sf(k, l, m), q%vf(eqn%adv)%sf(k, l, m), q%vf(eqn%cont)%sf(k, l, m), 1, blkmod(k, l, m))
        end do
      end do
    end do
  end subroutine
  subroutine by_element_const(q, n)        ! same with literal indices
    type(vector_field), intent(in) :: q
    integer, intent(in) :: n
    integer :: k, l, m
    !$acc parallel loop collapse(3) private(k, l, m)
    do m = 0, 0
      do l = 0, 0
        do k = 0, n
          call bulk_modulus(q%vf(1)%sf(k, l, m), q%vf(2)%sf(k, l, m), q%vf(3)%sf(k, l, m), 1, blkmod(k, l, m))
        end do
      end do
    end do
  end subroutine
  subroutine elements_in_scalar_out(q, n)
    type(vector_field), intent(in) :: q
    integer, intent(in) :: n
    integer :: k, l, m
    real(wp) :: b
    !$acc parallel loop collapse(3) private(k, l, m, b)
    do m = 0, 0
      do l = 0, 0
        do k = 0, n
          call bulk_modulus(q%vf(1)%sf(k, l, m), q%vf(2)%sf(k, l, m), q%vf(3)%sf(k, l, m), 1, b)
          blkmod(k, l, m) = b
        end do
      end do
    end do
  end subroutine
  subroutine scalars_in_element_out(q, n)
    type(vector_field), intent(in) :: q
    integer, intent(in) :: n
    integer :: k, l, m
    real(wp) :: p, a, ar
    !$acc parallel loop collapse(3) private(k, l, m, p, a, ar)
    do m = 0, 0
      do l = 0, 0
        do k = 0, n
          p = q%vf(1)%sf(k, l, m)
          a = q%vf(2)%sf(k, l, m)
          ar = q%vf(3)%sf(k, l, m)
          call bulk_modulus(p, a, ar, 1, blkmod(k, l, m))
        end do
      end do
    end do
  end subroutine
  subroutine by_scalar(q, n)               ! the fix: scalars in, scalar out
    type(vector_field), intent(in) :: q
    integer, intent(in) :: n
    integer :: k, l, m
    real(wp) :: p, a, ar, b
    !$acc parallel loop collapse(3) private(k, l, m, p, a, ar, b)
    do m = 0, 0
      do l = 0, 0
        do k = 0, n
          p = q%vf(eqn%e)%sf(k, l, m)
          a = q%vf(eqn%adv)%sf(k, l, m)
          ar = q%vf(eqn%cont)%sf(k, l, m)
          call bulk_modulus(p, a, ar, 1, b)
          blkmod(k, l, m) = b
        end do
      end do
    end do
  end subroutine
end module
main.f90
program repro
  use m_rhs
  use openacc
  implicit none
  integer, parameter :: n = 299
  type(vector_field) :: q
  integer :: i, k
  real(wp) :: ref(0:n), rho, gamma, pi_inf, dpi
  allocate(gam(2), pinf(2), rho0(2), c0(2), s(2), eos(2))
  !$acc enter data create(gam, pinf, rho0, c0, s, eos)
  gam = 2.5_wp; pinf = 0.0_wp; rho0 = 1.0_wp; c0 = 1.0_wp; s = 1.5_wp; eos = 1
  state_dependent = .true.   ! take the deep branch, as an MG case does
  !$acc update device(eqn, gam, pinf, rho0, c0, s, eos, state_dependent)
  allocate(q%vf(3))
  do i = 1, 3
    allocate(q%vf(i)%sf(0:n, 0:0, 0:0))
    q%vf(i)%sf(:, 0, 0) = [(0.9_wp + 0.1_wp*i + 1.0e-3_wp*k, k = 0, n)]
  end do
  !$acc enter data copyin(q)
  !$acc enter data copyin(q%vf)
  do i = 1, 3
    !$acc enter data copyin(q%vf(i))
    !$acc enter data copyin(q%vf(i)%sf)
  end do
  allocate(blkmod(0:n, 0:0, 0:0))
  !$acc enter data create(blkmod)
  do k = 0, n
    call bulk_modulus(q%vf(1)%sf(k, 0, 0), q%vf(2)%sf(k, 0, 0), q%vf(3)%sf(k, 0, 0), 1, ref(k))
  end do

  print '(a,i3)', 'devices: ', acc_get_num_devices(acc_get_device_type())
  blkmod = -1.0_wp
  !$acc update device(blkmod)
  !$acc parallel loop
  do k = 0, n
    blkmod(k, 0, 0) = real(k, wp)
  end do
  !$acc update host(blkmod)
  print '(a,i5,a,i5)', 'sanity store:      bad ', count(abs(blkmod(:, 0, 0) - [(real(k, wp), k = 0, n)]) > 0.0_wp), ' of ', n + 1

  call run(by_element_idx,   'element, type index: ')
  call run(by_element_const, 'element, const index:')
  call run(elements_in_scalar_out, 'elements in, scalar out:')
  call run(scalars_in_element_out, 'scalars in, element out:')
  call run(by_scalar,        'scalar:              ')
contains
  subroutine run(kernel, label)
    interface
      subroutine kernel(q, n)
        import :: vector_field
        type(vector_field), intent(in) :: q
        integer, intent(in) :: n
      end subroutine
    end interface
    character(*), intent(in) :: label
    blkmod = 0.0_wp
    !$acc update device(blkmod)
    call kernel(q, n)
    !$acc update host(blkmod)
    print '(a,a,i5,a,i5,2es14.6)', label, ' bad ', count(.not. (abs(blkmod(:, 0, 0) - ref) <= 1e-12_wp*abs(ref))), ' of ', n + 1, blkmod(0, 0, 0), ref(0)
  end subroutine
end program

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions