diff --git a/examples/p4est_2d_dgsem/elixir_euler_source_terms.jl b/examples/p4est_2d_dgsem/elixir_euler_source_terms.jl new file mode 100644 index 00000000000..f85658c5423 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_source_terms.jl @@ -0,0 +1,70 @@ +# The same setup as tree_2d_dgsem/elixir_euler_source_terms.jl +# to verify the P4estMesh implementation against TreeMesh + +using OrdinaryDiffEqLowStorageRK +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +initial_condition = initial_condition_convergence_test + +# Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of +# `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`. +# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`. +# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`. +# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`. +# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the +# `StepsizeCallback` (CFL-Condition) and less diffusion. +solver = DGSEM(polydeg = 3, surface_flux = FluxLaxFriedrichs(max_abs_speed_naive)) + +coordinates_min = (0.0, 0.0) +coordinates_max = (2.0, 2.0) + +trees_per_dimension = (16, 16) +mesh = P4estMesh(trees_per_dimension, + polydeg = 3, initial_refinement_level = 0, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + periodicity = true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; + source_terms = source_terms_convergence_test, + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +# Create ODE problem with time span from 0.0 to 1.0 +# Setting `real_type` allows to change the real number type, e.g., to `Float32`. +# This is particularly useful when changing the `storage_type` to a GPU array +# type such as `ROCArray` (AMD) or `CuArray` (NVIDIA CUDA). +ode = semidiscretize(semi, tspan; real_type = nothing, storage_type = nothing) + +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 = 1.0) + +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 + ode_default_options()..., callback = callbacks); diff --git a/examples/p4est_3d_dgsem/elixir_euler_source_terms.jl b/examples/p4est_3d_dgsem/elixir_euler_source_terms.jl new file mode 100644 index 00000000000..3fce3e54fb2 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_euler_source_terms.jl @@ -0,0 +1,72 @@ +# The same setup as tree_3d_dgsem/elixir_euler_source_terms.jl +# to verify the StructuredMesh implementation against TreeMesh + +using OrdinaryDiffEqLowStorageRK +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquations3D(gamma) + +initial_condition = initial_condition_convergence_test + +# Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of +# `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`. +# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`. +# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`. +# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`. +# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the +# `StepsizeCallback` (CFL-Condition) and less diffusion. +solver = DGSEM(polydeg = 3, surface_flux = FluxLaxFriedrichs(max_abs_speed_naive), + volume_integral = VolumeIntegralWeakForm()) + +coordinates_min = (0.0, 0.0, 0.0) +coordinates_max = (2.0, 2.0, 2.0) + +trees_per_dimension = (4, 4, 4) + +mesh = P4estMesh(trees_per_dimension, polydeg = 3, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + initial_refinement_level = 1, + periodicity = true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; + source_terms = source_terms_convergence_test, + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 5.0) +# Create ODE problem with time span from 0.0 to 1.0 +# Setting `real_type` allows to change the real number type, e.g., to `Float32`. +# This is particularly useful when changing the `storage_type` to a GPU array +# type such as `ROCArray` (AMD) or `CuArray` (NVIDIA CUDA). +ode = semidiscretize(semi, tspan; real_type = nothing, storage_type = nothing) + +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.6) + +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 + ode_default_options()..., callback = callbacks); diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index 34fe34f14dc..713cb491f03 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -93,6 +93,7 @@ include("dg_2d_gpu.jl") include("dg_3d.jl") include("dg_3d_parabolic.jl") include("dg_parallel.jl") +include("dg_3d_gpu.jl") # Subcell limiters include("subcell_limiters.jl") diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 47d9f649cf3..7bb0c928328 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -1245,7 +1245,7 @@ function rhs!(du, u, t, u_parent, semis, # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache) + calc_sources!(backend, du, u, t, source_terms, equations, dg, cache) end return nothing diff --git a/src/solvers/dgsem_p4est/dg_2d_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_gpu.jl index 925e7a809f8..17c5d3c2b9d 100644 --- a/src/solvers/dgsem_p4est/dg_2d_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_gpu.jl @@ -118,4 +118,33 @@ end apply_jacobian_per_quadrature_node!(du, MeshT, equations, dg, inverse_jacobian, i, j, element) end + +@kernel function calc_sources_KAkernel!(du, u, t, source_terms, + node_coordinates, + equations::AbstractEquations{2}, dg, cache) + i, j, element = @index(Global, NTuple) + u_local = get_node_vars(u, equations, dg, i, j, element) + x_local = get_node_coords(node_coordinates, equations, dg, i, j, element) + + du_local = source_terms(u_local, x_local, t, equations) + + add_to_node_vars!(du, du_local, equations, dg, i, j, element) +end + +function calc_sources!(backend::Backend, du, u, t, source_terms, + equations::AbstractEquations{2}, dg::DG, cache) + nelements(dg, cache) == 0 && return nothing + @unpack node_coordinates = cache.elements + kernel_cache = kernel_filter_cache(cache) + kernel! = calc_sources_KAkernel!(backend) + kernel!(du, u, t, source_terms, node_coordinates, equations, dg, kernel_cache, + ndrange = (nnodes(dg), nnodes(dg), nelements(dg, cache))) + + return nothing +end + +function calc_sources!(backend::Backend, du, u, t, source_terms::Nothing, + equations::AbstractEquations{2}, dg::DG, cache) + return nothing end +end #muladd diff --git a/src/solvers/dgsem_p4est/dg_3d.jl b/src/solvers/dgsem_p4est/dg_3d.jl index faac5fa46f1..ab3b910925c 100644 --- a/src/solvers/dgsem_p4est/dg_3d.jl +++ b/src/solvers/dgsem_p4est/dg_3d.jl @@ -106,27 +106,6 @@ function prolong2interfaces!(backend::Nothing, cache, u, return nothing end -function prolong2interfaces!(backend::Backend, cache, u, - mesh::Union{P4estMesh{3}, T8codeMesh{3}}, - equations, dg::DG) - @unpack interfaces = cache - @unpack neighbor_ids, node_indices = cache.interfaces - index_range = eachnode(dg) - - kernel! = prolong2interfaces_KAkernel!(backend) - kernel!(interfaces.u, u, typeof(mesh), equations, neighbor_ids, node_indices, - index_range, - ndrange = ninterfaces(interfaces)) - return nothing -end - -@kernel function prolong2interfaces_KAkernel!(interface_u, u, MeshT, equations, - neighbor_ids, node_indices, index_range) - interface = @index(Global) - prolong2interfaces_per_interface!(interface_u, u, MeshT, equations, neighbor_ids, - node_indices, index_range, interface) -end - @inline function prolong2interfaces_per_interface!(u_interface, u, ::Type{<:Union{P4estMesh{3}, T8codeMesh{3}}}, @@ -223,38 +202,6 @@ function calc_interface_flux!(backend::Nothing, surface_flux_values, return nothing end -function calc_interface_flux!(backend::Backend, surface_flux_values, - mesh::Union{P4estMesh{3}, T8codeMesh{3}}, - have_nonconservative_terms, - equations, surface_integral, dg::DG, cache) - @unpack neighbor_ids, node_indices = cache.interfaces - @unpack contravariant_vectors = cache.elements - index_range = eachnode(dg) - - kernel! = calc_interface_flux_KAkernel!(backend) - kernel!(surface_flux_values, typeof(mesh), have_nonconservative_terms, equations, - surface_integral, typeof(dg), cache.interfaces.u, - neighbor_ids, node_indices, contravariant_vectors, index_range, - ndrange = ninterfaces(cache.interfaces)) - return nothing -end - -@kernel function calc_interface_flux_KAkernel!(surface_flux_values, MeshT, - have_nonconservative_terms, equations, - surface_integral, SolverT, u_interface, - neighbor_ids, node_indices, - contravariant_vectors, index_range) - interface = @index(Global) - calc_interface_flux_per_interface!(surface_flux_values, - MeshT, - have_nonconservative_terms, - equations, surface_integral, SolverT, - u_interface, - neighbor_ids, node_indices, - contravariant_vectors, - index_range, interface) -end - @inline function calc_interface_flux_per_interface!(surface_flux_values, MeshT::Type{<:Union{P4estMesh{3}, T8codeMesh{3}}}, @@ -1026,29 +973,6 @@ function calc_surface_integral!(backend::Nothing, du, u, return nothing end -function calc_surface_integral!(backend::Backend, du, u, - mesh::Union{P4estMesh{3}, T8codeMesh{3}}, - equations, - surface_integral::SurfaceIntegralWeakForm, - dg::DGSEM, cache) - @unpack inverse_weights = dg.basis - @unpack surface_flux_values = cache.elements - - kernel! = calc_surface_integral_KAkernel!(backend) - kernel!(du, typeof(mesh), equations, surface_integral, dg, inverse_weights[1], - surface_flux_values, ndrange = nelements(cache.elements)) - return nothing -end - -@kernel function calc_surface_integral_KAkernel!(du, MeshT, equations, - surface_integral, dg, factor, - surface_flux_values) - element = @index(Global) - calc_surface_integral_per_element!(du, MeshT, - equations, surface_integral, dg, factor, - surface_flux_values, element) -end - @inline function calc_surface_integral_per_element!(du, ::Type{<:Union{P4estMesh{3}, T8codeMesh{3}}}, diff --git a/src/solvers/dgsem_p4est/dg_3d_gpu.jl b/src/solvers/dgsem_p4est/dg_3d_gpu.jl new file mode 100644 index 00000000000..046cf498dac --- /dev/null +++ b/src/solvers/dgsem_p4est/dg_3d_gpu.jl @@ -0,0 +1,112 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +function prolong2interfaces!(backend::Backend, cache, u, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, + equations, dg::DG) + @unpack interfaces = cache + @unpack neighbor_ids, node_indices = cache.interfaces + index_range = eachnode(dg) + + kernel! = prolong2interfaces_KAkernel!(backend) + kernel!(interfaces.u, u, typeof(mesh), equations, neighbor_ids, node_indices, + index_range, + ndrange = ninterfaces(interfaces)) + return nothing +end + +@kernel function prolong2interfaces_KAkernel!(interface_u, u, MeshT, equations, + neighbor_ids, node_indices, index_range) + interface = @index(Global) + prolong2interfaces_per_interface!(interface_u, u, MeshT, equations, neighbor_ids, + node_indices, index_range, interface) +end + +function calc_surface_integral!(backend::Backend, du, u, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, + equations, + surface_integral::SurfaceIntegralWeakForm, + dg::DGSEM, cache) + @unpack inverse_weights = dg.basis + @unpack surface_flux_values = cache.elements + + kernel! = calc_surface_integral_KAkernel!(backend) + kernel!(du, typeof(mesh), equations, surface_integral, dg, inverse_weights[1], + surface_flux_values, ndrange = nelements(cache.elements)) + return nothing +end + +@kernel function calc_surface_integral_KAkernel!(du, MeshT, equations, + surface_integral, dg, factor, + surface_flux_values) + element = @index(Global) + calc_surface_integral_per_element!(du, MeshT, + equations, surface_integral, dg, factor, + surface_flux_values, element) +end + +function calc_interface_flux!(backend::Backend, surface_flux_values, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, + have_nonconservative_terms, + equations, surface_integral, dg::DG, cache) + @unpack neighbor_ids, node_indices = cache.interfaces + @unpack contravariant_vectors = cache.elements + index_range = eachnode(dg) + + kernel! = calc_interface_flux_KAkernel!(backend) + kernel!(surface_flux_values, typeof(mesh), have_nonconservative_terms, equations, + surface_integral, typeof(dg), cache.interfaces.u, + neighbor_ids, node_indices, contravariant_vectors, index_range, + ndrange = ninterfaces(cache.interfaces)) + return nothing +end + +@kernel function calc_interface_flux_KAkernel!(surface_flux_values, MeshT, + have_nonconservative_terms, equations, + surface_integral, SolverT, u_interface, + neighbor_ids, node_indices, + contravariant_vectors, index_range) + interface = @index(Global) + calc_interface_flux_per_interface!(surface_flux_values, + MeshT, + have_nonconservative_terms, + equations, surface_integral, SolverT, + u_interface, + neighbor_ids, node_indices, + contravariant_vectors, + index_range, interface) +end + +@kernel function calc_sources_KAkernel!(du, u, t, source_terms, + node_coordinates, + equations::AbstractEquations{3}, dg, cache) + i, j, k, element = @index(Global, NTuple) + u_local = get_node_vars(u, equations, dg, i, j, k, element) + x_local = get_node_coords(node_coordinates, equations, dg, i, j, k, element) + + du_local = source_terms(u_local, x_local, t, equations) + + add_to_node_vars!(du, du_local, equations, dg, i, j, k, element) +end + +function calc_sources!(backend::Backend, du, u, t, source_terms, + equations::AbstractEquations{3}, dg::DG, cache) + nelements(dg, cache) == 0 && return nothing + @unpack node_coordinates = cache.elements + kernel_cache = kernel_filter_cache(cache) + kernel! = calc_sources_KAkernel!(backend) + kernel!(du, u, t, source_terms, node_coordinates, equations, dg, kernel_cache, + ndrange = (nnodes(dg), nnodes(dg), nnodes(dg), nelements(dg, cache))) + + return nothing +end + +function calc_sources!(backend::Backend, du, u, t, source_terms::Nothing, + equations::AbstractEquations{3}, dg::DG, cache) + return nothing +end +end #muladd diff --git a/src/solvers/dgsem_p4est/dg_3d_parallel.jl b/src/solvers/dgsem_p4est/dg_3d_parallel.jl index 91a700a847c..6fc4e2009b0 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parallel.jl @@ -107,7 +107,7 @@ function rhs!(du, u, t, # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache) + calc_sources!(backend, du, u, t, source_terms, equations, dg, cache) end # Finish to send MPI data diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index f3e8365b6f6..5ac27049151 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -88,7 +88,7 @@ function rhs!(du, u, t, # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache) + calc_sources!(backend, du, u, t, source_terms, equations, dg, cache) end return nothing diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 65383168835..8f195490767 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -111,7 +111,7 @@ function rhs!(du, u, t, # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache) + calc_sources!(backend, du, u, t, source_terms, equations, dg, cache) end return nothing @@ -765,12 +765,12 @@ function apply_jacobian!(backend::Nothing, du, mesh::TreeMesh{1}, end # Need dimension specific version to avoid error at dispatching -function calc_sources!(du, u, t, source_terms::Nothing, +function calc_sources!(backend::Nothing, du, u, t, source_terms::Nothing, equations::AbstractEquations{1}, dg::DG, cache) return nothing end -function calc_sources!(du, u, t, source_terms, +function calc_sources!(backend::Nothing, du, u, t, source_terms, equations::AbstractEquations{1}, dg::DG, cache) @unpack node_coordinates = cache.elements diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 815e0c0df13..29c881943ae 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -178,7 +178,7 @@ function rhs!(du, u, t, # Calculate source terms @trixi_timeit_ext backend timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache) + calc_sources!(backend, du, u, t, source_terms, equations, dg, cache) end return nothing @@ -1359,12 +1359,12 @@ function apply_jacobian!(backend::Nothing, du, mesh::TreeMesh{2}, end # Need dimension specific version to avoid error at dispatching -function calc_sources!(du, u, t, source_terms::Nothing, +function calc_sources!(backend::Nothing, du, u, t, source_terms::Nothing, equations::AbstractEquations{2}, dg::DG, cache) return nothing end -function calc_sources!(du, u, t, source_terms, +function calc_sources!(backend::Nothing, du, u, t, source_terms, equations::AbstractEquations{2}, dg::DG, cache) @unpack node_coordinates = cache.elements diff --git a/src/solvers/dgsem_tree/dg_2d_parallel.jl b/src/solvers/dgsem_tree/dg_2d_parallel.jl index d4b6192a4df..822ee9b7d26 100644 --- a/src/solvers/dgsem_tree/dg_2d_parallel.jl +++ b/src/solvers/dgsem_tree/dg_2d_parallel.jl @@ -554,7 +554,7 @@ function rhs!(du, u, t, # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache) + calc_sources!(backend, du, u, t, source_terms, equations, dg, cache) end # Finish to send MPI data diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl index 3fe1bc8b0a3..177dd5953a6 100644 --- a/src/solvers/dgsem_tree/dg_3d.jl +++ b/src/solvers/dgsem_tree/dg_3d.jl @@ -1422,12 +1422,12 @@ function apply_jacobian!(backend::Nothing, du, mesh::TreeMesh{3}, end # Need dimension specific version to avoid error at dispatching -function calc_sources!(du, u, t, source_terms::Nothing, +function calc_sources!(backend::Nothing, du, u, t, source_terms::Nothing, equations::AbstractEquations{3}, dg::DG, cache) return nothing end -function calc_sources!(du, u, t, source_terms, +function calc_sources!(backend::Nothing, du, u, t, source_terms, equations::AbstractEquations{3}, dg::DG, cache) @unpack node_coordinates = cache.elements diff --git a/test/test_amdgpu_2d.jl b/test/test_amdgpu_2d.jl index 13f98e1b958..145786eac24 100644 --- a/test/test_amdgpu_2d.jl +++ b/test/test_amdgpu_2d.jl @@ -67,6 +67,73 @@ end @test Trixi.storage_type(ode.p.cache.mortars) === ROCArray end +@trixi_testset "elixir_euler_source_terms.jl native" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), + # Expected errors are exactly the same as with TreeMesh! + l2=[9.321181254378498e-7, + 1.418121074369651e-6, + 1.4181210743821669e-6, + 4.824553091168877e-6], + linf=[9.577246532499473e-6, + 1.1707525985116263e-5, + 1.1707525982673772e-5, + 4.886961559069647e-5]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. + @test_allocations(Trixi.rhs!, semi, sol, 1000) + @test real(ode.p.solver) == Float64 + @test real(ode.p.solver.basis) == Float64 + @test real(ode.p.solver.mortar) == Float64 + # TODO: `mesh` is currently not `adapt`ed correctly + @test real(ode.p.mesh) == Float64 + @test typeof(equations.gamma) == Float64 + + @test ode.u0 isa Array + @test ode.p.solver.basis.derivative_matrix isa Array + + @test Trixi.storage_type(ode.p.cache.elements) === Array + @test Trixi.storage_type(ode.p.cache.interfaces) === Array + @test Trixi.storage_type(ode.p.cache.boundaries) === Array + @test Trixi.storage_type(ode.p.cache.mortars) === Array +end + +@trixi_testset "elixir_euler_source_terms.jl Float32 / AMDGPU" begin + # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules + using AMDGPU + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), + l2=Float32[2.4917018095933837e-6, + 2.7148269885239423e-6, + 2.695290306860358e-6, + 6.243861976167833e-6], + linf=Float32[1.6489475493930428e-5, + 1.7499923706143505e-5, + 1.893043518075288e-5, + 6.214141845717336e-5], + RealT_for_test_tolerances=Float32, + real_type=Float32, + storage_type=ROCArray, + gamma=Float32(1.4)) # TODO: This should not be required + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. + @test_allocations(Trixi.rhs!, semi, sol, 100_000) + @test real(ode.p.solver) == Float32 + @test real(ode.p.solver.basis) == Float32 + @test real(ode.p.solver.mortar) == Float32 + # TODO: `mesh` is currently not `adapt`ed correctly + @test real(ode.p.mesh) == Float64 + @test typeof(equations.gamma) == Float32 + + @test ode.u0 isa ROCArray + @test ode.p.solver.basis.derivative_matrix isa ROCArray + + @test Trixi.storage_type(ode.p.cache.elements) === ROCArray + @test Trixi.storage_type(ode.p.cache.interfaces) === ROCArray + @test Trixi.storage_type(ode.p.cache.boundaries) === ROCArray + @test Trixi.storage_type(ode.p.cache.mortars) === ROCArray +end + # Clean up afterwards: delete Trixi.jl output directory @test_nowarn isdir(outdir) && rm(outdir, recursive = true) end diff --git a/test/test_amdgpu_3d.jl b/test/test_amdgpu_3d.jl index b70ca55b962..328861732a2 100644 --- a/test/test_amdgpu_3d.jl +++ b/test/test_amdgpu_3d.jl @@ -67,6 +67,80 @@ end @test Trixi.storage_type(ode.p.cache.mortars) === ROCArray end +@trixi_testset "elixir_euler_source_terms.jl native" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), + l2=[ + 4.893619139889976e-5, + 5.3526950567182756e-5, + 5.35269505672133e-5, + 5.352695056735998e-5, + 0.00015172095200428318 + ], + linf=[ + 0.00031179856625374036, + 0.0003368725355339386, + 0.0003368725355383795, + 0.00033687253560787944, + 0.0013193387520935573 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. + @test_allocations(Trixi.rhs!, semi, sol, 1000) + @test real(ode.p.solver) == Float64 + @test real(ode.p.solver.basis) == Float64 + @test real(ode.p.solver.mortar) == Float64 + # TODO: `mesh` is currently not `adapt`ed correctly + @test real(ode.p.mesh) == Float64 + @test typeof(equations.gamma) == Float64 + + @test ode.u0 isa Array + @test ode.p.solver.basis.derivative_matrix isa Array + + @test Trixi.storage_type(ode.p.cache.elements) === Array + @test Trixi.storage_type(ode.p.cache.interfaces) === Array + @test Trixi.storage_type(ode.p.cache.boundaries) === Array + @test Trixi.storage_type(ode.p.cache.mortars) === Array +end + +@trixi_testset "elixir_euler_source_terms.jl Float32 / AMDGPU" begin + # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules + using AMDGPU + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), + l2=Float32[4.912578089985958e-5, + 5.3683407014580115e-5, + 5.368099834769191e-5, + 5.371664525206341e-5, + 0.00015186256300882088], + linf=Float32[0.00032772542853032327, + 0.00035144807715092874, + 0.0003549051465479014, + 0.00035573961157475686, + 0.0013591384887696734], + RealT_for_test_tolerances=Float32, + real_type=Float32, + storage_type=ROCArray, + gamma=Float32(1.4)) # TODO: This should not be required + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. + @test_allocations(Trixi.rhs!, semi, sol, 100_000) + @test real(ode.p.solver) == Float32 + @test real(ode.p.solver.basis) == Float32 + @test real(ode.p.solver.mortar) == Float32 + # TODO: `mesh` is currently not `adapt`ed correctly + @test real(ode.p.mesh) == Float64 + @test typeof(equations.gamma) == Float32 + + @test ode.u0 isa ROCArray + @test ode.p.solver.basis.derivative_matrix isa ROCArray + + @test Trixi.storage_type(ode.p.cache.elements) === ROCArray + @test Trixi.storage_type(ode.p.cache.interfaces) === ROCArray + @test Trixi.storage_type(ode.p.cache.boundaries) === ROCArray + @test Trixi.storage_type(ode.p.cache.mortars) === ROCArray +end + # Clean up afterwards: delete Trixi.jl output directory @test_nowarn isdir(outdir) && rm(outdir, recursive = true) end diff --git a/test/test_cuda_2d.jl b/test/test_cuda_2d.jl index 85fcb75139f..70c822814c6 100644 --- a/test/test_cuda_2d.jl +++ b/test/test_cuda_2d.jl @@ -67,6 +67,73 @@ end @test Trixi.storage_type(ode.p.cache.mortars) === CuArray end +@trixi_testset "elixir_euler_source_terms.jl native" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), + # Expected errors are exactly the same as with TreeMesh! + l2=[9.321181254378498e-7, + 1.418121074369651e-6, + 1.4181210743821669e-6, + 4.824553091168877e-6], + linf=[9.577246532499473e-6, + 1.1707525985116263e-5, + 1.1707525982673772e-5, + 4.886961559069647e-5]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. + @test_allocations(Trixi.rhs!, semi, sol, 1000) + @test real(ode.p.solver) == Float64 + @test real(ode.p.solver.basis) == Float64 + @test real(ode.p.solver.mortar) == Float64 + # TODO: `mesh` is currently not `adapt`ed correctly + @test real(ode.p.mesh) == Float64 + @test typeof(equations.gamma) == Float64 + + @test ode.u0 isa Array + @test ode.p.solver.basis.derivative_matrix isa Array + + @test Trixi.storage_type(ode.p.cache.elements) === Array + @test Trixi.storage_type(ode.p.cache.interfaces) === Array + @test Trixi.storage_type(ode.p.cache.boundaries) === Array + @test Trixi.storage_type(ode.p.cache.mortars) === Array +end + +@trixi_testset "elixir_euler_source_terms.jl Float32 / CUDA" begin + # Using CUDA inside the testset since otherwise the bindings are hiddend by the anonymous modules + using CUDA + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), + l2=Float32[2.4917018095933837e-6, + 2.7148269885239423e-6, + 2.695290306860358e-6, + 6.243861976167833e-6], + linf=Float32[1.6489475493930428e-5, + 1.7499923706143505e-5, + 1.893043518075288e-5, + 6.214141845717336e-5], + RealT_for_test_tolerances=Float32, + real_type=Float32, + storage_type=CuArray, + gamma=Float32(1.4)) # TODO: This should not be required + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. + @test_allocations(Trixi.rhs!, semi, sol, 100_000) + @test real(ode.p.solver) == Float32 + @test real(ode.p.solver.basis) == Float32 + @test real(ode.p.solver.mortar) == Float32 + # TODO: `mesh` is currently not `adapt`ed correctly + @test real(ode.p.mesh) == Float64 + @test typeof(equations.gamma) == Float32 + + @test ode.u0 isa CuArray + @test ode.p.solver.basis.derivative_matrix isa CuArray + + @test Trixi.storage_type(ode.p.cache.elements) === CuArray + @test Trixi.storage_type(ode.p.cache.interfaces) === CuArray + @test Trixi.storage_type(ode.p.cache.boundaries) === CuArray + @test Trixi.storage_type(ode.p.cache.mortars) === CuArray +end + # Clean up afterwards: delete Trixi.jl output directory @test_nowarn isdir(outdir) && rm(outdir, recursive = true) end diff --git a/test/test_cuda_3d.jl b/test/test_cuda_3d.jl index 6c590332555..95809d22146 100644 --- a/test/test_cuda_3d.jl +++ b/test/test_cuda_3d.jl @@ -67,6 +67,80 @@ end @test Trixi.storage_type(ode.p.cache.mortars) === CuArray end +@trixi_testset "elixir_euler_source_terms.jl native" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), + l2=[ + 4.893619139889976e-5, + 5.3526950567182756e-5, + 5.35269505672133e-5, + 5.352695056735998e-5, + 0.00015172095200428318 + ], + linf=[ + 0.00031179856625374036, + 0.0003368725355339386, + 0.0003368725355383795, + 0.00033687253560787944, + 0.0013193387520935573 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. + @test_allocations(Trixi.rhs!, semi, sol, 1000) + @test real(ode.p.solver) == Float64 + @test real(ode.p.solver.basis) == Float64 + @test real(ode.p.solver.mortar) == Float64 + # TODO: `mesh` is currently not `adapt`ed correctly + @test real(ode.p.mesh) == Float64 + @test typeof(equations.gamma) == Float64 + + @test ode.u0 isa Array + @test ode.p.solver.basis.derivative_matrix isa Array + + @test Trixi.storage_type(ode.p.cache.elements) === Array + @test Trixi.storage_type(ode.p.cache.interfaces) === Array + @test Trixi.storage_type(ode.p.cache.boundaries) === Array + @test Trixi.storage_type(ode.p.cache.mortars) === Array +end + +@trixi_testset "elixir_euler_source_terms.jl Float32 / CUDA" begin + # Using CUDA inside the testset since otherwise the bindings are hiddend by the anonymous modules + using CUDA + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), + l2=Float32[4.912578089985958e-5, + 5.3683407014580115e-5, + 5.368099834769191e-5, + 5.371664525206341e-5, + 0.00015186256300882088], + linf=Float32[0.00032772542853032327, + 0.00035144807715092874, + 0.0003549051465479014, + 0.00035573961157475686, + 0.0013591384887696734], + RealT_for_test_tolerances=Float32, + real_type=Float32, + storage_type=CuArray, + gamma=Float32(1.4)) # TODO: This should not be required + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. + @test_allocations(Trixi.rhs!, semi, sol, 100_000) + @test real(ode.p.solver) == Float32 + @test real(ode.p.solver.basis) == Float32 + @test real(ode.p.solver.mortar) == Float32 + # TODO: `mesh` is currently not `adapt`ed correctly + @test real(ode.p.mesh) == Float64 + @test typeof(equations.gamma) == Float64 + + @test ode.u0 isa CuArray + @test ode.p.solver.basis.derivative_matrix isa CuArray + + @test Trixi.storage_type(ode.p.cache.elements) === CuArray + @test Trixi.storage_type(ode.p.cache.interfaces) === CuArray + @test Trixi.storage_type(ode.p.cache.boundaries) === CuArray + @test Trixi.storage_type(ode.p.cache.mortars) === CuArray +end + # Clean up afterwards: delete Trixi.jl output directory @test_nowarn isdir(outdir) && rm(outdir, recursive = true) end diff --git a/test/test_kernelabstractions.jl b/test/test_kernelabstractions.jl index a1a771ee402..1e39e058131 100644 --- a/test/test_kernelabstractions.jl +++ b/test/test_kernelabstractions.jl @@ -44,6 +44,44 @@ end semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. @test_allocations(Trixi.rhs!, ode.p, sol, 60_000) end + +@trixi_testset "elixir_euler_source_terms.jl native" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", + "elixir_euler_source_terms.jl"), + # Expected errors are exactly the same as with TreeMesh! + l2=[9.321181254378498e-7, + 1.418121074369651e-6, + 1.4181210743821669e-6, + 4.824553091168877e-6], + linf=[9.577246532499473e-6, + 1.1707525985116263e-5, + 1.1707525982673772e-5, + 4.886961559069647e-5]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. + @test_allocations(Trixi.rhs!, semi, sol, 100_000) +end + +@trixi_testset "elixir_euler_source_terms.jl Float32" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", + "elixir_euler_source_terms.jl"), + l2=Float32[2.4917018095933837e-6, + 2.7148269885239423e-6, + 2.695290306860358e-6, + 6.243861976167833e-6], + linf=Float32[1.6489475493930428e-5, + 1.7499923706143505e-5, + 1.893043518075288e-5, + 6.214141845717336e-5], + RealT_for_test_tolerances=Float32, + real_type=Float32, + gamma=Float32(1.4)) # TODO: This should not be required + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. + @test_allocations(Trixi.rhs!, semi, sol, 100_000) +end end @testset "KernelAbstractions CPU 3D" begin @@ -74,6 +112,51 @@ end semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. @test_allocations(Trixi.rhs!, semi, sol, 370_000) end + +@trixi_testset "elixir_euler_source_terms.jl native" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem", + "elixir_euler_source_terms.jl"), + l2=[ + 4.893619139889976e-5, + 5.3526950567182756e-5, + 5.35269505672133e-5, + 5.352695056735998e-5, + 0.00015172095200428318 + ], + linf=[ + 0.00031179856625374036, + 0.0003368725355339386, + 0.0003368725355383795, + 0.00033687253560787944, + 0.0013193387520935573 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. + @test_allocations(Trixi.rhs!, semi, sol, 400_000) +end + +@trixi_testset "elixir_euler_source_terms.jl Float32" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem", + "elixir_euler_source_terms.jl"), + l2=Float32[4.912578089985958e-5, + 5.3683407014580115e-5, + 5.368099834769191e-5, + 5.371664525206341e-5, + 0.00015186256300882088], + linf=Float32[0.00032772542853032327, + 0.00035144807715092874, + 0.0003549051465479014, + 0.00035573961157475686, + 0.0013591384887696734], + RealT_for_test_tolerances=Float32, + real_type=Float32, + gamma=Float32(1.4)) # TODO: This should not be required + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. + @test_allocations(Trixi.rhs!, semi, sol, 400_000) +end end # Clean up afterwards: delete Trixi.jl output directory