different approach: use Ptr to Vector and native +

.github/workflows/ci.yml

@@ -116,7 +116,6 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
         env:
           PYTHON: ""
-      - run: julia --project=@. -e 'using Pkg; pkg"add MPI#vc/custom_ops"'
       - name: Run tests without coverage
         uses: julia-actions/julia-runtest@v1
         with:

src/auxiliary/mpi.jl

@@ -128,11 +128,3 @@ parallel execution of Trixi.jl.
 See the "Miscellaneous" section of the [documentation](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/).
 """
 ode_unstable_check(dt, u, semi, t) = isnan(dt)
-
-# Custom MPI operators to work around
-# https://github.com/trixi-framework/Trixi.jl/issues/1922
-function reduce_vector_plus(x, y)
-    x .+ y
-end
-MPI.@Op(reduce_vector_plus, SVector{4, Float64})
-MPI.@Op(reduce_vector_plus, SVector{5, Float64})

src/callbacks_step/analysis_dg2d_parallel.jl

@@ -162,10 +162,19 @@ function integrate_via_indices(func::Func, u,
                             normalize = normalize)
 
     # OBS! Global results are only calculated on MPI root, all other domains receive `nothing`
-    global_integral = MPI.Reduce!(Ref(local_integral), reduce_vector_plus, mpi_root(),
-                                  mpi_comm())
+    if local_integral isa Real
+        global_integral = MPI.Reduce!(Ref(local_integral), +, mpi_root(), mpi_comm())
+    else
+        global_integral = MPI.Reduce!(Base.unsafe_convert(Ptr{Float64}, Ref(local_integral)), +, mpi_root(), mpi_comm())
+    end
+
     if mpi_isroot()
-        integral = convert(typeof(local_integral), global_integral[])
+        if local_integral isa Real
+            integral = global_integral[]
+        else
+            global_wrapped = unsafe_wrap(Array, global_integral, length(local_integral))
+            integral = convert(typeof(local_integral), global_wrapped)
+        end
     else
         integral = convert(typeof(local_integral), NaN * local_integral)
     end