Skip to content

Stepsize callback returns unclear error on GPU when instability is detected #3182

Description

@MarcoArtiano

@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)

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ 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

#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:

Metadata

Metadata

Assignees

No one assigned

    Labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions