Stepsize callback returns unclear error on GPU when instability is detected

@amrueda discovered that some simulations were returning an unexpected error when running

`trixi_include("../Trixi.jl/examples/p4est_3d_dgsem/elixir_mhd_alfven_wave_combined_fluxes_nonperiodic.jl", polydeg = 6, storage_type= ROCArray)`

```julia
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Time integration │
│ ════════════════ │
│ Start time: ……………………………………………………………… 0.0 │
│ Final time: ……………………………………………………………… 1.0 │
│ time integrator: ………………………………………………… SSPRK43 │
│ adaptive: …………………………………………………………………… false │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Environment information │
│ ═══════════════════════ │
│ #threads: …………………………………………………………………… 1 │
│ threading backend: …………………………………………… polyester │
│ Backend: ……………………………………………………………………… KernelAbstractions AMDGPU │
│ AMDGPU devices: …………………………………………………… │
│ 1 (model): …………………………………………………………… AMD Instinct MI210 (gfx90a:sramecc+:xnack-) │
│ 1 (memory): ………………………………………………………… 61.070 GiB / 63.984 GiB available │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'IdealGlmMhdEquations3D' with DGSEM(polydeg=6)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 0 run time: 4.74000000e-07 s
Δt: 1.00000000e+00 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 0.00000000e+00 (0.000%) time/DOF/RHS: NaN s
PID: Inf s
#DOFs per field: 175616 alloc'd memory: 655.477 MiB
#elements: 512 device memory: 2984.250 MiB

Variable: rho rho_v1 rho_v2 rho_v3 rho_e_total B1 B2 B3 psi
L2 error: 3.32409319e-16 6.73218472e-08 3.36609236e-08 7.52681133e-08 7.47650072e-16 6.73218472e-08 3.36609236e-08 7.52681133e-08 0.00000000e+00
Linf error: 8.88178420e-16 1.49429967e-07 7.47149837e-08 1.67067783e-07 3.10862447e-15 1.49429968e-07 7.47149841e-08 1.67067783e-07 0.00000000e+00
∑S : 1.35975602e-16
────────────────────────────────────────────────────────────────────────────────────────────────────

#timesteps: 10 │ Δt: 6.8645e-03 │ sim. time: 6.8702e-02 (6.870%) │ run time: 5.1496e-01 s
#timesteps: 20 │ Δt: 6.8714e-03 │ sim. time: 1.3743e-01 (13.743%) │ run time: 8.3830e-01 s
#timesteps: 30 │ Δt: 6.8665e-03 │ sim. time: 2.0608e-01 (20.608%) │ run time: 1.1652e+00 s
ERROR: LoadError: Scalar indexing is disallowed.
Invocation of getindex resulted in scalar indexing of a GPU array.
This is typically caused by calling an iterating implementation of a method.
Such implementations *do not* execute on the GPU, but very slowly on the CPU,
and therefore should be avoided.

If you want to allow scalar iteration, use `allowscalar` or `@allowscalar`
to enable scalar iteration globally or for the operations in question.
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] errorscalar(op::String)
@ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:151
[3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
@ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:124
[4] assertscalar(op::String)
@ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:112
[5] getindex
@ ~/.julia/packages/GPUArrays/VNJ6B/src/host/indexing.jl:50 [inlined]
[6] log_numerical_instability(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}; jacobian_logging::Bool)
@ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/NVUZl/src/integrators/integrator_utils.jl:755
[7] log_numerical_instability
@ ~/.julia/packages/OrdinaryDiffEqCore/NVUZl/src/integrators/integrator_utils.jl:745 [inlined]
[8] check_error(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
@ SciMLBase ~/.julia/packages/SciMLBase/72tK7/src/integrator_interface.jl:1023
[9] check_error!
@ ~/.julia/packages/SciMLBase/72tK7/src/integrator_interface.jl:1069 [inlined]
[10] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
@ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/NVUZl/src/solve.jl:899
[11] #__solve#97
@ ~/.julia/packages/OrdinaryDiffEqCore/NVUZl/src/solve.jl:10 [inlined]
[12] __solve
@ ~/.julia/packages/OrdinaryDiffEqCore/NVUZl/src/solve.jl:1 [inlined]
[13] solve_call(_prob::ODEProblem{…}, args::SSPRK43{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/kdaLg/src/solve.jl:187
[14] solve_call
@ ~/.julia/packages/DiffEqBase/kdaLg/src/solve.jl:152 [inlined]
[15] #solve_up#40
@ ~/.julia/packages/DiffEqBase/kdaLg/src/solve.jl:661 [inlined]
[16] solve_up
@ ~/.julia/packages/DiffEqBase/kdaLg/src/solve.jl:634 [inlined]
[17] #solve#39
@ ~/.julia/packages/DiffEqBase/kdaLg/src/solve.jl:618 [inlined]
[18] top-level scope
@ ~/workspace/benchmarks/Trixi.jl/examples/p4est_3d_dgsem/elixir_mhd_alfven_wave_combined_fluxes_nonperiodic.jl:261
[19] include
@ ./Base.jl:496 [inlined]
[20] trixi_include(mapexpr::typeof(identity), mod::Module, elixir::String; enable_assignment_validation::Bool, replace_assignments_recursive::Bool, kwargs::@Kwargs{…})
@ TrixiBase ~/.julia/packages/TrixiBase/UX2bb/src/trixi_include.jl:84
[21] trixi_include
@ ~/.julia/packages/TrixiBase/UX2bb/src/trixi_include.jl:54 [inlined]
[22] #trixi_include#6
@ ~/.julia/packages/TrixiBase/UX2bb/src/trixi_include.jl:91 [inlined]
[23] trixi_include
@ ~/.julia/packages/TrixiBase/UX2bb/src/trixi_include.jl:90 [inlined]
[24] #trixi_include#7
@ ~/.julia/packages/TrixiBase/UX2bb/src/trixi_include.jl:95 [inlined]
[25] top-level scope
@ REPL[35]:1
in expression starting at /uni-mainz.de/homes/martiano/workspace/benchmarks/Trixi.jl/examples/p4est_3d_dgsem/elixir_mhd_alfven_wave_combined_fluxes_nonperiodic.jl:261
Some type information was truncated. Use `show(err)` to see complete types.
```
The reason seems to be the stepsize callback; the run indeed is fine when disabling it. The simulation crashes for stability issues and the error should be different.

