amdflang: array constructor assignment to private arrays in GPU_PARALLEL_LOOP produces wrong values (m_ibm.fpp, m_collisions.fpp)

[sbryngelson]

Summary

Array constructor assignments (var = [a, b, c]) to private local arrays inside GPU_PARALLEL_LOOP regions produce silently wrong values on AMD GPU when compiled with amdflang using both -fopenmp-assume-threads-oversubscription and -fopenmp-assume-teams-oversubscription. This is the flag combination MFC uses for all AMD GPU builds.

The bug is present in all AFAR drops we tested: 23.1.0, 23.2.0, 23.2.1 (gfx90a / MI250X).

Errors are large — O(10–100) in absolute value, not floating-point noise — affecting ~50% of loop iterations non-deterministically.

Trigger Conditions

The bug fires when all three conditions hold:

1. -fopenmp-assume-threads-oversubscription and -fopenmp-assume-teams-oversubscription are both present (either flag alone does not trigger it)
2. An array constructor [expr1, expr2, ...] is assigned to a private array inside a target teams distribute parallel do
3. The loop has enough iterations (≥ ~64)

Affected Code in MFC

src/simulation/m_ibm.fpp

s_ibm_correct_state (line ~195 GPU_PARALLEL_LOOP, lines 207/209):

real(wp), dimension(3) :: physical_loc   ! declared private
...
$:GPU_PARALLEL_LOOP(private='[..., physical_loc, ...]')
do i = 1, num_gps
    ...
    if (p > 0) then
        physical_loc = [x_cc(j), y_cc(k), z_cc(l)]   ! BUG: wrong values ~50% of iterations
    else
        physical_loc = [x_cc(j), y_cc(k), 0._wp]      ! BUG: same
    end if
    ...
end do

This subroutine is called every time step for all IBM cases. The wrong physical_loc values mean ghost point corrections are computed using corrupted spatial coordinates.

s_find_ghost_points (line ~397 GPU_PARALLEL_LOOP, lines 407/409): same pattern.

src/simulation/m_collisions.fpp

s_detect_ib_collisions_n2 (lines 337–338):

centroid_1 = [patch_ib(pid1)%x_centroid, patch_ib(pid1)%y_centroid, 0._wp]
centroid_2 = [patch_ib(pid2)%x_centroid, patch_ib(pid2)%y_centroid, 0._wp]

Inside a GPU_PARALLEL_LOOP. Only reachable when collision_model > 0 (not currently in CI).

Why CI Passes

The affected code paths (s_ibm_correct_state, s_find_ghost_points) are exercised by the GPU CI IBM test suite. However, the current regression tests appear to tolerate the corrupted ghost point coordinates — likely because the test geometries and flow conditions happen to be insensitive enough that wrong physical_loc values don't push the final solution norm outside the 1e-10 tolerance. This is a silent correctness bug.

Minimal Reproducer

program minimal_array_constructor
    implicit none
    integer, parameter :: N = 64
    real(8) :: x(N), y(N), z(N), out(N,3)
    real(8) :: loc(3)
    integer :: i, nerr

    do i = 1, N
        x(i) = real(i,8)*0.1d0; y(i) = real(i,8)*0.2d0; z(i) = real(i,8)*0.3d0
    end do

    !$omp target teams distribute parallel do map(to:x,y,z) map(from:out) private(loc)
    do i = 1, N
        loc = [x(i), y(i), z(i)]
        out(i,1) = loc(1); out(i,2) = loc(2); out(i,3) = loc(3)
    end do
    !$omp end target teams distribute parallel do

    nerr = 0
    do i = 1, N
        if (abs(out(i,1)-x(i)) > 1d-14 .or. &
            abs(out(i,2)-y(i)) > 1d-14 .or. &
            abs(out(i,3)-z(i)) > 1d-14) nerr = nerr + 1
    end do
    if (nerr == 0) then
        print *, "PASS"
    else
        print *, "FAIL:", nerr, "of", N, "(max abs error O(10-100))"
    end if
end program

Compile and run:

amdflang -fopenmp --offload-arch=gfx90a -fopenmp-target-fast \
  -fopenmp-assume-threads-oversubscription \
  -fopenmp-assume-teams-oversubscription \
  -o minimal minimal_array_constructor.f90
./minimal
# FAIL: 31 of 64 (max abs error O(10-100))

Workaround

Replace array constructor assignments with explicit element-by-element assignment:

! Instead of:
physical_loc = [x_cc(j), y_cc(k), z_cc(l)]

! Use:
physical_loc(1) = x_cc(j)
physical_loc(2) = y_cc(k)
physical_loc(3) = z_cc(l)

This is confirmed to produce correct results.

Upstream Compiler Bug

Being investigated in ROCm/llvm-project. A related fix for private VLA arrays is in ROCm/llvm-project#2422. The array constructor bug is a separate issue — the upstream report will follow.

cc @danieljvickers

[sbryngelson]

Updated scope — additional affected files and critical Riemann solver impact

Deeper analysis of the root cause has revealed two more affected files beyond those in the original report, and clarified that the bug is more severe than initially characterized.

What the bug actually does (corrected understanding)

