diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 54f23758c5..b8d3ebda79 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -86,8 +86,8 @@ jobs:
             arch: x64
             trixi_test: threaded_legacy
           - version: '1.10'
-            os: macos-13
-            arch: x64
+            os: macos-latest
+            arch: aarch64
             trixi_test: mpi
           - version: '1.10'
             os: macos-latest
@@ -109,6 +109,7 @@ jobs:
           arch: ${{ matrix.arch }}
       - run: julia -e 'using InteractiveUtils; versioninfo(verbose=true)'
       - uses: julia-actions/cache@v2
+      - run: julia --project=@. -e 'using Pkg; pkg"add MPI#1dc1740"'
       - uses: julia-actions/julia-buildpkg@v1
         env:
           PYTHON: ''
diff --git a/src/auxiliary/mpi.jl b/src/auxiliary/mpi.jl
index c85c23670b..66be89e233 100644
--- a/src/auxiliary/mpi.jl
+++ b/src/auxiliary/mpi.jl
@@ -128,3 +128,10 @@ 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, Any)
diff --git a/src/callbacks_step/analysis_dg2d_parallel.jl b/src/callbacks_step/analysis_dg2d_parallel.jl
index ba2989f721..783ea5a029 100644
--- a/src/callbacks_step/analysis_dg2d_parallel.jl
+++ b/src/callbacks_step/analysis_dg2d_parallel.jl
@@ -163,7 +163,8 @@ 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), +, mpi_root(), mpi_comm())
+    global_integral = MPI.Reduce!(Ref(local_integral), reduce_vector_plus, mpi_root(),
+                                  mpi_comm())
     if mpi_isroot()
         integral = convert(typeof(local_integral), global_integral[])
     else
@@ -196,7 +197,8 @@ function integrate_via_indices(func::Func, u,
         end
     end
 
-    global_integral = MPI.Reduce!(Ref(integral), +, mpi_root(), mpi_comm())
+    global_integral = MPI.Reduce!(Ref(integral), reduce_vector_plus, mpi_root(),
+                                  mpi_comm())
     total_volume = MPI.Reduce(volume, +, mpi_root(), mpi_comm())
     if mpi_isroot()
         integral = convert(typeof(integral), global_integral[])
diff --git a/src/callbacks_step/analysis_dg3d_parallel.jl b/src/callbacks_step/analysis_dg3d_parallel.jl
index 23aaa7319b..25fa3878ea 100644
--- a/src/callbacks_step/analysis_dg3d_parallel.jl
+++ b/src/callbacks_step/analysis_dg3d_parallel.jl
@@ -91,7 +91,8 @@ function integrate_via_indices(func::Func, u,
         end
     end
 
-    global_integral = MPI.Reduce!(Ref(integral), +, mpi_root(), mpi_comm())
+    global_integral = MPI.Reduce!(Ref(integral), reduce_vector_plus, mpi_root(),
+                                  mpi_comm())
     total_volume = MPI.Reduce(volume, +, mpi_root(), mpi_comm())
     if mpi_isroot()
         integral = convert(typeof(integral), global_integral[])