add storage_type, real_type to semidiscretize

Author: vchuravy

examples/p4est_2d_dgsem/elixir_advection_basic.jl

@@ -31,7 +31,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergen
 # ODE solvers, callbacks etc.
 
 # Create ODE problem with time span from 0.0 to 1.0
-ode = semidiscretize(semi, (0.0, 1.0))
+ode = semidiscretize(semi, (0.0, 1.0); real_type = nothing, storage_type = nothing)
 
 # At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
 # and resets the timers

src/semidiscretization/semidiscretization.jl

@@ -82,16 +82,35 @@ end
 
 Wrap the semidiscretization `semi` as an ODE problem in the time interval `tspan`
 that can be passed to `solve` from the [SciML ecosystem](https://diffeq.sciml.ai/latest/).
+
+The optional keyword arguments `storage_type` and `real_type` configure the underlying computational
+datastructures. `storage_type` changes the fundamental array type being used, allowing the
+experimental use of `CuArray` or other GPU array types. `real_type` changes the computational data type being used.
 """
 function semidiscretize(semi::AbstractSemidiscretization, tspan;
-                        reset_threads = true)
+                        reset_threads = true,
+                        storage_type = nothing,
+                        real_type = nothing)
     # Optionally reset Polyester.jl threads. See
     # https://github.com/trixi-framework/Trixi.jl/issues/1583
     # https://github.com/JuliaSIMD/Polyester.jl/issues/30
     if reset_threads
         Polyester.reset_threads!()
     end
 
+    if !(storage_type === nothing && real_type === nothing)
+        if storage_type === nothing
+            storage_type = Array
+        end
+        if real_type === nothing
+            real_type = Float64
+        end
+        semi = trixi_adapt(storage_type, real_type, semi)
+        if eltype(tspan) !== real_type
+            tspan = convert.(real_type, tspan)
+        end
+    end
+
     u0_ode = compute_coefficients(first(tspan), semi)
     # TODO: MPI, do we want to synchronize loading and print debug statements, e.g. using
     #       mpi_isparallel() && MPI.Barrier(mpi_comm())

test/test_p4est_2d.jl

@@ -35,6 +35,27 @@ isdir(outdir) && rm(outdir, recursive = true)
     @test real(semi32.mesh) == Float64
 end
 
+@trixi_testset "elixir_advection_basic.jl (Float32)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"),
+                        # Expected errors are exactly the same as with TreeMesh!
+                        l2=[8.311947673061856e-6],
+                        linf=[6.627000273229378e-5],
+                        real_type=Float32)
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+    @test real(ode.p.solver) == Float32
+    @test real(ode.p.solver.basis) == Float32
+    @test real(ode.p.solver.mortar) == Float32
+    # TODO: remake ignores the mesh itself as well
+    @test real(ode.p.mesh) == Float64
+end
+
 @trixi_testset "elixir_advection_nonconforming_flag.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_advection_nonconforming_flag.jl"),