Adapt DGSEM to enable caching of its data on GPU (#117)

* Start * Complete * Complete 1D * Fix * Fix * Complete 2D * Complete 3D * Test examples
trixi-gpu · Jan 20, 2025 · c64665d · c64665d
1 parent 3336c7c
commit c64665d
Show file tree

Hide file tree

Showing 84 changed files with 501 additions and 385 deletions.
diff --git a/README.md b/README.md
@@ -78,7 +78,8 @@ using OrdinaryDiffEq
 
 # Currently skip the issue of scalar indexing
 # See issues https://github.com/trixi-gpu/TrixiCUDA.jl/issues/59
-# and https://github.com/trixi-gpu/TrixiCUDA.jl/issues/113
+# https://github.com/trixi-gpu/TrixiCUDA.jl/issues/113
+# https://github.com/trixi-gpu/TrixiCUDA.jl/issues/118
 using CUDA
 CUDA.allowscalar(true)
 
@@ -88,7 +89,7 @@ CUDA.allowscalar(true)
 advection_velocity = 1.0
 equations = LinearScalarAdvectionEquation1D(advection_velocity)
 
-solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+solver = DGSEMGPU(polydeg = 3, surface_flux = flux_lax_friedrichs)
 
 coordinates_min = -1.0
 coordinates_max = 1.0 

diff --git a/examples/advection_basic_1d.jl b/examples/advection_basic_1d.jl
@@ -14,7 +14,7 @@ advection_velocity = 1.0
 equations = LinearScalarAdvectionEquation1D(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)
+solver = DGSEMGPU(polydeg = 3, surface_flux = flux_lax_friedrichs)
 
 coordinates_min = -1.0 # minimum coordinate
 coordinates_max = 1.0 # maximum coordinate

diff --git a/examples/advection_basic_2d.jl b/examples/advection_basic_2d.jl
@@ -14,7 +14,7 @@ advection_velocity = (0.2, -0.7)
 equations = LinearScalarAdvectionEquation2D(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)
+solver = DGSEMGPU(polydeg = 3, surface_flux = flux_lax_friedrichs)
 
 coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
 coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))

diff --git a/examples/advection_basic_3d.jl b/examples/advection_basic_3d.jl
@@ -14,7 +14,7 @@ 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)
+solver = DGSEMGPU(polydeg = 3, surface_flux = flux_lax_friedrichs)
 
 coordinates_min = (-1.0, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z))
 coordinates_max = (1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z))

diff --git a/examples/advection_mortar_2d.jl b/examples/advection_mortar_2d.jl
@@ -14,7 +14,7 @@ advection_velocity = (0.2, -0.7)
 equations = LinearScalarAdvectionEquation2D(advection_velocity)
 
 initial_condition = initial_condition_convergence_test
-solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+solver = DGSEMGPU(polydeg = 3, surface_flux = flux_lax_friedrichs)
 
 coordinates_min = (-1.0, -1.0)
 coordinates_max = (1.0, 1.0)

diff --git a/examples/advection_mortar_3d.jl b/examples/advection_mortar_3d.jl
@@ -15,7 +15,7 @@ advection_velocity = (0.2, -0.7, 0.5)
 equations = LinearScalarAdvectionEquation3D(advection_velocity)
 
 initial_condition = initial_condition_convergence_test
-solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+solver = DGSEMGPU(polydeg = 3, surface_flux = flux_lax_friedrichs)
 
 coordinates_min = (-1.0, -1.0, -1.0)
 coordinates_max = (1.0, 1.0, 1.0)

diff --git a/examples/euler_ec_1d.jl b/examples/euler_ec_1d.jl
@@ -14,8 +14,8 @@ equations = CompressibleEulerEquations1D(1.4)
 initial_condition = initial_condition_weak_blast_wave
 
 volume_flux = flux_ranocha
-solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha,
-               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+solver = DGSEMGPU(polydeg = 3, surface_flux = flux_ranocha,
+                  volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (-2.0,)
 coordinates_max = (2.0,)

diff --git a/examples/euler_ec_2d.jl b/examples/euler_ec_2d.jl
@@ -14,8 +14,8 @@ equations = CompressibleEulerEquations2D(1.4)
 initial_condition = initial_condition_weak_blast_wave
 
 volume_flux = flux_ranocha
-solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha,
-               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+solver = DGSEMGPU(polydeg = 3, surface_flux = flux_ranocha,
+                  volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (-2.0, -2.0)
 coordinates_max = (2.0, 2.0)

diff --git a/examples/euler_ec_3d.jl b/examples/euler_ec_3d.jl
@@ -15,8 +15,8 @@ equations = CompressibleEulerEquations3D(1.4)
 initial_condition = initial_condition_weak_blast_wave
 
 volume_flux = flux_ranocha
-solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha,
-               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+solver = DGSEMGPU(polydeg = 3, surface_flux = flux_ranocha,
+                  volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (-2.0, -2.0, -2.0)
 coordinates_max = (2.0, 2.0, 2.0)

diff --git a/examples/euler_shockcapturing_1d.jl b/examples/euler_shockcapturing_1d.jl
@@ -16,7 +16,7 @@ initial_condition = initial_condition_weak_blast_wave
 
 surface_flux = flux_lax_friedrichs
 volume_flux = flux_shima_etal
-basis = LobattoLegendreBasis(3)
+basis = LobattoLegendreBasisGPU(3)
 indicator_sc = IndicatorHennemannGassner(equations, basis,
                                          alpha_max = 0.5,
                                          alpha_min = 0.001,
@@ -25,7 +25,7 @@ indicator_sc = IndicatorHennemannGassner(equations, basis,
 volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                  volume_flux_dg = volume_flux,
                                                  volume_flux_fv = surface_flux)
-solver = DGSEM(basis, surface_flux, volume_integral)
+solver = DGSEMGPU(basis, surface_flux, volume_integral)
 
 coordinates_min = -2.0
 coordinates_max = 2.0

diff --git a/examples/euler_shockcapturing_2d.jl b/examples/euler_shockcapturing_2d.jl
@@ -16,7 +16,7 @@ initial_condition = initial_condition_weak_blast_wave
 
 surface_flux = flux_lax_friedrichs
 volume_flux = flux_shima_etal
-basis = LobattoLegendreBasis(3)
+basis = LobattoLegendreBasisGPU(3)
 indicator_sc = IndicatorHennemannGassner(equations, basis,
                                          alpha_max = 0.5,
                                          alpha_min = 0.001,
@@ -25,7 +25,7 @@ indicator_sc = IndicatorHennemannGassner(equations, basis,
 volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                  volume_flux_dg = volume_flux,
                                                  volume_flux_fv = surface_flux)
-solver = DGSEM(basis, surface_flux, volume_integral)
+solver = DGSEMGPU(basis, surface_flux, volume_integral)
 
 coordinates_min = (-2.0, -2.0)
 coordinates_max = (2.0, 2.0)

diff --git a/examples/euler_shockcapturing_3d.jl b/examples/euler_shockcapturing_3d.jl
@@ -18,7 +18,7 @@ surface_flux = flux_ranocha # OBS! Using a non-dissipative flux is only sensible
 # but not for real shock simulations
 volume_flux = flux_ranocha
 polydeg = 3
-basis = LobattoLegendreBasis(polydeg)
+basis = LobattoLegendreBasisGPU(polydeg)
 indicator_sc = IndicatorHennemannGassner(equations, basis,
                                          alpha_max = 0.5,
                                          alpha_min = 0.001,
@@ -27,7 +27,7 @@ indicator_sc = IndicatorHennemannGassner(equations, basis,
 volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                  volume_flux_dg = volume_flux,
                                                  volume_flux_fv = surface_flux)
-solver = DGSEM(basis, surface_flux, volume_integral)
+solver = DGSEMGPU(basis, surface_flux, volume_integral)
 
 coordinates_min = (-2.0, -2.0, -2.0)
 coordinates_max = (2.0, 2.0, 2.0)

diff --git a/examples/euler_source_terms_1d.jl b/examples/euler_source_terms_1d.jl
@@ -16,7 +16,7 @@ initial_condition = initial_condition_convergence_test
 
 # Note that the expected EOC of 5 is not reached with this flux.
 # Using flux_hll instead yields the expected EOC.
-solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs)
+solver = DGSEMGPU(polydeg = 4, surface_flux = flux_lax_friedrichs)
 
 coordinates_min = 0.0
 coordinates_max = 2.0

diff --git a/examples/euler_source_terms_2d.jl b/examples/euler_source_terms_2d.jl
@@ -13,7 +13,7 @@ CUDA.allowscalar(true)
 equations = CompressibleEulerEquations2D(1.4)
 
 initial_condition = initial_condition_convergence_test
-solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+solver = DGSEMGPU(polydeg = 3, surface_flux = flux_lax_friedrichs)
 
 coordinates_min = (0.0, 0.0)
 coordinates_max = (2.0, 2.0)

diff --git a/examples/euler_source_terms_3d.jl b/examples/euler_source_terms_3d.jl
@@ -14,8 +14,8 @@ equations = CompressibleEulerEquations3D(1.4)
 
 initial_condition = initial_condition_convergence_test
 
-solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs,
-               volume_integral = VolumeIntegralWeakForm())
+solver = DGSEMGPU(polydeg = 3, surface_flux = flux_lax_friedrichs,
+                  volume_integral = VolumeIntegralWeakForm())
 
 coordinates_min = (0.0, 0.0, 0.0)
 coordinates_max = (2.0, 2.0, 2.0)

diff --git a/examples/eulermulti_ec_1d.jl b/examples/eulermulti_ec_1d.jl
@@ -20,8 +20,8 @@ equations = CompressibleEulerMulticomponentEquations1D(gammas = (1.4, 1.4),
 initial_condition = initial_condition_weak_blast_wave
 
 volume_flux = flux_ranocha
-solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha,
-               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+solver = DGSEMGPU(polydeg = 3, surface_flux = flux_ranocha,
+                  volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (-2.0,)
 coordinates_max = (2.0,)

diff --git a/examples/eulermulti_ec_2d.jl b/examples/eulermulti_ec_2d.jl
@@ -15,8 +15,8 @@ equations = CompressibleEulerMulticomponentEquations2D(gammas = 1.4,
 initial_condition = initial_condition_weak_blast_wave
 
 volume_flux = flux_ranocha
-solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha,
-               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+solver = DGSEMGPU(polydeg = 3, surface_flux = flux_ranocha,
+                  volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (-2.0, -2.0)
 coordinates_max = (2.0, 2.0)

diff --git a/examples/hypdiff_nonperiodic_1d.jl b/examples/hypdiff_nonperiodic_1d.jl
@@ -16,7 +16,7 @@ initial_condition = initial_condition_poisson_nonperiodic
 
 boundary_conditions = boundary_condition_poisson_nonperiodic
 
-solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs)
+solver = DGSEMGPU(polydeg = 4, surface_flux = flux_lax_friedrichs)
 
 coordinates_min = 0.0
 coordinates_max = 1.0

diff --git a/examples/hypdiff_nonperiodic_2d.jl b/examples/hypdiff_nonperiodic_2d.jl
@@ -19,7 +19,7 @@ boundary_conditions = (x_neg = boundary_condition_poisson_nonperiodic,
                        y_neg = boundary_condition_periodic,
                        y_pos = boundary_condition_periodic)
 
-solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs)
+solver = DGSEMGPU(polydeg = 4, surface_flux = flux_lax_friedrichs)
 
 coordinates_min = (0.0, 0.0)
 coordinates_max = (1.0, 1.0)

diff --git a/examples/hypdiff_nonperiodic_3d.jl b/examples/hypdiff_nonperiodic_3d.jl
@@ -20,7 +20,7 @@ boundary_conditions = (x_neg = boundary_condition_poisson_nonperiodic,
                        z_neg = boundary_condition_periodic,
                        z_pos = boundary_condition_periodic)
 
-solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs)
+solver = DGSEMGPU(polydeg = 4, surface_flux = flux_lax_friedrichs)
 
 coordinates_min = (0.0, 0.0, 0.0)
 coordinates_max = (1.0, 1.0, 1.0)

diff --git a/examples/shallowwater_dirichlet_1d.jl b/examples/shallowwater_dirichlet_1d.jl
@@ -20,10 +20,10 @@ boundary_condition = BoundaryConditionDirichlet(initial_condition)
 # Get the DG approximation space
 
 volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-solver = DGSEM(polydeg = 4,
-               surface_flux = (flux_hll,
-                               flux_nonconservative_fjordholm_etal),
-               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+solver = DGSEMGPU(polydeg = 4,
+                  surface_flux = (flux_hll,
+                                  flux_nonconservative_fjordholm_etal),
+                  volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 ###############################################################################
 # Get the TreeMesh and setup a periodic mesh

diff --git a/examples/shallowwater_dirichlet_2d.jl b/examples/shallowwater_dirichlet_2d.jl
@@ -20,9 +20,9 @@ boundary_condition = BoundaryConditionDirichlet(initial_condition)
 # Get the DG approximation space
 
 volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-solver = DGSEM(polydeg = 3,
-               surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal),
-               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+solver = DGSEMGPU(polydeg = 3,
+                  surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal),
+                  volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 ###############################################################################
 # Get the TreeMesh and setup a periodic mesh

diff --git a/examples/shallowwater_ec_1d.jl b/examples/shallowwater_ec_1d.jl
@@ -49,9 +49,9 @@ initial_condition = initial_condition_ec_discontinuous_bottom
 # Get the DG approximation space
 
 volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-solver = DGSEM(polydeg = 4,
-               surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal),
-               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+solver = DGSEMGPU(polydeg = 4,
+                  surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal),
+                  volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 ###############################################################################
 # Get the TreeMesh and setup a periodic mesh

diff --git a/examples/shallowwater_ec_2d.jl b/examples/shallowwater_ec_2d.jl
@@ -21,9 +21,9 @@ initial_condition = initial_condition_weak_blast_wave
 # Get the DG approximation space
 
 volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-solver = DGSEM(polydeg = 4,
-               surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal),
-               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+solver = DGSEMGPU(polydeg = 4,
+                  surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal),
+                  volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 ###############################################################################
 # Get the TreeMesh and setup a periodic mesh

diff --git a/examples/shallowwater_source_terms_1d.jl b/examples/shallowwater_source_terms_1d.jl
@@ -18,9 +18,9 @@ initial_condition = initial_condition_convergence_test
 # Get the DG approximation space
 
 volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-solver = DGSEM(polydeg = 3,
-               surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal),
-               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+solver = DGSEMGPU(polydeg = 3,
+                  surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal),
+                  volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 ###############################################################################
 # Get the TreeMesh and setup a periodic mesh

diff --git a/examples/shallowwater_source_terms_2d.jl b/examples/shallowwater_source_terms_2d.jl
@@ -18,9 +18,9 @@ initial_condition = initial_condition_convergence_test # MMS EOC test
 # Get the DG approximation space
 
 volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-solver = DGSEM(polydeg = 3,
-               surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal),
-               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+solver = DGSEMGPU(polydeg = 3,
+                  surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal),
+                  volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 ###############################################################################
 # Get the TreeMesh and setup a periodic mesh

diff --git a/src/TrixiCUDA.jl b/src/TrixiCUDA.jl
@@ -7,21 +7,30 @@ using CUDA
 using CUDA: @cuda, CuArray, HostKernel,
             threadIdx, blockIdx, blockDim, reshape, similar, launch_configuration
 
-using Trixi: AbstractEquations, AbstractContainer, AbstractMesh, AbstractSemidiscretization,
-             True, False, TreeMesh, DGSEM, SemidiscretizationHyperbolic,
-             ElementContainer1D, ElementContainer2D, ElementContainer3D,
-             InterfaceContainer1D, InterfaceContainer2D, InterfaceContainer3D,
-             BoundaryContainer1D, BoundaryContainer2D, BoundaryContainer3D,
-             LobattoLegendreMortarL2, L2MortarContainer2D, L2MortarContainer3D,
-             BoundaryConditionPeriodic, BoundaryConditionDirichlet,
-             VolumeIntegralWeakForm, VolumeIntegralFluxDifferencing, VolumeIntegralShockCapturingHG,
-             allocate_coefficients, compute_coefficients, mesh_equations_solver_cache,
-             flux, ntuple, nvariables, nnodes, nelements, nmortars,
+# Trixi.jl methods
+using Trixi: allocate_coefficients, compute_coefficients, mesh_equations_solver_cache,
+             flux, ntuple, nnodes, nvariables, nelements,
              local_leaf_cells, init_elements, init_interfaces, init_boundaries, init_mortars,
              have_nonconservative_terms, boundary_condition_periodic,
              digest_boundary_conditions, check_periodicity_mesh_boundary_conditions,
+             gauss_lobatto_nodes_weights, vandermonde_legendre,
+             calc_dsplit, calc_dhat, calc_lhat, polynomial_derivative_matrix,
+             calc_forward_upper, calc_forward_lower, calc_reverse_upper, calc_reverse_lower,
              set_log_type!, set_sqrt_type!
 
+# Trixi.jl structs
+using Trixi: AbstractEquations, AbstractContainer, AbstractMesh, AbstractSemidiscretization,
+             AbstractSurfaceIntegral,
+             True, False, TreeMesh, DG, DGSEM, SemidiscretizationHyperbolic,
+             LobattoLegendreBasis, LobattoLegendreMortarL2,
+             ElementContainer1D, ElementContainer2D, ElementContainer3D,
+             InterfaceContainer1D, InterfaceContainer2D, InterfaceContainer3D,
+             BoundaryContainer1D, BoundaryContainer2D, BoundaryContainer3D,
+             L2MortarContainer2D, L2MortarContainer3D,
+             SurfaceIntegralWeakForm, VolumeIntegralWeakForm,
+             BoundaryConditionPeriodic,
+             VolumeIntegralFluxDifferencing, VolumeIntegralShockCapturingHG
+
 import Trixi: get_node_vars, get_node_coords, get_surface_node_vars,
               nelements, ninterfaces, nmortars, wrap_array, wrap_array_native
 
@@ -40,6 +49,8 @@ include("semidiscretization/semidiscretization.jl")
 include("solvers/solvers.jl")
 
 # Export the public APIs
+export LobattoLegendreBasisGPU
+export DGSEMGPU
 export SemidiscretizationHyperbolicGPU
 export semidiscretizeGPU
 

diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl
@@ -1,6 +1,9 @@
 # Everything specific about semidiscretization hyperbolic for PDE solvers.
 
-# Similar to `SemidiscretizationHyperbolic` in Trixi.jl but for GPU cache 
+# Similar to `SemidiscretizationHyperbolic` in Trixi.jl 
+# No need to adapt the `SemidiscretizationHyperbolic` struct as it is already GPU compatible
+
+# Outer constructor for GPU type
 function SemidiscretizationHyperbolicGPU(mesh, equations, initial_condition, solver;
                                          source_terms = nothing,
                                          boundary_conditions = boundary_condition_periodic,
@@ -16,7 +19,7 @@ function SemidiscretizationHyperbolicGPU(mesh, equations, initial_condition, sol
 
     check_periodicity_mesh_boundary_conditions(mesh, _boundary_conditions)
 
-    # Return the CPU type
+    # Return the CPU type (GPU compatible)
     SemidiscretizationHyperbolic{typeof(mesh), typeof(equations),
                                  typeof(initial_condition),
                                  typeof(_boundary_conditions), typeof(source_terms),