```
────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished. Final time: 1.0 Time steps: 166 (accepted), 173 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────

─────────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 10.2s / 98.0% 539MiB / 98.3%

Section ncalls time %tot avg alloc %tot avg
─────────────────────────────────────────────────────────────────────────────────────────────────────────
rhs_hyperbolic! 692 4.82s 48.3% 6.96ms 43.8MiB 8.3% 64.7KiB
volume integral 692 2.69s 26.9% 3.88ms 8.94MiB 1.7% 13.2KiB
boundary flux 692 969ms 9.7% 1.40ms 8.53MiB 1.6% 12.6KiB
prolong2boundaries 692 545ms 5.5% 787μs 7.54MiB 1.4% 11.2KiB
prolong2interfaces + flux 692 366ms 3.7% 529μs 7.70MiB 1.5% 11.4KiB
surface integral + Jacobian + source terms 692 146ms 1.5% 211μs 7.68MiB 1.4% 11.4KiB
reset ∂u/∂t 692 98.1ms 1.0% 142μs 3.36MiB 0.6% 4.97KiB
~rhs_hyperbolic!~ 692 6.27ms 0.1% 9.06μs 5.88KiB 0.0% 8.69B
prolong2mortars 692 1.98ms 0.0% 2.87μs 0.00B 0.0% 0.00B
mortar flux 692 1.30ms 0.0% 1.88μs 0.00B 0.0% 0.00B
analyze solution 3 4.68s 46.9% 1.56s 341MiB 64.4% 114MiB
I/O 4 473ms 4.7% 118ms 145MiB 27.3% 36.2MiB
save solution 3 434ms 4.4% 145ms 145MiB 27.3% 48.3MiB
~I/O~ 4 38.3ms 0.4% 9.58ms 11.1KiB 0.0% 2.79KiB
get element variables 3 38.3μs 0.0% 12.8μs 3.94KiB 0.0% 1.31KiB
get node variables 3 2.05μs 0.0% 682ns 0.00B 0.0% 0.00B
save mesh 3 1.27μs 0.0% 424ns 0.00B 0.0% 0.00B
─────────────────────────────────────────────────────────────────────────────────────────────────────────
```

