diff --git a/src/Trixi.jl b/src/Trixi.jl
index 6b20d9ec9b4..f05b55935e2 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -53,7 +53,7 @@ using RecipesBase: RecipesBase
 using Requires: @require
 using Static: Static, One
 @reexport using StaticArrays: SVector
-using StaticArrays: MVector, MArray, SMatrix, @SMatrix
+using StaticArrays: StaticArrays, MVector, MArray, SMatrix, @SMatrix
 using StrideArrays: PtrArray, StrideArray, StaticInt
 @reexport using StructArrays: StructArrays, StructArray
 using TimerOutputs: TimerOutputs, @notimeit, TimerOutput, print_timer, reset_timer!
@@ -203,7 +203,10 @@ export ControllerThreeLevel, ControllerThreeLevelCombined,
 
 export PositivityPreservingLimiterZhangShu
 
-export trixi_include, examples_dir, get_examples, default_example, default_example_unstructured
+export trixi_include, examples_dir, get_examples, default_example,
+       default_example_unstructured
+
+export ode_norm, ode_unstable_check
 
 export convergence_test, jacobian_fd, jacobian_ad_forward, linear_structure
 
diff --git a/src/auxiliary/mpi.jl b/src/auxiliary/mpi.jl
index 7fab358f55b..419342fbe1f 100644
--- a/src/auxiliary/mpi.jl
+++ b/src/auxiliary/mpi.jl
@@ -1,9 +1,3 @@
-# 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
-
 
 """
     init_mpi()
@@ -75,4 +69,67 @@ end
 end
 
 
-end # @muladd
+"""
+    ode_norm(u, t)
+
+Implementation of the weighted L2 norm of Hairer and Wanner used for error-based
+step size control in OrdinaryDiffEq.jl. This function is aware of MPI and uses
+global MPI communication when running in parallel.
+
+You must pass this function as a keyword argument
+`internalnorm=ode_norm`
+to OrdinaryDiffEq.jl's `solve` when using error-based step size control with MPI
+parallel execution of Trixi.jl.
+
+See [`https://diffeq.sciml.ai/stable/basics/common_solver_opts/#advanced_adaptive_stepsize_control`](https://diffeq.sciml.ai/stable/basics/common_solver_opts/#advanced_adaptive_stepsize_control)
+
+!!! warning "Experimental code"
+    This code is experimental and may be changed or removed in any future release.
+"""
+ode_norm(u::Number, t) = @fastmath abs(u)
+function ode_norm(u::AbstractArray, t)
+  local_sumabs2 = recursive_sum_abs2(u) # sum(abs2, u)
+  local_length  = recursive_length(u)   # length(u)
+  if mpi_isparallel()
+    global_sumabs2, global_length = MPI.Allreduce([local_sumabs2, local_length], +, mpi_comm())
+    return sqrt(global_sumabs2 / global_length)
+  else
+    return sqrt(local_sumabs2 / local_length)
+  end
+end
+
+# Recursive `sum(abs2, ...)` and `length(...)` are required when dealing with
+# arrays of arrays, e.g., when using `DGMulti` solvers with an array-of-structs
+# (`Array{SVector}`) or a structure-of-arrays (`StructArray`). We need to take
+# care of these situations when allowing to use `ode_norm` as default norm in
+# OrdinaryDiffEq.jl throughout all applications of Trixi.jl.
+recursive_sum_abs2(u::Number) = abs2(u)
+recursive_sum_abs2(u::AbstractArray) = sum(recursive_sum_abs2, u)
+
+recursive_length(u::Number) = length(u)
+recursive_length(u::AbstractArray{<:Number}) = length(u)
+recursive_length(u::AbstractArray{<:AbstractArray}) = sum(recursive_length, u)
+function recursive_length(u::AbstractArray{<:StaticArrays.StaticArray{S, <:Number}}) where {S}
+  prod(StaticArrays.Size(eltype(u))) * length(u)
+end
+
+
+"""
+    ode_unstable_check(dt, u, semi, t)
+
+Implementation of the basic check for instability used in OrdinaryDiffEq.jl.
+Instead of checking something like `any(isnan, u)`, this function just checks
+`isnan(dt)`. This helps when using MPI parallelization, since no additional
+global communication is required and all ranks will return the same result.
+
+You should pass this function as a keyword argument
+`unstable_check=ode_unstable_check`
+to OrdinaryDiffEq.jl's  `solve` when using error-based step size control with MPI
+parallel execution of Trixi.jl.
+
+See [`https://diffeq.sciml.ai/stable/basics/common_solver_opts/#Miscellaneous`](https://diffeq.sciml.ai/stable/basics/common_solver_opts/#Miscellaneous)
+
+!!! warning "Experimental code"
+    This code is experimental and may be changed or removed in any future release.
+"""
+ode_unstable_check(dt, u, semi, t) = isnan(dt)
diff --git a/test/test_mpi.jl b/test/test_mpi.jl
index 4a210f1ee6b..7f6cd9051a0 100644
--- a/test/test_mpi.jl
+++ b/test/test_mpi.jl
@@ -26,6 +26,22 @@ CI_ON_WINDOWS = (get(ENV, "GITHUB_ACTIONS", false) == "true") && Sys.iswindows()
   end
 end # MPI
 
+
+@trixi_testset "MPI supporting functionality" begin
+  using OrdinaryDiffEq
+
+  t = 0.5
+  let u = 1.0
+    @test ode_norm(u, t) ≈ OrdinaryDiffEq.ODE_DEFAULT_NORM(u, t)
+  end
+  let u = [1.0, -2.0]
+    @test ode_norm(u, t) ≈ OrdinaryDiffEq.ODE_DEFAULT_NORM(u, t)
+  end
+  let u = [SVector(1.0, -2.0), SVector(0.5, -0.1)]
+    @test ode_norm(u, t) ≈ OrdinaryDiffEq.ODE_DEFAULT_NORM(u, t)
+  end
+end # MPI supporting functionality
+
 # Clean up afterwards: delete Trixi output directory
 Trixi.mpi_isroot() && @test_nowarn rm(outdir, recursive=true)
 
diff --git a/test/test_mpi_p4est_2d.jl b/test/test_mpi_p4est_2d.jl
index 325767da068..9e6f1ebbeea 100644
--- a/test/test_mpi_p4est_2d.jl
+++ b/test/test_mpi_p4est_2d.jl
@@ -18,6 +18,21 @@ const EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "
       # Expected errors are exactly the same as with TreeMesh!
       l2   = [8.311947673061856e-6],
       linf = [6.627000273229378e-5])
+
+    @testset "error-based step size control" begin
+      Trixi.mpi_isroot() && println("-"^100)
+      Trixi.mpi_isroot() && println("elixir_advection_basic.jl with error-based step size control")
+
+      sol = solve(ode, RDPK3SpFSAL35(), abstol=1.0e-4, reltol=1.0e-4,
+                  save_everystep=false, callback=callbacks,
+                  internalnorm=ode_norm,
+                  unstable_check=ode_unstable_check); summary_callback()
+      errors = analysis_callback(sol)
+      if Trixi.mpi_isroot()
+        @test errors.l2 ≈ [3.3022040342579066e-5]    rtol=1.0e-4
+        @test errors.linf ≈ [0.00011787417954578494] rtol=1.0e-4
+      end
+    end
   end
 
   @trixi_testset "elixir_advection_nonconforming_flag.jl" begin
diff --git a/test/test_mpi_p4est_3d.jl b/test/test_mpi_p4est_3d.jl
index 3488e1bfbd0..f7eb1574ccf 100644
--- a/test/test_mpi_p4est_3d.jl
+++ b/test/test_mpi_p4est_3d.jl
@@ -18,6 +18,21 @@ const EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "
       # Expected errors are exactly the same as with TreeMesh!
       l2   = [0.00016263963870641478],
       linf = [0.0014537194925779984])
+
+    @testset "error-based step size control" begin
+      Trixi.mpi_isroot() && println("-"^100)
+      Trixi.mpi_isroot() && println("elixir_advection_basic.jl with error-based step size control")
+
+      sol = solve(ode, RDPK3SpFSAL35(), abstol=1.0e-4, reltol=1.0e-4,
+                  save_everystep=false, callback=callbacks,
+                  internalnorm=ode_norm,
+                  unstable_check=ode_unstable_check); summary_callback()
+      errors = analysis_callback(sol)
+      if Trixi.mpi_isroot()
+        @test errors.l2 ≈ [0.00016800412839949264]  rtol=1.0e-4
+        @test errors.linf ≈ [0.0014548839020096516] rtol=1.0e-4
+      end
+    end
   end
 
   @trixi_testset "elixir_advection_amr.jl" begin
diff --git a/test/test_mpi_tree.jl b/test/test_mpi_tree.jl
index d8f98f17c53..79ecd39287b 100644
--- a/test/test_mpi_tree.jl
+++ b/test/test_mpi_tree.jl
@@ -124,6 +124,21 @@ const EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
       l2   = [0.061751715597716854, 0.05018223615408711, 0.05018989446443463, 0.225871559730513],
       linf = [0.29347582879608825, 0.31081249232844693, 0.3107380389947736, 1.0540358049885143])
+
+    @testset "error-based step size control" begin
+      Trixi.mpi_isroot() && println("-"^100)
+      Trixi.mpi_isroot() && println("elixir_euler_ec.jl with error-based step size control")
+
+      sol = solve(ode, RDPK3SpFSAL35(), abstol=1.0e-4, reltol=1.0e-4,
+                  save_everystep=false, callback=callbacks,
+                  internalnorm=ode_norm,
+                  unstable_check=ode_unstable_check); summary_callback()
+      errors = analysis_callback(sol)
+      if Trixi.mpi_isroot()
+        @test errors.l2   ≈ [0.061653630426688116, 0.05006930431098764, 0.05007694316484242, 0.22550689872331683] rtol=1.0e-4
+        @test errors.linf ≈ [0.28516937484583693, 0.2983633696512788, 0.297812036335975, 1.027368795517512]       rtol=1.0e-4
+      end
+    end
   end
 
   @trixi_testset "elixir_euler_vortex.jl" begin