This is not a wrong-value bug — it is a missing-iteration bug. The affected iterations simply never execute; output cells for those iterations remain at their initial values (typically zero or uninitialized device memory). The failure is completely deterministic:

N_wrong = max(0,  N − (⌈N/32⌉ + 31))

For a loop of N=10000 cells: 9969 iterations never execute. For N=64: 31 skipped. Threshold is N>33.

Additional affected files

Beyond m_ibm.fpp and m_collisions.fpp previously reported, two more files contain array constructors inside GPU_PARALLEL_LOOP:

src/simulation/m_ib_patches.fpp — 8 patterns (NEW):

• xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] (~lines 234, 481, 719, 794)
• xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] (~lines 394, 881)
• xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center (~lines 600, 662)

All inside GPU_PARALLEL_LOOP over the computational domain. For any 3D simulation, N >> 33 and most cells will have wrong IB patch geometry.

src/simulation/m_riemann_solvers.fpp — 13 patterns (NEW and most critical):

Multispecies loops (three separate GPU_PARALLEL_LOOPs, ~lines 338, 1023, 2919):

Xs_L(:) = Ys_L(:)*MW_L/molecular_weights(:)
Xs_R(:) = Ys_R(:)*MW_R/molecular_weights(:)

HLLD MHD solver (~line 3459):

U_L     = [rho%L, rho%L*vel%L(1:3), B%L(2:3), E%L]
U_starL = [rhoL_star, rhoL_star*s_M, rhoL_star*vel%L(2:3), B%L(2:3), E_starL]
U_R     = [rho%R, rho%R*vel%R(1:3), B%R(2:3), E%R]
U_starR = [rhoR_star, rhoR_star*s_M, rhoR_star*vel%R(2:3), B%R(2:3), E_starR]
U_doubleL = [rhoL_star, rhoL_star*s_M, ...]
U_doubleR = [rhoR_star, rhoR_star*s_M, ...]

In production MFC runs over even modest grids (128³ cells → N ≈ millions), all of these loops will silently skip the vast majority of iterations, producing zero flux vectors for nearly all cells.

Why existing tests didn't catch this

The HLLD tests (test10, test11) in the original test suite ran with N=1 (single interface). By the formula, N=1 always gives 0 wrong values and passes. Production runs use N >> 33 and are fully exposed.

Upstream compiler fix

A fix has been filed at the compiler level:

• ROCm/llvm-project#2602 — one-line DeviceRTL fix in Workshare.cpp
• llvm/llvm-project#198621 — upstream tracking issue

Until that fix is deployed, all array constructor and whole-array slice assignments in GPU_PARALLEL_LOOP must be replaced with explicit element-by-element assignment.

Workaround pattern

! Instead of:
U_L = [rho%L, rho%L*vel%L(1:3), B%L(2:3), E%L]

! Use:
U_L(1)   = rho%L
U_L(2:4) = rho%L*vel%L(1:3)
U_L(5:6) = B%L(2:3)
U_L(7)   = E%L

[sbryngelson]

Why CI tests still pass

Three mechanisms interact to hide this bug from the regression suite.

---

1. Golden file comparison (primary reason)

MFC's test framework is kind: "golden" — each test compares its output against a pre-computed reference file. If those reference files were generated on AMD GPU hardware with the same buggy AFAR compiler (23.1.x or later), both the reference and the current test run have the same iterations silently skipped and produce the same wrong output. The comparison agrees and the test passes — not because the physics is correct, but because both runs are wrong in exactly the same way.

This is the most likely explanation for why IBM tests pass even with num_gps ≫ 33. The bug is consistent and deterministic (always the same iterations, always zero), so golden-file comparison will never catch it.

---

2. Grid sizes for IBM tests — right on the boundary

From toolchain/mfc/test/cases.py, the IBM test grid sizes are:

Dimension	Grid	Total cells	IB radius	Estimated num_gps	Bug fires?
1D	m=299	299	0.1	~2	No (N ≤ 33)
2D	m=49, n=39	1,911	0.1	~50–60	Yes (N_wrong ≈ 20–27)
3D	m=24³	13,824	0.1	~100–150	Yes (N_wrong ≈ 70–120)

The 1D case is safe. The 2D and 3D cases have num_gps above the threshold of 33, but the golden-file mechanism (above) prevents the tests from failing.

---

3. Multispecies and HLLD code paths not tested on AMD GPU, or covered by golden files

From cases.py, the MHD/HLLD tests use 2D grids of 50×50 = 2,500 cells. The formula gives N_wrong = 2,500 − (79 + 31) = 2,390 — the vast majority of cells would have zero fluxes. If these tests run on Frontier AMD at all, they are passing solely because the golden files were produced with the same broken compiler. If they run only on CPU (OpenMP host fallback), they are not exercising the GPU path.

---

What this means in practice

The MFC regression suite provides no signal that the physics is correct for AMD GPU runs involving:

• IBM with non-trivial geometries (any grid where num_gps > 33)
• Any multispecies or MHD simulation
• Any IB patch geometry computation

The tests pass, but they are comparing wrong-against-wrong. A simulation run on AMD GPU in production will produce silently incorrect results unless the golden files happen to have been generated with the fix in place.

The only way to verify correctness is to run the same case on CPU (where the bug does not exist) and compare against the AMD GPU result — they will differ significantly.