And reducing the polynomial degree, but reenabling the stepsize callback is also fine:
`trixi_include("../Trixi.jl/examples/p4est_3d_dgsem/elixir_mhd_alfven_wave_combined_fluxes_nonperiodic.jl", polydeg = 5, storage_type= ROCArray)`

```
────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished. Final time: 1.0 Time steps: 125 (accepted), 125 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────

─────────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 8.25s / 99.1% 450MiB / 99.0%

Section ncalls time %tot avg alloc %tot avg
─────────────────────────────────────────────────────────────────────────────────────────────────────────
analyze solution 3 4.62s 56.5% 1.54s 321MiB 72.1% 107MiB
rhs_hyperbolic! 500 2.93s 35.8% 5.86ms 31.8MiB 7.1% 65.0KiB
volume integral 500 1.65s 20.1% 3.29ms 6.40MiB 1.4% 13.1KiB
boundary flux 500 559ms 6.8% 1.12ms 6.00MiB 1.3% 12.3KiB
prolong2boundaries 500 316ms 3.9% 632μs 5.35MiB 1.2% 11.0KiB
prolong2interfaces + flux 500 259ms 3.2% 518μs 6.13MiB 1.4% 12.6KiB
surface integral + Jacobian + source terms 500 91.9ms 1.1% 184μs 5.44MiB 1.2% 11.1KiB
reset ∂u/∂t 500 52.9ms 0.6% 106μs 2.43MiB 0.5% 4.97KiB
~rhs_hyperbolic!~ 500 3.99ms 0.0% 7.98μs 5.88KiB 0.0% 12.0B
prolong2mortars 500 1.08ms 0.0% 2.16μs 0.00B 0.0% 0.00B
mortar flux 500 621μs 0.0% 1.24μs 0.00B 0.0% 0.00B
I/O 4 371ms 4.5% 92.7ms 91.2MiB 20.5% 22.8MiB
save solution 3 334ms 4.1% 111ms 91.2MiB 20.5% 30.4MiB
~I/O~ 4 37.0ms 0.5% 9.25ms 11.1KiB 0.0% 2.79KiB
get element variables 3 39.0μs 0.0% 13.0μs 4.59KiB 0.0% 1.53KiB
get node variables 3 9.42μs 0.0% 3.14μs 0.00B 0.0% 0.00B
save mesh 3 1.11μs 0.0% 370ns 0.00B 0.0% 0.00B
calculate dt 126 254ms 3.1% 2.01ms 1.40MiB 0.3% 11.4KiB
─────────────────────────────────────────────────────────────────────────────────────────────────────────
```

On the CPU indeed we have
```julia
#timesteps: 10 │ Δt: 6.8645e-03 │ sim. time: 6.8702e-02 (6.870%) │ run time: 4.0610e+00 s
#timesteps: 20 │ Δt: 6.8714e-03 │ sim. time: 1.3743e-01 (13.743%) │ run time: 7.8296e+00 s
#timesteps: 30 │ Δt: 6.8665e-03 │ sim. time: 2.0608e-01 (20.608%) │ run time: 1.1604e+01 s
┌ Warning: Verbosity toggle: instability
│ Instability detected. Aborting.
│
│ Diagnostics:
│
```

Start by reproducing the issue with examples/p4est_3d_dgsem/elixir_mhd_alfven_wave_combined_fluxes_nonperiodic.jl using polydeg = 6 and storage_type = ROCArray, then trace the stepsize callback through the log_numerical_instability and check_error stack entries. Done means GPU instability reports the intended stability error instead of a scalar-indexing error, while the polydeg = 5 and callback-disabled cases continue to work.

Written by the indexing model from the issue text.