trixi-framework · ranocha · Feb 6, 2026 · Jan 20, 2026 · Jan 20, 2026 · Jan 20, 2026
diff --git a/NEWS.md b/NEWS.md
@@ -9,6 +9,7 @@ for human readability.
 
 #### Changed
 
+- `convergence_test` now returns the complete convergence orders and the full errors matrix. To obtain the mean convergence rates, use `Trixi.calc_mean_convergence` on the convergence orders ([#2753]).
 - The serial and parallel mesh types have been renamed from `SerialTreeMesh`, `ParallelTreeMesh`, `SerialP4estMesh`, `ParallelP4estMesh`, `SerialT8codeMesh`, and `ParallelT8codeMesh` to `TreeMeshSerial`, `TreeMeshParallel`, `P4estMeshSerial`, `P4estMeshParallel`, `T8codeMeshSerial`, and `T8codeMeshParallel`, respectively ([#2787]).
 
 ## Changes in the v0.14 lifecycle
@@ -29,8 +30,8 @@ for human readability.
 
 #### Changed
 
-- van Leer's limiters changed to snake case: `vanLeer` => `vanleer`
-- A couple `struct`s have been made completely immutable, or only a couple fields thereof. Most notably, `save_solution.condition.save_initial_solution` where `save_solution isa SavesolutionCallback` can no longer be directly changed. Instead, the `@reset` macro from [Accessors.jl](https://github.com/JuliaObjects/Accessors.jl) is used in the elixirs.
+- van Leer's limiters changed to snake case: `vanLeer` => `vanleer` ([#2744])
+- A couple `struct`s have been made completely immutable, or only a couple fields thereof ([#2640]). Most notably, `save_solution.condition.save_initial_solution` where `save_solution isa SavesolutionCallback` can no longer be directly changed. Instead, the `@reset` macro from [Accessors.jl](https://github.com/JuliaObjects/Accessors.jl) is used in the elixirs.
 
 ## Changes in the v0.13 lifecycle
 

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -216,7 +216,7 @@ instance, in the example above the first execution of `trixi_include` takes abou
 To automatically determine the experimental order of convergence (EOC) for a
 given setup, execute
 ```julia
-julia> convergence_test(default_example(), 4)
+julia> convergence_test(default_example(), 4);
 ```
 This will run a convergence test with the elixir `default_example()`,
 using four iterations with different initial refinement levels. The initial
@@ -250,7 +250,7 @@ mean      3.99
 
 An example with multiple variables looks like this:
 ```julia
-julia> convergence_test(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_source_terms.jl"), 3)
+julia> eocs, errorsmatrix = convergence_test(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_source_terms.jl"), 3);
 ```
 ```
 [...]
@@ -273,6 +273,8 @@ error     EOC       error     EOC       error     EOC       error     EOC
 mean      3.94      mean      3.94      mean      3.94      mean      3.93
 --------------------------------------------------------------------------------
 ```
+The function `convergence_test` returns the experimental orders of convergence and the full errors
+matrix. To obtain the mean convergence rates, use `Trixi.calc_mean_convergence` on the convergence orders.
 
 ### Showcase of advanced features
 The presentation [From Mesh Generation to Adaptive Simulation: A Journey in Julia](https://youtu.be/_N4ozHr-t9E),

diff --git a/src/auxiliary/special_elixirs.jl b/src/auxiliary/special_elixirs.jl
@@ -36,7 +36,7 @@ function convergence_test(mod::Module, elixir::AbstractString, iterations,
 
         include_refined(mod, elixir, initial_resolution, iter; kwargs)
 
-        # @invokelatest is required for interactive use 
+        # @invokelatest is required for interactive use
         # due to world age issues on Julia 1.12 (and newer)
         l2_error, linf_error = @invokelatest mod.analysis_callback(@invokelatest mod.sol)
 
@@ -49,11 +49,26 @@ function convergence_test(mod::Module, elixir::AbstractString, iterations,
     end
 
     # Use raw error values to compute EOC
-    # @invokelatest is required for interactive use 
+    # @invokelatest is required for interactive use
     # due to world age issues on Julia 1.12 (and newer)
     return analyze_convergence(errors, iterations, (@invokelatest mod.semi))
 end
 
+"""
+    calc_mean_convergence(eocs)
+
+Calculate the mean convergence rates from the given experimental orders of convergence `eocs`.
+The `eocs` are expected to be in the format returned by [`convergence_test`](@ref), i.e., a `Dict` where
+the keys are the error types (e.g., `:l2`, `:linf`) and the values are matrices with the EOCs for each
+variable in the columns and the iterations in the rows.
+Returns a `Dict` with the same keys as `eocs` and the mean convergence rates for all variables as values.
+"""
+function calc_mean_convergence(eocs)
+    return Dict(kind => [sum(eocs[kind][:, v]) / length(eocs[kind][:, v])
+                         for v in 1:size(eocs[kind], 2)]
+                for kind in keys(eocs))
+end
+
 # Analyze convergence for any semidiscretization
 # Note: this intermediate method is to allow dispatching on the semidiscretization
 function analyze_convergence(errors, iterations, semi::AbstractSemidiscretization)
@@ -77,8 +92,7 @@ function analyze_convergence(errors, iterations,
     eocs = Dict(kind => log.(error[2:end, :] ./ error[1:(end - 1), :]) ./ log(1 / 2)
                 for (kind, error) in errorsmatrix)
 
-    eoc_mean_values = Dict{Symbol, Any}()
-    eoc_mean_values[:variables] = variablenames
+    eoc_mean_values = calc_mean_convergence(eocs)
 
     for (kind, error) in errorsmatrix
         println(kind)
@@ -112,18 +126,15 @@ function analyze_convergence(errors, iterations,
         println("")
 
         # Print mean EOCs
-        mean_values = zeros(eltype(errors[:l2]), nvariables)
         for v in 1:nvariables
-            mean_values[v] = sum(eocs[kind][:, v]) ./ length(eocs[kind][:, v])
             @printf("%-10s", "mean")
-            @printf("%-10.2f", mean_values[v])
+            @printf("%-10.2f", eoc_mean_values[kind][v])
         end
-        eoc_mean_values[kind] = mean_values
         println("")
         println("-"^100)
     end
 
-    return eoc_mean_values
+    return eocs, errorsmatrix
 end
 
 function convergence_test(elixir::AbstractString, iterations, RealT = Float64;

diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl
@@ -777,7 +777,8 @@ function analyze_convergence(errors_coupled, iterations,
         append!(errors[i][:linf], errors_coupled[:linf][first:last])
     end
 
-    eoc_mean_values = Vector{Dict{Symbol, Any}}(undef, nsystems(semi_coupled))
+    eocs = Vector{Dict{Symbol, Any}}(undef, nsystems(semi_coupled))
+    errorsmatrix = Vector{Dict{Symbol, Matrix{Float64}}}(undef, nsystems(semi_coupled))
     for i in eachsystem(semi_coupled)
         # Use visual cues to separate output from multiple systems
         println()
@@ -790,9 +791,10 @@ function analyze_convergence(errors_coupled, iterations,
         _, equations, _, _ = mesh_equations_solver_cache(semi)
         variablenames = varnames(cons2cons, equations)
 
-        eoc_mean_values[i] = analyze_convergence(errors[i], iterations, variablenames)
+        eocs[i], errorsmatrix[i] = analyze_convergence(errors[i], iterations,
+                                                       variablenames)
     end
 
-    return eoc_mean_values
+    return eocs, errorsmatrix
 end
 end # @muladd
diff --git a/test/test_dgmulti_1d.jl b/test/test_dgmulti_1d.jl
@@ -66,9 +66,10 @@ end
 
 @trixi_testset "elixir_euler_flux_diff.jl (convergence)" begin
     using Trixi: convergence_test
-    mean_convergence = convergence_test(@__MODULE__,
-                                        joinpath(EXAMPLES_DIR,
-                                                 "elixir_euler_flux_diff.jl"), 3)
+    eocs, _ = convergence_test(@__MODULE__,
+                               joinpath(EXAMPLES_DIR,
+                                        "elixir_euler_flux_diff.jl"), 3)
+    mean_convergence = Trixi.calc_mean_convergence(eocs)
     @test isapprox(mean_convergence[:l2],
                    [4.1558759698638434, 3.977911306037128, 4.041421206468769],
                    rtol = 0.05)

diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl
@@ -262,9 +262,10 @@ end
 
 @trixi_testset "elixir_euler_weakform.jl (convergence)" begin
     using Trixi: convergence_test
-    mean_convergence = convergence_test(@__MODULE__,
-                                        joinpath(EXAMPLES_DIR,
-                                                 "elixir_euler_weakform.jl"), 2)
+    eocs, _ = convergence_test(@__MODULE__,
+                               joinpath(EXAMPLES_DIR,
+                                        "elixir_euler_weakform.jl"), 2)
+    mean_convergence = Trixi.calc_mean_convergence(eocs)
     @test isapprox(mean_convergence[:l2],
                    [
                        4.243843382379403,

diff --git a/test/test_special_elixirs.jl b/test/test_special_elixirs.jl
@@ -21,68 +21,77 @@ EXAMPLES_DIR = examples_dir()
 
 @testset "Convergence test" begin
     @timed_testset "tree_2d_dgsem" begin
-        mean_convergence = convergence_test(@__MODULE__,
-                                            joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
-                                                     "elixir_advection_extended.jl"),
-                                            3, initial_refinement_level = 2)
+        eocs, _ = convergence_test(@__MODULE__,
+                                   joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
+                                            "elixir_advection_extended.jl"),
+                                   3, initial_refinement_level = 2)
+        mean_convergence = Trixi.calc_mean_convergence(eocs)
         @test isapprox(mean_convergence[:l2], [4.0], rtol = 0.05)
     end
 
     @timed_testset "structured_2d_dgsem" begin
-        mean_convergence = convergence_test(@__MODULE__,
-                                            joinpath(EXAMPLES_DIR,
-                                                     "structured_2d_dgsem",
-                                                     "elixir_advection_extended.jl"),
-                                            3, cells_per_dimension = (5, 9))
+        eocs, _ = convergence_test(@__MODULE__,
+                                   joinpath(EXAMPLES_DIR,
+                                            "structured_2d_dgsem",
+                                            "elixir_advection_extended.jl"),
+                                   3, cells_per_dimension = (5, 9))
+        mean_convergence = Trixi.calc_mean_convergence(eocs)
         @test isapprox(mean_convergence[:l2], [4.0], rtol = 0.05)
     end
 
     @timed_testset "structured_2d_dgsem coupled" begin
-        mean_convergence = convergence_test(@__MODULE__,
-                                            joinpath(EXAMPLES_DIR,
-                                                     "structured_2d_dgsem",
-                                                     "elixir_advection_coupled.jl"),
-                                            3)
-        @test isapprox(mean_convergence[1][:l2], [4.0], rtol = 0.05)
-        @test isapprox(mean_convergence[2][:l2], [4.0], rtol = 0.05)
+        eocs, _ = convergence_test(@__MODULE__,
+                                   joinpath(EXAMPLES_DIR,
+                                            "structured_2d_dgsem",
+                                            "elixir_advection_coupled.jl"),
+                                   3)
+        for i in Trixi.eachsystem(semi)
+            mean_convergence = Trixi.calc_mean_convergence(eocs[i])
+            @test isapprox(mean_convergence[:l2], [4.0], rtol = 0.05)
+            @test isapprox(mean_convergence[:l2], [4.0], rtol = 0.05)
+        end
     end
 
     @timed_testset "p4est_2d_dgsem" begin
         # Run convergence test on unrefined mesh
         no_refine = @cfunction((p4est, which_tree, quadrant)->Cint(0), Cint,
                                (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t},
                                 Ptr{Trixi.p4est_quadrant_t}))
-        mean_convergence = convergence_test(@__MODULE__,
-                                            joinpath(EXAMPLES_DIR, "p4est_2d_dgsem",
-                                                     "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"),
-                                            2, refine_fn_c = no_refine)
+        eocs, _ = convergence_test(@__MODULE__,
+                                   joinpath(EXAMPLES_DIR, "p4est_2d_dgsem",
+                                            "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"),
+                                   2, refine_fn_c = no_refine)
+        mean_convergence = Trixi.calc_mean_convergence(eocs)
         @test isapprox(mean_convergence[:linf], [3.2, 3.2, 4.0, 3.7], rtol = 0.05)
     end
 
     @timed_testset "structured_3d_dgsem" begin
-        mean_convergence = convergence_test(@__MODULE__,
-                                            joinpath(EXAMPLES_DIR,
-                                                     "structured_3d_dgsem",
-                                                     "elixir_advection_basic.jl"),
-                                            2, cells_per_dimension = (7, 4, 5))
+        eocs, _ = convergence_test(@__MODULE__,
+                                   joinpath(EXAMPLES_DIR,
+                                            "structured_3d_dgsem",
+                                            "elixir_advection_basic.jl"),
+                                   2, cells_per_dimension = (7, 4, 5))
+        mean_convergence = Trixi.calc_mean_convergence(eocs)
         @test isapprox(mean_convergence[:l2], [4.0], rtol = 0.05)
     end
 
     @timed_testset "p4est_3d_dgsem" begin
-        mean_convergence = convergence_test(@__MODULE__,
-                                            joinpath(EXAMPLES_DIR, "p4est_3d_dgsem",
-                                                     "elixir_advection_unstructured_curved.jl"),
-                                            2, initial_refinement_level = 0)
+        eocs, _ = convergence_test(@__MODULE__,
+                                   joinpath(EXAMPLES_DIR, "p4est_3d_dgsem",
+                                            "elixir_advection_unstructured_curved.jl"),
+                                   2, initial_refinement_level = 0)
+        mean_convergence = Trixi.calc_mean_convergence(eocs)
         @test isapprox(mean_convergence[:l2], [2.7], rtol = 0.05)
     end
 
     @timed_testset "paper_self_gravitating_gas_dynamics" begin
-        mean_convergence = convergence_test(@__MODULE__,
-                                            joinpath(EXAMPLES_DIR,
-                                                     "paper_self_gravitating_gas_dynamics",
-                                                     "elixir_eulergravity_convergence.jl"),
-                                            2, tspan = (0.0, 0.25),
-                                            initial_refinement_level = 1)
+        eocs, _ = convergence_test(@__MODULE__,
+                                   joinpath(EXAMPLES_DIR,
+                                            "paper_self_gravitating_gas_dynamics",
+                                            "elixir_eulergravity_convergence.jl"),
+                                   2, tspan = (0.0, 0.25),
+                                   initial_refinement_level = 1)
+        mean_convergence = Trixi.calc_mean_convergence(eocs)
         @test isapprox(mean_convergence[:l2], 4 * ones(4), atol = 0.4)
     end
 end