GPU offloading

Project.toml

@@ -4,6 +4,7 @@ authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor
 version = "0.9.12-DEV"
 
 [deps]
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
@@ -17,6 +18,7 @@ FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
@@ -30,7 +32,6 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
-RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -64,21 +65,23 @@ TrixiMakieExt = "Makie"
 TrixiNLsolveExt = "NLsolve"
 
 [compat]
+Adapt = "3.7, 4.0"
 Accessors = "0.1.12"
 CodeTracking = "1.0.5"
 ConstructionBase = "1.3"
 Convex = "0.16"
 DataStructures = "0.18.15"
 DelimitedFiles = "1"
-DiffEqBase = "6 - 6.143"
-DiffEqCallbacks = "2.25"
+DiffEqBase = "6"
+DiffEqCallbacks = "2.25, 3, 4"
 Downloads = "1.6"
 ECOS = "1.1.2"
 EllipsisNotation = "1.0"
 FillArrays = "0.13.2, 1"
 ForwardDiff = "0.10.24"
 HDF5 = "0.16.10, 0.17"
 IfElse = "0.1"
+KernelAbstractions = "0.9"
 LinearAlgebra = "1"
 LinearMaps = "2.7, 3.0"
 LoopVectorization = "0.12.151"
@@ -94,15 +97,14 @@ PrecompileTools = "1.1"
 Preferences = "1.3"
 Printf = "1"
 RecipesBase = "1.1"
-RecursiveArrayTools = "2.38.10"
 Reexport = "1.0"
 Requires = "1.1"
 SciMLBase = "1.90, 2"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
 StableRNGs = "1.0.2"
 StartUpDG = "0.17.7, 1.1.5"
-Static = "0.8.7"
+Static = "0.8.7, 1"
 StaticArrayInterface = "1.4"
 StaticArrays = "1.5"
 StrideArrays = "0.1.26"

examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl

@@ -1,6 +1,7 @@
 
 using OrdinaryDiffEq
 using Trixi
+using CUDA
 
 ###############################################################################
 # semidiscretization of the compressible Euler equations
@@ -37,7 +38,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
 # ODE solvers, callbacks etc.
 
 tspan = (0.0, 5.0)
-ode = semidiscretize(semi, tspan)
+ode = semidiscretize(semi, tspan; adapt_to=CuArray)
 
 summary_callback = SummaryCallback()
 

examples/p4est_3d_dgsem/elixir_euler_taylor_green_vortex.jl

@@ -0,0 +1,79 @@
+using OrdinaryDiffEq
+using Trixi
+using CUDA
+CUDA.allowscalar(false)
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations3D(1.4)
+
+"""
+    initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D)
+
+The classical inviscid Taylor-Green vortex.
+"""
+function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D)
+  A  = 1.0 # magnitude of speed
+  Ms = 0.1 # maximum Mach number
+
+  rho = 1.0
+  v1  =  A * sin(x[1]) * cos(x[2]) * cos(x[3])
+  v2  = -A * cos(x[1]) * sin(x[2]) * cos(x[3])
+  v3  = 0.0
+  p   = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms
+  p   = p + 1.0/16.0 * A^2 * rho * (cos(2*x[1])*cos(2*x[3]) + 2*cos(2*x[2]) + 2*cos(2*x[1]) + cos(2*x[2])*cos(2*x[3]))
+
+  return prim2cons(SVector(rho, v1, v2, v3, p), equations)
+end
+
+initial_condition = initial_condition_taylor_green_vortex
+
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
+               volume_integral=VolumeIntegralFluxDifferencing(flux_lax_friedrichs))
+
+coordinates_min = (-1.0, -1.0, -1.0) .* pi
+coordinates_max = ( 1.0,  1.0,  1.0) .* pi
+
+# Create P4estMesh with 8 x 8 x 8 elements (note `refinement_level=1`)
+trees_per_dimension = (4, 4, 4)
+mesh = P4estMesh(trees_per_dimension, polydeg=1,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 initial_refinement_level=2)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)#5.0)
+ode = semidiscretize(semi, tspan; adapt_to=CuArray)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl=0.9)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, 
+                        alive_callback,
+                        #save_solution,
+                        stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

examples/structured_2d_dgsem/elixir_advection_coupled.jl

@@ -1,5 +1,6 @@
 using OrdinaryDiffEq
 using Trixi
+import Trixi.Indexing
 
 ###############################################################################
 # Coupled semidiscretization of four linear advection systems using converter functions such that
@@ -61,13 +62,13 @@ coupling_function12 = (x, u, equations_other, equations_own) -> u
 coupling_function13 = (x, u, equations_other, equations_own) -> u
 
 # Define the coupling boundary conditions and the system it is coupled to.
