Scalar indexing on GPU arrays in integration process

[huiyuxie]

This issue was found in #112. When running the examples, a scalar indexing issue was encountered during integration. Specifically, this can be traced back to the integrate function (refer to analysis_dg1d.jl, analysis_dg2d.jl, and analysis_dg3d.jl in the callbacks_step directory).

An MWE can be easily reproduced by setting CUDA.allowscalar(false) at the top of any example script.

Also, it can be further traced down to this example (updated):

using CUDA
CUDA.allowscalar(false)

using Trixi
using TrixiCUDA

# Set the precision
RealT = Float64

advection_velocity = 1.0
equations = LinearScalarAdvectionEquation1D(advection_velocity)

solver_gpu = DGSEMGPU(polydeg = 3, surface_flux = flux_lax_friedrichs, RealT = RealT)

coordinates_min = -1.0
coordinates_max = 1.0

mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level = 4,
                n_cells_max = 30_000)

semi_gpu = SemidiscretizationHyperbolicGPU(mesh, equations, initial_condition_convergence_test,
                                           solver_gpu)

tspan_gpu = (0.0, 1.0)
t_gpu = 0.0

# Semi on GPU
equations_gpu, mesh_gpu, solver_gpu = semi_gpu.equations, semi_gpu.mesh, semi_gpu.solver
cache_gpu, cache_cpu = semi_gpu.cache_gpu, semi_gpu.cache_cpu
boundary_conditions_gpu, source_terms_gpu = semi_gpu.boundary_conditions, semi_gpu.source_terms

# ODE on GPU
ode_gpu = semidiscretizeGPU(semi_gpu, tspan_gpu)
u_gpu_ = copy(ode_gpu.u0)
du_gpu_ = similar(u_gpu_)
u_gpu = TrixiCUDA.wrap_array(u_gpu_, mesh_gpu, equations_gpu, solver_gpu, cache_gpu)
du_gpu = TrixiCUDA.wrap_array(du_gpu_, mesh_gpu, equations_gpu, solver_gpu, cache_gpu)

# This is the root cause
Trixi.get_node_vars(u_gpu, equations_gpu, solver_gpu, 1, 1)

Currently, this issue is skipped by setting CUDA.allowscalar(true).

[LeeLizuoLiu]

Hi, it seems the integrate function in Trixi/src/solvers/dgsem/basis_lobatto_legendre.jl may also need to be rewritten. Below is the code

using CUDA, Trixi, TrixiCUDA
CUDA.allowscalar(false)

advection_velocity = 1.0
equations = LinearScalarAdvectionEquation1D(advection_velocity)

basis_gpu = LobattoLegendreBasisGPU(4) # Use GPU basis

solver = DGSEMGPU(basis_gpu, flux_lax_friedrichs) # Change to DGSEMGPU 

coordinates_min = -1.0
coordinates_max = 1.0

mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level = 4,
                n_cells_max = 30_000)

semi_gpu = SemidiscretizationHyperbolicGPU(mesh, equations, initial_condition_convergence_test,
                                           solver)

tspan_gpu = (0.0, 1.0)
t_gpu = 0.0

# Semi on GPU
equations_gpu = semi_gpu.equations
mesh_gpu, solver_gpu, cache_gpu = semi_gpu.mesh, semi_gpu.solver, semi_gpu.cache
boundary_conditions_gpu = semi_gpu.boundary_conditions
source_terms_gpu = semi_gpu.source_terms

# ODE on GPU
ode_gpu = semidiscretizeGPU(semi_gpu, tspan_gpu)
u_ode_gpu = copy(ode_gpu.u0)
du_ode_gpu = similar(u_ode_gpu)
u_gpu = Trixi.wrap_array(u_ode_gpu, mesh_gpu, equations_gpu, solver_gpu, cache_gpu)
du_gpu = Trixi.wrap_array(du_ode_gpu, mesh_gpu, equations_gpu, solver_gpu, cache_gpu)

# This is the root cause
Trixi.get_node_vars(u_gpu, equations_gpu, solver_gpu, 1, 1) 

with the same project.toml

name = "TrixiCUDA"
uuid = "5056661c-942a-472c-b3bc-01507f717a60"
authors = ["Huiyu Xie"]
version = "1.0.0-DEV"

[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"

[compat]
CUDA = "5.7"

[extras]
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Test"]

ERROR: 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/uiVyU/src/host/indexing.jl:50 [inlined]
  [6] first
    @ ./abstractarray.jl:452 [inlined]
  [7] integrate(f::ComposedFunction{…}, u::CuArray{…}, basis::LobattoLegendreBasis{…})
    @ Trixi ~/.julia/packages/Trixi/8wTPc/src/solvers/dgsem/basis_lobatto_legendre.jl:134
  [8] init_elements!(elements::Trixi.ElementContainer1D{…}, cell_ids::Vector{…}, mesh::TreeMesh{…}, basis::LobattoLegendreBasis{…})
    @ Trixi ~/.julia/packages/Trixi/8wTPc/src/solvers/dgsem_tree/containers_1d.jl:105
  [9] init_elements(cell_ids::Vector{…}, mesh::TreeMesh{…}, equations::LinearScalarAdvectionEquation1D{…}, basis::LobattoLegendreBasis{…}, ::Type{…}, ::Type{…})
    @ Trixi ~/.julia/packages/Trixi/8wTPc/src/solvers/dgsem_tree/containers_1d.jl:96
 [10] create_cache_gpu(mesh::TreeMesh{…}, equations::LinearScalarAdvectionEquation1D{…}, dg::DGSEM{…}, RealT::Type, uEltype::Type)
    @ TrixiCUDA ~/DGCFN.jl/TrixiCUDA.jl/src/solvers/cache.jl:19
 [11] SemidiscretizationHyperbolicGPU(mesh::TreeMesh{…}, equations::LinearScalarAdvectionEquation1D{…}, initial_condition::typeof(initial_condition_convergence_test), solver::DGSEM{…}; source_terms::Nothing, boundary_conditions::Trixi.BoundaryConditionPeriodic, RealT::Type, uEltype::Type, initial_cache_gpu::@NamedTuple{}, initial_cache_cpu::@NamedTuple{})
    @ TrixiCUDA ~/DGCFN.jl/TrixiCUDA.jl/src/semidiscretization/semidiscretization_hyperbolic.jl:54
 [12] SemidiscretizationHyperbolicGPU(mesh::TreeMesh{…}, equations::LinearScalarAdvectionEquation1D{…}, initial_condition::Function, solver::DGSEM{…})
    @ TrixiCUDA ~/DGCFN.jl/TrixiCUDA.jl/src/semidiscretization/semidiscretization_hyperbolic.jl:45
 [13] top-level scope
    @ ~/DGCFN.jl/TrixiCUDA.jl/debug_integrate.jl:18
Some type information was truncated. Use `show(err)` to see complete types.

[huiyuxie]

Exactly thanks for looking into this issue! I will resolve it as soon as possible.

[LeeLizuoLiu]

You are so welcome! Just take your time. I will also try to fix it.