From 3d24fb889d5061912a97aae40653055260e8ba7e Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 25 Jun 2026 11:51:00 +0200 Subject: [PATCH 01/46] adapt_structure --- src/solvers/dgsem_p4est/dg_parallel.jl | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_parallel.jl b/src/solvers/dgsem_p4est/dg_parallel.jl index aea544026ef..934350f046a 100644 --- a/src/solvers/dgsem_p4est/dg_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_parallel.jl @@ -50,8 +50,27 @@ end @inline Base.eltype(::P4estMPICache{BufferType}) where {BufferType} = eltype(BufferType) -# @eval due to @muladd -@eval Adapt.@adapt_structure(P4estMPICache) +function Adapt.adapt_structure(to, mpi_cache::P4estMPICache) + mpi_neighbor_ranks = mpi_cache.mpi_neighbor_ranks + mpi_neighbor_interfaces = Adapt.adapt_structure(to, mpi_cache.mpi_neighbor_interfaces) + mpi_neighbor_mortars = Adapt.adapt_structure(to, mpi_cache.mpi_neighbor_mortars) + mpi_send_buffers = Adapt.adapt_structure(to, mpi_cache.mpi_send_buffers) + mpi_recv_buffers = Adapt.adapt_structure(to, mpi_cache.mpi_recv_buffers) + mpi_send_requests = mpi_cache.mpi_send_requests + mpi_recv_requests = mpi_cache.mpi_recv_requests + n_elements_by_rank = mpi_cache.n_elements_by_rank + n_elements_global = mpi_cache.n_elements_global + first_element_global_id = mpi_cache.first_element_global_id + + @assert eltype(mpi_send_buffers) == eltype(mpi_recv_buffers) + BufferType = eltype(mpi_send_buffers) + VecInt = eltype(mpi_neighbor_interfaces) + return P4estMPICache{BufferType, VecInt}(mpi_neighbor_ranks, mpi_neighbor_interfaces, + mpi_neighbor_mortars, mpi_send_buffers, + mpi_recv_buffers, mpi_send_requests, + mpi_recv_requests, n_elements_by_rank, + n_elements_global, first_element_global_id) +end ## # Note that the code in `start_mpi_send`/`finish_mpi_receive!` is sensitive to inference on (at least) Julia 1.10. From 072b754b5bd60084bcb43c220b1a1b09175e58d9 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Fri, 26 Jun 2026 12:27:37 +0200 Subject: [PATCH 02/46] prolong2mpiinterfaces --- src/solvers/dgsem_p4est/dg_2d_parallel.jl | 72 +++++++++++++++-------- src/solvers/dgsem_tree/dg_2d_parallel.jl | 2 +- 2 files changed, 48 insertions(+), 26 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index 2cefabd9539..7be060c89d8 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -5,41 +5,63 @@ @muladd begin #! format: noindent -function prolong2mpiinterfaces!(cache, u, +function prolong2mpiinterfaces!(backend::Nothing, cache, u, mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, equations, surface_integral, dg::DG) - @unpack mpi_interfaces = cache + @unpack local_sides, local_neighbor_ids, node_indices = cache.mpi_interfaces index_range = eachnode(dg) + variables_range = eachvariable(equations) @threaded for interface in eachmpiinterface(dg, cache) - # Copy solution data from the local element using "delayed indexing" with - # a start value and a step size to get the correct face and orientation. - # Note that in the current implementation, the interface will be - # "aligned at the primary element", i.e., the index of the primary side - # will always run forwards. - local_side = mpi_interfaces.local_sides[interface] - local_element = mpi_interfaces.local_neighbor_ids[interface] - local_indices = mpi_interfaces.node_indices[interface] + prolong2mpiinterfaces_per_interface!(cache.mpi_interfaces.u, interface, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u) + end - i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], - index_range) - j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], - index_range) + return nothing +end - i_element = i_element_start - j_element = j_element_start - for i in eachnode(dg) - for v in eachvariable(equations) - mpi_interfaces.u[local_side, v, i, interface] = u[v, i_element, - j_element, - local_element] - end - i_element += i_element_step - j_element += j_element_step +function prolong2mpiinterfaces!(backend::Backend, cache, u, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, + equations, surface_integral, dg::DG) + nmpiinterfaces(dg, cache) == 0 && return nothing + @unpack local_sides, local_neighbor_ids, node_indices = cache.mpi_interfaces + index_range = eachnode(dg) + variables_range = eachvariable(equations) + + kernel! = prolong2mpiinterfaces_kernel!(backend) + kernel!(cache.mpi_interfaces.u, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u, + ndrange = nmpiinterfaces(dg, cache)) + return nothing +end + +@kernel function prolong2mpiinterfaces_kernel!(mpi_interfaces_u, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u) + interface = @index(Global) + prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, interface, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u) +end + +@inline function prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, interface, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u) + # Copy solution data from the local element using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + # Note that in the current implementation, the interface will be + # "aligned at the primary element", i.e., the index of the primary side + # will always run forwards. + local_side = local_sides[interface] + local_element = local_neighbor_ids[interface] + local_indices = node_indices[interface] + + i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], index_range) + j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], index_range) + + i_element = i_element_start + j_element = j_element_start + for i in index_range + for v in variables_range + mpi_interfaces_u[local_side, v, i, interface] = u[v, i_element, j_element, local_element] end + i_element += i_element_step + j_element += j_element_step end - return nothing end diff --git a/src/solvers/dgsem_tree/dg_2d_parallel.jl b/src/solvers/dgsem_tree/dg_2d_parallel.jl index 57ed0033804..9ec05180295 100644 --- a/src/solvers/dgsem_tree/dg_2d_parallel.jl +++ b/src/solvers/dgsem_tree/dg_2d_parallel.jl @@ -462,7 +462,7 @@ function rhs!(du, u, t, # Prolong solution to MPI interfaces @trixi_timeit timer() "prolong2mpiinterfaces" begin - prolong2mpiinterfaces!(cache, u, mesh, equations, dg.surface_integral, dg) + prolong2mpiinterfaces!(backend, cache, u, mesh, equations, dg.surface_integral, dg) end # Prolong solution to MPI mortars From ed8fc4724f98dd3511bdafd34bd72092600804fc Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Wed, 2 Sep 2026 16:44:37 +0200 Subject: [PATCH 03/46] integrate_via_indices 2D --- src/callbacks_step/analysis_dg2d.jl | 4 +- src/callbacks_step/analysis_dg2d_parallel.jl | 57 +++++++++++++++++++- 2 files changed, 59 insertions(+), 2 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 22331ff4c3e..da789d54187 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -339,7 +339,9 @@ function integrate_via_indices(func::Func, u, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, - T8codeMesh{2}}, + T8codeMesh{2}, + P4estMeshParallel{2}, + T8codeMeshParallel{2}}, equations, dg::DGSEM, cache, args...; normalize = true) where {Func} return integrate_via_indices(func, trixi_backend(u), u, mesh, equations, dg, cache, diff --git a/src/callbacks_step/analysis_dg2d_parallel.jl b/src/callbacks_step/analysis_dg2d_parallel.jl index fa2325ae453..5c5dd6a103f 100644 --- a/src/callbacks_step/analysis_dg2d_parallel.jl +++ b/src/callbacks_step/analysis_dg2d_parallel.jl @@ -173,7 +173,7 @@ function integrate_via_indices(func::Func, u, return integral end -function integrate_via_indices(func::Func, u, +function integrate_via_indices(func::Func, backend::Nothing, u, mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, equations, dg::DGSEM, cache, args...; normalize = true) where {Func} @@ -215,4 +215,59 @@ function integrate_via_indices(func::Func, u, return integral end + + +function integrate_via_indices(func::Func, backend::Backend, u, + mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, + equations, + dg::DGSEM, cache, args...; normalize = true) where {Func} + @unpack weights = dg.basis + @unpack inverse_jacobian = cache.elements + + function local_plus((integral, total_volume), (integral_element, volume_element)) + return (integral + integral_element, total_volume + volume_element) + end + + # `func(u, 1,1,1, equations, dg, args...)` might access GPU memory + # so we have to rely on the compiler to correctly infer the type of the integral here. + # TODO: Technically we need device_promote_op here that "infers" the function within the context of the GPU. + integral0 = zero(Base.promote_op(func, typeof(u), Int, Int, Int, typeof(equations), + typeof(dg), map(typeof, args)...)) + volume0 = zero(eltype(weights)) * zero(eltype(inverse_jacobian)) + init = neutral = (integral0, volume0) + + # Use quadrature to numerically integrate over entire domain + num_elements = nelements(dg, cache) + integral, volume = AcceleratedKernels.mapreduce(local_plus, 1:num_elements, + backend; init, + neutral) do element + # Initialize integral with zeros of the right shapeu + local_integral, local_total_volume = neutral + + for j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(inverse_jacobian[i, j, element])) + local_integral += volume_jacobian * weights[i] * weights[j] * + func(u, i, j, element, equations, dg, args...) + local_total_volume += volume_jacobian * weights[i] * weights[j] + end + + return (local_integral, local_total_volume) + end + + global_integral = MPI.Reduce!(Ref(integral), +, mpi_root(), mpi_comm()) + total_volume = MPI.Reduce(volume, +, mpi_root(), mpi_comm()) + if mpi_isroot() + integral = convert(typeof(integral), global_integral[]) + else + integral = convert(typeof(integral), NaN * integral) + total_volume = volume # non-root processes receive nothing from reduce -> overwrite + end + + # Normalize with total volume + if normalize + integral = integral / total_volume + end + + return integral +end end # @muladd From 78ccbda86c1e34aa5fe8587c6ee6a0fc092494c9 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 14 Sep 2026 12:03:22 +0200 Subject: [PATCH 04/46] move to CPU as in non-parallel case --- src/callbacks_step/analysis_dg2d_parallel.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/callbacks_step/analysis_dg2d_parallel.jl b/src/callbacks_step/analysis_dg2d_parallel.jl index 47f3c4def0c..41c9ba06e32 100644 --- a/src/callbacks_step/analysis_dg2d_parallel.jl +++ b/src/callbacks_step/analysis_dg2d_parallel.jl @@ -98,6 +98,13 @@ function calc_error_norms(func, u, t, analyzer, @unpack node_coordinates, inverse_jacobian = cache.elements @unpack u_local, u_tmp1, x_local, x_tmp1, jacobian_local, jacobian_tmp1 = cache_analysis + # Calculate error norms on the CPU, to ensure the order of summation is the same. + if trixi_backend(u) !== nothing + node_coordinates = Array(node_coordinates) + inverse_jacobian = Array(inverse_jacobian) + u = Array(u) + end + # Set up data structures l2_error = zero(func(get_node_vars(u, equations, dg, 1, 1, 1), equations)) linf_error = copy(l2_error) From a737166d02e3b4352f2be5d24047606162c57d74 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 14 Sep 2026 22:30:16 +0200 Subject: [PATCH 05/46] P4estMeshParallel{2} GPU support - prolong2mpiinterfaces - calc_mpi_interface_flux - start_mpi_send - finish_mpi_receive --- src/solvers/dgsem_p4est/dg_2d_gpu.jl | 4 +- src/solvers/dgsem_p4est/dg_2d_parallel.jl | 165 +++++++++++++--------- src/solvers/dgsem_p4est/dg_parallel.jl | 89 +++++++++++- src/solvers/dgsem_tree/dg_2d.jl | 2 +- src/solvers/dgsem_tree/dg_2d_parallel.jl | 11 +- 5 files changed, 192 insertions(+), 79 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_gpu.jl index 81e43a8849c..4e8ccf8d137 100644 --- a/src/solvers/dgsem_p4est/dg_2d_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_gpu.jl @@ -7,12 +7,12 @@ function rhs_hyperbolic!(backend::Backend, du, u, t, - mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}, + mesh::Union{P4estMeshSerial{2}, P4estMeshView{2}, + T8codeMesh{2}, P4estMesh{3}, T8codeMesh{3}}, equations, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} - # Calculate volume integral @trixi_timeit_ext backend timer() "volume integral" begin calc_volume_integral!(backend, du, u, mesh, diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index 7be060c89d8..6e4ea8c5fb7 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -14,7 +14,9 @@ function prolong2mpiinterfaces!(backend::Nothing, cache, u, variables_range = eachvariable(equations) @threaded for interface in eachmpiinterface(dg, cache) - prolong2mpiinterfaces_per_interface!(cache.mpi_interfaces.u, interface, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u) + prolong2mpiinterfaces_per_interface!(cache.mpi_interfaces.u, interface, + local_sides, local_neighbor_ids, node_indices, + index_range, variables_range, u) end return nothing @@ -30,17 +32,25 @@ function prolong2mpiinterfaces!(backend::Backend, cache, u, variables_range = eachvariable(equations) kernel! = prolong2mpiinterfaces_kernel!(backend) - kernel!(cache.mpi_interfaces.u, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u, + kernel!(cache.mpi_interfaces.u, local_sides, local_neighbor_ids, node_indices, + index_range, variables_range, u, ndrange = nmpiinterfaces(dg, cache)) return nothing end -@kernel function prolong2mpiinterfaces_kernel!(mpi_interfaces_u, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u) +@kernel function prolong2mpiinterfaces_kernel!(mpi_interfaces_u, local_sides, + local_neighbor_ids, node_indices, + index_range, variables_range, u) interface = @index(Global) - prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, interface, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u) + prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, interface, local_sides, + local_neighbor_ids, node_indices, index_range, + variables_range, u) end -@inline function prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, interface, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u) +@inline function prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, interface, + local_sides, local_neighbor_ids, + node_indices, index_range, + variables_range, u) # Copy solution data from the local element using "delayed indexing" with # a start value and a step size to get the correct face and orientation. # Note that in the current implementation, the interface will be @@ -50,14 +60,17 @@ end local_element = local_neighbor_ids[interface] local_indices = node_indices[interface] - i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], index_range) - j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], index_range) + i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], + index_range) + j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], + index_range) i_element = i_element_start j_element = j_element_start for i in index_range for v in variables_range - mpi_interfaces_u[local_side, v, i, interface] = u[v, i_element, j_element, local_element] + mpi_interfaces_u[local_side, v, i, interface] = u[v, i_element, j_element, + local_element] end i_element += i_element_step j_element += j_element_step @@ -65,84 +78,102 @@ end return nothing end -function calc_mpi_interface_flux!(surface_flux_values, +function calc_mpi_interface_flux!(backend::Nothing, surface_flux_values, mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, have_nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces @unpack contravariant_vectors = cache.elements + @unpack u = cache.mpi_interfaces index_range = eachnode(dg) - index_end = last(index_range) @threaded for interface in eachmpiinterface(dg, cache) - # Get element and side index information on the local element - local_element = local_neighbor_ids[interface] - local_indices = node_indices[interface] - local_direction = indices2direction(local_indices) - local_side = local_sides[interface] - - # Create the local i,j indexing on the local element used to pull normal direction information - i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], - index_range) - j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], - index_range) - - i_element = i_element_start - j_element = j_element_start - - # Initiate the node index to be used in the surface for loop, - # the surface flux storage must be indexed in alignment with the local element indexing - if :i_backward in local_indices - surface_node = index_end - surface_node_step = -1 - else - surface_node = 1 - surface_node_step = 1 - end + calc_mpi_interface_flux_per_interface!(surface_flux_values, typeof(mesh), + have_nonconservative_terms, equations, + surface_integral, typeof(dg), + local_neighbor_ids, node_indices, + local_sides, contravariant_vectors, u, + index_range, interface) + end - for node in eachnode(dg) - # Get the normal direction on the local element - # Contravariant vectors at interfaces in negative coordinate direction - # are pointing inwards. This is handled by `get_normal_direction`. - normal_direction = get_normal_direction(local_direction, - contravariant_vectors, - i_element, j_element, local_element) - - calc_mpi_interface_flux!(surface_flux_values, mesh, - have_nonconservative_terms, - equations, - surface_integral, dg, cache, - interface, normal_direction, - node, local_side, - surface_node, local_direction, local_element) - - # Increment local element indices to pull the normal direction - i_element += i_element_step - j_element += j_element_step - - # Increment the surface node index along the local element - surface_node += surface_node_step - end + return nothing +end + +@inline function calc_mpi_interface_flux_per_interface!(surface_flux_values, + MeshT::Type{<:Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}}, + have_nonconservative_terms, equations, + surface_integral, SolverT::Type{<:DG}, + local_neighbor_ids, node_indices, local_sides, + contravariant_vectors, u_mpi_interfaces, + index_range, interface) + + # Get element and side index information on the local element + local_element = local_neighbor_ids[interface] + local_indices = node_indices[interface] + local_direction = indices2direction(local_indices) + local_side = local_sides[interface] + + # Create the local i,j indexing on the local element used to pull normal direction information + i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], + index_range) + j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], + index_range) + + i_element = i_element_start + j_element = j_element_start + + # Initiate the node index to be used in the surface for loop, + # the surface flux storage must be indexed in alignment with the local element indexing + if :i_backward in local_indices + surface_node = last(index_range) + surface_node_step = -1 + else + surface_node = 1 + surface_node_step = 1 end + for node in index_range + # Get the normal direction on the local element + # Contravariant vectors at interfaces in negative coordinate direction + # are pointing inwards. This is handled by `get_normal_direction`. + normal_direction = get_normal_direction(local_direction, + contravariant_vectors, + i_element, j_element, local_element) + + calc_mpi_interface_flux!(surface_flux_values, MeshT, + have_nonconservative_terms, equations, + surface_integral, SolverT, + u_mpi_interfaces, interface, normal_direction, + node, local_side, surface_node, local_direction, + local_element) + + # Increment local element indices to pull the normal direction + i_element += i_element_step + j_element += j_element_step + + # Increment the surface node index along the local element + surface_node += surface_node_step + end return nothing end # Inlined version of the interface flux computation for conservation laws @inline function calc_mpi_interface_flux!(surface_flux_values, - mesh::Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}, + MeshT::Type{<:Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}}, have_nonconservative_terms::False, equations, - surface_integral, dg::DG, cache, + surface_integral, SolverT::Type{<:DG}, + u_mpi_interfaces, interface_index, normal_direction, interface_node_index, local_side, surface_node_index, local_direction_index, local_element_index) - @unpack u = cache.mpi_interfaces @unpack surface_flux = surface_integral - u_ll, u_rr = get_surface_node_vars(u, equations, dg, interface_node_index, + u_ll, u_rr = get_surface_node_vars(u_mpi_interfaces, equations, SolverT, + interface_node_index, interface_index) if local_side == 1 @@ -160,18 +191,18 @@ end # Inlined version of the interface flux computation for non-conservative equations @inline function calc_mpi_interface_flux!(surface_flux_values, - mesh::Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}, + MeshT::Type{<:Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}}, have_nonconservative_terms::True, equations, - surface_integral, dg::DG, cache, + surface_integral, SolverT::Type{<:DG}, + u_mpi_interfaces, interface_index, normal_direction, interface_node_index, local_side, surface_node_index, local_direction_index, local_element_index) - @unpack u = cache.mpi_interfaces surface_flux, nonconservative_flux = surface_integral.surface_flux - u_ll, u_rr = get_surface_node_vars(u, equations, dg, + u_ll, u_rr = get_surface_node_vars(u_mpi_interfaces, equations, SolverT, interface_node_index, interface_index) diff --git a/src/solvers/dgsem_p4est/dg_parallel.jl b/src/solvers/dgsem_p4est/dg_parallel.jl index 934350f046a..87f8f29a0da 100644 --- a/src/solvers/dgsem_p4est/dg_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_parallel.jl @@ -52,7 +52,8 @@ end function Adapt.adapt_structure(to, mpi_cache::P4estMPICache) mpi_neighbor_ranks = mpi_cache.mpi_neighbor_ranks - mpi_neighbor_interfaces = Adapt.adapt_structure(to, mpi_cache.mpi_neighbor_interfaces) + mpi_neighbor_interfaces = Adapt.adapt_structure(to, + mpi_cache.mpi_neighbor_interfaces) mpi_neighbor_mortars = Adapt.adapt_structure(to, mpi_cache.mpi_neighbor_mortars) mpi_send_buffers = Adapt.adapt_structure(to, mpi_cache.mpi_send_buffers) mpi_recv_buffers = Adapt.adapt_structure(to, mpi_cache.mpi_recv_buffers) @@ -65,7 +66,8 @@ function Adapt.adapt_structure(to, mpi_cache::P4estMPICache) @assert eltype(mpi_send_buffers) == eltype(mpi_recv_buffers) BufferType = eltype(mpi_send_buffers) VecInt = eltype(mpi_neighbor_interfaces) - return P4estMPICache{BufferType, VecInt}(mpi_neighbor_ranks, mpi_neighbor_interfaces, + return P4estMPICache{BufferType, VecInt}(mpi_neighbor_ranks, + mpi_neighbor_interfaces, mpi_neighbor_mortars, mpi_send_buffers, mpi_recv_buffers, mpi_send_requests, mpi_recv_requests, n_elements_by_rank, @@ -82,7 +84,8 @@ precompile(Base.reindex, (Tuple{Base.Slice{Base.OneTo{Int64}}, Int64, Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Int64}, Tuple{Int64, Int64, Int64})) -function start_mpi_send!(mpi_cache::P4estMPICache, mesh, equations, dg, cache) +function start_mpi_send!(backend::Nothing, mpi_cache::P4estMPICache, mesh, equations, + dg, cache) data_size = nvariables(equations) * nnodes(dg)^(ndims(mesh) - 1) n_small_elements = 2^(ndims(mesh) - 1) @@ -137,6 +140,45 @@ function start_mpi_send!(mpi_cache::P4estMPICache, mesh, equations, dg, cache) return nothing end +# TODO GPU: MPI mortars +function start_mpi_send!(backend::Backend, mpi_cache::P4estMPICache, + mesh::P4estMeshParallel{2}, equations, dg, cache) + @unpack mpi_neighbor_ranks, mpi_neighbor_interfaces = mpi_cache + @unpack mpi_send_buffers, mpi_send_requests = mpi_cache + @unpack local_sides, u = cache.mpi_interfaces + + kernel! = start_mpi_send_KAkernel!(backend) + + for (rank_index, neighbor_rank) in enumerate(mpi_neighbor_ranks) + send_buffer = mpi_send_buffers[rank_index] + neighbor_interfaces = mpi_neighbor_interfaces[rank_index] + kernel!(send_buffer, neighbor_interfaces, local_sides, u, + Val(nvariables(equations)), Val(ndims(mesh)), + ndrange = (nnodes(dg), length(neighbor_interfaces))) + + # wait for the kernel to return before sending the buffer + KernelAbstractions.synchronize(backend) + mpi_send_requests[rank_index] = MPI.Isend(send_buffer, neighbor_rank, + mpi_rank(), mpi_comm()) + end +end + +@kernel function start_mpi_send_KAkernel!(send_buffer, neighbor_interfaces, local_sides, + u_mpi_interfaces, ::Val{NVARS}, + ::Val{2}) where {NVARS} + index_node, index_interface = @index(Global, NTuple) + index_linear = @index(Global, Linear) + + buffer_offset = (index_linear - 1) * NVARS + interface = neighbor_interfaces[index_interface] + local_side = local_sides[interface] + + for v in 1:NVARS + send_buffer[buffer_offset + v] = u_mpi_interfaces[local_side, v, index_node, + index_interface] + end +end + function start_mpi_receive!(mpi_cache::P4estMPICache) for (index, rank) in enumerate(mpi_cache.mpi_neighbor_ranks) mpi_cache.mpi_recv_requests[index] = MPI.Irecv!(mpi_cache.mpi_recv_buffers[index], @@ -150,7 +192,8 @@ function finish_mpi_send!(mpi_cache::P4estMPICache) return MPI.Waitall(mpi_cache.mpi_send_requests, MPI.Status) end -function finish_mpi_receive!(mpi_cache::P4estMPICache, mesh, equations, dg, cache) +function finish_mpi_receive!(backend::Nothing, mpi_cache::P4estMPICache, mesh, + equations, dg, cache) data_size = nvariables(equations) * nnodes(dg)^(ndims(mesh) - 1) n_small_elements = 2^(ndims(mesh) - 1) n_positions = n_small_elements + 1 @@ -200,6 +243,44 @@ function finish_mpi_receive!(mpi_cache::P4estMPICache, mesh, equations, dg, cach return nothing end +# TODO GPU: MPI mortars +function finish_mpi_receive!(backend::Backend, mpi_cache::P4estMPICache, + mesh::P4estMeshParallel{2}, equations, dg, cache) + @unpack mpi_neighbor_interfaces = mpi_cache + @unpack mpi_recv_buffers, mpi_recv_requests = mpi_cache + @unpack local_sides, u = cache.mpi_interfaces + + kernel! = finish_mpi_receive_KAkernel!(backend) + + # Start receiving and unpack received data until all communication is finished + data = MPI.Waitany(mpi_recv_requests) + while data !== nothing + recv_buffer = mpi_recv_buffers[data] + neighbor_interfaces = mpi_neighbor_interfaces[data] + kernel!(recv_buffer, neighbor_interfaces, local_sides, u, + Val(nvariables(equations)), Val(ndims(mesh)), + ndrange = (nnodes(dg), length(neighbor_interfaces))) + + data = MPI.Waitany(mpi_recv_requests) + end + # Wait for the last kernel to return ? + KernelAbstractions.synchronize(backend) +end + +@kernel function finish_mpi_receive_KAkernel!(recv_buffer, neighbor_interfaces, + local_sides, + u_mpi_interfaces, ::Val{NVARS}, + ::Val{3}) where {NVARS} + index_node, index_interface = @index(Global, NTuple) + index_linear = @index(Global, Linear) + buffer_offset = (index_linear - 1) * NVARS + interface = neighbor_interfaces[index_interface] + remote_side = local_sides[interface] == 1 ? 2 : 1 + for v in 1:NVARS + u_mpi_interfaces[remote_side, v, index_node, interface] = recv_buffer[buffer_offset + v] + end +end + # Return a tuple `indices` where indices[position] is a `(first, last)` tuple for accessing the # data corresponding to the `position` part of a mortar in an MPI buffer. The mortar data must begin # at `index_base`+1 in the MPI buffer. `data_size` is the data size associated with each small diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 9022f9100a8..14af40e6a4c 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -112,7 +112,7 @@ end # 2D/3D `TreeMesh`es, `P4estMesh`es, and `T8codeMesh`es. function rhs_hyperbolic!(backend::Nothing, du, u, t, - mesh::Union{TreeMesh{2}, P4estMesh{2}, P4estMeshView{2}, + mesh::Union{TreeMesh{2}, P4estMeshSerial{2}, P4estMeshView{2}, T8codeMesh{2}, TreeMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations, diff --git a/src/solvers/dgsem_tree/dg_2d_parallel.jl b/src/solvers/dgsem_tree/dg_2d_parallel.jl index 17406e9036c..c1e5655e417 100644 --- a/src/solvers/dgsem_tree/dg_2d_parallel.jl +++ b/src/solvers/dgsem_tree/dg_2d_parallel.jl @@ -450,7 +450,7 @@ function init_mpi_neighbor_connectivity(elements, mpi_interfaces, mpi_mortars, return mpi_neighbor_ranks, mpi_neighbor_interfaces, mpi_neighbor_mortars end -function rhs_hyperbolic!(backend::Nothing, +function rhs_hyperbolic!(backend, du, u, t, mesh::Union{TreeMeshParallel{2}, P4estMeshParallel{2}, T8codeMeshParallel{2}}, equations, @@ -461,7 +461,8 @@ function rhs_hyperbolic!(backend::Nothing, # Prolong solution to MPI interfaces @trixi_timeit timer() "prolong2mpiinterfaces" begin - prolong2mpiinterfaces!(backend, cache, u, mesh, equations, dg.surface_integral, dg) + prolong2mpiinterfaces!(backend, cache, u, mesh, equations, dg.surface_integral, + dg) end # Prolong solution to MPI mortars @@ -472,7 +473,7 @@ function rhs_hyperbolic!(backend::Nothing, # Start to send MPI data @trixi_timeit timer() "start MPI send" begin - start_mpi_send!(cache.mpi_cache, mesh, equations, dg, cache) + start_mpi_send!(backend, cache.mpi_cache, mesh, equations, dg, cache) end # Reset du @@ -524,12 +525,12 @@ function rhs_hyperbolic!(backend::Nothing, # Finish to receive MPI data @trixi_timeit timer() "finish MPI receive" begin - finish_mpi_receive!(cache.mpi_cache, mesh, equations, dg, cache) + finish_mpi_receive!(backend, cache.mpi_cache, mesh, equations, dg, cache) end # Calculate MPI interface fluxes @trixi_timeit timer() "MPI interface flux" begin - calc_mpi_interface_flux!(cache.elements.surface_flux_values, mesh, + calc_mpi_interface_flux!(backend, cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, dg.surface_integral, dg, cache) end From b935cedfacee0acbc043b6efa352a861d7cfebd4 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 15 Sep 2026 09:53:07 +0200 Subject: [PATCH 06/46] rhs_hyperbolic for parallel GPU --- src/solvers/dgsem_p4est/dg_2d_gpu.jl | 2 +- src/solvers/dgsem_p4est/dg_2d_parallel.jl | 40 +-- src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl | 254 ++++++++++++++++++ src/solvers/dgsem_p4est/dg_parallel.jl | 78 +----- src/solvers/dgsem_tree/dg_2d.jl | 2 +- src/solvers/dgsem_tree/dg_2d_parallel.jl | 2 +- 6 files changed, 268 insertions(+), 110 deletions(-) create mode 100644 src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl diff --git a/src/solvers/dgsem_p4est/dg_2d_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_gpu.jl index 4e8ccf8d137..528efbbcd7f 100644 --- a/src/solvers/dgsem_p4est/dg_2d_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_gpu.jl @@ -7,7 +7,7 @@ function rhs_hyperbolic!(backend::Backend, du, u, t, - mesh::Union{P4estMeshSerial{2}, P4estMeshView{2}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}, P4estMesh{3}, T8codeMesh{3}}, equations, diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index 6e4ea8c5fb7..2f033058cc6 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -15,38 +15,14 @@ function prolong2mpiinterfaces!(backend::Nothing, cache, u, @threaded for interface in eachmpiinterface(dg, cache) prolong2mpiinterfaces_per_interface!(cache.mpi_interfaces.u, interface, - local_sides, local_neighbor_ids, node_indices, + local_sides, local_neighbor_ids, + node_indices, index_range, variables_range, u) end return nothing end -function prolong2mpiinterfaces!(backend::Backend, cache, u, - mesh::Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}, - equations, surface_integral, dg::DG) - nmpiinterfaces(dg, cache) == 0 && return nothing - @unpack local_sides, local_neighbor_ids, node_indices = cache.mpi_interfaces - index_range = eachnode(dg) - variables_range = eachvariable(equations) - - kernel! = prolong2mpiinterfaces_kernel!(backend) - kernel!(cache.mpi_interfaces.u, local_sides, local_neighbor_ids, node_indices, - index_range, variables_range, u, - ndrange = nmpiinterfaces(dg, cache)) - return nothing -end - -@kernel function prolong2mpiinterfaces_kernel!(mpi_interfaces_u, local_sides, - local_neighbor_ids, node_indices, - index_range, variables_range, u) - interface = @index(Global) - prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, interface, local_sides, - local_neighbor_ids, node_indices, index_range, - variables_range, u) -end - @inline function prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, interface, local_sides, local_neighbor_ids, node_indices, index_range, @@ -103,10 +79,14 @@ end @inline function calc_mpi_interface_flux_per_interface!(surface_flux_values, MeshT::Type{<:Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}}, - have_nonconservative_terms, equations, - surface_integral, SolverT::Type{<:DG}, - local_neighbor_ids, node_indices, local_sides, - contravariant_vectors, u_mpi_interfaces, + have_nonconservative_terms, + equations, + surface_integral, + SolverT::Type{<:DG}, + local_neighbor_ids, + node_indices, local_sides, + contravariant_vectors, + u_mpi_interfaces, index_range, interface) # Get element and side index information on the local element diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl new file mode 100644 index 00000000000..bc6060c3ed9 --- /dev/null +++ b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl @@ -0,0 +1,254 @@ +# 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 rhs_hyperbolic!(backend::Backend, + du, u, t, + mesh::Union{TreeMeshParallel{2}, P4estMeshParallel{2}, + T8codeMeshParallel{2}}, equations, + boundary_conditions, source_terms::Source, + dg::DG, cache) where {Source} + # Start to receive MPI data + @trixi_timeit timer() "start MPI receive" start_mpi_receive!(cache.mpi_cache) + + # Prolong solution to MPI interfaces + @trixi_timeit timer() "prolong2mpiinterfaces" begin + prolong2mpiinterfaces!(backend, cache, u, mesh, equations, dg.surface_integral, + dg) + end + + # Prolong solution to MPI mortars + @trixi_timeit timer() "prolong2mpimortars" begin + prolong2mpimortars!(cache, u, mesh, equations, + dg.mortar, dg) + end + + # Start to send MPI data + @trixi_timeit timer() "start MPI send" begin + start_mpi_send!(backend, cache.mpi_cache, mesh, equations, dg, cache) + end + + # Reset du + @trixi_timeit timer() "reset ∂u/∂t" set_zero!(du, dg, cache) + + # Calculate volume integral + @trixi_timeit timer() "volume integral" begin + calc_volume_integral!(backend, du, u, mesh, + have_nonconservative_terms(equations), equations, + dg.volume_integral, dg, cache) + end + + # Calculate interface fluxes + @trixi_timeit_ext backend timer() "prolong2interfaces + flux" begin + prolong2interfaces_and_calc_interface_flux!(backend, + cache.elements.surface_flux_values, + u, mesh, + have_nonconservative_terms(equations), + equations, + dg.surface_integral, dg, cache) + end + + # Prolong solution to boundaries + @trixi_timeit timer() "prolong2boundaries" begin + prolong2boundaries!(backend, cache, u, mesh, equations, dg) + end + + # Calculate boundary fluxes + @trixi_timeit timer() "boundary flux" begin + calc_boundary_flux!(backend, cache, t, boundary_conditions, mesh, equations, + dg.surface_integral, dg) + end + + # Prolong solution to mortars + @trixi_timeit timer() "prolong2mortars" begin + prolong2mortars!(cache, u, mesh, equations, + dg.mortar, dg) + end + + # Calculate mortar fluxes + @trixi_timeit timer() "mortar flux" begin + calc_mortar_flux!(cache.elements.surface_flux_values, mesh, + have_nonconservative_terms(equations), equations, + dg.mortar, dg.surface_integral, dg, cache) + end + + # Finish to receive MPI data + @trixi_timeit timer() "finish MPI receive" begin + finish_mpi_receive!(backend, cache.mpi_cache, mesh, equations, dg, cache) + end + + # Calculate MPI interface fluxes + @trixi_timeit timer() "MPI interface flux" begin + calc_mpi_interface_flux!(backend, cache.elements.surface_flux_values, mesh, + have_nonconservative_terms(equations), equations, + dg.surface_integral, dg, cache) + end + + # Calculate MPI mortar fluxes + @trixi_timeit timer() "MPI mortar flux" begin + calc_mpi_mortar_flux!(cache.elements.surface_flux_values, mesh, + have_nonconservative_terms(equations), equations, + dg.mortar, dg.surface_integral, dg, cache) + end + + # Calculate surface integrals, apply Jacobian from mapping to reference element + # and calculate source terms + @trixi_timeit_ext backend timer() "surface, Jacobian + source terms" begin + calc_surface_integral_and_apply_jacobian_and_calc_sources!(backend, du, u, t, + source_terms, mesh, + equations, + dg.surface_integral, + dg, cache) + end + + # Finish to send MPI data + @trixi_timeit timer() "finish MPI send" finish_mpi_send!(cache.mpi_cache) + + return nothing +end + +function prolong2mpiinterfaces!(backend::Backend, cache, u, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, + equations, surface_integral, dg::DG) + nmpiinterfaces(dg, cache) == 0 && return nothing + @unpack local_sides, local_neighbor_ids, node_indices = cache.mpi_interfaces + index_range = eachnode(dg) + variables_range = eachvariable(equations) + + kernel! = prolong2mpiinterfaces_KAkernel!(backend) + kernel!(cache.mpi_interfaces.u, local_sides, local_neighbor_ids, node_indices, + index_range, variables_range, u, + ndrange = nmpiinterfaces(dg, cache)) + return nothing +end + +@kernel function prolong2mpiinterfaces_KAkernel!(mpi_interfaces_u, local_sides, + local_neighbor_ids, node_indices, + index_range, variables_range, u) + interface = @index(Global) + prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, interface, local_sides, + local_neighbor_ids, node_indices, index_range, + variables_range, u) +end + +function calc_mpi_interface_flux!(backend::Backend, surface_flux_values, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, + have_nonconservative_terms, + equations, surface_integral, dg::DG, cache) + @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces + @unpack contravariant_vectors = cache.elements + @unpack u = cache.mpi_interfaces + index_range = eachnode(dg) + + kernel! = calc_mpi_interface_flux_KAkernel!(backend) + kernel!(surface_flux_values, typeof(mesh), + have_nonconservative_terms, equations, + surface_integral, typeof(dg), + local_neighbor_ids, node_indices, local_sides, contravariant_vectors, + u, index_range, + ndrange = nmpiinterfaces(dg, cache)) + return nothing +end + +@kernel function calc_mpi_interface_flux_KAkernel!(surface_flux_values, + MeshT::Type{<:Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}}, + have_nonconservative_terms, + equations, + surface_integral, + SolverT::Type{<:DG}, + local_neighbor_ids, node_indices, + local_sides, contravariant_vectors, + u, + index_range) + interface = @index(Global) + calc_mpi_interface_flux_per_interface!(surface_flux_values, MeshT, + have_nonconservative_terms, equations, + surface_integral, SolverT, + local_neighbor_ids, node_indices, + local_sides, contravariant_vectors, u, + index_range, interface) +end + +# TODO GPU: MPI mortars +function start_mpi_send!(backend::Backend, mpi_cache::P4estMPICache, + mesh::P4estMeshParallel{2}, equations, dg, cache) + @unpack mpi_neighbor_ranks, mpi_neighbor_interfaces = mpi_cache + @unpack mpi_send_buffers, mpi_send_requests = mpi_cache + @unpack local_sides, u = cache.mpi_interfaces + + kernel! = start_mpi_send_KAkernel!(backend) + + for (rank_index, neighbor_rank) in enumerate(mpi_neighbor_ranks) + send_buffer = mpi_send_buffers[rank_index] + neighbor_interfaces = mpi_neighbor_interfaces[rank_index] + kernel!(send_buffer, neighbor_interfaces, local_sides, u, + Val(nvariables(equations)), Val(ndims(mesh)), + ndrange = (nnodes(dg), length(neighbor_interfaces))) + + # wait for the kernel to return before sending the buffer + KernelAbstractions.synchronize(backend) + mpi_send_requests[rank_index] = MPI.Isend(send_buffer, neighbor_rank, + mpi_rank(), mpi_comm()) + end +end + +@kernel function start_mpi_send_KAkernel!(send_buffer, neighbor_interfaces, local_sides, + u_mpi_interfaces, ::Val{NVARS}, + ::Val{2}) where {NVARS} + index_node, index_interface = @index(Global, NTuple) + index_linear = @index(Global, Linear) + + buffer_offset = (index_linear - 1) * NVARS + interface = neighbor_interfaces[index_interface] + local_side = local_sides[interface] + + for v in 1:NVARS + send_buffer[buffer_offset + v] = u_mpi_interfaces[local_side, v, index_node, + index_interface] + end +end + +# TODO GPU: MPI mortars +function finish_mpi_receive!(backend::Backend, mpi_cache::P4estMPICache, + mesh::P4estMeshParallel{2}, equations, dg, cache) + @unpack mpi_neighbor_interfaces = mpi_cache + @unpack mpi_recv_buffers, mpi_recv_requests = mpi_cache + @unpack local_sides, u = cache.mpi_interfaces + + kernel! = finish_mpi_receive_KAkernel!(backend) + + # Start receiving and unpack received data until all communication is finished + data = MPI.Waitany(mpi_recv_requests) + while data !== nothing + recv_buffer = mpi_recv_buffers[data] + neighbor_interfaces = mpi_neighbor_interfaces[data] + kernel!(recv_buffer, neighbor_interfaces, local_sides, u, + Val(nvariables(equations)), Val(ndims(mesh)), + ndrange = (nnodes(dg), length(neighbor_interfaces))) + + data = MPI.Waitany(mpi_recv_requests) + end + # Wait for the last kernel to return ? + KernelAbstractions.synchronize(backend) +end + +@kernel function finish_mpi_receive_KAkernel!(recv_buffer, neighbor_interfaces, + local_sides, + u_mpi_interfaces, ::Val{NVARS}, + ::Val{2}) where {NVARS} + index_node, index_interface = @index(Global, NTuple) + index_linear = @index(Global, Linear) + buffer_offset = (index_linear - 1) * NVARS + interface = neighbor_interfaces[index_interface] + remote_side = local_sides[interface] == 1 ? 2 : 1 + for v in 1:NVARS + u_mpi_interfaces[remote_side, v, index_node, interface] = recv_buffer[buffer_offset + v] + end +end +end diff --git a/src/solvers/dgsem_p4est/dg_parallel.jl b/src/solvers/dgsem_p4est/dg_parallel.jl index 87f8f29a0da..f456576bf5a 100644 --- a/src/solvers/dgsem_p4est/dg_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_parallel.jl @@ -140,45 +140,6 @@ function start_mpi_send!(backend::Nothing, mpi_cache::P4estMPICache, mesh, equat return nothing end -# TODO GPU: MPI mortars -function start_mpi_send!(backend::Backend, mpi_cache::P4estMPICache, - mesh::P4estMeshParallel{2}, equations, dg, cache) - @unpack mpi_neighbor_ranks, mpi_neighbor_interfaces = mpi_cache - @unpack mpi_send_buffers, mpi_send_requests = mpi_cache - @unpack local_sides, u = cache.mpi_interfaces - - kernel! = start_mpi_send_KAkernel!(backend) - - for (rank_index, neighbor_rank) in enumerate(mpi_neighbor_ranks) - send_buffer = mpi_send_buffers[rank_index] - neighbor_interfaces = mpi_neighbor_interfaces[rank_index] - kernel!(send_buffer, neighbor_interfaces, local_sides, u, - Val(nvariables(equations)), Val(ndims(mesh)), - ndrange = (nnodes(dg), length(neighbor_interfaces))) - - # wait for the kernel to return before sending the buffer - KernelAbstractions.synchronize(backend) - mpi_send_requests[rank_index] = MPI.Isend(send_buffer, neighbor_rank, - mpi_rank(), mpi_comm()) - end -end - -@kernel function start_mpi_send_KAkernel!(send_buffer, neighbor_interfaces, local_sides, - u_mpi_interfaces, ::Val{NVARS}, - ::Val{2}) where {NVARS} - index_node, index_interface = @index(Global, NTuple) - index_linear = @index(Global, Linear) - - buffer_offset = (index_linear - 1) * NVARS - interface = neighbor_interfaces[index_interface] - local_side = local_sides[interface] - - for v in 1:NVARS - send_buffer[buffer_offset + v] = u_mpi_interfaces[local_side, v, index_node, - index_interface] - end -end - function start_mpi_receive!(mpi_cache::P4estMPICache) for (index, rank) in enumerate(mpi_cache.mpi_neighbor_ranks) mpi_cache.mpi_recv_requests[index] = MPI.Irecv!(mpi_cache.mpi_recv_buffers[index], @@ -243,44 +204,6 @@ function finish_mpi_receive!(backend::Nothing, mpi_cache::P4estMPICache, mesh, return nothing end -# TODO GPU: MPI mortars -function finish_mpi_receive!(backend::Backend, mpi_cache::P4estMPICache, - mesh::P4estMeshParallel{2}, equations, dg, cache) - @unpack mpi_neighbor_interfaces = mpi_cache - @unpack mpi_recv_buffers, mpi_recv_requests = mpi_cache - @unpack local_sides, u = cache.mpi_interfaces - - kernel! = finish_mpi_receive_KAkernel!(backend) - - # Start receiving and unpack received data until all communication is finished - data = MPI.Waitany(mpi_recv_requests) - while data !== nothing - recv_buffer = mpi_recv_buffers[data] - neighbor_interfaces = mpi_neighbor_interfaces[data] - kernel!(recv_buffer, neighbor_interfaces, local_sides, u, - Val(nvariables(equations)), Val(ndims(mesh)), - ndrange = (nnodes(dg), length(neighbor_interfaces))) - - data = MPI.Waitany(mpi_recv_requests) - end - # Wait for the last kernel to return ? - KernelAbstractions.synchronize(backend) -end - -@kernel function finish_mpi_receive_KAkernel!(recv_buffer, neighbor_interfaces, - local_sides, - u_mpi_interfaces, ::Val{NVARS}, - ::Val{3}) where {NVARS} - index_node, index_interface = @index(Global, NTuple) - index_linear = @index(Global, Linear) - buffer_offset = (index_linear - 1) * NVARS - interface = neighbor_interfaces[index_interface] - remote_side = local_sides[interface] == 1 ? 2 : 1 - for v in 1:NVARS - u_mpi_interfaces[remote_side, v, index_node, interface] = recv_buffer[buffer_offset + v] - end -end - # Return a tuple `indices` where indices[position] is a `(first, last)` tuple for accessing the # data corresponding to the `position` part of a mortar in an MPI buffer. The mortar data must begin # at `index_base`+1 in the MPI buffer. `data_size` is the data size associated with each small @@ -707,6 +630,7 @@ end end include("dg_2d_parallel.jl") +include("dg_2d_parallel_gpu.jl") include("dg_3d_parallel.jl") include("dg_2d_parabolic_parallel.jl") end # muladd diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 14af40e6a4c..9022f9100a8 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -112,7 +112,7 @@ end # 2D/3D `TreeMesh`es, `P4estMesh`es, and `T8codeMesh`es. function rhs_hyperbolic!(backend::Nothing, du, u, t, - mesh::Union{TreeMesh{2}, P4estMeshSerial{2}, P4estMeshView{2}, + mesh::Union{TreeMesh{2}, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}, TreeMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations, diff --git a/src/solvers/dgsem_tree/dg_2d_parallel.jl b/src/solvers/dgsem_tree/dg_2d_parallel.jl index c1e5655e417..e1b5a61d099 100644 --- a/src/solvers/dgsem_tree/dg_2d_parallel.jl +++ b/src/solvers/dgsem_tree/dg_2d_parallel.jl @@ -450,7 +450,7 @@ function init_mpi_neighbor_connectivity(elements, mpi_interfaces, mpi_mortars, return mpi_neighbor_ranks, mpi_neighbor_interfaces, mpi_neighbor_mortars end -function rhs_hyperbolic!(backend, +function rhs_hyperbolic!(backend::Nothing, du, u, t, mesh::Union{TreeMeshParallel{2}, P4estMeshParallel{2}, T8codeMeshParallel{2}}, equations, From 1f9d981ef303f044fc5e47eda29710583827e8bb Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 15 Sep 2026 10:05:11 +0200 Subject: [PATCH 07/46] remove reset du --- src/solvers/dgsem_p4est/dg_2d_gpu.jl | 4 ++-- src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl | 3 --- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_gpu.jl index 528efbbcd7f..81e43a8849c 100644 --- a/src/solvers/dgsem_p4est/dg_2d_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_gpu.jl @@ -7,12 +7,12 @@ function rhs_hyperbolic!(backend::Backend, du, u, t, - mesh::Union{P4estMesh{2}, P4estMeshView{2}, - T8codeMesh{2}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}, P4estMesh{3}, T8codeMesh{3}}, equations, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} + # Calculate volume integral @trixi_timeit_ext backend timer() "volume integral" begin calc_volume_integral!(backend, du, u, mesh, diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl index bc6060c3ed9..7568a2e1d5e 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl @@ -31,9 +31,6 @@ function rhs_hyperbolic!(backend::Backend, start_mpi_send!(backend, cache.mpi_cache, mesh, equations, dg, cache) end - # Reset du - @trixi_timeit timer() "reset ∂u/∂t" set_zero!(du, dg, cache) - # Calculate volume integral @trixi_timeit timer() "volume integral" begin calc_volume_integral!(backend, du, u, mesh, From 7f1ba04efe6e82e85bf3c9a9a89bda8fd5e37230 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 15 Sep 2026 12:48:29 +0200 Subject: [PATCH 08/46] fmt --- src/callbacks_step/analysis_dg2d_parallel.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d_parallel.jl b/src/callbacks_step/analysis_dg2d_parallel.jl index 41c9ba06e32..b2bbdbd72ff 100644 --- a/src/callbacks_step/analysis_dg2d_parallel.jl +++ b/src/callbacks_step/analysis_dg2d_parallel.jl @@ -223,7 +223,6 @@ function integrate_via_indices(func::Func, backend::Nothing, u, return integral end - function integrate_via_indices(func::Func, backend::Backend, u, mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, equations, @@ -246,8 +245,8 @@ function integrate_via_indices(func::Func, backend::Backend, u, # Use quadrature to numerically integrate over entire domain num_elements = nelements(dg, cache) integral, volume = AcceleratedKernels.mapreduce(local_plus, 1:num_elements, - backend; init, - neutral) do element + backend; init, + neutral) do element # Initialize integral with zeros of the right shapeu local_integral, local_total_volume = neutral From a48f151d839c303eedb667a2f4fbe01d769d37b2 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 15 Sep 2026 15:16:16 +0200 Subject: [PATCH 09/46] fix bug --- src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl index 7568a2e1d5e..f1c7e032e68 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl @@ -207,7 +207,7 @@ end for v in 1:NVARS send_buffer[buffer_offset + v] = u_mpi_interfaces[local_side, v, index_node, - index_interface] + interface] end end From 5fcbd10787376601f0561f3b7cc7ebb9cd2964e6 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 15 Sep 2026 15:16:56 +0200 Subject: [PATCH 10/46] use trixi_timeit_ext --- src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl | 34 +++++++++---------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl index f1c7e032e68..6c59316676a 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl @@ -7,32 +7,32 @@ function rhs_hyperbolic!(backend::Backend, du, u, t, - mesh::Union{TreeMeshParallel{2}, P4estMeshParallel{2}, - T8codeMeshParallel{2}}, equations, + mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, + equations, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} # Start to receive MPI data @trixi_timeit timer() "start MPI receive" start_mpi_receive!(cache.mpi_cache) # Prolong solution to MPI interfaces - @trixi_timeit timer() "prolong2mpiinterfaces" begin + @trixi_timeit_ext backend timer() "prolong2mpiinterfaces" begin prolong2mpiinterfaces!(backend, cache, u, mesh, equations, dg.surface_integral, dg) end # Prolong solution to MPI mortars - @trixi_timeit timer() "prolong2mpimortars" begin - prolong2mpimortars!(cache, u, mesh, equations, + @trixi_timeit_ext backend timer() "prolong2mpimortars" begin + prolong2mpimortars!(backend, cache, u, mesh, equations, dg.mortar, dg) end # Start to send MPI data - @trixi_timeit timer() "start MPI send" begin + @trixi_timeit_ext backend timer() "start MPI send" begin start_mpi_send!(backend, cache.mpi_cache, mesh, equations, dg, cache) end # Calculate volume integral - @trixi_timeit timer() "volume integral" begin + @trixi_timeit_ext backend timer() "volume integral" begin calc_volume_integral!(backend, du, u, mesh, have_nonconservative_terms(equations), equations, dg.volume_integral, dg, cache) @@ -49,44 +49,44 @@ function rhs_hyperbolic!(backend::Backend, end # Prolong solution to boundaries - @trixi_timeit timer() "prolong2boundaries" begin + @trixi_timeit_ext backend timer() "prolong2boundaries" begin prolong2boundaries!(backend, cache, u, mesh, equations, dg) end # Calculate boundary fluxes - @trixi_timeit timer() "boundary flux" begin + @trixi_timeit_ext backend timer() "boundary flux" begin calc_boundary_flux!(backend, cache, t, boundary_conditions, mesh, equations, dg.surface_integral, dg) end # Prolong solution to mortars - @trixi_timeit timer() "prolong2mortars" begin - prolong2mortars!(cache, u, mesh, equations, + @trixi_timeit_ext backend timer() "prolong2mortars" begin + prolong2mortars!(backend, cache, u, mesh, equations, dg.mortar, dg) end # Calculate mortar fluxes - @trixi_timeit timer() "mortar flux" begin - calc_mortar_flux!(cache.elements.surface_flux_values, mesh, + @trixi_timeit_ext backend timer() "mortar flux" begin + calc_mortar_flux!(backend, cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, dg.mortar, dg.surface_integral, dg, cache) end # Finish to receive MPI data - @trixi_timeit timer() "finish MPI receive" begin + @trixi_timeit_ext backend timer() "finish MPI receive" begin finish_mpi_receive!(backend, cache.mpi_cache, mesh, equations, dg, cache) end # Calculate MPI interface fluxes - @trixi_timeit timer() "MPI interface flux" begin + @trixi_timeit_ext backend timer() "MPI interface flux" begin calc_mpi_interface_flux!(backend, cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, dg.surface_integral, dg, cache) end # Calculate MPI mortar fluxes - @trixi_timeit timer() "MPI mortar flux" begin - calc_mpi_mortar_flux!(cache.elements.surface_flux_values, mesh, + @trixi_timeit_ext backend timer() "MPI mortar flux" begin + calc_mpi_mortar_flux!(backend, cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, dg.mortar, dg.surface_integral, dg, cache) end From 18df67e92642edc093087c88d11e28a31df757d9 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 15 Sep 2026 15:17:47 +0200 Subject: [PATCH 11/46] avoid kernel launches with zero size --- src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl | 58 ++++++++++++++----- 1 file changed, 45 insertions(+), 13 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl index 6c59316676a..75f9623896b 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl @@ -137,6 +137,7 @@ function calc_mpi_interface_flux!(backend::Backend, surface_flux_values, T8codeMeshParallel{2}}, have_nonconservative_terms, equations, surface_integral, dg::DG, cache) + nmpiinterfaces(dg, cache) == 0 && return nothing @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces @unpack contravariant_vectors = cache.elements @unpack u = cache.mpi_interfaces @@ -172,24 +173,50 @@ end index_range, interface) end -# TODO GPU: MPI mortars +# For GPU backends MPI mortars are not yet implemented +function prolong2mpimortars!(backend::Backend, cache, u, + mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, + equations, mortar, dg::DG) + @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" + return nothing +end + +# For GPU backends MPI mortars are not yet implemented +function calc_mpi_mortar_flux!(backend::Backend, surface_flux_values, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, + have_nonconservative_terms, equations, mortar, + surface_integral, dg::DG, cache) + @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" + return nothing +end + function start_mpi_send!(backend::Backend, mpi_cache::P4estMPICache, - mesh::P4estMeshParallel{2}, equations, dg, cache) + mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, + equations, dg, cache) + # TODO GPU: MPI mortars are not yet implemented + @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" + @unpack mpi_neighbor_ranks, mpi_neighbor_interfaces = mpi_cache @unpack mpi_send_buffers, mpi_send_requests = mpi_cache @unpack local_sides, u = cache.mpi_interfaces kernel! = start_mpi_send_KAkernel!(backend) - for (rank_index, neighbor_rank) in enumerate(mpi_neighbor_ranks) + for (rank_index, _) in enumerate(mpi_neighbor_ranks) send_buffer = mpi_send_buffers[rank_index] neighbor_interfaces = mpi_neighbor_interfaces[rank_index] - kernel!(send_buffer, neighbor_interfaces, local_sides, u, - Val(nvariables(equations)), Val(ndims(mesh)), - ndrange = (nnodes(dg), length(neighbor_interfaces))) + if !isempty(neighbor_interfaces) + kernel!(send_buffer, neighbor_interfaces, local_sides, u, + Val(nvariables(equations)), Val(ndims(mesh)), + ndrange = (nnodes(dg), length(neighbor_interfaces))) + end + end - # wait for the kernel to return before sending the buffer - KernelAbstractions.synchronize(backend) + # wait for the kernels to return before sending the buffers + KernelAbstractions.synchronize(backend) + + for (rank_index, neighbor_rank) in enumerate(mpi_neighbor_ranks) mpi_send_requests[rank_index] = MPI.Isend(send_buffer, neighbor_rank, mpi_rank(), mpi_comm()) end @@ -211,9 +238,12 @@ end end end -# TODO GPU: MPI mortars function finish_mpi_receive!(backend::Backend, mpi_cache::P4estMPICache, - mesh::P4estMeshParallel{2}, equations, dg, cache) + mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, + equations, dg, cache) + # TODO GPU: MPI mortars are not yet implemented + @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" + @unpack mpi_neighbor_interfaces = mpi_cache @unpack mpi_recv_buffers, mpi_recv_requests = mpi_cache @unpack local_sides, u = cache.mpi_interfaces @@ -225,9 +255,11 @@ function finish_mpi_receive!(backend::Backend, mpi_cache::P4estMPICache, while data !== nothing recv_buffer = mpi_recv_buffers[data] neighbor_interfaces = mpi_neighbor_interfaces[data] - kernel!(recv_buffer, neighbor_interfaces, local_sides, u, - Val(nvariables(equations)), Val(ndims(mesh)), - ndrange = (nnodes(dg), length(neighbor_interfaces))) + if !isempty(neighbor_interfaces) + kernel!(recv_buffer, neighbor_interfaces, local_sides, u, + Val(nvariables(equations)), Val(ndims(mesh)), + ndrange = (nnodes(dg), length(neighbor_interfaces))) + end data = MPI.Waitany(mpi_recv_requests) end From cc61eba56b9ec413c80075fb65aefc7e93769104 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 15 Sep 2026 15:24:03 +0200 Subject: [PATCH 12/46] missed line --- src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl index 75f9623896b..b99f36de141 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl @@ -217,6 +217,7 @@ function start_mpi_send!(backend::Backend, mpi_cache::P4estMPICache, KernelAbstractions.synchronize(backend) for (rank_index, neighbor_rank) in enumerate(mpi_neighbor_ranks) + send_buffer = mpi_send_buffers[rank_index] mpi_send_requests[rank_index] = MPI.Isend(send_buffer, neighbor_rank, mpi_rank(), mpi_comm()) end From 84abfff059ef985aaab272456e6122bc3f51be16 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 15 Sep 2026 16:47:41 +0200 Subject: [PATCH 13/46] removed parallel meshes should already be included --- src/callbacks_step/analysis_dg2d.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index da789d54187..22331ff4c3e 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -339,9 +339,7 @@ function integrate_via_indices(func::Func, u, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, - T8codeMesh{2}, - P4estMeshParallel{2}, - T8codeMeshParallel{2}}, + T8codeMesh{2}}, equations, dg::DGSEM, cache, args...; normalize = true) where {Func} return integrate_via_indices(func, trixi_backend(u), u, mesh, equations, dg, cache, From 8ea605b051e13369cace4d1435eda4bfd031d66a Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 15 Sep 2026 16:48:34 +0200 Subject: [PATCH 14/46] try to set up test with mpi and kernelabstractions --- .github/workflows/ci.yml | 1 + test/runtests.jl | 16 +++++++++++ test/test_kernelabstractions_mpi.jl | 42 +++++++++++++++++++++++++++++ 3 files changed, 59 insertions(+) create mode 100644 test/test_kernelabstractions_mpi.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 409c6b19fee..af1665a28b0 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -85,6 +85,7 @@ jobs: - mpi - threaded - kernelabstractions + - kernelabstractions_mpi include: - version: '1.11' os: ubuntu-latest diff --git a/test/runtests.jl b/test/runtests.jl index 0fef08091d8..1229501e28b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -70,6 +70,10 @@ end if "kernelabstractions" in SUITES_TO_DISPATCH import Trixi end +if "kernelabstractions_mpi" in SUITES_TO_DISPATCH + import MPI + import Trixi +end # The GPU suites run their items only when the respective backend is `functional()`, # so requesting `CUDA`/`AMDGPU` on a machine without that GPU warns and runs nothing @@ -115,6 +119,18 @@ function dispatch_special_suite(suite) finally Trixi.set_threading_backend!(Symbol(previous_backend)) end + elseif suite == "kernelabstractions_mpi" + # The threading backend is selected via a preference that is read on Julia + # startup, so we set it here (in the parent) and restore it afterwards; + # the relaunched worker picks it up. + previous_backend = Trixi._PREFERENCE_THREADING + Trixi.set_threading_backend!(:kernelabstractions) + try + cmd = `$(MPI.mpiexec()) -n $TRIXI_MPI_NPROCS $julia --threads=1 --check-bounds=yes --heap-size-hint=0.5G --project=$project $(@__FILE__)` + run_worker(cmd, suite) + finally + Trixi.set_threading_backend!(Symbol(previous_backend)) + end else error("We should not reach this branch; something is wrong.") end diff --git a/test/test_kernelabstractions_mpi.jl b/test/test_kernelabstractions_mpi.jl new file mode 100644 index 00000000000..da5bb1cdf83 --- /dev/null +++ b/test/test_kernelabstractions_mpi.jl @@ -0,0 +1,42 @@ +@testsnippet KernelAbstractionsMPIExamples begin + EXAMPLES_DIR = examples_dir() +end + +@testitem "KernelAbstractions MPI backend preference" setup=[Setup] tags=[:kernelabstractions_mpi] begin + @test Trixi._PREFERENCE_THREADING == :kernelabstractions +end + +@testitem "KernelAbstractions MPI GPU 2D: elixir_advection_basic.jl" setup=[ + Setup, + KernelAbstractionsMPIExamples +] tags=[:kernelabstractions_mpi] begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", + "elixir_advection_basic.jl"), + # Expected errors are exactly the same as with TreeMesh! + l2=8.311947673061856e-6, + linf=6.627000273229378e-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 real(semi.solver) == Float64 + @test real(semi.solver.basis) == Float64 + @test real(semi.solver.mortar) == Float64 + # TODO: `mesh` is currently not `adapt`ed correctly + @test real(semi.mesh) == Float64 + @test typeof(semi.equations.advection_velocity) == SVector{2, Float64} + + @test ode.u0 isa Array + @test semi.solver.basis.derivative_matrix isa Array + + @test Trixi.storage_type(semi.cache.elements) === Array + @test Trixi.storage_type(semi.cache.interfaces) === Array + @test Trixi.storage_type(semi.cache.boundaries) === Array + @test Trixi.storage_type(semi.cache.mortars) === Array + + # Ensure that the RHS computation overwrites existing data in `du` correctly. + u_ode = copy(ode.u0) + du_ode = similar(u_ode) + fill!(du_ode, convert(eltype(du_ode), NaN)) + Trixi.rhs_hyperbolic!(du_ode, u_ode, ode.p, first(ode.tspan)) + @test all(isfinite, du_ode) +end From d317322aac19bac7e032f116eb973ad6ec0797c4 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Fri, 18 Sep 2026 08:11:13 +0200 Subject: [PATCH 15/46] add kernelabstractions_mpi to SPECIAL_PROCESS_SUITES --- test/runtests.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 1229501e28b..fa0991861e7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -28,7 +28,8 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) # also works for every MPI rank). The `TRIXI_TEST` value of such a suite selects # the items via the equally-named tag; `threaded_legacy` reuses the `threaded` # items but is launched on a different Julia version by CI. -const SPECIAL_PROCESS_SUITES = ("mpi", "threaded", "kernelabstractions") +const SPECIAL_PROCESS_SUITES = ("mpi", "threaded", "kernelabstractions", + "kernelabstractions_mp") const IN_WORKER = haskey(ENV, "TRIXI_TEST_RUN_ITEMS") # Remove Trixi's output directory `out`, where examples write solution/restart/mesh From 0d31cc5091e1e7e3799722c2c114a5d61864bcf0 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Fri, 18 Sep 2026 08:12:28 +0200 Subject: [PATCH 16/46] backend arguments --- src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl | 12 +++++++----- src/solvers/dgsem_p4est/dg_2d_parallel.jl | 2 +- src/solvers/dgsem_p4est/dg_3d_parallel.jl | 4 ++-- src/solvers/dgsem_tree/dg_2d_parallel.jl | 12 +++++++----- 4 files changed, 17 insertions(+), 13 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl index eac9b2f0565..dc8825bb74c 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -28,7 +28,7 @@ function rhs_parabolic!(backend::Nothing, du, u, t, # Prolong transformed variables to MPI mortars @trixi_timeit timer() "prolong2mpimortars gradient" begin - prolong2mpimortars!(cache, u_transformed, mesh, equations_parabolic, + prolong2mpimortars!(backend, cache, u_transformed, mesh, equations_parabolic, dg.mortar, dg) end @@ -40,7 +40,7 @@ function rhs_parabolic!(backend::Nothing, du, u, t, # Start gradient MPI send @trixi_timeit timer() "start MPI send gradient" begin - start_mpi_send!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) + start_mpi_send!(backend, cache.mpi_cache, mesh, equations_parabolic, dg, cache) end # Local gradient computation @@ -52,7 +52,8 @@ function rhs_parabolic!(backend::Nothing, du, u, t, # Finish gradient MPI receive @trixi_timeit timer() "finish MPI receive gradient" begin - finish_mpi_receive!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) + finish_mpi_receive!(backend, cache.mpi_cache, mesh, equations_parabolic, dg, + cache) end # MPI interface fluxes for gradients @@ -122,7 +123,7 @@ function rhs_parabolic!(backend::Nothing, du, u, t, # Start divergence MPI send @trixi_timeit timer() "start MPI send divergence" begin - start_mpi_send!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) + start_mpi_send!(backend, cache.mpi_cache, mesh, equations_parabolic, dg, cache) end # Local, i.e., non-MPI interface/boundary/mortar fluxes @@ -137,7 +138,8 @@ function rhs_parabolic!(backend::Nothing, du, u, t, # Finish divergence MPI receive @trixi_timeit timer() "finish MPI receive divergence" begin - finish_mpi_receive!(cache.mpi_cache, mesh, equations_parabolic, dg, cache) + finish_mpi_receive!(backend, cache.mpi_cache, mesh, equations_parabolic, dg, + cache) end # MPI interface fluxes for divergence diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index 2f033058cc6..dd6a3bcb9a0 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -204,7 +204,7 @@ end return nothing end -function prolong2mpimortars!(cache, u, +function prolong2mpimortars!(backend::Nothing, cache, u, mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, equations, mortar_l2::LobattoLegendreMortarL2, diff --git a/src/solvers/dgsem_p4est/dg_3d_parallel.jl b/src/solvers/dgsem_p4est/dg_3d_parallel.jl index 54667502a3b..d97c476b06e 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parallel.jl @@ -27,7 +27,7 @@ function rhs_hyperbolic!(backend::Nothing, # Start to send MPI data @trixi_timeit timer() "start MPI send" begin - start_mpi_send!(cache.mpi_cache, mesh, equations, dg, cache) + start_mpi_send!(backend, cache.mpi_cache, mesh, equations, dg, cache) end # Reset du @@ -78,7 +78,7 @@ function rhs_hyperbolic!(backend::Nothing, # Finish to receive MPI data @trixi_timeit timer() "finish MPI receive" begin - finish_mpi_receive!(cache.mpi_cache, mesh, equations, dg, cache) + finish_mpi_receive!(backend, cache.mpi_cache, mesh, equations, dg, cache) end # Calculate MPI interface fluxes diff --git a/src/solvers/dgsem_tree/dg_2d_parallel.jl b/src/solvers/dgsem_tree/dg_2d_parallel.jl index e1b5a61d099..dcd208efd7d 100644 --- a/src/solvers/dgsem_tree/dg_2d_parallel.jl +++ b/src/solvers/dgsem_tree/dg_2d_parallel.jl @@ -58,7 +58,8 @@ function start_mpi_receive!(mpi_cache::MPICache) end # TODO: MPI dimension agnostic -function start_mpi_send!(mpi_cache::MPICache, mesh, equations, dg, cache) +function start_mpi_send!(backend::Nothing, mpi_cache::MPICache, mesh, equations, dg, + cache) data_size = nvariables(equations) * nnodes(dg)^(ndims(mesh) - 1) for rank in 1:length(mpi_cache.mpi_neighbor_ranks) @@ -159,7 +160,8 @@ function finish_mpi_send!(mpi_cache::MPICache) end # TODO: MPI dimension agnostic -function finish_mpi_receive!(mpi_cache::MPICache, mesh, equations, dg, cache) +function finish_mpi_receive!(backend::Nothing, mpi_cache::MPICache, mesh, equations, dg, + cache) data_size = nvariables(equations) * nnodes(dg)^(ndims(mesh) - 1) # Start receiving and unpack received data until all communication is finished @@ -467,7 +469,7 @@ function rhs_hyperbolic!(backend::Nothing, # Prolong solution to MPI mortars @trixi_timeit timer() "prolong2mpimortars" begin - prolong2mpimortars!(cache, u, mesh, equations, + prolong2mpimortars!(backend, cache, u, mesh, equations, dg.mortar, dg) end @@ -599,7 +601,7 @@ function prolong2mpiinterfaces!(cache, u, return nothing end -function prolong2mpimortars!(cache, u, +function prolong2mpimortars!(backend::Nothing, cache, u, mesh::TreeMeshParallel{2}, equations, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM) @@ -725,7 +727,7 @@ function prolong2mpimortars!(cache, u, return nothing end -function calc_mpi_interface_flux!(surface_flux_values, +function calc_mpi_interface_flux!(backend::Nothing, surface_flux_values, mesh::TreeMeshParallel{2}, have_nonconservative_terms::False, equations, surface_integral, dg::DG, cache) From 429de3d07346f215c2522f73a9c52d05e405e81b Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 21 Sep 2026 12:50:57 +0200 Subject: [PATCH 17/46] typo --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index fa0991861e7..74090300e81 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -29,7 +29,7 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) # the items via the equally-named tag; `threaded_legacy` reuses the `threaded` # items but is launched on a different Julia version by CI. const SPECIAL_PROCESS_SUITES = ("mpi", "threaded", "kernelabstractions", - "kernelabstractions_mp") + "kernelabstractions_mpi") const IN_WORKER = haskey(ENV, "TRIXI_TEST_RUN_ITEMS") # Remove Trixi's output directory `out`, where examples write solution/restart/mesh From cbccef79f417df9873827c48043dfe5178851be3 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 21 Sep 2026 13:20:24 +0200 Subject: [PATCH 18/46] another backend argument --- src/solvers/dgsem_tree/dg_2d_parallel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_tree/dg_2d_parallel.jl b/src/solvers/dgsem_tree/dg_2d_parallel.jl index dcd208efd7d..985f257d70b 100644 --- a/src/solvers/dgsem_tree/dg_2d_parallel.jl +++ b/src/solvers/dgsem_tree/dg_2d_parallel.jl @@ -565,7 +565,7 @@ function rhs_hyperbolic!(backend::Nothing, return nothing end -function prolong2mpiinterfaces!(cache, u, +function prolong2mpiinterfaces!(backend::Nothing, cache, u, mesh::TreeMeshParallel{2}, equations, surface_integral, dg::DG) @unpack mpi_interfaces = cache From daaa8983cb58d8a52d0a73bfb528e5bcf1854d20 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 21 Sep 2026 16:36:13 +0200 Subject: [PATCH 19/46] another one --- src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl index dc8825bb74c..76023f25edf 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic_parallel.jl @@ -34,7 +34,7 @@ function rhs_parabolic!(backend::Nothing, du, u, t, # Prolong transformed variables to MPI interfaces @trixi_timeit timer() "prolong2mpiinterfaces gradient" begin - prolong2mpiinterfaces!(cache, u_transformed, mesh, equations_parabolic, + prolong2mpiinterfaces!(backend, cache, u_transformed, mesh, equations_parabolic, dg.surface_integral, dg) end From 36b923540a0db210766e0fde992577dd8f764adf Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 22 Sep 2026 22:35:10 +0200 Subject: [PATCH 20/46] added NEWS.md --- NEWS.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/NEWS.md b/NEWS.md index dbfa419fe3a..701aa764c96 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,9 @@ for human readability. ## Changes in the v0.17 lifecycle #### Added +- Added multi-GPU support for 2D P4estMesh. A local CUDA-aware MPI installation is required. + For configuration see [Parallelization](https://trixi-framework.org/TrixiDocumentation/stable/parallelization/#parallel_system_MPI) + in the docs. - TimerOutputs.jl v1 is now supported in addition to v0.5 ([#3172]). When TimerOutputs.jl v1 is used, the new preference `Trixi.set_timer_bars!` toggles the bars visualizing the fraction of time and allocations spent in From 5200c4011c9a19573bd7d1a296652c7c7e58976a Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Wed, 23 Sep 2026 10:26:21 +0200 Subject: [PATCH 21/46] unpack --- src/solvers/dgsem_p4est/dg_2d_parallel.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index 2f033058cc6..b326c557595 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -10,11 +10,12 @@ function prolong2mpiinterfaces!(backend::Nothing, cache, u, T8codeMeshParallel{2}}, equations, surface_integral, dg::DG) @unpack local_sides, local_neighbor_ids, node_indices = cache.mpi_interfaces + mpi_interfaces_u = cache.mpi_interfaces.u index_range = eachnode(dg) variables_range = eachvariable(equations) @threaded for interface in eachmpiinterface(dg, cache) - prolong2mpiinterfaces_per_interface!(cache.mpi_interfaces.u, interface, + prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, interface, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u) From 596ad550d264ba7825de9d82e43687802a4ad6d4 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Fri, 2 Oct 2026 11:52:54 +0200 Subject: [PATCH 22/46] fmt --- src/solvers/dgsem_p4est/dg_2d_parallel.jl | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index 98d7ccd8b38..4a60eab76ee 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -45,16 +45,17 @@ end local_indices = node_indices[interface] i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], - index_range) + index_range) j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], - index_range) + index_range) i_element = i_element_start j_element = j_element_start for i in index_range for v in variables_range - mpi_interfaces_u[local_side, v, i, interface] = u[v, i_element, j_element, - local_element] + mpi_interfaces_u[local_side, v, i, interface] = u[v, i_element, + j_element, + local_element] end i_element += i_element_step j_element += j_element_step @@ -113,9 +114,9 @@ end # Create the local i,j indexing on the local element used to pull normal direction information i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], - index_range) + index_range) j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], - index_range) + index_range) i_element = i_element_start j_element = j_element_start @@ -139,11 +140,11 @@ end i_element, j_element, local_element) calc_mpi_interface_flux!(surface_flux_values, MeshT, - have_nonconservative_terms, equations, - surface_integral, SolverT, - u_mpi_interfaces, interface, normal_direction, - node, local_side, surface_node, local_direction, - local_element) + have_nonconservative_terms, equations, + surface_integral, SolverT, + u_mpi_interfaces, interface, normal_direction, + node, local_side, surface_node, local_direction, + local_element) # Increment local element indices to pull the normal direction i_element += i_element_step From ae83c141672a6504170255957658f759750537ab Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 12 Oct 2026 15:27:42 +0200 Subject: [PATCH 23/46] missed in merge --- NEWS.md | 3 --- 1 file changed, 3 deletions(-) diff --git a/NEWS.md b/NEWS.md index 6567bdcfdd8..c1d9befb338 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,11 +8,9 @@ for human readability. ## Changes in the v0.17 lifecycle #### Added -<<<<<<< HEAD - Added multi-GPU support for 2D P4estMesh. A local CUDA-aware MPI installation is required. For configuration see [Parallelization](https://trixi-framework.org/TrixiDocumentation/stable/parallelization/#parallel_system_MPI) in the docs. -======= - `semidiscretize` of a `SemidiscretizationHyperbolicParabolic` supports the keyword arguments `storage_type`, `real_type`, and `flux_differencing_kernel` for the `P4estMesh` ([#3313]). @@ -20,7 +18,6 @@ for human readability. `max_abs_speed` are now available for `PassiveTracerEquations` ([#3319]). - New entropy-stable surface flux `FluxTracerEquationsUpwind` for `PassiveTracerEquations`, which upwinds the tracers with the mass flux of an entropy-stable flow flux ([#3320]). ->>>>>>> main - Contour plotting using Makie.jl is now supported for `PlotData2DCartesian` data ([#3238]). - TimerOutputs.jl v1 is now supported in addition to v0.5 ([#3172]). When TimerOutputs.jl v1 is used, the new preference `Trixi.set_timer_bars!` From 0a6378b1b7168307bb9913c299e064de6c569bb5 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 12 Oct 2026 15:28:24 +0200 Subject: [PATCH 24/46] add inbounds=true to kernels --- src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl index b99f36de141..b6e7115bb7f 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl @@ -123,7 +123,7 @@ function prolong2mpiinterfaces!(backend::Backend, cache, u, return nothing end -@kernel function prolong2mpiinterfaces_KAkernel!(mpi_interfaces_u, local_sides, +@kernel inbounds=true function prolong2mpiinterfaces_KAkernel!(mpi_interfaces_u, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u) interface = @index(Global) @@ -153,7 +153,7 @@ function calc_mpi_interface_flux!(backend::Backend, surface_flux_values, return nothing end -@kernel function calc_mpi_interface_flux_KAkernel!(surface_flux_values, +@kernel inbounds=true function calc_mpi_interface_flux_KAkernel!(surface_flux_values, MeshT::Type{<:Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}}, have_nonconservative_terms, @@ -223,7 +223,7 @@ function start_mpi_send!(backend::Backend, mpi_cache::P4estMPICache, end end -@kernel function start_mpi_send_KAkernel!(send_buffer, neighbor_interfaces, local_sides, +@kernel inbounds=true function start_mpi_send_KAkernel!(send_buffer, neighbor_interfaces, local_sides, u_mpi_interfaces, ::Val{NVARS}, ::Val{2}) where {NVARS} index_node, index_interface = @index(Global, NTuple) @@ -268,7 +268,7 @@ function finish_mpi_receive!(backend::Backend, mpi_cache::P4estMPICache, KernelAbstractions.synchronize(backend) end -@kernel function finish_mpi_receive_KAkernel!(recv_buffer, neighbor_interfaces, +@kernel inbounds=true function finish_mpi_receive_KAkernel!(recv_buffer, neighbor_interfaces, local_sides, u_mpi_interfaces, ::Val{NVARS}, ::Val{2}) where {NVARS} From 35e2c516eff7ff37f10f02db977ecda2dbe3a0d7 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 12 Oct 2026 15:28:42 +0200 Subject: [PATCH 25/46] to the docs as well --- docs/src/heterogeneous.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/heterogeneous.md b/docs/src/heterogeneous.md index 7b3acacd180..c6254cde788 100644 --- a/docs/src/heterogeneous.md +++ b/docs/src/heterogeneous.md @@ -190,7 +190,7 @@ end return nothing end - @kernel function rhs_fct_kernel!(unpacked_args, args) + @kernel inbounds=true function rhs_fct_kernel!(unpacked_args, args) element = @index(Global) rhs_fct_per_element(element, unpacked_args, args) end From ada4452931a0e62204f5c3520e45f6d0e4fafd32 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 12 Oct 2026 15:29:32 +0200 Subject: [PATCH 26/46] use Base.@propagate_inbounds --- src/solvers/dgsem_p4est/dg_2d_parallel.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index 4a60eab76ee..0d306ac79fe 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -158,7 +158,7 @@ end end # Inlined version of the interface flux computation for conservation laws -@inline function calc_mpi_interface_flux!(surface_flux_values, +Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, MeshT::Type{<:Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}}, have_nonconservative_terms::False, equations, @@ -188,7 +188,7 @@ end end # Inlined version of the interface flux computation for non-conservative equations -@inline function calc_mpi_interface_flux!(surface_flux_values, +Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, MeshT::Type{<:Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}}, have_nonconservative_terms::True, equations, From 1fbb6e7ae3c95c0d7f9234af45a5115f915398ed Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 12 Oct 2026 15:29:50 +0200 Subject: [PATCH 27/46] fix --- src/solvers/dgsem_p4est/dg_2d_parallel.jl | 2 +- test/test_kernelabstractions_mpi.jl | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index 0d306ac79fe..7686402ea8b 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -17,7 +17,7 @@ function prolong2mpiinterfaces!(backend::Nothing, cache, u, # Explicit bounds check, which allows us to assume inbounds access below @boundscheck begin check_axes(u, mesh, equations, dg, cache) - check_axes(mpi_interfaces, equations, dg, cache) + check_axes(cache.mpi_interfaces, equations, dg, cache) end @threaded for interface in eachmpiinterface(dg, cache) diff --git a/test/test_kernelabstractions_mpi.jl b/test/test_kernelabstractions_mpi.jl index da5bb1cdf83..f958c4ace92 100644 --- a/test/test_kernelabstractions_mpi.jl +++ b/test/test_kernelabstractions_mpi.jl @@ -15,8 +15,6 @@ end # Expected errors are exactly the same as with TreeMesh! l2=8.311947673061856e-6, linf=6.627000273229378e-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 real(semi.solver) == Float64 @test real(semi.solver.basis) == Float64 From 5c67317bf9172fdd3891825d95349c678ab1087f Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 12 Oct 2026 15:30:34 +0200 Subject: [PATCH 28/46] fmt --- src/solvers/dgsem_p4est/dg_2d_parallel.jl | 46 ++++++++++------- src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl | 49 +++++++++++-------- 2 files changed, 57 insertions(+), 38 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index 7686402ea8b..26113c20cd6 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -159,15 +159,20 @@ end # Inlined version of the interface flux computation for conservation laws Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, - MeshT::Type{<:Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}}, - have_nonconservative_terms::False, equations, - surface_integral, SolverT::Type{<:DG}, - u_mpi_interfaces, - interface_index, normal_direction, - interface_node_index, local_side, - surface_node_index, local_direction_index, - local_element_index) + MeshT::Type{<:Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}}, + have_nonconservative_terms::False, + equations, + surface_integral, + SolverT::Type{<:DG}, + u_mpi_interfaces, + interface_index, + normal_direction, + interface_node_index, + local_side, + surface_node_index, + local_direction_index, + local_element_index) @unpack surface_flux = surface_integral u_ll, u_rr = get_surface_node_vars(u_mpi_interfaces, equations, SolverT, @@ -189,15 +194,20 @@ end # Inlined version of the interface flux computation for non-conservative equations Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, - MeshT::Type{<:Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}}, - have_nonconservative_terms::True, equations, - surface_integral, SolverT::Type{<:DG}, - u_mpi_interfaces, - interface_index, normal_direction, - interface_node_index, local_side, - surface_node_index, local_direction_index, - local_element_index) + MeshT::Type{<:Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}}, + have_nonconservative_terms::True, + equations, + surface_integral, + SolverT::Type{<:DG}, + u_mpi_interfaces, + interface_index, + normal_direction, + interface_node_index, + local_side, + surface_node_index, + local_direction_index, + local_element_index) surface_flux, nonconservative_flux = surface_integral.surface_flux u_ll, u_rr = get_surface_node_vars(u_mpi_interfaces, equations, SolverT, diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl index b6e7115bb7f..4bf8f289e6b 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl @@ -123,9 +123,12 @@ function prolong2mpiinterfaces!(backend::Backend, cache, u, return nothing end -@kernel inbounds=true function prolong2mpiinterfaces_KAkernel!(mpi_interfaces_u, local_sides, - local_neighbor_ids, node_indices, - index_range, variables_range, u) +@kernel inbounds=true function prolong2mpiinterfaces_KAkernel!(mpi_interfaces_u, + local_sides, + local_neighbor_ids, + node_indices, + index_range, + variables_range, u) interface = @index(Global) prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, interface, local_sides, local_neighbor_ids, node_indices, index_range, @@ -154,16 +157,18 @@ function calc_mpi_interface_flux!(backend::Backend, surface_flux_values, end @kernel inbounds=true function calc_mpi_interface_flux_KAkernel!(surface_flux_values, - MeshT::Type{<:Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}}, - have_nonconservative_terms, - equations, - surface_integral, - SolverT::Type{<:DG}, - local_neighbor_ids, node_indices, - local_sides, contravariant_vectors, - u, - index_range) + MeshT::Type{<:Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}}, + have_nonconservative_terms, + equations, + surface_integral, + SolverT::Type{<:DG}, + local_neighbor_ids, + node_indices, + local_sides, + contravariant_vectors, + u, + index_range) interface = @index(Global) calc_mpi_interface_flux_per_interface!(surface_flux_values, MeshT, have_nonconservative_terms, equations, @@ -223,9 +228,11 @@ function start_mpi_send!(backend::Backend, mpi_cache::P4estMPICache, end end -@kernel inbounds=true function start_mpi_send_KAkernel!(send_buffer, neighbor_interfaces, local_sides, - u_mpi_interfaces, ::Val{NVARS}, - ::Val{2}) where {NVARS} +@kernel inbounds=true function start_mpi_send_KAkernel!(send_buffer, + neighbor_interfaces, + local_sides, + u_mpi_interfaces, ::Val{NVARS}, + ::Val{2}) where {NVARS} index_node, index_interface = @index(Global, NTuple) index_linear = @index(Global, Linear) @@ -268,10 +275,12 @@ function finish_mpi_receive!(backend::Backend, mpi_cache::P4estMPICache, KernelAbstractions.synchronize(backend) end -@kernel inbounds=true function finish_mpi_receive_KAkernel!(recv_buffer, neighbor_interfaces, - local_sides, - u_mpi_interfaces, ::Val{NVARS}, - ::Val{2}) where {NVARS} +@kernel inbounds=true function finish_mpi_receive_KAkernel!(recv_buffer, + neighbor_interfaces, + local_sides, + u_mpi_interfaces, + ::Val{NVARS}, + ::Val{2}) where {NVARS} index_node, index_interface = @index(Global, NTuple) index_linear = @index(Global, Linear) buffer_offset = (index_linear - 1) * NVARS From 62796e18abb375f727bae269bfa5d8c221c19cc6 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 13 Oct 2026 22:13:56 +0200 Subject: [PATCH 29/46] 3D --- NEWS.md | 2 +- src/callbacks_step/analysis_dg3d_parallel.jl | 64 +++- src/solvers/dgsem_p4est/dg_2d_parallel.jl | 7 +- src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl | 12 +- src/solvers/dgsem_p4est/dg_3d_parallel.jl | 300 ++++++++++-------- src/solvers/dgsem_p4est/dg_parallel.jl | 1 + test/test_kernelabstractions_mpi.jl | 33 ++ 7 files changed, 280 insertions(+), 139 deletions(-) diff --git a/NEWS.md b/NEWS.md index 77c9966ffe2..597663f29b6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,7 +8,7 @@ for human readability. ## Changes in the v0.17 lifecycle #### Added -- Added multi-GPU support for 2D P4estMesh. A local CUDA-aware MPI installation is required. +- Added multi-GPU support for 2D and 3D P4estMesh. A local CUDA-aware MPI installation is required. For configuration see [Parallelization](https://trixi-framework.org/TrixiDocumentation/stable/parallelization/#parallel_system_MPI) in the docs. - The new method `trixi_adapt(semi; storage_type, real_type, flux_differencing_kernel)` diff --git a/src/callbacks_step/analysis_dg3d_parallel.jl b/src/callbacks_step/analysis_dg3d_parallel.jl index 0870fe215eb..56805372874 100644 --- a/src/callbacks_step/analysis_dg3d_parallel.jl +++ b/src/callbacks_step/analysis_dg3d_parallel.jl @@ -13,6 +13,13 @@ function calc_error_norms(func, u, t, analyzer, @unpack node_coordinates, inverse_jacobian = cache.elements @unpack u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2, jacobian_local, jacobian_tmp1, jacobian_tmp2 = cache_analysis + # Calculate error norms on the CPU, to ensure the order of summation is the same. + if trixi_backend(u) !== nothing + node_coordinates = Array(node_coordinates) + inverse_jacobian = Array(inverse_jacobian) + u = Array(u) + end + # Set up data structures l2_error = zero(func(get_node_vars(u, equations, dg, 1, 1, 1, 1), equations)) linf_error = copy(l2_error) @@ -66,7 +73,7 @@ function calc_error_norms(func, u, t, analyzer, return l2_error, linf_error end -function integrate_via_indices(func::Func, u, +function integrate_via_indices(func::Func, backend::Nothing, u, mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, equations, dg::DGSEM, cache, args...; normalize = true) where {Func} @@ -108,4 +115,59 @@ function integrate_via_indices(func::Func, u, return integral end + +function integrate_via_indices(func::Func, backend::Backend, u, + mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, + equations, + dg::DGSEM, cache, args...; normalize = true) where {Func} + @unpack weights = dg.basis + @unpack inverse_jacobian = cache.elements + + function local_plus((integral, total_volume), (integral_element, volume_element)) + return (integral + integral_element, total_volume + volume_element) + end + + # `func(u, 1,1,1,1, equations, dg, args...)` might access GPU memory + # so we have to rely on the compiler to correctly infer the type of the integral here. + # TODO: Technically we need device_promote_op here that "infers" the function within the context of the GPU. + integral0 = zero(Base.promote_op(func, typeof(u), Int, Int, Int, Int, + typeof(equations), typeof(dg), + map(typeof, args)...)) + volume0 = zero(eltype(weights)) * zero(eltype(inverse_jacobian)) + init = neutral = (integral0, volume0) + + # Use quadrature to numerically integrate over entire domain + num_elements = nelements(dg, cache) + integral, volume = AcceleratedKernels.mapreduce(local_plus, 1:num_elements, + backend; init, + neutral) do element + # Initialize integral with zeros of the right shape + local_integral, local_total_volume = neutral + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(inverse_jacobian[i, j, k, element])) + local_integral += volume_jacobian * weights[i] * weights[j] * weights[k] * + func(u, i, j, k, element, equations, dg, args...) + local_total_volume += volume_jacobian * weights[i] * weights[j] * weights[k] + end + + return (local_integral, local_total_volume) + end + + global_integral = MPI.Reduce!(Ref(integral), +, mpi_root(), mpi_comm()) + total_volume = MPI.Reduce(volume, +, mpi_root(), mpi_comm()) + if mpi_isroot() + integral = convert(typeof(integral), global_integral[]) + else + integral = convert(typeof(integral), NaN * integral) + total_volume = volume # non-root processes receive nothing from reduce -> overwrite + end + + # Normalize with total volume + if normalize + integral = integral / total_volume + end + + return integral +end end # @muladd diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index 26113c20cd6..7232b8f7863 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -21,7 +21,7 @@ function prolong2mpiinterfaces!(backend::Nothing, cache, u, end @threaded for interface in eachmpiinterface(dg, cache) - prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, interface, + prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, typeof(mesh), interface, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u) @@ -30,7 +30,10 @@ function prolong2mpiinterfaces!(backend::Nothing, cache, u, return nothing end -@inline function prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, interface, +@inline function prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, + ::Type{<:Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}}, + interface, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl index 4bf8f289e6b..c8d37407761 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl @@ -117,22 +117,24 @@ function prolong2mpiinterfaces!(backend::Backend, cache, u, variables_range = eachvariable(equations) kernel! = prolong2mpiinterfaces_KAkernel!(backend) - kernel!(cache.mpi_interfaces.u, local_sides, local_neighbor_ids, node_indices, - index_range, variables_range, u, + kernel!(cache.mpi_interfaces.u, typeof(mesh), local_sides, local_neighbor_ids, + node_indices, index_range, variables_range, u, ndrange = nmpiinterfaces(dg, cache)) return nothing end @kernel inbounds=true function prolong2mpiinterfaces_KAkernel!(mpi_interfaces_u, + MeshT::Type{<:Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}}, local_sides, local_neighbor_ids, node_indices, index_range, variables_range, u) interface = @index(Global) - prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, interface, local_sides, - local_neighbor_ids, node_indices, index_range, - variables_range, u) + prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, MeshT, interface, + local_sides, local_neighbor_ids, node_indices, + index_range, variables_range, u) end function calc_mpi_interface_flux!(backend::Backend, surface_flux_values, diff --git a/src/solvers/dgsem_p4est/dg_3d_parallel.jl b/src/solvers/dgsem_p4est/dg_3d_parallel.jl index 1ac3d9b28c7..694a8c12392 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parallel.jl @@ -16,12 +16,13 @@ function rhs_hyperbolic!(backend::Nothing, # Prolong solution to MPI interfaces @trixi_timeit timer() "prolong2mpiinterfaces" begin - prolong2mpiinterfaces!(cache, u, mesh, equations, dg.surface_integral, dg) + prolong2mpiinterfaces!(backend, cache, u, mesh, equations, dg.surface_integral, + dg) end # Prolong solution to MPI mortars @trixi_timeit timer() "prolong2mpimortars" begin - prolong2mpimortars!(cache, u, mesh, equations, + prolong2mpimortars!(backend, cache, u, mesh, equations, dg.mortar, dg) end @@ -83,7 +84,7 @@ function rhs_hyperbolic!(backend::Nothing, # Calculate MPI interface fluxes @trixi_timeit timer() "MPI interface flux" begin - calc_mpi_interface_flux!(cache.elements.surface_flux_values, mesh, + calc_mpi_interface_flux!(backend, cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, dg.surface_integral, dg, cache) end @@ -116,70 +117,87 @@ function rhs_hyperbolic!(backend::Nothing, return nothing end -function prolong2mpiinterfaces!(cache, u, +function prolong2mpiinterfaces!(backend::Nothing, cache, u, mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, equations, surface_integral, dg::DG) - @unpack mpi_interfaces = cache + @unpack local_sides, local_neighbor_ids, node_indices = cache.mpi_interfaces + mpi_interfaces_u = cache.mpi_interfaces.u index_range = eachnode(dg) + variables_range = eachvariable(equations) # Explicit bounds check, which allows us to assume inbounds access below @boundscheck begin check_axes(u, mesh, equations, dg, cache) - check_axes(mpi_interfaces, equations, dg, cache) + check_axes(cache.mpi_interfaces, equations, dg, cache) end @threaded for interface in eachmpiinterface(dg, cache) - @inbounds begin - # Copy solution data from the local element using "delayed indexing" with - # a start value and a step size to get the correct face and orientation. - # Note that in the current implementation, the interface will be - # "aligned at the primary element", i.e., the index of the primary side - # will always run forwards. - local_side = mpi_interfaces.local_sides[interface] - local_element = mpi_interfaces.local_neighbor_ids[interface] - local_indices = mpi_interfaces.node_indices[interface] - - i_element_start, i_element_step_i, i_element_step_j = index_to_start_step_3d(local_indices[1], - index_range) - j_element_start, j_element_step_i, j_element_step_j = index_to_start_step_3d(local_indices[2], - index_range) - k_element_start, k_element_step_i, k_element_step_j = index_to_start_step_3d(local_indices[3], - index_range) - - i_element = i_element_start - j_element = j_element_start - k_element = k_element_start - for j in eachnode(dg) - for i in eachnode(dg) - for v in eachvariable(equations) - mpi_interfaces.u[local_side, v, i, j, interface] = u[v, - i_element, - j_element, - k_element, - local_element] - end - i_element += i_element_step_i - j_element += j_element_step_i - k_element += k_element_step_i + prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, typeof(mesh), interface, + local_sides, local_neighbor_ids, + node_indices, + index_range, variables_range, u) + end + + return nothing +end + +@inline function prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, + ::Type{<:Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}}, + interface, + local_sides, local_neighbor_ids, + node_indices, index_range, + variables_range, u) + @inbounds begin + # Copy solution data from the local element using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + # Note that in the current implementation, the interface will be + # "aligned at the primary element", i.e., the index of the primary side + # will always run forwards. + local_side = local_sides[interface] + local_element = local_neighbor_ids[interface] + local_indices = node_indices[interface] + + i_element_start, i_element_step_i, i_element_step_j = index_to_start_step_3d(local_indices[1], + index_range) + j_element_start, j_element_step_i, j_element_step_j = index_to_start_step_3d(local_indices[2], + index_range) + k_element_start, k_element_step_i, k_element_step_j = index_to_start_step_3d(local_indices[3], + index_range) + + i_element = i_element_start + j_element = j_element_start + k_element = k_element_start + for j in index_range + for i in index_range + for v in variables_range + mpi_interfaces_u[local_side, v, i, j, interface] = u[v, + i_element, + j_element, + k_element, + local_element] end - i_element += i_element_step_j - j_element += j_element_step_j - k_element += k_element_step_j + i_element += i_element_step_i + j_element += j_element_step_i + k_element += k_element_step_i end + i_element += i_element_step_j + j_element += j_element_step_j + k_element += k_element_step_j end end - return nothing end -function calc_mpi_interface_flux!(surface_flux_values, +function calc_mpi_interface_flux!(backend::Nothing, surface_flux_values, mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, have_nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces @unpack contravariant_vectors = cache.elements + @unpack u = cache.mpi_interfaces index_range = eachnode(dg) # Explicit bounds check, which allows us to assume inbounds access below @@ -190,85 +208,107 @@ function calc_mpi_interface_flux!(surface_flux_values, end @threaded for interface in eachmpiinterface(dg, cache) - @inbounds begin - # Get element and side index information on the local element - local_element = local_neighbor_ids[interface] - local_indices = node_indices[interface] - local_direction = indices2direction(local_indices) - local_side = local_sides[interface] - - # Create the local i,j,k indexing on the local element used to pull normal direction information - i_element_start, i_element_step_i, i_element_step_j = index_to_start_step_3d(local_indices[1], - index_range) - j_element_start, j_element_step_i, j_element_step_j = index_to_start_step_3d(local_indices[2], - index_range) - k_element_start, k_element_step_i, k_element_step_j = index_to_start_step_3d(local_indices[3], - index_range) - - i_element = i_element_start - j_element = j_element_start - k_element = k_element_start - - # Initiate the node indices to be used in the surface for loop, - # the surface flux storage must be indexed in alignment with the local element indexing - local_surface_indices = surface_indices(local_indices) - i_surface_start, i_surface_step_i, i_surface_step_j = index_to_start_step_3d(local_surface_indices[1], - index_range) - j_surface_start, j_surface_step_i, j_surface_step_j = index_to_start_step_3d(local_surface_indices[2], - index_range) - i_surface = i_surface_start - j_surface = j_surface_start + calc_mpi_interface_flux_per_interface!(surface_flux_values, typeof(mesh), + have_nonconservative_terms, equations, + surface_integral, typeof(dg), + local_neighbor_ids, node_indices, + local_sides, contravariant_vectors, u, + index_range, interface) + end + + return nothing +end + +@inline function calc_mpi_interface_flux_per_interface!(surface_flux_values, + MeshT::Type{<:Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}}, + have_nonconservative_terms, + equations, + surface_integral, + SolverT::Type{<:DG}, + local_neighbor_ids, + node_indices, local_sides, + contravariant_vectors, + u_mpi_interfaces, + index_range, interface) + @inbounds begin + # Get element and side index information on the local element + local_element = local_neighbor_ids[interface] + local_indices = node_indices[interface] + local_direction = indices2direction(local_indices) + local_side = local_sides[interface] + + # Create the local i,j,k indexing on the local element used to pull normal direction information + i_element_start, i_element_step_i, i_element_step_j = index_to_start_step_3d(local_indices[1], + index_range) + j_element_start, j_element_step_i, j_element_step_j = index_to_start_step_3d(local_indices[2], + index_range) + k_element_start, k_element_step_i, k_element_step_j = index_to_start_step_3d(local_indices[3], + index_range) + + i_element = i_element_start + j_element = j_element_start + k_element = k_element_start + + # Initiate the node indices to be used in the surface for loop, + # the surface flux storage must be indexed in alignment with the local element indexing + local_surface_indices = surface_indices(local_indices) + i_surface_start, i_surface_step_i, i_surface_step_j = index_to_start_step_3d(local_surface_indices[1], + index_range) + j_surface_start, j_surface_step_i, j_surface_step_j = index_to_start_step_3d(local_surface_indices[2], + index_range) + i_surface = i_surface_start + j_surface = j_surface_start + + for j in index_range + for i in index_range + # Get the normal direction on the local element + # Contravariant vectors at interfaces in negative coordinate direction + # are pointing inwards. This is handled by `get_normal_direction`. + normal_direction = get_normal_direction(local_direction, + contravariant_vectors, + i_element, + j_element, + k_element, + local_element) + + calc_mpi_interface_flux!(surface_flux_values, MeshT, + have_nonconservative_terms, equations, + surface_integral, SolverT, + u_mpi_interfaces, interface, normal_direction, + i, j, local_side, + i_surface, j_surface, local_direction, + local_element) - for j in eachnode(dg) - for i in eachnode(dg) - # Get the normal direction on the local element - # Contravariant vectors at interfaces in negative coordinate direction - # are pointing inwards. This is handled by `get_normal_direction`. - normal_direction = get_normal_direction(local_direction, - contravariant_vectors, - i_element, - j_element, - k_element, - local_element) - - calc_mpi_interface_flux!(surface_flux_values, mesh, - have_nonconservative_terms, equations, - surface_integral, dg, cache, - interface, normal_direction, - i, j, local_side, - i_surface, j_surface, local_direction, - local_element) - - # Increment local element indices to pull the normal direction - i_element += i_element_step_i - j_element += j_element_step_i - k_element += k_element_step_i - # Increment the surface node indices along the local element - i_surface += i_surface_step_i - j_surface += j_surface_step_i - end # Increment local element indices to pull the normal direction - i_element += i_element_step_j - j_element += j_element_step_j - k_element += k_element_step_j + i_element += i_element_step_i + j_element += j_element_step_i + k_element += k_element_step_i # Increment the surface node indices along the local element - i_surface += i_surface_step_j - j_surface += j_surface_step_j + i_surface += i_surface_step_i + j_surface += j_surface_step_i end + # Increment local element indices to pull the normal direction + i_element += i_element_step_j + j_element += j_element_step_j + k_element += k_element_step_j + # Increment the surface node indices along the local element + i_surface += i_surface_step_j + j_surface += j_surface_step_j end end - return nothing end # Inlined version of the interface flux computation for conservation laws Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, - mesh::Union{P4estMeshParallel{3}, - T8codeMeshParallel{3}}, + MeshT::Type{<:Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}}, have_nonconservative_terms::False, equations, - surface_integral, dg::DG, - cache, + surface_integral, + SolverT::Type{<:DG}, + u_mpi_interfaces, interface_index, normal_direction, interface_i_node_index, @@ -278,10 +318,9 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, surface_j_node_index, local_direction_index, local_element_index) - @unpack u = cache.mpi_interfaces @unpack surface_flux = surface_integral - u_ll, u_rr = get_surface_node_vars(u, equations, dg, + u_ll, u_rr = get_surface_node_vars(u_mpi_interfaces, equations, SolverT, interface_i_node_index, interface_j_node_index, interface_index) @@ -301,12 +340,13 @@ end # Inlined version of the interface flux computation for non-conservative equations Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, - mesh::Union{P4estMeshParallel{3}, - T8codeMeshParallel{3}}, + MeshT::Type{<:Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}}, have_nonconservative_terms::True, equations, - surface_integral, dg::DG, - cache, + surface_integral, + SolverT::Type{<:DG}, + u_mpi_interfaces, interface_index, normal_direction, interface_i_node_index, @@ -317,12 +357,12 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, local_direction_index, local_element_index) calc_mpi_interface_flux!(surface_flux_values, - mesh, + MeshT, have_nonconservative_terms, combine_conservative_and_nonconservative_fluxes(surface_integral.surface_flux, equations), equations, - surface_integral, dg, cache, + surface_integral, SolverT, u_mpi_interfaces, interface_index, normal_direction, interface_i_node_index, interface_j_node_index, local_side, @@ -332,13 +372,14 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, end Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, - mesh::Union{P4estMeshParallel{3}, - T8codeMeshParallel{3}}, + MeshT::Type{<:Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}}, have_nonconservative_terms::True, combine_conservative_and_nonconservative_fluxes::False, equations, - surface_integral, dg::DG, - cache, + surface_integral, + SolverT::Type{<:DG}, + u_mpi_interfaces, interface_index, normal_direction, interface_i_node_index, @@ -348,10 +389,9 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, surface_j_node_index, local_direction_index, local_element_index) - @unpack u = cache.mpi_interfaces surface_flux, nonconservative_flux = surface_integral.surface_flux - u_ll, u_rr = get_surface_node_vars(u, equations, dg, + u_ll, u_rr = get_surface_node_vars(u_mpi_interfaces, equations, SolverT, interface_i_node_index, interface_j_node_index, interface_index) @@ -374,13 +414,14 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, end Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, - mesh::Union{P4estMeshParallel{3}, - T8codeMeshParallel{3}}, + MeshT::Type{<:Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}}, have_nonconservative_terms::True, combine_conservative_and_nonconservative_fluxes::True, equations, - surface_integral, dg::DG, - cache, + surface_integral, + SolverT::Type{<:DG}, + u_mpi_interfaces, interface_index, normal_direction, interface_i_node_index, @@ -390,10 +431,9 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, surface_j_node_index, local_direction_index, local_element_index) - @unpack u = cache.mpi_interfaces surface_flux = surface_integral.surface_flux - u_ll, u_rr = get_surface_node_vars(u, equations, dg, + u_ll, u_rr = get_surface_node_vars(u_mpi_interfaces, equations, SolverT, interface_i_node_index, interface_j_node_index, interface_index) @@ -414,7 +454,7 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, return nothing end -function prolong2mpimortars!(cache, u, +function prolong2mpimortars!(backend::Nothing, cache, u, mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, equations, mortar_l2::LobattoLegendreMortarL2, diff --git a/src/solvers/dgsem_p4est/dg_parallel.jl b/src/solvers/dgsem_p4est/dg_parallel.jl index d84e33810ca..76d662bedc3 100644 --- a/src/solvers/dgsem_p4est/dg_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_parallel.jl @@ -637,6 +637,7 @@ end include("dg_2d_parallel.jl") include("dg_2d_parallel_gpu.jl") include("dg_3d_parallel.jl") +include("dg_3d_parallel_gpu.jl") include("dg_2d_parabolic_parallel.jl") include("dg_3d_parabolic_parallel.jl") end # muladd diff --git a/test/test_kernelabstractions_mpi.jl b/test/test_kernelabstractions_mpi.jl index f958c4ace92..da147b897f8 100644 --- a/test/test_kernelabstractions_mpi.jl +++ b/test/test_kernelabstractions_mpi.jl @@ -38,3 +38,36 @@ end Trixi.rhs_hyperbolic!(du_ode, u_ode, ode.p, first(ode.tspan)) @test all(isfinite, du_ode) end + +@testitem "KernelAbstractions MPI GPU 3D: elixir_advection_basic.jl" setup=[ + Setup, + KernelAbstractionsMPIExamples +] tags=[:kernelabstractions_mpi] begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem", + "elixir_advection_basic.jl"), + # Expected errors are exactly the same as with TreeMesh! + l2=[0.00016263963870641478], + linf=[0.0014537194925779984]) + semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem. + @test real(semi.solver) == Float64 + @test real(semi.solver.basis) == Float64 + @test real(semi.solver.mortar) == Float64 + # TODO: `mesh` is currently not `adapt`ed correctly + @test real(semi.mesh) == Float64 + @test typeof(semi.equations.advection_velocity) == SVector{3, Float64} + + @test ode.u0 isa Array + @test semi.solver.basis.derivative_matrix isa Array + + @test Trixi.storage_type(semi.cache.elements) === Array + @test Trixi.storage_type(semi.cache.interfaces) === Array + @test Trixi.storage_type(semi.cache.boundaries) === Array + @test Trixi.storage_type(semi.cache.mortars) === Array + + # Ensure that the RHS computation overwrites existing data in `du` correctly. + u_ode = copy(ode.u0) + du_ode = similar(u_ode) + fill!(du_ode, convert(eltype(du_ode), NaN)) + Trixi.rhs_hyperbolic!(du_ode, u_ode, ode.p, first(ode.tspan)) + @test all(isfinite, du_ode) +end From 3d3fb3140008e69e30bb50fde2baad17c024bf88 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Tue, 13 Oct 2026 22:21:32 +0200 Subject: [PATCH 30/46] missed file --- src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl | 297 ++++++++++++++++++ 1 file changed, 297 insertions(+) create mode 100644 src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl diff --git a/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl new file mode 100644 index 00000000000..cf4415b2cc8 --- /dev/null +++ b/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl @@ -0,0 +1,297 @@ +# 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 rhs_hyperbolic!(backend::Backend, + du, u, t, + mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, + equations, + boundary_conditions, source_terms::Source, + dg::DG, cache) where {Source} + # Start to receive MPI data + @trixi_timeit timer() "start MPI receive" start_mpi_receive!(cache.mpi_cache) + + # Prolong solution to MPI interfaces + @trixi_timeit_ext backend timer() "prolong2mpiinterfaces" begin + prolong2mpiinterfaces!(backend, cache, u, mesh, equations, dg.surface_integral, + dg) + end + + # Prolong solution to MPI mortars + @trixi_timeit_ext backend timer() "prolong2mpimortars" begin + prolong2mpimortars!(backend, cache, u, mesh, equations, + dg.mortar, dg) + end + + # Start to send MPI data + @trixi_timeit_ext backend timer() "start MPI send" begin + start_mpi_send!(backend, cache.mpi_cache, mesh, equations, dg, cache) + end + + # Calculate volume integral + @trixi_timeit_ext backend timer() "volume integral" begin + calc_volume_integral!(backend, du, u, mesh, + have_nonconservative_terms(equations), equations, + dg.volume_integral, dg, cache) + end + + # Calculate interface fluxes + @trixi_timeit_ext backend timer() "prolong2interfaces + flux" begin + prolong2interfaces_and_calc_interface_flux!(backend, + cache.elements.surface_flux_values, + u, mesh, + have_nonconservative_terms(equations), + equations, + dg.surface_integral, dg, cache) + end + + # Prolong solution to boundaries + @trixi_timeit_ext backend timer() "prolong2boundaries" begin + prolong2boundaries!(backend, cache, u, mesh, equations, dg) + end + + # Calculate boundary fluxes + @trixi_timeit_ext backend timer() "boundary flux" begin + calc_boundary_flux!(backend, cache, t, boundary_conditions, mesh, equations, + dg.surface_integral, dg) + end + + # Prolong solution to mortars + @trixi_timeit_ext backend timer() "prolong2mortars" begin + prolong2mortars!(backend, cache, u, mesh, equations, + dg.mortar, dg) + end + + # Calculate mortar fluxes + @trixi_timeit_ext backend timer() "mortar flux" begin + calc_mortar_flux!(backend, cache.elements.surface_flux_values, mesh, + have_nonconservative_terms(equations), equations, + dg.mortar, dg.surface_integral, dg, cache) + end + + # Finish to receive MPI data + @trixi_timeit_ext backend timer() "finish MPI receive" begin + finish_mpi_receive!(backend, cache.mpi_cache, mesh, equations, dg, cache) + end + + # Calculate MPI interface fluxes + @trixi_timeit_ext backend timer() "MPI interface flux" begin + calc_mpi_interface_flux!(backend, cache.elements.surface_flux_values, mesh, + have_nonconservative_terms(equations), equations, + dg.surface_integral, dg, cache) + end + + # Calculate MPI mortar fluxes + @trixi_timeit_ext backend timer() "MPI mortar flux" begin + calc_mpi_mortar_flux!(backend, cache.elements.surface_flux_values, mesh, + have_nonconservative_terms(equations), equations, + dg.mortar, dg.surface_integral, dg, cache) + end + + # Calculate surface integrals, apply Jacobian from mapping to reference element + # and calculate source terms + @trixi_timeit_ext backend timer() "surface, Jacobian + source terms" begin + calc_surface_integral_and_apply_jacobian_and_calc_sources!(backend, du, u, t, + source_terms, mesh, + equations, + dg.surface_integral, + dg, cache) + end + + # Finish to send MPI data + @trixi_timeit timer() "finish MPI send" finish_mpi_send!(cache.mpi_cache) + + return nothing +end + +function prolong2mpiinterfaces!(backend::Backend, cache, u, + mesh::Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}, + equations, surface_integral, dg::DG) + nmpiinterfaces(dg, cache) == 0 && return nothing + @unpack local_sides, local_neighbor_ids, node_indices = cache.mpi_interfaces + index_range = eachnode(dg) + variables_range = eachvariable(equations) + + kernel! = prolong2mpiinterfaces_KAkernel!(backend) + kernel!(cache.mpi_interfaces.u, typeof(mesh), local_sides, local_neighbor_ids, + node_indices, index_range, variables_range, u, + ndrange = nmpiinterfaces(dg, cache)) + return nothing +end + +@kernel inbounds=true function prolong2mpiinterfaces_KAkernel!(mpi_interfaces_u, + MeshT::Type{<:Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}}, + local_sides, + local_neighbor_ids, + node_indices, + index_range, + variables_range, u) + interface = @index(Global) + prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, MeshT, interface, + local_sides, local_neighbor_ids, node_indices, + index_range, variables_range, u) +end + +function calc_mpi_interface_flux!(backend::Backend, surface_flux_values, + mesh::Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}, + have_nonconservative_terms, + equations, surface_integral, dg::DG, cache) + nmpiinterfaces(dg, cache) == 0 && return nothing + @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces + @unpack contravariant_vectors = cache.elements + @unpack u = cache.mpi_interfaces + index_range = eachnode(dg) + + kernel! = calc_mpi_interface_flux_KAkernel!(backend) + kernel!(surface_flux_values, typeof(mesh), + have_nonconservative_terms, equations, + surface_integral, typeof(dg), + local_neighbor_ids, node_indices, local_sides, contravariant_vectors, + u, index_range, + ndrange = nmpiinterfaces(dg, cache)) + return nothing +end + +@kernel inbounds=true function calc_mpi_interface_flux_KAkernel!(surface_flux_values, + MeshT::Type{<:Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}}, + have_nonconservative_terms, + equations, + surface_integral, + SolverT::Type{<:DG}, + local_neighbor_ids, + node_indices, + local_sides, + contravariant_vectors, + u, + index_range) + interface = @index(Global) + calc_mpi_interface_flux_per_interface!(surface_flux_values, MeshT, + have_nonconservative_terms, equations, + surface_integral, SolverT, + local_neighbor_ids, node_indices, + local_sides, contravariant_vectors, u, + index_range, interface) +end + +# For GPU backends MPI mortars are not yet implemented +function prolong2mpimortars!(backend::Backend, cache, u, + mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, + equations, mortar, dg::DG) + @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" + return nothing +end + +# For GPU backends MPI mortars are not yet implemented +function calc_mpi_mortar_flux!(backend::Backend, surface_flux_values, + mesh::Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}, + have_nonconservative_terms, equations, mortar, + surface_integral, dg::DG, cache) + @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" + return nothing +end + +function start_mpi_send!(backend::Backend, mpi_cache::P4estMPICache, + mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, + equations, dg, cache) + # TODO GPU: MPI mortars are not yet implemented + @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" + + @unpack mpi_neighbor_ranks, mpi_neighbor_interfaces = mpi_cache + @unpack mpi_send_buffers, mpi_send_requests = mpi_cache + @unpack local_sides, u = cache.mpi_interfaces + + kernel! = start_mpi_send_KAkernel!(backend) + + for (rank_index, _) in enumerate(mpi_neighbor_ranks) + send_buffer = mpi_send_buffers[rank_index] + neighbor_interfaces = mpi_neighbor_interfaces[rank_index] + if !isempty(neighbor_interfaces) + kernel!(send_buffer, neighbor_interfaces, local_sides, u, + Val(nvariables(equations)), Val(ndims(mesh)), + ndrange = (nnodes(dg), nnodes(dg), length(neighbor_interfaces))) + end + end + + # wait for the kernels to return before sending the buffers + KernelAbstractions.synchronize(backend) + + for (rank_index, neighbor_rank) in enumerate(mpi_neighbor_ranks) + send_buffer = mpi_send_buffers[rank_index] + mpi_send_requests[rank_index] = MPI.Isend(send_buffer, neighbor_rank, + mpi_rank(), mpi_comm()) + end +end + +@kernel inbounds=true function start_mpi_send_KAkernel!(send_buffer, + neighbor_interfaces, + local_sides, + u_mpi_interfaces, ::Val{NVARS}, + ::Val{3}) where {NVARS} + i, j, index_interface = @index(Global, NTuple) + index_linear = @index(Global, Linear) + + # The data of each interface is stored contiguously in the buffer, with the + # variables running fastest, followed by the surface node indices `i` and `j`. + buffer_offset = (index_linear - 1) * NVARS + interface = neighbor_interfaces[index_interface] + local_side = local_sides[interface] + + for v in 1:NVARS + send_buffer[buffer_offset + v] = u_mpi_interfaces[local_side, v, i, j, + interface] + end +end + +function finish_mpi_receive!(backend::Backend, mpi_cache::P4estMPICache, + mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, + equations, dg, cache) + # TODO GPU: MPI mortars are not yet implemented + @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" + + @unpack mpi_neighbor_interfaces = mpi_cache + @unpack mpi_recv_buffers, mpi_recv_requests = mpi_cache + @unpack local_sides, u = cache.mpi_interfaces + + kernel! = finish_mpi_receive_KAkernel!(backend) + + # Start receiving and unpack received data until all communication is finished + data = MPI.Waitany(mpi_recv_requests) + while data !== nothing + recv_buffer = mpi_recv_buffers[data] + neighbor_interfaces = mpi_neighbor_interfaces[data] + if !isempty(neighbor_interfaces) + kernel!(recv_buffer, neighbor_interfaces, local_sides, u, + Val(nvariables(equations)), Val(ndims(mesh)), + ndrange = (nnodes(dg), nnodes(dg), length(neighbor_interfaces))) + end + + data = MPI.Waitany(mpi_recv_requests) + end + # Wait for the last kernel to return ? + KernelAbstractions.synchronize(backend) +end + +@kernel inbounds=true function finish_mpi_receive_KAkernel!(recv_buffer, + neighbor_interfaces, + local_sides, + u_mpi_interfaces, + ::Val{NVARS}, + ::Val{3}) where {NVARS} + i, j, index_interface = @index(Global, NTuple) + index_linear = @index(Global, Linear) + buffer_offset = (index_linear - 1) * NVARS + interface = neighbor_interfaces[index_interface] + remote_side = local_sides[interface] == 1 ? 2 : 1 + for v in 1:NVARS + u_mpi_interfaces[remote_side, v, i, j, interface] = recv_buffer[buffer_offset + v] + end +end +end From 3979a79e0f365d0cfcb926f81ae3b35347edb8af Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Sat, 17 Oct 2026 23:28:52 +0200 Subject: [PATCH 31/46] fix merge [skip ci] --- src/callbacks_step/analysis_dg2d_parallel.jl | 54 ------- src/solvers/dgsem_p4est/dg_2d_parallel.jl | 162 +++++-------------- 2 files changed, 41 insertions(+), 175 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d_parallel.jl b/src/callbacks_step/analysis_dg2d_parallel.jl index d20b0d7776b..7119016d75b 100644 --- a/src/callbacks_step/analysis_dg2d_parallel.jl +++ b/src/callbacks_step/analysis_dg2d_parallel.jl @@ -287,58 +287,4 @@ function integrate_via_indices(func::Func, backend::Backend, u, return reduce_integral_and_volume_mpi(integral, volume, normalize) end - -function integrate_via_indices(func::Func, backend::Backend, u, - mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, - equations, - dg::DGSEM, cache, args...; normalize = true) where {Func} - @unpack weights = dg.basis - @unpack inverse_jacobian = cache.elements - - function local_plus((integral, total_volume), (integral_element, volume_element)) - return (integral + integral_element, total_volume + volume_element) - end - - # `func(u, 1,1,1, equations, dg, args...)` might access GPU memory - # so we have to rely on the compiler to correctly infer the type of the integral here. - # TODO: Technically we need device_promote_op here that "infers" the function within the context of the GPU. - integral0 = zero(Base.promote_op(func, typeof(u), Int, Int, Int, typeof(equations), - typeof(dg), map(typeof, args)...)) - volume0 = zero(eltype(weights)) * zero(eltype(inverse_jacobian)) - init = neutral = (integral0, volume0) - - # Use quadrature to numerically integrate over entire domain - num_elements = nelements(dg, cache) - integral, volume = AcceleratedKernels.mapreduce(local_plus, 1:num_elements, - backend; init, - neutral) do element - # Initialize integral with zeros of the right shapeu - local_integral, local_total_volume = neutral - - for j in eachnode(dg), i in eachnode(dg) - volume_jacobian = abs(inv(inverse_jacobian[i, j, element])) - local_integral += volume_jacobian * weights[i] * weights[j] * - func(u, i, j, element, equations, dg, args...) - local_total_volume += volume_jacobian * weights[i] * weights[j] - end - - return (local_integral, local_total_volume) - end - - global_integral = MPI.Reduce!(Ref(integral), +, mpi_root(), mpi_comm()) - total_volume = MPI.Reduce(volume, +, mpi_root(), mpi_comm()) - if mpi_isroot() - integral = convert(typeof(integral), global_integral[]) - else - integral = convert(typeof(integral), NaN * integral) - total_volume = volume # non-root processes receive nothing from reduce -> overwrite - end - - # Normalize with total volume - if normalize - integral = integral / total_volume - end - - return integral -end end # @muladd diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index 062438698cb..87f48b7decc 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -9,61 +9,46 @@ function prolong2mpiinterfaces!(backend::Nothing, cache, u, mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, equations, surface_integral, dg::DG) - @unpack local_sides, local_neighbor_ids, node_indices = cache.mpi_interfaces - mpi_interfaces_u = cache.mpi_interfaces.u + @unpack mpi_interfaces = cache index_range = eachnode(dg) - variables_range = eachvariable(equations) # Explicit bounds check, which allows us to assume inbounds access below @boundscheck begin check_axes(u, mesh, equations, dg, cache) - check_axes(cache.mpi_interfaces, equations, dg, cache) + check_axes(mpi_interfaces, equations, dg, cache) end @threaded for interface in eachmpiinterface(dg, cache) - prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, typeof(mesh), interface, - local_sides, local_neighbor_ids, - node_indices, - index_range, variables_range, u) - end - - return nothing -end - -@inline function prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, - ::Type{<:Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}}, - interface, - local_sides, local_neighbor_ids, - node_indices, index_range, - variables_range, u) - @inbounds begin - # Copy solution data from the local element using "delayed indexing" with - # a start value and a step size to get the correct face and orientation. - # Note that in the current implementation, the interface will be - # "aligned at the primary element", i.e., the index of the primary side - # will always run forwards. - local_side = local_sides[interface] - local_element = local_neighbor_ids[interface] - local_indices = node_indices[interface] - - i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], - index_range) - j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], - index_range) - - i_element = i_element_start - j_element = j_element_start - for i in index_range - for v in variables_range - mpi_interfaces_u[local_side, v, i, interface] = u[v, i_element, - j_element, - local_element] + @inbounds begin + # Copy solution data from the local element using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + # Note that in the current implementation, the interface will be + # "aligned at the primary element", i.e., the index of the primary side + # will always run forwards. + local_side = mpi_interfaces.local_sides[interface] + local_element = mpi_interfaces.local_neighbor_ids[interface] + local_indices = mpi_interfaces.node_indices[interface] + + i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], + index_range) + j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], + index_range) + + i_element = i_element_start + j_element = j_element_start + for i in eachnode(dg) + for v in eachvariable(equations) + mpi_interfaces.u[local_side, v, i, interface] = u[v, + i_element, + j_element, + local_element] + end + i_element += i_element_step + j_element += j_element_step end - i_element += i_element_step - j_element += j_element_step end end + return nothing end @@ -74,6 +59,7 @@ function calc_mpi_interface_flux!(backend::Nothing, surface_flux_values, equations, surface_integral, dg::DG, cache) @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces index_range = eachnode(dg) + index_end = last(index_range) # Explicit bounds check, which allows us to assume inbounds access below @boundscheck begin @@ -123,79 +109,14 @@ function calc_mpi_interface_flux!(backend::Nothing, surface_flux_values, return nothing end -@inline function calc_mpi_interface_flux_per_interface!(surface_flux_values, - MeshT::Type{<:Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}}, - have_nonconservative_terms, - equations, - surface_integral, - SolverT::Type{<:DG}, - local_neighbor_ids, - node_indices, local_sides, - contravariant_vectors, - u_mpi_interfaces, - index_range, interface) - @inbounds begin - # Get element and side index information on the local element - local_element = local_neighbor_ids[interface] - local_indices = node_indices[interface] - local_direction = indices2direction(local_indices) - local_side = local_sides[interface] - - # Create the local i,j indexing on the local element used to pull normal direction information - i_element_start, i_element_step = index_to_start_step_2d(local_indices[1], - index_range) - j_element_start, j_element_step = index_to_start_step_2d(local_indices[2], - index_range) - - i_element = i_element_start - j_element = j_element_start - - # Initiate the node index to be used in the surface for loop, - # the surface flux storage must be indexed in alignment with the local element indexing - if :i_backward in local_indices - surface_node = last(index_range) - surface_node_step = -1 - else - surface_node = 1 - surface_node_step = 1 - end - - for node in index_range - # Get the normal direction on the local element - # Contravariant vectors at interfaces in negative coordinate direction - # are pointing inwards. This is handled by `get_normal_direction`. - normal_direction = get_normal_direction(local_direction, - contravariant_vectors, - i_element, j_element, local_element) - - calc_mpi_interface_flux!(surface_flux_values, MeshT, - have_nonconservative_terms, equations, - surface_integral, SolverT, - u_mpi_interfaces, interface, normal_direction, - node, local_side, surface_node, local_direction, - local_element) - - # Increment local element indices to pull the normal direction - i_element += i_element_step - j_element += j_element_step - - # Increment the surface node index along the local element - surface_node += surface_node_step - end - end - return nothing -end - # Inlined version of the interface flux computation for conservation laws Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, - MeshT::Type{<:Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}}, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, have_nonconservative_terms::False, equations, - surface_integral, - SolverT::Type{<:DG}, - u_mpi_interfaces, + surface_integral, dg::DG, + cache, interface_index, normal_direction, interface_node_index, @@ -203,10 +124,10 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, surface_node_index, local_direction_index, local_element_index) + @unpack u = cache.mpi_interfaces @unpack surface_flux = surface_integral - u_ll, u_rr = get_surface_node_vars(u_mpi_interfaces, equations, SolverT, - interface_node_index, + u_ll, u_rr = get_surface_node_vars(u, equations, dg, interface_node_index, interface_index) # `normal_direction` is the normal direction of the primary element. @@ -227,13 +148,12 @@ end # Inlined version of the interface flux computation for non-conservative equations Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, - MeshT::Type{<:Union{P4estMeshParallel{2}, - T8codeMeshParallel{2}}}, + mesh::Union{P4estMeshParallel{2}, + T8codeMeshParallel{2}}, have_nonconservative_terms::True, equations, - surface_integral, - SolverT::Type{<:DG}, - u_mpi_interfaces, + surface_integral, dg::DG, + cache, interface_index, normal_direction, interface_node_index, @@ -273,7 +193,7 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, @unpack u = cache.mpi_interfaces surface_flux, nonconservative_flux = surface_integral.surface_flux - u_ll, u_rr = get_surface_node_vars(u_mpi_interfaces, equations, SolverT, + u_ll, u_rr = get_surface_node_vars(u, equations, dg, interface_node_index, interface_index) From 81aa553706ea9c5f5f68d1458e37d5a534c422df Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Sun, 18 Oct 2026 00:08:53 +0200 Subject: [PATCH 32/46] revert to main --- src/solvers/dgsem_p4est/dg_3d_parallel.jl | 178 ++++++++-------------- 1 file changed, 64 insertions(+), 114 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_3d_parallel.jl b/src/solvers/dgsem_p4est/dg_3d_parallel.jl index 1995fe8bc52..33e3fde3cf3 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parallel.jl @@ -121,22 +121,54 @@ function prolong2mpiinterfaces!(backend::Nothing, cache, u, mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, equations, surface_integral, dg::DG) - @unpack local_sides, local_neighbor_ids, node_indices = cache.mpi_interfaces - mpi_interfaces_u = cache.mpi_interfaces.u + @unpack mpi_interfaces = cache index_range = eachnode(dg) - variables_range = eachvariable(equations) # Explicit bounds check, which allows us to assume inbounds access below @boundscheck begin check_axes(u, mesh, equations, dg, cache) - check_axes(cache.mpi_interfaces, equations, dg, cache) + check_axes(mpi_interfaces, equations, dg, cache) end @threaded for interface in eachmpiinterface(dg, cache) - prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, typeof(mesh), interface, - local_sides, local_neighbor_ids, - node_indices, - index_range, variables_range, u) + @inbounds begin + # Copy solution data from the local element using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + # Note that in the current implementation, the interface will be + # "aligned at the primary element", i.e., the index of the primary side + # will always run forwards. + local_side = mpi_interfaces.local_sides[interface] + local_element = mpi_interfaces.local_neighbor_ids[interface] + local_indices = mpi_interfaces.node_indices[interface] + + i_element_start, i_element_step_i, i_element_step_j = index_to_start_step_3d(local_indices[1], + index_range) + j_element_start, j_element_step_i, j_element_step_j = index_to_start_step_3d(local_indices[2], + index_range) + k_element_start, k_element_step_i, k_element_step_j = index_to_start_step_3d(local_indices[3], + index_range) + + i_element = i_element_start + j_element = j_element_start + k_element = k_element_start + for j in eachnode(dg) + for i in eachnode(dg) + for v in eachvariable(equations) + mpi_interfaces.u[local_side, v, i, j, interface] = u[v, + i_element, + j_element, + k_element, + local_element] + end + i_element += i_element_step_i + j_element += j_element_step_i + k_element += k_element_step_i + end + i_element += i_element_step_j + j_element += j_element_step_j + k_element += k_element_step_j + end + end end return nothing @@ -204,96 +236,14 @@ function calc_mpi_interface_flux!(backend::Nothing, surface_flux_values, return nothing end -@inline function calc_mpi_interface_flux_per_interface!(surface_flux_values, - MeshT::Type{<:Union{P4estMeshParallel{3}, - T8codeMeshParallel{3}}}, - have_nonconservative_terms, - equations, - surface_integral, - SolverT::Type{<:DG}, - local_neighbor_ids, - node_indices, local_sides, - contravariant_vectors, - u_mpi_interfaces, - index_range, interface) - @inbounds begin - # Get element and side index information on the local element - local_element = local_neighbor_ids[interface] - local_indices = node_indices[interface] - local_direction = indices2direction(local_indices) - local_side = local_sides[interface] - - # Create the local i,j,k indexing on the local element used to pull normal direction information - i_element_start, i_element_step_i, i_element_step_j = index_to_start_step_3d(local_indices[1], - index_range) - j_element_start, j_element_step_i, j_element_step_j = index_to_start_step_3d(local_indices[2], - index_range) - k_element_start, k_element_step_i, k_element_step_j = index_to_start_step_3d(local_indices[3], - index_range) - - i_element = i_element_start - j_element = j_element_start - k_element = k_element_start - - # Initiate the node indices to be used in the surface for loop, - # the surface flux storage must be indexed in alignment with the local element indexing - local_surface_indices = surface_indices(local_indices) - i_surface_start, i_surface_step_i, i_surface_step_j = index_to_start_step_3d(local_surface_indices[1], - index_range) - j_surface_start, j_surface_step_i, j_surface_step_j = index_to_start_step_3d(local_surface_indices[2], - index_range) - i_surface = i_surface_start - j_surface = j_surface_start - - for j in index_range - for i in index_range - # Get the normal direction on the local element - # Contravariant vectors at interfaces in negative coordinate direction - # are pointing inwards. This is handled by `get_normal_direction`. - normal_direction = get_normal_direction(local_direction, - contravariant_vectors, - i_element, - j_element, - k_element, - local_element) - - calc_mpi_interface_flux!(surface_flux_values, MeshT, - have_nonconservative_terms, equations, - surface_integral, SolverT, - u_mpi_interfaces, interface, normal_direction, - i, j, local_side, - i_surface, j_surface, local_direction, - local_element) - - # Increment local element indices to pull the normal direction - i_element += i_element_step_i - j_element += j_element_step_i - k_element += k_element_step_i - # Increment the surface node indices along the local element - i_surface += i_surface_step_i - j_surface += j_surface_step_i - end - # Increment local element indices to pull the normal direction - i_element += i_element_step_j - j_element += j_element_step_j - k_element += k_element_step_j - # Increment the surface node indices along the local element - i_surface += i_surface_step_j - j_surface += j_surface_step_j - end - end - return nothing -end - # Inlined version of the interface flux computation for conservation laws Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, - MeshT::Type{<:Union{P4estMeshParallel{3}, - T8codeMeshParallel{3}}}, + mesh::Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}, have_nonconservative_terms::False, equations, - surface_integral, - SolverT::Type{<:DG}, - u_mpi_interfaces, + surface_integral, dg::DG, + cache, interface_index, normal_direction, interface_i_node_index, @@ -303,9 +253,10 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, surface_j_node_index, local_direction_index, local_element_index) + @unpack u = cache.mpi_interfaces @unpack surface_flux = surface_integral - u_ll, u_rr = get_surface_node_vars(u_mpi_interfaces, equations, SolverT, + u_ll, u_rr = get_surface_node_vars(u, equations, dg, interface_i_node_index, interface_j_node_index, interface_index) @@ -328,13 +279,12 @@ end # Inlined version of the interface flux computation for non-conservative equations Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, - MeshT::Type{<:Union{P4estMeshParallel{3}, - T8codeMeshParallel{3}}}, + mesh::Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}, have_nonconservative_terms::True, equations, - surface_integral, - SolverT::Type{<:DG}, - u_mpi_interfaces, + surface_integral, dg::DG, + cache, interface_index, normal_direction, interface_i_node_index, @@ -345,12 +295,12 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, local_direction_index, local_element_index) calc_mpi_interface_flux!(surface_flux_values, - MeshT, + mesh, have_nonconservative_terms, combine_conservative_and_nonconservative_fluxes(surface_integral.surface_flux, equations), equations, - surface_integral, SolverT, u_mpi_interfaces, + surface_integral, dg, cache, interface_index, normal_direction, interface_i_node_index, interface_j_node_index, local_side, @@ -360,14 +310,13 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, end Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, - MeshT::Type{<:Union{P4estMeshParallel{3}, - T8codeMeshParallel{3}}}, + mesh::Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}, have_nonconservative_terms::True, combine_conservative_and_nonconservative_fluxes::False, equations, - surface_integral, - SolverT::Type{<:DG}, - u_mpi_interfaces, + surface_integral, dg::DG, + cache, interface_index, normal_direction, interface_i_node_index, @@ -377,9 +326,10 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, surface_j_node_index, local_direction_index, local_element_index) + @unpack u = cache.mpi_interfaces surface_flux, nonconservative_flux = surface_integral.surface_flux - u_ll, u_rr = get_surface_node_vars(u_mpi_interfaces, equations, SolverT, + u_ll, u_rr = get_surface_node_vars(u, equations, dg, interface_i_node_index, interface_j_node_index, interface_index) @@ -405,14 +355,13 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, end Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, - MeshT::Type{<:Union{P4estMeshParallel{3}, - T8codeMeshParallel{3}}}, + mesh::Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}, have_nonconservative_terms::True, combine_conservative_and_nonconservative_fluxes::True, equations, - surface_integral, - SolverT::Type{<:DG}, - u_mpi_interfaces, + surface_integral, dg::DG, + cache, interface_index, normal_direction, interface_i_node_index, @@ -422,9 +371,10 @@ Base.@propagate_inbounds function calc_mpi_interface_flux!(surface_flux_values, surface_j_node_index, local_direction_index, local_element_index) + @unpack u = cache.mpi_interfaces surface_flux = surface_integral.surface_flux - u_ll, u_rr = get_surface_node_vars(u_mpi_interfaces, equations, SolverT, + u_ll, u_rr = get_surface_node_vars(u, equations, dg, interface_i_node_index, interface_j_node_index, interface_index) From 297eabe1c46359883cd64aa4799b3f42175a56dc Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 19 Oct 2026 10:35:37 +0200 Subject: [PATCH 33/46] docs --- NEWS.md | 3 ++- docs/src/heterogeneous.md | 10 +++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/NEWS.md b/NEWS.md index 97203a460ae..6c96bacfa41 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,7 +9,8 @@ for human readability. #### Added - Added multi-GPU support for 2D and 3D `P4estMesh`es, i.e., GPU computations with MPI - ([#3240]). A local CUDA/ROCm-aware MPI installation is required. For configuration, see + ([#3240], [#3364]). A local CUDA/ROCm-aware MPI installation is required. + For configuration, see [Parallelization](https://trixi-framework.org/TrixiDocumentation/stable/parallelization/#parallel_system_MPI) in the docs. Currently, only hyperbolic problems (no `SemidiscretizationHyperbolicParabolic`) on conforming meshes (no mortars, thus also no AMR) are supported. diff --git a/docs/src/heterogeneous.md b/docs/src/heterogeneous.md index 5249d07addb..baea379ab5a 100644 --- a/docs/src/heterogeneous.md +++ b/docs/src/heterogeneous.md @@ -57,12 +57,12 @@ are supported on conforming (non-adaptive) `P4estMesh{2}`s and `P4estMesh{3}`s a see, e.g., `examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl` and `examples/p4est_3d_dgsem/elixir_navierstokes_convergence.jl`. -Hyperbolic problems on conforming `P4estMesh`es in 2D can also be run on multiple GPUs -with MPI, using one MPI rank per GPU. The setup is the same as above, and the elixir is -started with MPI as described in [Parallelization](@ref parallel_system_MPI). +Hyperbolic problems on conforming `P4estMesh`es in 2D and 3D can also be run on +multiple GPUs with MPI, using one MPI rank per GPU. The setup is the same as above, and +the elixir is started with MPI as described in [Parallelization](@ref parallel_system_MPI). This requires a local CUDA/ROCm-aware MPI installation, since the MPI buffers are stored -on the GPU and passed directly to MPI. Mixed hyperbolic-parabolic problems, 3D meshes, -and non-conforming meshes (mortars, AMR) are not supported with multiple GPUs yet. +on the GPU and passed directly to MPI. Mixed hyperbolic-parabolic problems and +non-conforming meshes (mortars, AMR) are not supported with multiple GPUs yet. See the section on [GPU-aware MPI support](https://juliaparallel.org/MPI.jl/stable/usage/#GPU-aware-MPI-support) in the MPI.jl documentation for details, including example scripts such as From 26e6ebd3a6c094fb2c3f22f3c0e6881f963b6cb7 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 19 Oct 2026 10:36:22 +0200 Subject: [PATCH 34/46] dimension agnostic start_mpi_send! and finish_mpi_receive! --- src/solvers/dgsem_p4est/dg_parallel.jl | 1 + src/solvers/dgsem_p4est/dg_parallel_gpu.jl | 36 ++++++++++++++++++++++ 2 files changed, 37 insertions(+) create mode 100644 src/solvers/dgsem_p4est/dg_parallel_gpu.jl diff --git a/src/solvers/dgsem_p4est/dg_parallel.jl b/src/solvers/dgsem_p4est/dg_parallel.jl index 196d45fad4a..7d9bef03035 100644 --- a/src/solvers/dgsem_p4est/dg_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_parallel.jl @@ -791,6 +791,7 @@ Base.@propagate_inbounds function get_normal_direction(container::Union{P4estMPI Val(ndims(container)))) end +include("dg_parallel_gpu.jl") include("dg_2d_parallel.jl") include("dg_2d_parallel_gpu.jl") include("dg_3d_parallel.jl") diff --git a/src/solvers/dgsem_p4est/dg_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_parallel_gpu.jl new file mode 100644 index 00000000000..f55de4fbff5 --- /dev/null +++ b/src/solvers/dgsem_p4est/dg_parallel_gpu.jl @@ -0,0 +1,36 @@ +# 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 + +# The MPI send buffer is filled in `prolong2mpiinterfaces!`, and the kernels have +# finished when this is called, see `rhs_hyperbolic!`. +# This function is dimension-agnostic and also used in 3D. +function start_mpi_send!(backend::Backend, mpi_cache::P4estMPICache, + mesh::Union{P4estMeshParallel, T8codeMeshParallel}, + equations, dg, cache) + # TODO GPU: MPI mortars are not yet implemented + @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" + + for (index, rank) in enumerate(mpi_cache.mpi_neighbor_ranks) + send_buffer = mpi_buffer(get_mpi_send_buffer(mpi_cache, index)) + mpi_cache.mpi_send_requests[index] = MPI.Isend(send_buffer, rank, mpi_rank(), + mpi_comm()) + end + + return nothing +end + +# The received data is not unpacked but read directly from the MPI receive buffer +# by `calc_mpi_interface_flux!`. Thus, we can wait for all data at once. +# This function is dimension-agnostic and also used in 3D. +function finish_mpi_receive!(backend::Backend, mpi_cache::P4estMPICache, + mesh::Union{P4estMeshParallel, T8codeMeshParallel}, + equations, dg, cache) + MPI.Waitall(mpi_cache.mpi_recv_requests) + + return nothing +end +end # muladd \ No newline at end of file From dbd634a544e2ea4acdfb1414ab91e0c906c5106e Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 19 Oct 2026 10:40:13 +0200 Subject: [PATCH 35/46] fix merge --- src/callbacks_step/analysis_dg3d_parallel.jl | 16 +--------------- 1 file changed, 1 insertion(+), 15 deletions(-) diff --git a/src/callbacks_step/analysis_dg3d_parallel.jl b/src/callbacks_step/analysis_dg3d_parallel.jl index bdc7ad0661b..0f41120f0d2 100644 --- a/src/callbacks_step/analysis_dg3d_parallel.jl +++ b/src/callbacks_step/analysis_dg3d_parallel.jl @@ -123,20 +123,6 @@ function integrate_via_indices(func::Func, backend::Backend, u, return (local_integral, local_total_volume) end - global_integral = MPI.Reduce!(Ref(integral), +, mpi_root(), mpi_comm()) - total_volume = MPI.Reduce(volume, +, mpi_root(), mpi_comm()) - if mpi_isroot() - integral = convert(typeof(integral), global_integral[]) - else - integral = convert(typeof(integral), NaN * integral) - total_volume = volume # non-root processes receive nothing from reduce -> overwrite - end - - # Normalize with total volume - if normalize - integral = integral / total_volume - end - - return integral + return reduce_integral_and_volume_mpi(integral, volume, normalize) end end # @muladd From 9d079c3da7cbb8c8ade4c4de94df2a6fc26ae941 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 19 Oct 2026 10:57:03 +0200 Subject: [PATCH 36/46] tests --- test/test_kernelabstractions_mpi.jl | 77 +++++++++++++++++++++++++++++ 1 file changed, 77 insertions(+) diff --git a/test/test_kernelabstractions_mpi.jl b/test/test_kernelabstractions_mpi.jl index 0629809e9bc..5018b7453d9 100644 --- a/test/test_kernelabstractions_mpi.jl +++ b/test/test_kernelabstractions_mpi.jl @@ -83,3 +83,80 @@ end ], tspan=(0.0, 0.1)) end + +@testitem "KernelAbstractions MPI CPU 3D: elixir_advection_basic.jl" setup=[ + Setup, + KernelAbstractionsMPIExamples +] tags=[:kernelabstractions_mpi] begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem", + "elixir_advection_basic.jl"), + # Expected errors are exactly the same as with TreeMesh! + l2=[0.00016263963870641478], + linf=[0.0014537194925779984]) + # The elixir adapts `semi` before creating the ODE problem + @test semi === ode.p + @test real(semi.solver) == Float64 + @test real(semi.solver.basis) == Float64 + @test real(semi.solver.mortar) == Float64 + # TODO: `mesh` is currently not `adapt`ed correctly + @test real(semi.mesh) == Float64 + @test typeof(semi.equations.advection_velocity) == SVector{3, Float64} + + @test ode.u0 isa Array + @test semi.solver.basis.derivative_matrix isa Array + + @test Trixi.storage_type(semi.cache.elements) === Array + @test Trixi.storage_type(semi.cache.interfaces) === Array + @test Trixi.storage_type(semi.cache.boundaries) === Array + @test Trixi.storage_type(semi.cache.mortars) === Array + + # Ensure that the RHS computation overwrites existing data in `du` correctly. + u_ode = copy(ode.u0) + du_ode = similar(u_ode) + fill!(du_ode, convert(eltype(du_ode), NaN)) + Trixi.rhs_hyperbolic!(du_ode, u_ode, semi, first(ode.tspan)) + @test all(isfinite, du_ode) +end + +# The unstructured mesh contains MPI interfaces with flipped orientations +@testitem "KernelAbstractions MPI CPU 3D: elixir_advection_unstructured_curved.jl" setup=[Setup, KernelAbstractionsMPIExamples] tags=[:kernelabstractions_mpi] begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem", + "elixir_advection_unstructured_curved.jl"), + # Expected errors agree with the non-MPI version up to round-off + l2=[0.0004750004258546538], + linf=[0.026527551737137167]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000) +end + +@testitem "KernelAbstractions MPI CPU 3D: elixir_mhd_alfven_wave_combined_fluxes_nonperiodic.jl native" setup=[ + Setup, + KernelAbstractionsMPIExamples +] tags=[:kernelabstractions_mpi] begin + @test_trixi_include(joinpath(EXAMPLES_DIR,"p4est_3d_dgsem", + "elixir_mhd_alfven_wave_combined_fluxes_nonperiodic.jl"), + l2=[ + 0.00020523607157908708, + 0.000636037375360857, + 0.0002768722721598528, + 0.000780404043245244, + 0.000709669882361381, + 0.0006335012491207207, + 0.0003174622818755188, + 0.0007723389125485991, + 0.00023680042414966782 + ], + linf=[ + 0.0011634542712232188, + 0.004143013256670491, + 0.0016758874736150725, + 0.004705809053315052, + 0.0043527912850840345, + 0.00389837136357718, + 0.0014571305070021356, + 0.00445327133732587, + 0.0011999511085407923 + ], + tspan=(0.0, 0.1)) +end From d0c3b7b047248d83affc9a1c9ff40bd7665262cc Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 19 Oct 2026 10:59:01 +0200 Subject: [PATCH 37/46] revert --- src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl index 8c9ff4bd843..772e15765d9 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl @@ -93,7 +93,7 @@ function rhs_hyperbolic!(backend::Backend, end # Finish to receive MPI data - @trixi_timeit_ext backend timer() "finish MPI receive" begin + @trixi_timeit timer() "finish MPI receive" begin finish_mpi_receive!(backend, cache.mpi_cache, mesh, equations, dg, cache) end From 084085eee76944a166e2507bde8ec15b8f3ed980 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 19 Oct 2026 10:59:11 +0200 Subject: [PATCH 38/46] moved --- src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl | 27 ------------------- 1 file changed, 27 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl index 772e15765d9..e3da66beba4 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl @@ -206,33 +206,6 @@ Base.@propagate_inbounds function get_mpi_buffer_node_vars(mpi_buffer, equations Val(nvariables(equations)))) end -# The MPI send buffer is filled in `prolong2mpiinterfaces!`, and the kernels have -# finished when this is called, see `rhs_hyperbolic!`. -function start_mpi_send!(backend::Backend, mpi_cache::P4estMPICache, - mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, - equations, dg, cache) - # TODO GPU: MPI mortars are not yet implemented - @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" - - for (index, rank) in enumerate(mpi_cache.mpi_neighbor_ranks) - send_buffer = mpi_buffer(get_mpi_send_buffer(mpi_cache, index)) - mpi_cache.mpi_send_requests[index] = MPI.Isend(send_buffer, rank, mpi_rank(), - mpi_comm()) - end - - return nothing -end - -# The received data is not unpacked but read directly from the MPI receive buffer -# by `calc_mpi_interface_flux!`. Thus, we can wait for all data at once. -function finish_mpi_receive!(backend::Backend, mpi_cache::P4estMPICache, - mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, - equations, dg, cache) - MPI.Waitall(mpi_cache.mpi_recv_requests) - - return nothing -end - # Compute the surface fluxes at all MPI interfaces. In contrast to the CPU version, # the local data is taken directly from `u` and the remote data directly from the # MPI receive buffer. From cbf6f6896b26341110eb0890246752809e7a12b5 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 19 Oct 2026 11:59:00 +0200 Subject: [PATCH 39/46] align with 2D --- src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl | 9 - src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl | 417 ++++++++++++------ src/solvers/dgsem_p4est/dg_parallel_gpu.jl | 11 +- 3 files changed, 286 insertions(+), 151 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl index e3da66beba4..7529d026a90 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl @@ -197,15 +197,6 @@ Base.@propagate_inbounds function mpi_buffer_node_offset(mpi_interfaces_buffer_o return mpi_interfaces_buffer_offsets[interface] + (i - 1) * nvariables(equations) end -Base.@propagate_inbounds function get_mpi_buffer_node_vars(mpi_buffer, equations, - offset) - # Explicit bounds check, which can be removed by calling this function with `@inbounds` - @boundscheck checkbounds(mpi_buffer, offset .+ eachvariable(equations)) - # Assume inbounds access now - return SVector(ntuple(@inline(v->@inbounds mpi_buffer[offset + v]), - Val(nvariables(equations)))) -end - # Compute the surface fluxes at all MPI interfaces. In contrast to the CPU version, # the local data is taken directly from `u` and the remote data directly from the # MPI receive buffer. diff --git a/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl index cf4415b2cc8..7f06beea679 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl @@ -5,19 +5,22 @@ @muladd begin #! format: noindent +# In contrast to the CPU version, the GPU version does not use the MPI interface +# container `cache.mpi_interfaces.u`. Instead, `prolong2mpiinterfaces!` writes the +# local data directly into the MPI send buffer and `calc_mpi_interface_flux!` reads +# the remote data directly from the MPI receive buffer, while the local data is read +# from `u`. Both kernels handle the MPI interfaces of all neighbor ranks at once and +# parallelize over all nodes of these interfaces. function rhs_hyperbolic!(backend::Backend, du, u, t, mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, equations, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} - # Start to receive MPI data - @trixi_timeit timer() "start MPI receive" start_mpi_receive!(cache.mpi_cache) - - # Prolong solution to MPI interfaces + # Prolong solution to MPI interfaces, i.e., fill the MPI send buffer @trixi_timeit_ext backend timer() "prolong2mpiinterfaces" begin - prolong2mpiinterfaces!(backend, cache, u, mesh, equations, dg.surface_integral, - dg) + prolong2mpiinterfaces!(backend, cache, u, mesh, equations, + dg.surface_integral, dg) end # Prolong solution to MPI mortars @@ -26,8 +29,25 @@ function rhs_hyperbolic!(backend::Backend, dg.mortar, dg) end + # Wait for the kernels above to finish before sending their results via MPI. + # This also guarantees that the MPI interface flux kernel of the previous call, + # which reads from the MPI receive buffer, has finished. Thus, we can only start + # to receive MPI data afterwards, in contrast to the CPU version. + # MPI is not ordered with respect to the kernels queued on the GPU, which run + # asynchronously. Thus, all kernels writing to the MPI send buffer must be launched + # before this synchronization (and after `finish_mpi_send!` of the previous call), + # and all kernels reading from the MPI receive buffer must be launched after + # `finish_mpi_receive!`. Keep this in mind when adding kernels using the MPI + # buffers (e.g., for MPI mortars). Missing synchronization cannot be detected by + # tests with the CPU backend of KernelAbstractions.jl, which runs kernels + # synchronously. + @trixi_timeit timer() "synchronize" KernelAbstractions.synchronize(backend) + + # Start to receive MPI data + @trixi_timeit timer() "start MPI receive" start_mpi_receive!(cache.mpi_cache) + # Start to send MPI data - @trixi_timeit_ext backend timer() "start MPI send" begin + @trixi_timeit timer() "start MPI send" begin start_mpi_send!(backend, cache.mpi_cache, mesh, equations, dg, cache) end @@ -73,13 +93,13 @@ function rhs_hyperbolic!(backend::Backend, end # Finish to receive MPI data - @trixi_timeit_ext backend timer() "finish MPI receive" begin + @trixi_timeit timer() "finish MPI receive" begin finish_mpi_receive!(backend, cache.mpi_cache, mesh, equations, dg, cache) end # Calculate MPI interface fluxes @trixi_timeit_ext backend timer() "MPI interface flux" begin - calc_mpi_interface_flux!(backend, cache.elements.surface_flux_values, mesh, + calc_mpi_interface_flux!(backend, cache.elements.surface_flux_values, u, mesh, have_nonconservative_terms(equations), equations, dg.surface_integral, dg, cache) end @@ -107,191 +127,306 @@ function rhs_hyperbolic!(backend::Backend, return nothing end +# Copy the local solution data at all MPI interfaces directly into the MPI send buffer. +# The data of each interface node is stored at the same position as in the CPU version, +# which copies `cache.mpi_interfaces.u[local_side, :, i, j, interface]` into the buffer. function prolong2mpiinterfaces!(backend::Backend, cache, u, mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, equations, surface_integral, dg::DG) nmpiinterfaces(dg, cache) == 0 && return nothing - @unpack local_sides, local_neighbor_ids, node_indices = cache.mpi_interfaces + @unpack local_neighbor_ids, node_indices = cache.mpi_interfaces + @unpack mpi_send_buffer, mpi_interfaces_buffer_offsets = cache.mpi_cache index_range = eachnode(dg) - variables_range = eachvariable(equations) + + # Explicit bounds check, which allows us to assume inbounds access in the kernel. + # The offsets are constructed in `init_mpi_buffers` such that the data of all + # MPI interfaces fits into the MPI buffers. + @boundscheck begin + check_axes(u, mesh, equations, dg, cache) + check_axes(cache.mpi_interfaces, equations, dg, cache) + check_axes(mpi_interfaces_buffer_offsets, (eachmpiinterface(dg, cache),)) + checkbounds(mpi_send_buffer, + 1:(nvariables(equations) * nnodes(dg)^2 * + nmpiinterfaces(dg, cache))) + end kernel! = prolong2mpiinterfaces_KAkernel!(backend) - kernel!(cache.mpi_interfaces.u, typeof(mesh), local_sides, local_neighbor_ids, - node_indices, index_range, variables_range, u, - ndrange = nmpiinterfaces(dg, cache)) + kernel!(mpi_send_buffer, u, typeof(mesh), equations, dg, local_neighbor_ids, + node_indices, mpi_interfaces_buffer_offsets, index_range, + ndrange = (nnodes(dg), nnodes(dg), nmpiinterfaces(dg, cache))) return nothing end -@kernel inbounds=true function prolong2mpiinterfaces_KAkernel!(mpi_interfaces_u, +@kernel inbounds=true function prolong2mpiinterfaces_KAkernel!(mpi_send_buffer, u, MeshT::Type{<:Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}}, - local_sides, + equations, dg, local_neighbor_ids, node_indices, - index_range, - variables_range, u) - interface = @index(Global) - prolong2mpiinterfaces_per_interface!(mpi_interfaces_u, MeshT, interface, - local_sides, local_neighbor_ids, node_indices, - index_range, variables_range, u) + mpi_interfaces_buffer_offsets, + index_range) + i, j, interface = @index(Global, NTuple) + n = length(index_range) + + # Copy solution data from the local element using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + # Note that in the current implementation, the interface will be + # "aligned at the primary element", i.e., the indices of the primary side + # will always run forwards. + local_element = local_neighbor_ids[interface] + local_indices = node_indices[interface] + + i_element_start, i_element_step_i, i_element_step_j = index_to_start_step_3d(local_indices[1], + index_range) + j_element_start, j_element_step_i, j_element_step_j = index_to_start_step_3d(local_indices[2], + index_range) + k_element_start, k_element_step_i, k_element_step_j = index_to_start_step_3d(local_indices[3], + index_range) + + i_element = delayed_index_3d(i_element_start, i_element_step_i, i_element_step_j, + i, j, n) + j_element = delayed_index_3d(j_element_start, j_element_step_i, j_element_step_j, + i, j, n) + k_element = delayed_index_3d(k_element_start, k_element_step_i, k_element_step_j, + i, j, n) + + u_local = get_node_vars(u, equations, dg, i_element, j_element, k_element, + local_element) + offset = mpi_buffer_node_offset(mpi_interfaces_buffer_offsets, equations, dg, + i, j, interface) + for v in eachvariable(equations) + mpi_send_buffer[offset + v] = u_local[v] + end +end + +# Offset of the data of node `(i, j)` of an MPI interface in the MPI send/receive buffer +Base.@propagate_inbounds function mpi_buffer_node_offset(mpi_interfaces_buffer_offsets, + equations, dg, i, j, interface) + return mpi_interfaces_buffer_offsets[interface] + (i + (j - 1) * nnodes(dg) - 1) * nvariables(equations) end -function calc_mpi_interface_flux!(backend::Backend, surface_flux_values, +# Compute the surface fluxes at all MPI interfaces. In contrast to the CPU version, +# the local data is taken directly from `u` and the remote data directly from the +# MPI receive buffer. +function calc_mpi_interface_flux!(backend::Backend, surface_flux_values, u, mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, have_nonconservative_terms, equations, surface_integral, dg::DG, cache) nmpiinterfaces(dg, cache) == 0 && return nothing - @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces - @unpack contravariant_vectors = cache.elements - @unpack u = cache.mpi_interfaces + @unpack local_neighbor_ids, node_indices, local_sides, normal_directions = cache.mpi_interfaces + @unpack mpi_recv_buffer, mpi_interfaces_buffer_offsets = cache.mpi_cache index_range = eachnode(dg) + # Explicit bounds check, which allows us to assume inbounds access in the kernel. + # The offsets are constructed in `init_mpi_buffers` such that the data of all + # MPI interfaces fits into the MPI buffers. + @boundscheck begin + check_axes(u, mesh, equations, dg, cache) + check_axes(cache.mpi_interfaces, equations, dg, cache) + check_axes_surface_flux_values(surface_flux_values, mesh, equations, dg, cache) + check_axes(mpi_interfaces_buffer_offsets, (eachmpiinterface(dg, cache),)) + checkbounds(mpi_recv_buffer, + 1:(nvariables(equations) * nnodes(dg)^2 * + nmpiinterfaces(dg, cache))) + end + kernel! = calc_mpi_interface_flux_KAkernel!(backend) - kernel!(surface_flux_values, typeof(mesh), - have_nonconservative_terms, equations, - surface_integral, typeof(dg), - local_neighbor_ids, node_indices, local_sides, contravariant_vectors, - u, index_range, - ndrange = nmpiinterfaces(dg, cache)) + kernel!(surface_flux_values, u, mpi_recv_buffer, typeof(mesh), + have_nonconservative_terms, equations, surface_integral, dg, + local_neighbor_ids, node_indices, local_sides, + mpi_interfaces_buffer_offsets, normal_directions, index_range, + ndrange = (nnodes(dg), nnodes(dg), nmpiinterfaces(dg, cache))) return nothing end -@kernel inbounds=true function calc_mpi_interface_flux_KAkernel!(surface_flux_values, +@kernel inbounds=true function calc_mpi_interface_flux_KAkernel!(surface_flux_values, u, + mpi_recv_buffer, MeshT::Type{<:Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}}, have_nonconservative_terms, equations, surface_integral, - SolverT::Type{<:DG}, + dg, local_neighbor_ids, node_indices, local_sides, - contravariant_vectors, - u, + mpi_interfaces_buffer_offsets, + normal_directions, index_range) - interface = @index(Global) - calc_mpi_interface_flux_per_interface!(surface_flux_values, MeshT, - have_nonconservative_terms, equations, - surface_integral, SolverT, - local_neighbor_ids, node_indices, - local_sides, contravariant_vectors, u, - index_range, interface) -end + i, j, interface = @index(Global, NTuple) + n = length(index_range) -# For GPU backends MPI mortars are not yet implemented -function prolong2mpimortars!(backend::Backend, cache, u, - mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, - equations, mortar, dg::DG) - @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" - return nothing -end - -# For GPU backends MPI mortars are not yet implemented -function calc_mpi_mortar_flux!(backend::Backend, surface_flux_values, - mesh::Union{P4estMeshParallel{3}, - T8codeMeshParallel{3}}, - have_nonconservative_terms, equations, mortar, - surface_integral, dg::DG, cache) - @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" - return nothing -end + # Get element and side index information on the local element + local_element = local_neighbor_ids[interface] + local_indices = node_indices[interface] + local_direction = indices2direction(local_indices) + local_side = local_sides[interface] -function start_mpi_send!(backend::Backend, mpi_cache::P4estMPICache, - mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, - equations, dg, cache) - # TODO GPU: MPI mortars are not yet implemented - @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" + # Create the local i,j,k indexing on the local element used to pull the local data + i_element_start, i_element_step_i, i_element_step_j = index_to_start_step_3d(local_indices[1], + index_range) + j_element_start, j_element_step_i, j_element_step_j = index_to_start_step_3d(local_indices[2], + index_range) + k_element_start, k_element_step_i, k_element_step_j = index_to_start_step_3d(local_indices[3], + index_range) + + i_element = delayed_index_3d(i_element_start, i_element_step_i, i_element_step_j, + i, j, n) + j_element = delayed_index_3d(j_element_start, j_element_step_i, j_element_step_j, + i, j, n) + k_element = delayed_index_3d(k_element_start, k_element_step_i, k_element_step_j, + i, j, n) + + # The surface flux storage must be indexed in alignment with the local element + # indexing + local_surface_indices = surface_indices(local_indices) + i_surface_start, i_surface_step_i, i_surface_step_j = index_to_start_step_3d(local_surface_indices[1], + index_range) + j_surface_start, j_surface_step_i, j_surface_step_j = index_to_start_step_3d(local_surface_indices[2], + index_range) + i_surface = delayed_index_3d(i_surface_start, i_surface_step_i, i_surface_step_j, + i, j, n) + j_surface = delayed_index_3d(j_surface_start, j_surface_step_i, j_surface_step_j, + i, j, n) + + # The interface is "aligned at the primary element", i.e., `u_ll` is the solution + # on the primary side and `u_rr` the solution on the secondary side + u_local = get_node_vars(u, equations, dg, i_element, j_element, k_element, + local_element) + offset = mpi_buffer_node_offset(mpi_interfaces_buffer_offsets, equations, dg, + i, j, interface) + u_remote = get_mpi_buffer_node_vars(mpi_recv_buffer, equations, offset) + if local_side == 1 + u_ll, u_rr = u_local, u_remote + else # local_side == 2 + u_ll, u_rr = u_remote, u_local + end - @unpack mpi_neighbor_ranks, mpi_neighbor_interfaces = mpi_cache - @unpack mpi_send_buffers, mpi_send_requests = mpi_cache - @unpack local_sides, u = cache.mpi_interfaces + # Use the normal direction of the primary element, which is available on both + # MPI ranks. Thus, both ranks compute the same numerical flux, as in the serial + # `calc_interface_flux!`, see `exchange_normal_directions!`. + normal_direction = SVector(normal_directions[1, i, j, interface], + normal_directions[2, i, j, interface], + normal_directions[3, i, j, interface]) + + calc_mpi_interface_flux_per_node!(surface_flux_values, MeshT, + have_nonconservative_terms, equations, + surface_integral, u_ll, u_rr, normal_direction, + local_side, i_surface, j_surface, local_direction, + local_element) +end - kernel! = start_mpi_send_KAkernel!(backend) +# Inlined version of the interface flux computation for conservation laws +Base.@propagate_inbounds function calc_mpi_interface_flux_per_node!(surface_flux_values, + MeshT::Type{<:Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}}, + have_nonconservative_terms::False, + equations, + surface_integral, + u_ll, u_rr, + normal_direction, + local_side, + surface_i_node_index, + surface_j_node_index, + local_direction_index, + local_element_index) + @unpack surface_flux = surface_integral + + # `normal_direction` is the normal direction of the primary element. + # The flux entering the primary element leaves the secondary element + # (note the `-` for the `else` case). + if local_side == 1 + flux_ = surface_flux(u_ll, u_rr, normal_direction, equations) + else # local_side == 2 + flux_ = -surface_flux(u_ll, u_rr, normal_direction, equations) + end - for (rank_index, _) in enumerate(mpi_neighbor_ranks) - send_buffer = mpi_send_buffers[rank_index] - neighbor_interfaces = mpi_neighbor_interfaces[rank_index] - if !isempty(neighbor_interfaces) - kernel!(send_buffer, neighbor_interfaces, local_sides, u, - Val(nvariables(equations)), Val(ndims(mesh)), - ndrange = (nnodes(dg), nnodes(dg), length(neighbor_interfaces))) - end + for v in eachvariable(equations) + surface_flux_values[v, surface_i_node_index, surface_j_node_index, + local_direction_index, local_element_index] = flux_[v] end - # wait for the kernels to return before sending the buffers - KernelAbstractions.synchronize(backend) + return nothing +end - for (rank_index, neighbor_rank) in enumerate(mpi_neighbor_ranks) - send_buffer = mpi_send_buffers[rank_index] - mpi_send_requests[rank_index] = MPI.Isend(send_buffer, neighbor_rank, - mpi_rank(), mpi_comm()) - end +# Inlined version of the interface flux computation for non-conservative equations. +# As for the other surface fluxes on GPUs, see `prolong2interfaces_and_calc_interface_flux_per_node!`, +# only surface fluxes combining the conservative and nonconservative terms are supported. +Base.@propagate_inbounds function calc_mpi_interface_flux_per_node!(surface_flux_values, + MeshT::Type{<:Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}}, + have_nonconservative_terms::True, + equations, + surface_integral, + u_ll, u_rr, + normal_direction, + local_side, + surface_i_node_index, + surface_j_node_index, + local_direction_index, + local_element_index) + calc_mpi_interface_flux_per_node!(surface_flux_values, MeshT, + have_nonconservative_terms, + combine_conservative_and_nonconservative_fluxes(surface_integral.surface_flux, + equations), + equations, surface_integral, u_ll, u_rr, + normal_direction, local_side, + surface_i_node_index, surface_j_node_index, + local_direction_index, local_element_index) + return nothing end -@kernel inbounds=true function start_mpi_send_KAkernel!(send_buffer, - neighbor_interfaces, - local_sides, - u_mpi_interfaces, ::Val{NVARS}, - ::Val{3}) where {NVARS} - i, j, index_interface = @index(Global, NTuple) - index_linear = @index(Global, Linear) - - # The data of each interface is stored contiguously in the buffer, with the - # variables running fastest, followed by the surface node indices `i` and `j`. - buffer_offset = (index_linear - 1) * NVARS - interface = neighbor_interfaces[index_interface] - local_side = local_sides[interface] +Base.@propagate_inbounds function calc_mpi_interface_flux_per_node!(surface_flux_values, + MeshT::Type{<:Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}}, + have_nonconservative_terms::True, + combine_conservative_and_nonconservative_fluxes::True, + equations, + surface_integral, + u_ll, u_rr, + normal_direction, + local_side, + surface_i_node_index, + surface_j_node_index, + local_direction_index, + local_element_index) + @unpack surface_flux = surface_integral + + # Compute the flux including the non-conservative term for this side of the interface. + # `normal_direction` is the normal direction of the primary element + # (note the `-` for the `else` case). + flux_left, flux_right = surface_flux(u_ll, u_rr, normal_direction, equations) + if local_side == 1 + flux_ = flux_left + else # local_side == 2 + flux_ = -flux_right + end - for v in 1:NVARS - send_buffer[buffer_offset + v] = u_mpi_interfaces[local_side, v, i, j, - interface] + for v in eachvariable(equations) + surface_flux_values[v, surface_i_node_index, surface_j_node_index, + local_direction_index, local_element_index] = flux_[v] end + + return nothing end -function finish_mpi_receive!(backend::Backend, mpi_cache::P4estMPICache, +# For GPU backends MPI mortars are not yet implemented +function prolong2mpimortars!(backend::Backend, cache, u, mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, - equations, dg, cache) - # TODO GPU: MPI mortars are not yet implemented + equations, mortar, dg::DG) @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" - - @unpack mpi_neighbor_interfaces = mpi_cache - @unpack mpi_recv_buffers, mpi_recv_requests = mpi_cache - @unpack local_sides, u = cache.mpi_interfaces - - kernel! = finish_mpi_receive_KAkernel!(backend) - - # Start receiving and unpack received data until all communication is finished - data = MPI.Waitany(mpi_recv_requests) - while data !== nothing - recv_buffer = mpi_recv_buffers[data] - neighbor_interfaces = mpi_neighbor_interfaces[data] - if !isempty(neighbor_interfaces) - kernel!(recv_buffer, neighbor_interfaces, local_sides, u, - Val(nvariables(equations)), Val(ndims(mesh)), - ndrange = (nnodes(dg), nnodes(dg), length(neighbor_interfaces))) - end - - data = MPI.Waitany(mpi_recv_requests) - end - # Wait for the last kernel to return ? - KernelAbstractions.synchronize(backend) + return nothing end -@kernel inbounds=true function finish_mpi_receive_KAkernel!(recv_buffer, - neighbor_interfaces, - local_sides, - u_mpi_interfaces, - ::Val{NVARS}, - ::Val{3}) where {NVARS} - i, j, index_interface = @index(Global, NTuple) - index_linear = @index(Global, Linear) - buffer_offset = (index_linear - 1) * NVARS - interface = neighbor_interfaces[index_interface] - remote_side = local_sides[interface] == 1 ? 2 : 1 - for v in 1:NVARS - u_mpi_interfaces[remote_side, v, i, j, interface] = recv_buffer[buffer_offset + v] - end -end +# For GPU backends MPI mortars are not yet implemented +function calc_mpi_mortar_flux!(backend::Backend, surface_flux_values, + mesh::Union{P4estMeshParallel{3}, + T8codeMeshParallel{3}}, + have_nonconservative_terms, equations, mortar, + surface_integral, dg::DG, cache) + @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" + return nothing end +end # @muladd diff --git a/src/solvers/dgsem_p4est/dg_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_parallel_gpu.jl index f55de4fbff5..7f8482fb204 100644 --- a/src/solvers/dgsem_p4est/dg_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_parallel_gpu.jl @@ -33,4 +33,13 @@ function finish_mpi_receive!(backend::Backend, mpi_cache::P4estMPICache, return nothing end -end # muladd \ No newline at end of file + +Base.@propagate_inbounds function get_mpi_buffer_node_vars(mpi_buffer, equations, + offset) + # Explicit bounds check, which can be removed by calling this function with `@inbounds` + @boundscheck checkbounds(mpi_buffer, offset .+ eachvariable(equations)) + # Assume inbounds access now + return SVector(ntuple(@inline(v->@inbounds mpi_buffer[offset + v]), + Val(nvariables(equations)))) +end +end # muladd From 8b002a99f2caf1c78a33abeb2e50c2cbc1f65081 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 19 Oct 2026 11:59:23 +0200 Subject: [PATCH 40/46] fmt --- src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl | 5 +++-- test/test_kernelabstractions_mpi.jl | 7 +++++-- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl index 7f06beea679..2f8a18d75bf 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl @@ -203,7 +203,8 @@ end # Offset of the data of node `(i, j)` of an MPI interface in the MPI send/receive buffer Base.@propagate_inbounds function mpi_buffer_node_offset(mpi_interfaces_buffer_offsets, equations, dg, i, j, interface) - return mpi_interfaces_buffer_offsets[interface] + (i + (j - 1) * nnodes(dg) - 1) * nvariables(equations) + return mpi_interfaces_buffer_offsets[interface] + + (i + (j - 1) * nnodes(dg) - 1) * nvariables(equations) end # Compute the surface fluxes at all MPI interfaces. In contrast to the CPU version, @@ -295,7 +296,7 @@ end # on the primary side and `u_rr` the solution on the secondary side u_local = get_node_vars(u, equations, dg, i_element, j_element, k_element, local_element) - offset = mpi_buffer_node_offset(mpi_interfaces_buffer_offsets, equations, dg, + offset = mpi_buffer_node_offset(mpi_interfaces_buffer_offsets, equations, dg, i, j, interface) u_remote = get_mpi_buffer_node_vars(mpi_recv_buffer, equations, offset) if local_side == 1 diff --git a/test/test_kernelabstractions_mpi.jl b/test/test_kernelabstractions_mpi.jl index 5018b7453d9..1055b6bfde8 100644 --- a/test/test_kernelabstractions_mpi.jl +++ b/test/test_kernelabstractions_mpi.jl @@ -119,7 +119,10 @@ end end # The unstructured mesh contains MPI interfaces with flipped orientations -@testitem "KernelAbstractions MPI CPU 3D: elixir_advection_unstructured_curved.jl" setup=[Setup, KernelAbstractionsMPIExamples] tags=[:kernelabstractions_mpi] begin +@testitem "KernelAbstractions MPI CPU 3D: elixir_advection_unstructured_curved.jl" setup=[ + Setup, + KernelAbstractionsMPIExamples +] tags=[:kernelabstractions_mpi] begin @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem", "elixir_advection_unstructured_curved.jl"), # Expected errors agree with the non-MPI version up to round-off @@ -134,7 +137,7 @@ end Setup, KernelAbstractionsMPIExamples ] tags=[:kernelabstractions_mpi] begin - @test_trixi_include(joinpath(EXAMPLES_DIR,"p4est_3d_dgsem", + @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem", "elixir_mhd_alfven_wave_combined_fluxes_nonperiodic.jl"), l2=[ 0.00020523607157908708, From 0bed854e35c84cf5348f4bb9483f767018241576 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 19 Oct 2026 13:31:49 +0200 Subject: [PATCH 41/46] adapt reference values --- test/test_kernelabstractions_mpi.jl | 36 ++++++++++++++--------------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/test/test_kernelabstractions_mpi.jl b/test/test_kernelabstractions_mpi.jl index 1055b6bfde8..50d1aad51c1 100644 --- a/test/test_kernelabstractions_mpi.jl +++ b/test/test_kernelabstractions_mpi.jl @@ -140,26 +140,26 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem", "elixir_mhd_alfven_wave_combined_fluxes_nonperiodic.jl"), l2=[ - 0.00020523607157908708, - 0.000636037375360857, - 0.0002768722721598528, - 0.000780404043245244, - 0.000709669882361381, - 0.0006335012491207207, - 0.0003174622818755188, - 0.0007723389125485991, - 0.00023680042414966782 + 0.00016680585980070498, + 0.0005174249857687416, + 0.00027338410141732186, + 0.0006282227969118499, + 0.0006628677315607643, + 0.0005583953631900831, + 0.000181088291202701, + 0.0006112889578653593, + 0.00023674944227366533 ], linf=[ - 0.0011634542712232188, - 0.004143013256670491, - 0.0016758874736150725, - 0.004705809053315052, - 0.0043527912850840345, - 0.00389837136357718, - 0.0014571305070021356, - 0.00445327133732587, - 0.0011999511085407923 + 0.0009684118857454793, + 0.0032142734886018848, + 0.001806828008193248, + 0.0030945203920691444, + 0.0037003192368261217, + 0.0031660143839563615, + 0.0011578323492360543, + 0.003133093207396312, + 0.0008290351027180755 ], tspan=(0.0, 0.1)) end From 454a90eb44ed68d4fdc928e8cda4b93a014960e0 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 19 Oct 2026 15:38:09 +0200 Subject: [PATCH 42/46] no allocation test --- test/test_kernelabstractions_mpi.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/test_kernelabstractions_mpi.jl b/test/test_kernelabstractions_mpi.jl index 50d1aad51c1..7194cecfbc6 100644 --- a/test/test_kernelabstractions_mpi.jl +++ b/test/test_kernelabstractions_mpi.jl @@ -128,9 +128,6 @@ end # Expected errors agree with the non-MPI version up to round-off l2=[0.0004750004258546538], linf=[0.026527551737137167]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - @test_allocations(Trixi.rhs_hyperbolic!, semi, sol, 1000) end @testitem "KernelAbstractions MPI CPU 3D: elixir_mhd_alfven_wave_combined_fluxes_nonperiodic.jl native" setup=[ From e1609c1fac07cb2a8a7566b320cb0cbd8e3494a9 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Wed, 4 Nov 2026 10:19:56 +0100 Subject: [PATCH 43/46] address review comments [skip ci] --- src/solvers/dgsem_p4est/dg_parallel_gpu.jl | 10 ++++++---- test/test_kernelabstractions_mpi.jl | 2 +- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_parallel_gpu.jl index 7f8482fb204..b65c751670d 100644 --- a/src/solvers/dgsem_p4est/dg_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_parallel_gpu.jl @@ -7,9 +7,10 @@ # The MPI send buffer is filled in `prolong2mpiinterfaces!`, and the kernels have # finished when this is called, see `rhs_hyperbolic!`. -# This function is dimension-agnostic and also used in 3D. +# This function is dimension-agnostic and used in both 2D and 3D. function start_mpi_send!(backend::Backend, mpi_cache::P4estMPICache, - mesh::Union{P4estMeshParallel, T8codeMeshParallel}, + mesh::Union{P4estMeshParallel{2}, P4estMeshParallel{3}, + T8codeMeshParallel{2}, T8codeMeshParallel{3}}, equations, dg, cache) # TODO GPU: MPI mortars are not yet implemented @assert isempty(eachmpimortar(dg, cache)) "MPI mortars are not yet supported on GPU backends" @@ -25,9 +26,10 @@ end # The received data is not unpacked but read directly from the MPI receive buffer # by `calc_mpi_interface_flux!`. Thus, we can wait for all data at once. -# This function is dimension-agnostic and also used in 3D. +# This function is dimension-agnostic and used in both 2D and 3D. function finish_mpi_receive!(backend::Backend, mpi_cache::P4estMPICache, - mesh::Union{P4estMeshParallel, T8codeMeshParallel}, + mesh::Union{P4estMeshParallel{2}, P4estMeshParallel{3}, + T8codeMeshParallel{2}, T8codeMeshParallel{3}}, equations, dg, cache) MPI.Waitall(mpi_cache.mpi_recv_requests) diff --git a/test/test_kernelabstractions_mpi.jl b/test/test_kernelabstractions_mpi.jl index 7194cecfbc6..0ab1b719d77 100644 --- a/test/test_kernelabstractions_mpi.jl +++ b/test/test_kernelabstractions_mpi.jl @@ -130,7 +130,7 @@ end linf=[0.026527551737137167]) end -@testitem "KernelAbstractions MPI CPU 3D: elixir_mhd_alfven_wave_combined_fluxes_nonperiodic.jl native" setup=[ +@testitem "KernelAbstractions MPI CPU 3D: elixir_mhd_alfven_wave_combined_fluxes_nonperiodic.jl" setup=[ Setup, KernelAbstractionsMPIExamples ] tags=[:kernelabstractions_mpi] begin From 4a401cce0133d849c6415013b824924d2d537e6c Mon Sep 17 00:00:00 2001 From: Benedict <135045760+benegee@users.noreply.github.com> Date: Thu, 5 Nov 2026 12:44:35 +0100 Subject: [PATCH 44/46] Update src/callbacks_step/analysis_dg3d_parallel.jl Co-authored-by: Hendrik Ranocha --- src/callbacks_step/analysis_dg3d_parallel.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/callbacks_step/analysis_dg3d_parallel.jl b/src/callbacks_step/analysis_dg3d_parallel.jl index f79cb49a2fa..d06d54d9af2 100644 --- a/src/callbacks_step/analysis_dg3d_parallel.jl +++ b/src/callbacks_step/analysis_dg3d_parallel.jl @@ -96,9 +96,9 @@ function integrate_via_indices(func::Func, backend::Backend, u, # `func(u, 1,1,1,1, equations, dg, args...)` might access GPU memory # so we have to rely on the compiler to correctly infer the type of the integral here. # TODO: Technically we need device_promote_op here that "infers" the function within the context of the GPU. - integral0 = zero(Base.promote_op(func, typeof(u), Int, Int, Int, Int, - typeof(equations), typeof(dg), - map(typeof, args)...)) + integral0 = zero_return_type(func, typeof(u), Int, Int, Int, Int, + typeof(equations), typeof(dg), + map(typeof, args)...)) volume0 = zero(eltype(weights)) * zero(eltype(inverse_jacobian)) init = neutral = (integral0, volume0) From 5783cbea8994617d0a213f1db7f8ad1af9ff882f Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 5 Nov 2026 13:21:16 +0100 Subject: [PATCH 45/46] fixup --- src/callbacks_step/analysis_dg3d_parallel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/analysis_dg3d_parallel.jl b/src/callbacks_step/analysis_dg3d_parallel.jl index d06d54d9af2..a3d6c482715 100644 --- a/src/callbacks_step/analysis_dg3d_parallel.jl +++ b/src/callbacks_step/analysis_dg3d_parallel.jl @@ -98,7 +98,7 @@ function integrate_via_indices(func::Func, backend::Backend, u, # TODO: Technically we need device_promote_op here that "infers" the function within the context of the GPU. integral0 = zero_return_type(func, typeof(u), Int, Int, Int, Int, typeof(equations), typeof(dg), - map(typeof, args)...)) + map(typeof, args)...) volume0 = zero(eltype(weights)) * zero(eltype(inverse_jacobian)) init = neutral = (integral0, volume0) From 61e74e1b05f479ee182a6f45b93a4966bc946b00 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 5 Nov 2026 13:21:32 +0100 Subject: [PATCH 46/46] dimension agnostic rhs_hyperbolic --- src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl | 122 ----------------- src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl | 122 ----------------- src/solvers/dgsem_p4est/dg_parallel_gpu.jl | 123 ++++++++++++++++++ 3 files changed, 123 insertions(+), 244 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl index 7529d026a90..cf6bdb640db 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel_gpu.jl @@ -5,128 +5,6 @@ @muladd begin #! format: noindent -# In contrast to the CPU version, the GPU version does not use the MPI interface -# container `cache.mpi_interfaces.u`. Instead, `prolong2mpiinterfaces!` writes the -# local data directly into the MPI send buffer and `calc_mpi_interface_flux!` reads -# the remote data directly from the MPI receive buffer, while the local data is read -# from `u`. Both kernels handle the MPI interfaces of all neighbor ranks at once and -# parallelize over all nodes of these interfaces. -function rhs_hyperbolic!(backend::Backend, - du, u, t, - mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}}, - equations, - boundary_conditions, source_terms::Source, - dg::DG, cache) where {Source} - # Prolong solution to MPI interfaces, i.e., fill the MPI send buffer - @trixi_timeit_ext backend timer() "prolong2mpiinterfaces" begin - prolong2mpiinterfaces!(backend, cache, u, mesh, equations, - dg.surface_integral, dg) - end - - # Prolong solution to MPI mortars - @trixi_timeit_ext backend timer() "prolong2mpimortars" begin - prolong2mpimortars!(backend, cache, u, mesh, equations, - dg.mortar, dg) - end - - # Wait for the kernels above to finish before sending their results via MPI. - # This also guarantees that the MPI interface flux kernel of the previous call, - # which reads from the MPI receive buffer, has finished. Thus, we can only start - # to receive MPI data afterwards, in contrast to the CPU version. - # MPI is not ordered with respect to the kernels queued on the GPU, which run - # asynchronously. Thus, all kernels writing to the MPI send buffer must be launched - # before this synchronization (and after `finish_mpi_send!` of the previous call), - # and all kernels reading from the MPI receive buffer must be launched after - # `finish_mpi_receive!`. Keep this in mind when adding kernels using the MPI - # buffers (e.g., for MPI mortars). Missing synchronization cannot be detected by - # tests with the CPU backend of KernelAbstractions.jl, which runs kernels - # synchronously. - @trixi_timeit timer() "synchronize" KernelAbstractions.synchronize(backend) - - # Start to receive MPI data - @trixi_timeit timer() "start MPI receive" start_mpi_receive!(cache.mpi_cache) - - # Start to send MPI data - @trixi_timeit timer() "start MPI send" begin - start_mpi_send!(backend, cache.mpi_cache, mesh, equations, dg, cache) - end - - # Calculate volume integral - @trixi_timeit_ext backend timer() "volume integral" begin - calc_volume_integral!(backend, du, u, mesh, - have_nonconservative_terms(equations), equations, - dg.volume_integral, dg, cache) - end - - # Calculate interface fluxes - @trixi_timeit_ext backend timer() "prolong2interfaces + flux" begin - prolong2interfaces_and_calc_interface_flux!(backend, - cache.elements.surface_flux_values, - u, mesh, - have_nonconservative_terms(equations), - equations, - dg.surface_integral, dg, cache) - end - - # Prolong solution to boundaries - @trixi_timeit_ext backend timer() "prolong2boundaries" begin - prolong2boundaries!(backend, cache, u, mesh, equations, dg) - end - - # Calculate boundary fluxes - @trixi_timeit_ext backend timer() "boundary flux" begin - calc_boundary_flux!(backend, cache, t, boundary_conditions, mesh, equations, - dg.surface_integral, dg) - end - - # Prolong solution to mortars - @trixi_timeit_ext backend timer() "prolong2mortars" begin - prolong2mortars!(backend, cache, u, mesh, equations, - dg.mortar, dg) - end - - # Calculate mortar fluxes - @trixi_timeit_ext backend timer() "mortar flux" begin - calc_mortar_flux!(backend, cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, - dg.mortar, dg.surface_integral, dg, cache) - end - - # Finish to receive MPI data - @trixi_timeit timer() "finish MPI receive" begin - finish_mpi_receive!(backend, cache.mpi_cache, mesh, equations, dg, cache) - end - - # Calculate MPI interface fluxes - @trixi_timeit_ext backend timer() "MPI interface flux" begin - calc_mpi_interface_flux!(backend, cache.elements.surface_flux_values, u, mesh, - have_nonconservative_terms(equations), equations, - dg.surface_integral, dg, cache) - end - - # Calculate MPI mortar fluxes - @trixi_timeit_ext backend timer() "MPI mortar flux" begin - calc_mpi_mortar_flux!(backend, cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, - dg.mortar, dg.surface_integral, dg, cache) - end - - # Calculate surface integrals, apply Jacobian from mapping to reference element - # and calculate source terms - @trixi_timeit_ext backend timer() "surface, Jacobian + source terms" begin - calc_surface_integral_and_apply_jacobian_and_calc_sources!(backend, du, u, t, - source_terms, mesh, - equations, - dg.surface_integral, - dg, cache) - end - - # Finish to send MPI data - @trixi_timeit timer() "finish MPI send" finish_mpi_send!(cache.mpi_cache) - - return nothing -end - # Copy the local solution data at all MPI interfaces directly into the MPI send buffer. # The data of each interface node is stored at the same position as in the CPU version, # which copies `cache.mpi_interfaces.u[local_side, :, i, interface]` into the buffer. diff --git a/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl index 2f8a18d75bf..20617ee7cd9 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parallel_gpu.jl @@ -5,128 +5,6 @@ @muladd begin #! format: noindent -# In contrast to the CPU version, the GPU version does not use the MPI interface -# container `cache.mpi_interfaces.u`. Instead, `prolong2mpiinterfaces!` writes the -# local data directly into the MPI send buffer and `calc_mpi_interface_flux!` reads -# the remote data directly from the MPI receive buffer, while the local data is read -# from `u`. Both kernels handle the MPI interfaces of all neighbor ranks at once and -# parallelize over all nodes of these interfaces. -function rhs_hyperbolic!(backend::Backend, - du, u, t, - mesh::Union{P4estMeshParallel{3}, T8codeMeshParallel{3}}, - equations, - boundary_conditions, source_terms::Source, - dg::DG, cache) where {Source} - # Prolong solution to MPI interfaces, i.e., fill the MPI send buffer - @trixi_timeit_ext backend timer() "prolong2mpiinterfaces" begin - prolong2mpiinterfaces!(backend, cache, u, mesh, equations, - dg.surface_integral, dg) - end - - # Prolong solution to MPI mortars - @trixi_timeit_ext backend timer() "prolong2mpimortars" begin - prolong2mpimortars!(backend, cache, u, mesh, equations, - dg.mortar, dg) - end - - # Wait for the kernels above to finish before sending their results via MPI. - # This also guarantees that the MPI interface flux kernel of the previous call, - # which reads from the MPI receive buffer, has finished. Thus, we can only start - # to receive MPI data afterwards, in contrast to the CPU version. - # MPI is not ordered with respect to the kernels queued on the GPU, which run - # asynchronously. Thus, all kernels writing to the MPI send buffer must be launched - # before this synchronization (and after `finish_mpi_send!` of the previous call), - # and all kernels reading from the MPI receive buffer must be launched after - # `finish_mpi_receive!`. Keep this in mind when adding kernels using the MPI - # buffers (e.g., for MPI mortars). Missing synchronization cannot be detected by - # tests with the CPU backend of KernelAbstractions.jl, which runs kernels - # synchronously. - @trixi_timeit timer() "synchronize" KernelAbstractions.synchronize(backend) - - # Start to receive MPI data - @trixi_timeit timer() "start MPI receive" start_mpi_receive!(cache.mpi_cache) - - # Start to send MPI data - @trixi_timeit timer() "start MPI send" begin - start_mpi_send!(backend, cache.mpi_cache, mesh, equations, dg, cache) - end - - # Calculate volume integral - @trixi_timeit_ext backend timer() "volume integral" begin - calc_volume_integral!(backend, du, u, mesh, - have_nonconservative_terms(equations), equations, - dg.volume_integral, dg, cache) - end - - # Calculate interface fluxes - @trixi_timeit_ext backend timer() "prolong2interfaces + flux" begin - prolong2interfaces_and_calc_interface_flux!(backend, - cache.elements.surface_flux_values, - u, mesh, - have_nonconservative_terms(equations), - equations, - dg.surface_integral, dg, cache) - end - - # Prolong solution to boundaries - @trixi_timeit_ext backend timer() "prolong2boundaries" begin - prolong2boundaries!(backend, cache, u, mesh, equations, dg) - end - - # Calculate boundary fluxes - @trixi_timeit_ext backend timer() "boundary flux" begin - calc_boundary_flux!(backend, cache, t, boundary_conditions, mesh, equations, - dg.surface_integral, dg) - end - - # Prolong solution to mortars - @trixi_timeit_ext backend timer() "prolong2mortars" begin - prolong2mortars!(backend, cache, u, mesh, equations, - dg.mortar, dg) - end - - # Calculate mortar fluxes - @trixi_timeit_ext backend timer() "mortar flux" begin - calc_mortar_flux!(backend, cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, - dg.mortar, dg.surface_integral, dg, cache) - end - - # Finish to receive MPI data - @trixi_timeit timer() "finish MPI receive" begin - finish_mpi_receive!(backend, cache.mpi_cache, mesh, equations, dg, cache) - end - - # Calculate MPI interface fluxes - @trixi_timeit_ext backend timer() "MPI interface flux" begin - calc_mpi_interface_flux!(backend, cache.elements.surface_flux_values, u, mesh, - have_nonconservative_terms(equations), equations, - dg.surface_integral, dg, cache) - end - - # Calculate MPI mortar fluxes - @trixi_timeit_ext backend timer() "MPI mortar flux" begin - calc_mpi_mortar_flux!(backend, cache.elements.surface_flux_values, mesh, - have_nonconservative_terms(equations), equations, - dg.mortar, dg.surface_integral, dg, cache) - end - - # Calculate surface integrals, apply Jacobian from mapping to reference element - # and calculate source terms - @trixi_timeit_ext backend timer() "surface, Jacobian + source terms" begin - calc_surface_integral_and_apply_jacobian_and_calc_sources!(backend, du, u, t, - source_terms, mesh, - equations, - dg.surface_integral, - dg, cache) - end - - # Finish to send MPI data - @trixi_timeit timer() "finish MPI send" finish_mpi_send!(cache.mpi_cache) - - return nothing -end - # Copy the local solution data at all MPI interfaces directly into the MPI send buffer. # The data of each interface node is stored at the same position as in the CPU version, # which copies `cache.mpi_interfaces.u[local_side, :, i, j, interface]` into the buffer. diff --git a/src/solvers/dgsem_p4est/dg_parallel_gpu.jl b/src/solvers/dgsem_p4est/dg_parallel_gpu.jl index b65c751670d..09de17c1621 100644 --- a/src/solvers/dgsem_p4est/dg_parallel_gpu.jl +++ b/src/solvers/dgsem_p4est/dg_parallel_gpu.jl @@ -5,6 +5,129 @@ @muladd begin #! format: noindent +# In contrast to the CPU version, the GPU version does not use the MPI interface +# container `cache.mpi_interfaces.u`. Instead, `prolong2mpiinterfaces!` writes the +# local data directly into the MPI send buffer and `calc_mpi_interface_flux!` reads +# the remote data directly from the MPI receive buffer, while the local data is read +# from `u`. Both kernels handle the MPI interfaces of all neighbor ranks at once and +# parallelize over all nodes of these interfaces. +function rhs_hyperbolic!(backend::Backend, + du, u, t, + mesh::Union{P4estMeshParallel{2}, T8codeMeshParallel{2}, + P4estMeshParallel{3}, T8codeMeshParallel{3}}, + equations, + boundary_conditions, source_terms::Source, + dg::DG, cache) where {Source} + # Prolong solution to MPI interfaces, i.e., fill the MPI send buffer + @trixi_timeit_ext backend timer() "prolong2mpiinterfaces" begin + prolong2mpiinterfaces!(backend, cache, u, mesh, equations, + dg.surface_integral, dg) + end + + # Prolong solution to MPI mortars + @trixi_timeit_ext backend timer() "prolong2mpimortars" begin + prolong2mpimortars!(backend, cache, u, mesh, equations, + dg.mortar, dg) + end + + # Wait for the kernels above to finish before sending their results via MPI. + # This also guarantees that the MPI interface flux kernel of the previous call, + # which reads from the MPI receive buffer, has finished. Thus, we can only start + # to receive MPI data afterwards, in contrast to the CPU version. + # MPI is not ordered with respect to the kernels queued on the GPU, which run + # asynchronously. Thus, all kernels writing to the MPI send buffer must be launched + # before this synchronization (and after `finish_mpi_send!` of the previous call), + # and all kernels reading from the MPI receive buffer must be launched after + # `finish_mpi_receive!`. Keep this in mind when adding kernels using the MPI + # buffers (e.g., for MPI mortars). Missing synchronization cannot be detected by + # tests with the CPU backend of KernelAbstractions.jl, which runs kernels + # synchronously. + @trixi_timeit timer() "synchronize" KernelAbstractions.synchronize(backend) + + # Start to receive MPI data + @trixi_timeit timer() "start MPI receive" start_mpi_receive!(cache.mpi_cache) + + # Start to send MPI data + @trixi_timeit timer() "start MPI send" begin + start_mpi_send!(backend, cache.mpi_cache, mesh, equations, dg, cache) + end + + # Calculate volume integral + @trixi_timeit_ext backend timer() "volume integral" begin + calc_volume_integral!(backend, du, u, mesh, + have_nonconservative_terms(equations), equations, + dg.volume_integral, dg, cache) + end + + # Calculate interface fluxes + @trixi_timeit_ext backend timer() "prolong2interfaces + flux" begin + prolong2interfaces_and_calc_interface_flux!(backend, + cache.elements.surface_flux_values, + u, mesh, + have_nonconservative_terms(equations), + equations, + dg.surface_integral, dg, cache) + end + + # Prolong solution to boundaries + @trixi_timeit_ext backend timer() "prolong2boundaries" begin + prolong2boundaries!(backend, cache, u, mesh, equations, dg) + end + + # Calculate boundary fluxes + @trixi_timeit_ext backend timer() "boundary flux" begin + calc_boundary_flux!(backend, cache, t, boundary_conditions, mesh, equations, + dg.surface_integral, dg) + end + + # Prolong solution to mortars + @trixi_timeit_ext backend timer() "prolong2mortars" begin + prolong2mortars!(backend, cache, u, mesh, equations, + dg.mortar, dg) + end + + # Calculate mortar fluxes + @trixi_timeit_ext backend timer() "mortar flux" begin + calc_mortar_flux!(backend, cache.elements.surface_flux_values, mesh, + have_nonconservative_terms(equations), equations, + dg.mortar, dg.surface_integral, dg, cache) + end + + # Finish to receive MPI data + @trixi_timeit timer() "finish MPI receive" begin + finish_mpi_receive!(backend, cache.mpi_cache, mesh, equations, dg, cache) + end + + # Calculate MPI interface fluxes + @trixi_timeit_ext backend timer() "MPI interface flux" begin + calc_mpi_interface_flux!(backend, cache.elements.surface_flux_values, u, mesh, + have_nonconservative_terms(equations), equations, + dg.surface_integral, dg, cache) + end + + # Calculate MPI mortar fluxes + @trixi_timeit_ext backend timer() "MPI mortar flux" begin + calc_mpi_mortar_flux!(backend, cache.elements.surface_flux_values, mesh, + have_nonconservative_terms(equations), equations, + dg.mortar, dg.surface_integral, dg, cache) + end + + # Calculate surface integrals, apply Jacobian from mapping to reference element + # and calculate source terms + @trixi_timeit_ext backend timer() "surface, Jacobian + source terms" begin + calc_surface_integral_and_apply_jacobian_and_calc_sources!(backend, du, u, t, + source_terms, mesh, + equations, + dg.surface_integral, + dg, cache) + end + + # Finish to send MPI data + @trixi_timeit timer() "finish MPI send" finish_mpi_send!(cache.mpi_cache) + + return nothing +end + # The MPI send buffer is filled in `prolong2mpiinterfaces!`, and the kernels have # finished when this is called, see `rhs_hyperbolic!`. # This function is dimension-agnostic and used in both 2D and 3D.