-boundary_conditions_x_neg1 = BoundaryConditionCoupled(2, (:end, :i_forward), Float64,
+boundary_conditions_x_neg1 = BoundaryConditionCoupled(2, (Indexing.last, Indexing.i_forward), Float64,
                                                       coupling_function12)
-boundary_conditions_x_pos1 = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64,
+boundary_conditions_x_pos1 = BoundaryConditionCoupled(2, (Indexing.first, Indexing.i_forward), Float64,
                                                       coupling_function12)
-boundary_conditions_y_neg1 = BoundaryConditionCoupled(3, (:i_forward, :end), Float64,
+boundary_conditions_y_neg1 = BoundaryConditionCoupled(3, (Indexing.i_forward, Indexing.last), Float64,
                                                       coupling_function13)
-boundary_conditions_y_pos1 = BoundaryConditionCoupled(3, (:i_forward, :begin), Float64,
+boundary_conditions_y_pos1 = BoundaryConditionCoupled(3, (Indexing.i_forward, Indexing.first), Float64,
                                                       coupling_function13)
 
 # A semidiscretization collects data structures and functions for the spatial discretization
@@ -93,13 +94,13 @@ coupling_function21 = (x, u, equations_other, equations_own) -> u
 coupling_function24 = (x, u, equations_other, equations_own) -> u
 
 # Define the coupling boundary conditions and the system it is coupled to.
-boundary_conditions_x_neg2 = BoundaryConditionCoupled(1, (:end, :i_forward), Float64,
+boundary_conditions_x_neg2 = BoundaryConditionCoupled(1, (Indexing.last, Indexing.i_forward), Float64,
                                                       coupling_function21)
-boundary_conditions_x_pos2 = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64,
+boundary_conditions_x_pos2 = BoundaryConditionCoupled(1, (Indexing.first, Indexing.i_forward), Float64,
                                                       coupling_function21)
-boundary_conditions_y_neg2 = BoundaryConditionCoupled(4, (:i_forward, :end), Float64,
+boundary_conditions_y_neg2 = BoundaryConditionCoupled(4, (Indexing.i_forward, Indexing.last), Float64,
                                                       coupling_function24)
-boundary_conditions_y_pos2 = BoundaryConditionCoupled(4, (:i_forward, :begin), Float64,
+boundary_conditions_y_pos2 = BoundaryConditionCoupled(4, (Indexing.i_forward, Indexing.first), Float64,
                                                       coupling_function24)
 
 # A semidiscretization collects data structures and functions for the spatial discretization
@@ -125,13 +126,13 @@ coupling_function34 = (x, u, equations_other, equations_own) -> u
 coupling_function31 = (x, u, equations_other, equations_own) -> u
 
 # Define the coupling boundary conditions and the system it is coupled to.
-boundary_conditions_x_neg3 = BoundaryConditionCoupled(4, (:end, :i_forward), Float64,
+boundary_conditions_x_neg3 = BoundaryConditionCoupled(4, (Indexing.last, Indexing.i_forward), Float64,
                                                       coupling_function34)
-boundary_conditions_x_pos3 = BoundaryConditionCoupled(4, (:begin, :i_forward), Float64,
+boundary_conditions_x_pos3 = BoundaryConditionCoupled(4, (Indexing.first, Indexing.i_forward), Float64,
                                                       coupling_function34)
-boundary_conditions_y_neg3 = BoundaryConditionCoupled(1, (:i_forward, :end), Float64,
+boundary_conditions_y_neg3 = BoundaryConditionCoupled(1, (Indexing.i_forward, Indexing.last), Float64,
                                                       coupling_function31)
-boundary_conditions_y_pos3 = BoundaryConditionCoupled(1, (:i_forward, :begin), Float64,
+boundary_conditions_y_pos3 = BoundaryConditionCoupled(1, (Indexing.i_forward, Indexing.first), Float64,
                                                       coupling_function31)
 
 # A semidiscretization collects data structures and functions for the spatial discretization
@@ -157,13 +158,13 @@ coupling_function43 = (x, u, equations_other, equations_own) -> u
 coupling_function42 = (x, u, equations_other, equations_own) -> u
 
 # Define the coupling boundary conditions and the system it is coupled to.
-boundary_conditions_x_neg4 = BoundaryConditionCoupled(3, (:end, :i_forward), Float64,
+boundary_conditions_x_neg4 = BoundaryConditionCoupled(3, (Indexing.last, Indexing.i_forward), Float64,
                                                       coupling_function43)
-boundary_conditions_x_pos4 = BoundaryConditionCoupled(3, (:begin, :i_forward), Float64,
+boundary_conditions_x_pos4 = BoundaryConditionCoupled(3, (Indexing.first, Indexing.i_forward), Float64,
                                                       coupling_function43)
-boundary_conditions_y_neg4 = BoundaryConditionCoupled(2, (:i_forward, :end), Float64,
+boundary_conditions_y_neg4 = BoundaryConditionCoupled(2, (Indexing.i_forward, Indexing.last), Float64,
                                                       coupling_function42)
-boundary_conditions_y_pos4 = BoundaryConditionCoupled(2, (:i_forward, :begin), Float64,
+boundary_conditions_y_pos4 = BoundaryConditionCoupled(2, (Indexing.i_forward, Indexing.first), Float64,
                                                       coupling_function42)
 
 # A semidiscretization collects data structures and functions for the spatial discretization

examples/structured_2d_dgsem/elixir_advection_meshview.jl

@@ -51,21 +51,21 @@ mesh2 = StructuredMeshView(parent_mesh; indices_min = (9, 1), indices_max = (16,
 coupling_function = (x, u, equations_other, equations_own) -> u
 
 # Define the coupled boundary conditions
-# The indices (:end, :i_forward) and (:begin, :i_forward) denote the interface indexing.
+# The indices (Indexing.last, Indexing.i_forward) and (Indexing.first, Indexing.i_forward) denote the interface indexing.
 # For a system with coupling in x and y see examples/structured_2d_dgsem/elixir_advection_coupled.jl.
 boundary_conditions1 = (
                         # Connect left boundary with right boundary of left mesh
-                        x_neg = BoundaryConditionCoupled(2, (:end, :i_forward), Float64,
+                        x_neg = BoundaryConditionCoupled(2, (Indexing.last, Indexing.i_forward), Float64,
                                                          coupling_function),
-                        x_pos = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64,
+                        x_pos = BoundaryConditionCoupled(2, (Indexing.first, Indexing.i_forward), Float64,
                                                          coupling_function),
                         y_neg = boundary_condition_periodic,
                         y_pos = boundary_condition_periodic)
 boundary_conditions2 = (
                         # Connect left boundary with right boundary of left mesh
-                        x_neg = BoundaryConditionCoupled(1, (:end, :i_forward), Float64,
+                        x_neg = BoundaryConditionCoupled(1, (Indexing.last, Indexing.i_forward), Float64,
                                                          coupling_function),
-                        x_pos = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64,
+                        x_pos = BoundaryConditionCoupled(1, (Indexing.first, Indexing.i_forward), Float64,
                                                          coupling_function),
                         y_neg = boundary_condition_periodic,
                         y_pos = boundary_condition_periodic)

examples/structured_2d_dgsem/elixir_mhd_coupled.jl

@@ -49,9 +49,9 @@ mesh1 = StructuredMesh(cells_per_dimension,
                        periodicity = (false, true))
 
 coupling_function1 = (x, u, equations_other, equations_own) -> u
-boundary_conditions1 = (x_neg = BoundaryConditionCoupled(2, (:end, :i_forward), Float64,
+boundary_conditions1 = (x_neg = BoundaryConditionCoupled(2, (Indexing.last, Indexing.i_forward), Float64,
                                                          coupling_function1),
-                        x_pos = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64,
+                        x_pos = BoundaryConditionCoupled(2, (Indexing.first, Indexing.i_forward), Float64,
                                                          coupling_function1),
                         y_neg = boundary_condition_periodic,
                         y_pos = boundary_condition_periodic)
@@ -72,9 +72,9 @@ mesh2 = StructuredMesh(cells_per_dimension,
                        periodicity = (false, true))
 
 coupling_function2 = (x, u, equations_other, equations_own) -> u
-boundary_conditions2 = (x_neg = BoundaryConditionCoupled(1, (:end, :i_forward), Float64,
+boundary_conditions2 = (x_neg = BoundaryConditionCoupled(1, (Indexing.last, Indexing.i_forward), Float64,
                                                          coupling_function2),
-                        x_pos = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64,
+                        x_pos = BoundaryConditionCoupled(1, (Indexing.first, Indexing.i_forward), Float64,
                                                          coupling_function2),
                         y_neg = boundary_condition_periodic,
                         y_pos = boundary_condition_periodic)

examples/t8code_3d_dgsem/elixir_advection_cubed_sphere.jl

@@ -0,0 +1,61 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.2, -0.7, 0.5)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+initial_condition = initial_condition_convergence_test
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict(:inside => boundary_condition,
+                           :outside => boundary_condition)
+
+# Note that the first argument refers to the level of refinement, unlike in for p4est
+mesh = Trixi.T8codeMeshCubedSphere(2, 3, 0.5, 0.5;
+                                   polydeg = 3, initial_refinement_level = 0)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval = 100)
+
+# # The SaveSolutionCallback allows to save the solution to a file in regular intervals
+# save_solution = SaveSolutionCallback(interval = 100,
+#                                      solution_variables = cons2prim)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 1.2)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, # save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()