From 00d1e519dac603c0d3747ee240271dcd042e8453 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Tue, 6 Feb 2024 16:46:38 +0000
Subject: [PATCH 01/58] Added capability to the glm speed calculations to cope
 with coupled semidiscretizations.

---
 src/callbacks_step/glm_speed.jl | 40 +++++++++++++++++++++------------
 1 file changed, 26 insertions(+), 14 deletions(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 036f61a522b..37b1555126a 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -6,7 +6,7 @@
 #! format: noindent
 
 """
-    GlmSpeedCallback(; glm_scale=0.5, cfl)
+    GlmSpeedCallback(; glm_scale=0.5, cfl, semi_indices=())
 
 Update the divergence cleaning wave speed `c_h` according to the time step
 computed in [`StepsizeCallback`](@ref) for the ideal GLM-MHD equations.
@@ -14,18 +14,21 @@ The `cfl` number should be set to the same value as for the time step size calcu
 `glm_scale` ensures that the GLM wave speed is lower than the fastest physical waves in the MHD
 solution and should thus be set to a value within the interval [0,1]. Note that `glm_scale = 0`
 deactivates the divergence cleaning.
+In case of a couplings semidiscretization specify for which semi_index the divergence
+cleaning should be applied to.
 """
 struct GlmSpeedCallback{RealT <: Real}
     glm_scale::RealT
     cfl::RealT
+    semi_indices::Tuple
 end
 
 function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:GlmSpeedCallback})
     @nospecialize cb # reduce precompilation time
 
     glm_speed_callback = cb.affect!
-    @unpack glm_scale, cfl = glm_speed_callback
-    print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ", cfl=", cfl, ")")
+    @unpack glm_scale, cfl, semi_infices = glm_speed_callback
+    print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ", cfl=", cfl, "semi_indices=", semi_infices, ")")
 end
 
 function Base.show(io::IO, ::MIME"text/plain",
@@ -40,15 +43,16 @@ function Base.show(io::IO, ::MIME"text/plain",
         setup = [
             "GLM wave speed scaling" => glm_speed_callback.glm_scale,
             "Expected CFL number" => glm_speed_callback.cfl,
+#            "Coupled semidiscretization indices" => glm_speed_callback.semi_indices,
         ]
         summary_box(io, "GlmSpeedCallback", setup)
     end
 end
 
-function GlmSpeedCallback(; glm_scale = 0.5, cfl)
+function GlmSpeedCallback(; glm_scale = 0.5, cfl, semi_indices = ())
     @assert 0<=glm_scale<=1 "glm_scale must be between 0 and 1"
 
-    glm_speed_callback = GlmSpeedCallback(glm_scale, cfl)
+    glm_speed_callback = GlmSpeedCallback(glm_scale, cfl, semi_indices)
 
     DiscreteCallback(glm_speed_callback, glm_speed_callback, # the first one is the condition, the second the affect!
                      save_positions = (false, false),
@@ -69,17 +73,25 @@ end
 @inline function (glm_speed_callback::GlmSpeedCallback)(integrator)
     dt = get_proposed_dt(integrator)
     semi = integrator.p
-    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
-    @unpack glm_scale, cfl = glm_speed_callback
+    
+    @unpack glm_scale, cfl, semi_indices = glm_speed_callback
 
-    # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR)
-    c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache)
-
-    # c_h is proportional to its own time step divided by the complete MHD time step
-    equations.c_h = glm_scale * c_h_deltat / dt
+    if length(semi_indices) == 0
+        semi = tuple(semi)
+    end
 
-    # avoid re-evaluating possible FSAL stages
-    u_modified!(integrator, false)
+    for semi_index in semi_indices
+        mesh, equations, solver, cache = mesh_equations_solver_cache(semi.semis[semi_index])
+
+        # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR)
+        c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache)
+    
+        # c_h is proportional to its own time step divided by the complete MHD time step
+        equations.c_h = glm_scale * c_h_deltat / dt
+    
+        # avoid re-evaluating possible FSAL stages
+        u_modified!(integrator, false)
+    end
 
     return nothing
 end

From d37ab1891ac5d8a243724d80f022b2998f037517 Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Wed, 7 Feb 2024 15:11:02 +0000
Subject: [PATCH 02/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 src/callbacks_step/glm_speed.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 37b1555126a..c8ef2327c4a 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -27,7 +27,7 @@ function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:GlmSpeedCallback})
     @nospecialize cb # reduce precompilation time
 
     glm_speed_callback = cb.affect!
-    @unpack glm_scale, cfl, semi_infices = glm_speed_callback
+    @unpack glm_scale, cfl, semi_indices = glm_speed_callback
     print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ", cfl=", cfl, "semi_indices=", semi_infices, ")")
 end
 

From 02afd0942606d9056216cfcd0d3642cb2d3ab170 Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Wed, 7 Feb 2024 15:11:16 +0000
Subject: [PATCH 03/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 src/callbacks_step/glm_speed.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index c8ef2327c4a..40c43393c91 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -28,7 +28,7 @@ function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:GlmSpeedCallback})
 
     glm_speed_callback = cb.affect!
     @unpack glm_scale, cfl, semi_indices = glm_speed_callback
-    print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ", cfl=", cfl, "semi_indices=", semi_infices, ")")
+    print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ", cfl=", cfl, "semi_indices=", semi_indices, ")")
 end
 
 function Base.show(io::IO, ::MIME"text/plain",

From d00e1c72f728742c634d2fd50037c6d4b6fed8bf Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Wed, 7 Feb 2024 16:06:16 +0000
Subject: [PATCH 04/58] Corrected glm_speed issue for non-coupling simulations.
 Just had to make sure it can handle the situation where there is only one
 semidiscretization.

---
 src/callbacks_step/glm_speed.jl | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 37b1555126a..ef27926d52e 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -77,11 +77,13 @@ end
     @unpack glm_scale, cfl, semi_indices = glm_speed_callback
 
     if length(semi_indices) == 0
-        semi = tuple(semi)
+        semis = tuple(semi)
+    else
+        semis = semi.semis
     end
 
-    for semi_index in semi_indices
-        mesh, equations, solver, cache = mesh_equations_solver_cache(semi.semis[semi_index])
+    for semi in semis
+        mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
 
         # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR)
         c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache)

From 19e271d67263a0f7ed214b5eb426c58728036069 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Thu, 8 Feb 2024 13:57:11 +0000
Subject: [PATCH 05/58] Applied autoformatter.

---
 src/callbacks_step/glm_speed.jl | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 21b9b3a7e89..9a5de905b49 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -28,7 +28,8 @@ function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:GlmSpeedCallback})
 
     glm_speed_callback = cb.affect!
     @unpack glm_scale, cfl, semi_indices = glm_speed_callback
-    print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ", cfl=", cfl, "semi_indices=", semi_indices, ")")
+    print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ", cfl=", cfl, "semi_indices=",
+          semi_indices, ")")
 end
 
 function Base.show(io::IO, ::MIME"text/plain",
@@ -43,7 +44,7 @@ function Base.show(io::IO, ::MIME"text/plain",
         setup = [
             "GLM wave speed scaling" => glm_speed_callback.glm_scale,
             "Expected CFL number" => glm_speed_callback.cfl,
-#            "Coupled semidiscretization indices" => glm_speed_callback.semi_indices,
+            #            "Coupled semidiscretization indices" => glm_speed_callback.semi_indices,
         ]
         summary_box(io, "GlmSpeedCallback", setup)
     end
@@ -73,7 +74,7 @@ end
 @inline function (glm_speed_callback::GlmSpeedCallback)(integrator)
     dt = get_proposed_dt(integrator)
     semi = integrator.p
-    
+
     @unpack glm_scale, cfl, semi_indices = glm_speed_callback
 
     if length(semi_indices) == 0
@@ -87,10 +88,10 @@ end
 
         # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR)
         c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache)
-    
+
         # c_h is proportional to its own time step divided by the complete MHD time step
         equations.c_h = glm_scale * c_h_deltat / dt
-    
+
         # avoid re-evaluating possible FSAL stages
         u_modified!(integrator, false)
     end

From 52c674fa2d4d70a7243c2625bcfe185f45b9860e Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Thu, 8 Feb 2024 14:06:33 +0000
Subject: [PATCH 06/58] Added comment to glm speed.

---
 src/callbacks_step/glm_speed.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 9a5de905b49..50f9100b9e7 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -77,6 +77,7 @@ end
 
     @unpack glm_scale, cfl, semi_indices = glm_speed_callback
 
+    # Make sure we can handle multiple semidiscretization in coupled simulations.
     if length(semi_indices) == 0
         semis = tuple(semi)
     else

From 8136c5fff49e7449e640e48bc57eb96bd02bc5d8 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Mon, 12 Feb 2024 19:43:13 +0000
Subject: [PATCH 07/58] Added temptative example elixir for coupled GlmMhd
 equations. This requires the modified callback for the GLM cleaning speed.

---
 .../structured_2d_dgsem/elixir_mhd_coupled.jl | 138 ++++++++++++++++++
 1 file changed, 138 insertions(+)
 create mode 100644 examples/structured_2d_dgsem/elixir_mhd_coupled.jl

diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
new file mode 100644
index 00000000000..84a762d9a59
--- /dev/null
+++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
@@ -0,0 +1,138 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Coupled semidiscretization of two ideal glmMhd systems using converter functions such that
+# they are also coupled across the domain boundaries to generate a periodic system.
+#
+# In this elixir, we have a square domain that is divided into a left and right half.
+# On each half of the domain, a completely independent SemidiscretizationHyperbolic is created for the
+# linear ideal glmMhd equations. The four systems are coupled in the x and y-direction.
+# For a high-level overview, see also the figure below:
+#
+# (-2,  2)                                   ( 2,  2)
+#     ┌────────────────────┬────────────────────┐
+#     │    ↑ periodic ↑    │    ↑ periodic ↑    │
+#     │                    │                    │
+#     │     =========      │     =========      │
+#     │     system #1      │     system #2      │
+#     │     =========      │     =========      │
+#     │                    │                    │
+#     │<-- coupled         │<-- coupled         │
+#     │         coupled -->│         coupled -->│
+#     │                    │                    │
+#     │    ↓ periodic ↓    │    ↓ periodic ↓    │
+#     └────────────────────┴────────────────────┘
+# (-2, -2)                                   ( 2, -2)
+
+
+equations = IdealGlmMhdEquations2D(1.4)
+
+cells_per_dimension = (32, 64)
+
+volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+function initial_condition_constant(x, t, equations::IdealGlmMhdEquations2D)
+  rho = 1.0
+  v1 = 0.0
+  v2 = 0.0
+  v3 = 0.0
+  p = rho^equations.gamma
+  B1 = 0.0
+  B2 = 0.0
+  B3 = 0.0
+  psi = 0.0
+  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations)
+end
+
+###########
+# system #1
+###########
+
+initial_condition1 = initial_condition_constant
+coordinates_min1 = (-2.0, -2.0)
+coordinates_max1 = ( 0.0,  2.0)
+mesh1 = StructuredMesh(cells_per_dimension,
+                       coordinates_min1,
+                       coordinates_max1)
+
+coupling_function1 = (x, u, equations_other, equations_own) -> u
+boundary_conditions1 = (
+                       x_neg=BoundaryConditionCoupled(2, (:end, :i_forward), Float64, coupling_function1),
+                       x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward), Float64, coupling_function1),
+                       y_neg=boundary_condition_periodic,
+                       y_pos=boundary_condition_periodic,
+                      )
+
+semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition1, solver,
+                                     boundary_conditions=boundary_conditions1)
+
+###########
+# system #2
+###########
+
+initial_condition2 = initial_condition_constant
+coordinates_min2 = ( 0.0, -2.0)
+coordinates_max2 = ( 2.0,  2.0)
+mesh2 = StructuredMesh(cells_per_dimension,
+                       coordinates_min2,
+                       coordinates_max2)
+
+coupling_function2 = (x, u, equations_other, equations_own) -> u
+boundary_conditions2 = (
+                       x_neg=BoundaryConditionCoupled(1, (:end, :i_forward), Float64, coupling_function2),
+                       x_pos=BoundaryConditionCoupled(1, (:begin, :i_forward), Float64, coupling_function2),
+                       y_neg=boundary_condition_periodic,
+                       y_pos=boundary_condition_periodic,
+                      )
+
+semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition2, solver,
+                                     boundary_conditions=boundary_conditions2)
+
+# Create a semidiscretization that bundles all the semidiscretizations.
+semi = SemidiscretizationCoupled(semi1, semi2)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.1)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+
+analysis_callback1 = AnalysisCallback(semi1, interval=100)
+analysis_callback2 = AnalysisCallback(semi2, interval=100)
+analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=50,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+cfl = 1.0
+
+stepsize_callback = StepsizeCallback(cfl=cfl)
+
+glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl, semi_indices=tuple(1))
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback,
+                        glm_speed_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=0.01, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

From 06f3a967806c284642eb132fcc7801cea792fe62 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Tue, 13 Feb 2024 14:27:43 +0000
Subject: [PATCH 08/58] Added ability to deal with both, single
 semidiscretizations and multiple one when it comes to GlmMhd.

---
 src/callbacks_step/glm_speed.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 50f9100b9e7..a3295737e93 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -77,14 +77,16 @@ end
 
     @unpack glm_scale, cfl, semi_indices = glm_speed_callback
 
-    # Make sure we can handle multiple semidiscretization in coupled simulations.
+    # Make sure we can handle multiple semidiscretizationis in coupled simulations.
     if length(semi_indices) == 0
         semis = tuple(semi)
+	semi_indices = (1)
     else
         semis = semi.semis
     end
 
-    for semi in semis
+    for semi_index in semi_indices
+	semi = semis[semi_index]
         mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
 
         # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR)

From 8dea565fc2e9859e4145819181c5d4079f8c60f5 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Thu, 15 Feb 2024 15:35:15 +0000
Subject: [PATCH 09/58] Added temporary fix for flux_nonconservative_powell
 function for IdealGlmMhdEquations2D.

---
 src/equations/ideal_glm_mhd_2d.jl | 35 +++++++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl
index 4366cd32f08..91e4477bf62 100644
--- a/src/equations/ideal_glm_mhd_2d.jl
+++ b/src/equations/ideal_glm_mhd_2d.jl
@@ -281,6 +281,41 @@ end
     return f
 end
 
+@inline function flux_nonconservative_powell(u_ll, u_rr,
+                                             normal_direction_ll::AbstractVector,
+                                             equations::IdealGlmMhdEquations2D)
+    rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr
+
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+    v3_ll = rho_v3_ll / rho_ll
+    v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
+
+    # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector
+    # at the same node location) while `B_dot_n_rr` uses the averaged normal
+    # direction. The reason for this is that `v_dot_n_ll` depends only on the left
+    # state and multiplies some gradient while `B_dot_n_rr` is used to compute
+    # the divergence of B.
+    v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2]
+    B_dot_n_rr = B1_rr * normal_direction_ll[1] +
+                 B2_rr * normal_direction_ll[2]
+
+    # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
+    # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
+    f = SVector(0,
+                B1_ll * B_dot_n_rr,
+                B2_ll * B_dot_n_rr,
+                B3_ll * B_dot_n_rr,
+                v_dot_B_ll * B_dot_n_rr + v_dot_n_ll * psi_ll * psi_rr,
+                v1_ll * B_dot_n_rr,
+                v2_ll * B_dot_n_rr,
+                v3_ll * B_dot_n_rr,
+                v_dot_n_ll * psi_rr)
+
+    return f
+end
+
 """
     flux_nonconservative_powell_local_symmetric(u_ll, u_rr,
                                                 orientation::Integer,

From 7f7a0ff3aa96230309660f7ba7ceb6558ff61854 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Thu, 15 Feb 2024 15:35:58 +0000
Subject: [PATCH 10/58] Fix in cumputing boundary fluxes in coupled domains for
 the case of conservative and non-conservative fluxes.

---
 .../semidiscretization_coupled.jl              | 18 ++++++++++++++----
 1 file changed, 14 insertions(+), 4 deletions(-)

diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl
index dc21dbe9a1e..1edaa23a574 100644
--- a/src/semidiscretization/semidiscretization_coupled.jl
+++ b/src/semidiscretization/semidiscretization_coupled.jl
@@ -435,10 +435,20 @@ function (boundary_condition::BoundaryConditionCoupled)(u_inner, orientation, di
                                 Val(nvariables(equations))))
 
     # Calculate boundary flux
-    if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
-        flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
-    else # u_boundary is "left" of boundary, u_inner is "right" of boundary
-        flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
+    if surface_flux_function isa Tuple
+        if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
+            flux = surface_flux_function[1](u_inner, u_boundary, orientation, equations) +
+                   0.5 * surface_flux_function[2](u_inner, u_boundary, orientation, equations)
+        else # u_boundary is "left" of boundary, u_inner is "right" of boundary
+            flux = surface_flux_function[1](u_boundary, u_inner, orientation, equations) +
+                   0.5 * surface_flux_function[2](u_boundary, u_inner, orientation, equations)
+        end
+    else
+        if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
+            flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
+        else # u_boundary is "left" of boundary, u_inner is "right" of boundary
+            flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
+        end
     end
 
     return flux

From 21717ba8d6eacc0a94bebfed49b3c9161224f1a3 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Thu, 15 Feb 2024 15:42:19 +0000
Subject: [PATCH 11/58] Applied autoformatter.

---
 src/callbacks_step/glm_speed.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index a3295737e93..075706ce4b8 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -80,13 +80,13 @@ end
     # Make sure we can handle multiple semidiscretizationis in coupled simulations.
     if length(semi_indices) == 0
         semis = tuple(semi)
-	semi_indices = (1)
+        semi_indices = (1)
     else
         semis = semi.semis
     end
 
     for semi_index in semi_indices
-	semi = semis[semi_index]
+        semi = semis[semi_index]
         mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
 
         # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR)

From 0a346c542b33e95ad26ee0f1eeb89bf8159d5319 Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Fri, 16 Feb 2024 18:54:37 +0000
Subject: [PATCH 12/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>
---
 src/callbacks_step/glm_speed.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 075706ce4b8..f266e96e333 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -14,8 +14,8 @@ The `cfl` number should be set to the same value as for the time step size calcu
 `glm_scale` ensures that the GLM wave speed is lower than the fastest physical waves in the MHD
 solution and should thus be set to a value within the interval [0,1]. Note that `glm_scale = 0`
 deactivates the divergence cleaning.
-In case of a couplings semidiscretization specify for which semi_index the divergence
-cleaning should be applied to.
+In case of coupling semidiscretizations, specify for which semi_index the divergence
+cleaning should be applied.
 """
 struct GlmSpeedCallback{RealT <: Real}
     glm_scale::RealT

From 199106bf805814dcac203096ac4a0afd9ed4af60 Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Fri, 16 Feb 2024 18:58:29 +0000
Subject: [PATCH 13/58] Update
 examples/structured_2d_dgsem/elixir_mhd_coupled.jl

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>
---
 examples/structured_2d_dgsem/elixir_mhd_coupled.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
index 84a762d9a59..289d4d18d44 100644
--- a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
+++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
@@ -2,12 +2,12 @@ using OrdinaryDiffEq
 using Trixi
 
 ###############################################################################
-# Coupled semidiscretization of two ideal glmMhd systems using converter functions such that
-# they are also coupled across the domain boundaries to generate a periodic system.
+# Two semidiscretizations of the ideal GLM-MHD systems using converter functions such that
+# they are coupled across the domain boundaries to generate a periodic system.
 #
 # In this elixir, we have a square domain that is divided into a left and right half.
-# On each half of the domain, a completely independent SemidiscretizationHyperbolic is created for the
-# linear ideal glmMhd equations. The four systems are coupled in the x and y-direction.
+# On each half of the domain, an independent SemidiscretizationHyperbolic is created for
+# each set of ideal GLM-MHD equations. The two systems are coupled in the x and y-direction.
 # For a high-level overview, see also the figure below:
 #
 # (-2,  2)                                   ( 2,  2)

From cf51f9089af8d494df06235990a5a4592440e1c0 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Wed, 28 Feb 2024 15:24:45 +0000
Subject: [PATCH 14/58] Applied autoformatter.

---
 .../structured_2d_dgsem/elixir_mhd_coupled.jl | 80 +++++++++----------
 .../semidiscretization_coupled.jl             | 16 ++--
 2 files changed, 50 insertions(+), 46 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
index 289d4d18d44..3ef5c19074f 100644
--- a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
+++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
@@ -25,7 +25,6 @@ using Trixi
 #     └────────────────────┴────────────────────┘
 # (-2, -2)                                   ( 2, -2)
 
-
 equations = IdealGlmMhdEquations2D(1.4)
 
 cells_per_dimension = (32, 64)
@@ -36,16 +35,16 @@ solver = DGSEM(polydeg = 3,
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 function initial_condition_constant(x, t, equations::IdealGlmMhdEquations2D)
-  rho = 1.0
-  v1 = 0.0
-  v2 = 0.0
-  v3 = 0.0
-  p = rho^equations.gamma
-  B1 = 0.0
-  B2 = 0.0
-  B3 = 0.0
-  psi = 0.0
-  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations)
+    rho = 1.0
+    v1 = 0.0
+    v2 = 0.0
+    v3 = 0.0
+    p = rho^equations.gamma
+    B1 = 0.0
+    B2 = 0.0
+    B3 = 0.0
+    psi = 0.0
+    return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations)
 end
 
 ###########
@@ -54,43 +53,43 @@ end
 
 initial_condition1 = initial_condition_constant
 coordinates_min1 = (-2.0, -2.0)
-coordinates_max1 = ( 0.0,  2.0)
+coordinates_max1 = (0.0, 2.0)
 mesh1 = StructuredMesh(cells_per_dimension,
                        coordinates_min1,
                        coordinates_max1)
 
 coupling_function1 = (x, u, equations_other, equations_own) -> u
-boundary_conditions1 = (
-                       x_neg=BoundaryConditionCoupled(2, (:end, :i_forward), Float64, coupling_function1),
-                       x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward), Float64, coupling_function1),
-                       y_neg=boundary_condition_periodic,
-                       y_pos=boundary_condition_periodic,
-                      )
+boundary_conditions1 = (x_neg = BoundaryConditionCoupled(2, (:end, :i_forward), Float64,
+                                                         coupling_function1),
+                        x_pos = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64,
+                                                         coupling_function1),
+                        y_neg = boundary_condition_periodic,
+                        y_pos = boundary_condition_periodic)
 
 semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition1, solver,
-                                     boundary_conditions=boundary_conditions1)
+                                     boundary_conditions = boundary_conditions1)
 
 ###########
 # system #2
 ###########
 
 initial_condition2 = initial_condition_constant
-coordinates_min2 = ( 0.0, -2.0)
-coordinates_max2 = ( 2.0,  2.0)
+coordinates_min2 = (0.0, -2.0)
+coordinates_max2 = (2.0, 2.0)
 mesh2 = StructuredMesh(cells_per_dimension,
                        coordinates_min2,
                        coordinates_max2)
 
 coupling_function2 = (x, u, equations_other, equations_own) -> u
-boundary_conditions2 = (
-                       x_neg=BoundaryConditionCoupled(1, (:end, :i_forward), Float64, coupling_function2),
-                       x_pos=BoundaryConditionCoupled(1, (:begin, :i_forward), Float64, coupling_function2),
-                       y_neg=boundary_condition_periodic,
-                       y_pos=boundary_condition_periodic,
-                      )
+boundary_conditions2 = (x_neg = BoundaryConditionCoupled(1, (:end, :i_forward), Float64,
+                                                         coupling_function2),
+                        x_pos = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64,
+                                                         coupling_function2),
+                        y_neg = boundary_condition_periodic,
+                        y_pos = boundary_condition_periodic)
 
 semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition2, solver,
-                                     boundary_conditions=boundary_conditions2)
+                                     boundary_conditions = boundary_conditions2)
 
 # Create a semidiscretization that bundles all the semidiscretizations.
 semi = SemidiscretizationCoupled(semi1, semi2)
@@ -105,22 +104,22 @@ summary_callback = SummaryCallback()
 
 analysis_interval = 100
 
-analysis_callback1 = AnalysisCallback(semi1, interval=100)
-analysis_callback2 = AnalysisCallback(semi2, interval=100)
+analysis_callback1 = AnalysisCallback(semi1, interval = 100)
+analysis_callback2 = AnalysisCallback(semi2, interval = 100)
 analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2)
 
-alive_callback = AliveCallback(analysis_interval=analysis_interval)
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
-save_solution = SaveSolutionCallback(interval=50,
-                                     save_initial_solution=true,
-                                     save_final_solution=true,
-                                     solution_variables=cons2prim)
+save_solution = SaveSolutionCallback(interval = 50,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
 
 cfl = 1.0
 
-stepsize_callback = StepsizeCallback(cfl=cfl)
+stepsize_callback = StepsizeCallback(cfl = cfl)
 
-glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl, semi_indices=tuple(1))
+glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl, semi_indices = tuple(1))
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
@@ -128,11 +127,10 @@ callbacks = CallbackSet(summary_callback,
                         stepsize_callback,
                         glm_speed_callback)
 
-
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=0.01, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 0.01, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
 summary_callback() # print the timer summary
diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl
index 1edaa23a574..69993338d9e 100644
--- a/src/semidiscretization/semidiscretization_coupled.jl
+++ b/src/semidiscretization/semidiscretization_coupled.jl
@@ -437,11 +437,15 @@ function (boundary_condition::BoundaryConditionCoupled)(u_inner, orientation, di
     # Calculate boundary flux
     if surface_flux_function isa Tuple
         if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
-            flux = surface_flux_function[1](u_inner, u_boundary, orientation, equations) +
-                   0.5 * surface_flux_function[2](u_inner, u_boundary, orientation, equations)
+            flux = surface_flux_function[1](u_inner, u_boundary, orientation,
+                                            equations) +
+                   0.5 *
+                   surface_flux_function[2](u_inner, u_boundary, orientation, equations)
         else # u_boundary is "left" of boundary, u_inner is "right" of boundary
-            flux = surface_flux_function[1](u_boundary, u_inner, orientation, equations) +
-                   0.5 * surface_flux_function[2](u_boundary, u_inner, orientation, equations)
+            flux = surface_flux_function[1](u_boundary, u_inner, orientation,
+                                            equations) +
+                   0.5 *
+                   surface_flux_function[2](u_boundary, u_inner, orientation, equations)
         end
     else
         if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
@@ -475,7 +479,9 @@ function allocate_coupled_boundary_condition(boundary_condition, direction, mesh
 end
 
 # In 2D
-function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2},
+function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{
+                                                                                          2
+                                                                                          },
                                              direction, mesh, equations, dg::DGSEM)
     if direction in (1, 2)
         cell_size = size(mesh, 2)

From 5ac5bde7028a7d5b3e83c3d73db25f7001aa8599 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Wed, 28 Feb 2024 16:18:26 +0000
Subject: [PATCH 15/58] Added coupled 2d GLM-MHD simulations on structured grid
 in the tests.

---
 test/test_structured_2d.jl | 55 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 55 insertions(+)

diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl
index 64a1faf05b8..2a494342c84 100644
--- a/test/test_structured_2d.jl
+++ b/test/test_structured_2d.jl
@@ -917,6 +917,61 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_mhd_coupled.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_coupled.jl"),
+                        l2=[
+                            1.0743426980507015e-7, 0.030901698521864966,
+                            0.030901698662039206, 0.04370160129981656,
+                            8.259193827852516e-8, 0.03090169908364623,
+                            0.030901699039770684, 0.04370160128147447,
+                            8.735923402748945e-9, 1.0743426996067106e-7,
+                            0.03090169852186498, 0.030901698662039206,
+                            0.04370160129981657, 8.259193829690747e-8,
+                            0.03090169908364624, 0.030901699039770726,
+                            0.04370160128147445, 8.73592340076897e-9,
+                        ],
+                        linf=[
+                            9.021023431587949e-7, 0.043701454182710486,
+                            0.043701458294527366, 0.061803146322536154,
+                            9.487023335807976e-7, 0.043701561010342616,
+                            0.04370147392153734, 0.06180318786081025,
+                            3.430673132525334e-8, 9.02102342825728e-7,
+                            0.043701454182710764, 0.043701458294525895,
+                            0.06180314632253597, 9.487023254761695e-7,
+                            0.04370156101034084, 0.04370147392153745,
+                            0.06180318786081015, 3.430672973680963e-8,
+                        ],
+                        coverage_override=(maxiters = 10^5,))
+
+    @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin
+        errors = analysis_callback(sol)
+        @test errors.l2≈[
+            1.0743426980507015e-7, 0.030901698521864966, 0.030901698662039206,
+            0.04370160129981656, 8.259193827852516e-8, 0.03090169908364623,
+            0.030901699039770684, 0.04370160128147447, 8.735923402748945e-9,
+            1.0743426996067106e-7, 0.03090169852186498, 0.030901698662039206,
+            0.04370160129981657, 8.259193829690747e-8, 0.03090169908364624,
+            0.030901699039770726, 0.04370160128147445, 8.73592340076897e-9,
+        ] rtol=1.0e-4
+        @test errors.linf≈[
+            9.021023431587949e-7, 0.043701454182710486, 0.043701458294527366,
+            0.061803146322536154, 9.487023335807976e-7, 0.043701561010342616,
+            0.04370147392153734, 0.06180318786081025, 3.430673132525334e-8,
+            9.02102342825728e-7, 0.043701454182710764, 0.043701458294525895,
+            0.06180314632253597, 9.487023254761695e-7, 0.04370156101034084,
+            0.04370147392153745, 0.06180318786081015, 3.430672973680963e-8,
+        ] rtol=1.0e-4
+        # Ensure that we do not have excessive memory allocations
+        # (e.g., from type instabilities)
+        let
+            t = sol.t[end]
+            u_ode = sol.u[end]
+            du_ode = similar(u_ode)
+            @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        end
+    end
+end
 end
 
 # Clean up afterwards: delete Trixi.jl output directory

From dc81243f25980bde15dc172e7bac8b9538125a73 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Wed, 28 Feb 2024 16:23:36 +0000
Subject: [PATCH 16/58] Added comment on adding non-conservative fluxes to the
 conservative fluxes.

---
 src/semidiscretization/semidiscretization_coupled.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl
index 69993338d9e..2da840be1dc 100644
--- a/src/semidiscretization/semidiscretization_coupled.jl
+++ b/src/semidiscretization/semidiscretization_coupled.jl
@@ -436,6 +436,7 @@ function (boundary_condition::BoundaryConditionCoupled)(u_inner, orientation, di
 
     # Calculate boundary flux
     if surface_flux_function isa Tuple
+	# In case of conservative and non-conservative fluxes add the non-conservative one with a factor of 1/2.
         if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
             flux = surface_flux_function[1](u_inner, u_boundary, orientation,
                                             equations) +

From cda53ad0e7c260a55af7f4546c66d7eb32acc6d2 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Wed, 28 Feb 2024 16:24:29 +0000
Subject: [PATCH 17/58] Changed the coupled GLM-MHD 2d example to an Alfven
 wave system.

---
 .../structured_2d_dgsem/elixir_mhd_coupled.jl | 31 ++++++-------------
 1 file changed, 10 insertions(+), 21 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
index 3ef5c19074f..2d66be5e948 100644
--- a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
+++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
@@ -7,7 +7,8 @@ using Trixi
 #
 # In this elixir, we have a square domain that is divided into a left and right half.
 # On each half of the domain, an independent SemidiscretizationHyperbolic is created for
-# each set of ideal GLM-MHD equations. The two systems are coupled in the x and y-direction.
+# each set of ideal GLM-MHD equations. The two systems are coupled in the x-direction
+# and are periodic in the y-direction.
 # For a high-level overview, see also the figure below:
 #
 # (-2,  2)                                   ( 2,  2)
@@ -25,7 +26,8 @@ using Trixi
 #     └────────────────────┴────────────────────┘
 # (-2, -2)                                   ( 2, -2)
 
-equations = IdealGlmMhdEquations2D(1.4)
+gamma = 5 / 3
+equations = IdealGlmMhdEquations2D(gamma)
 
 cells_per_dimension = (32, 64)
 
@@ -34,26 +36,13 @@ solver = DGSEM(polydeg = 3,
                surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell),
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
-function initial_condition_constant(x, t, equations::IdealGlmMhdEquations2D)
-    rho = 1.0
-    v1 = 0.0
-    v2 = 0.0
-    v3 = 0.0
-    p = rho^equations.gamma
-    B1 = 0.0
-    B2 = 0.0
-    B3 = 0.0
-    psi = 0.0
-    return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations)
-end
-
 ###########
 # system #1
 ###########
 
-initial_condition1 = initial_condition_constant
-coordinates_min1 = (-2.0, -2.0)
-coordinates_max1 = (0.0, 2.0)
+initial_condition1 = initial_condition_convergence_test
+coordinates_min1 = (-1/sin(pi/4), -1/sin(pi/4))
+coordinates_max1 = (0.0, 1/sin(pi/4))
 mesh1 = StructuredMesh(cells_per_dimension,
                        coordinates_min1,
                        coordinates_max1)
@@ -73,9 +62,9 @@ semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition1, solve
 # system #2
 ###########
 
-initial_condition2 = initial_condition_constant
-coordinates_min2 = (0.0, -2.0)
-coordinates_max2 = (2.0, 2.0)
+initial_condition2 = initial_condition_convergence_test
+coordinates_min2 = (0.0, -1/sin(pi/4))
+coordinates_max2 = (1/sin(pi/4), 1/sin(pi/4))
 mesh2 = StructuredMesh(cells_per_dimension,
                        coordinates_min2,
                        coordinates_max2)

From 73551c55e682766b0028965e995c80bc4c0d6fac Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Wed, 28 Feb 2024 16:58:11 +0000
Subject: [PATCH 18/58] Removed redefined non-conservative Powell flux.

---
 src/equations/ideal_glm_mhd_2d.jl | 35 -------------------------------
 1 file changed, 35 deletions(-)

diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl
index 91e4477bf62..4366cd32f08 100644
--- a/src/equations/ideal_glm_mhd_2d.jl
+++ b/src/equations/ideal_glm_mhd_2d.jl
@@ -281,41 +281,6 @@ end
     return f
 end
 
-@inline function flux_nonconservative_powell(u_ll, u_rr,
-                                             normal_direction_ll::AbstractVector,
-                                             equations::IdealGlmMhdEquations2D)
-    rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
-    rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr
-
-    v1_ll = rho_v1_ll / rho_ll
-    v2_ll = rho_v2_ll / rho_ll
-    v3_ll = rho_v3_ll / rho_ll
-    v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
-
-    # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector
-    # at the same node location) while `B_dot_n_rr` uses the averaged normal
-    # direction. The reason for this is that `v_dot_n_ll` depends only on the left
-    # state and multiplies some gradient while `B_dot_n_rr` is used to compute
-    # the divergence of B.
-    v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2]
-    B_dot_n_rr = B1_rr * normal_direction_ll[1] +
-                 B2_rr * normal_direction_ll[2]
-
-    # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
-    # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
-    f = SVector(0,
-                B1_ll * B_dot_n_rr,
-                B2_ll * B_dot_n_rr,
-                B3_ll * B_dot_n_rr,
-                v_dot_B_ll * B_dot_n_rr + v_dot_n_ll * psi_ll * psi_rr,
-                v1_ll * B_dot_n_rr,
-                v2_ll * B_dot_n_rr,
-                v3_ll * B_dot_n_rr,
-                v_dot_n_ll * psi_rr)
-
-    return f
-end
-
 """
     flux_nonconservative_powell_local_symmetric(u_ll, u_rr,
                                                 orientation::Integer,

From 7a584a54d9894fecfca7d9cccd96c57ac624b0de Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Wed, 28 Feb 2024 16:58:31 +0000
Subject: [PATCH 19/58] Redefined non-conservative Powell flux so it can be
 used with more arguments.

---
 .../structured_2d_dgsem/elixir_mhd_coupled.jl | 20 +++++++++++++------
 1 file changed, 14 insertions(+), 6 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
index 2d66be5e948..16c3277184e 100644
--- a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
+++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
@@ -31,9 +31,17 @@ equations = IdealGlmMhdEquations2D(gamma)
 
 cells_per_dimension = (32, 64)
 
-volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
+# Define a non-conservative Powell flux functions that uses the normal direction for the average.
+function flux_nonconservative_powell_reduced(u_ll, u_rr,
+                                             normal_direction_ll::AbstractVector,
+					     normal_direction_average::AbstractVector,
+					     equations::IdealGlmMhdEquations2D)
+    flux_nonconservative_powell(u_ll, u_rr, normal_direction_ll, normal_direction_ll,
+                                equations)
+end
+volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell_reduced)
 solver = DGSEM(polydeg = 3,
-               surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell),
+               surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell_reduced),
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 ###########
@@ -41,8 +49,8 @@ solver = DGSEM(polydeg = 3,
 ###########
 
 initial_condition1 = initial_condition_convergence_test
-coordinates_min1 = (-1/sin(pi/4), -1/sin(pi/4))
-coordinates_max1 = (0.0, 1/sin(pi/4))
+coordinates_min1 = (-1 / sin(pi / 4), -1 / sin(pi / 4))
+coordinates_max1 = (0.0, 1 / sin(pi / 4))
 mesh1 = StructuredMesh(cells_per_dimension,
                        coordinates_min1,
                        coordinates_max1)
@@ -63,8 +71,8 @@ semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition1, solve
 ###########
 
 initial_condition2 = initial_condition_convergence_test
-coordinates_min2 = (0.0, -1/sin(pi/4))
-coordinates_max2 = (1/sin(pi/4), 1/sin(pi/4))
+coordinates_min2 = (0.0, -1 / sin(pi / 4))
+coordinates_max2 = (1 / sin(pi / 4), 1 / sin(pi / 4))
 mesh2 = StructuredMesh(cells_per_dimension,
                        coordinates_min2,
                        coordinates_max2)

From 28e7e7b6344319de4636219f9e108946e2a67f74 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Wed, 28 Feb 2024 17:25:05 +0000
Subject: [PATCH 20/58] Debugged the extension of flux_nonconservative_powell.

---
 examples/structured_2d_dgsem/elixir_mhd_coupled.jl | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
index 16c3277184e..ec111c5fe64 100644
--- a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
+++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
@@ -31,17 +31,17 @@ equations = IdealGlmMhdEquations2D(gamma)
 
 cells_per_dimension = (32, 64)
 
-# Define a non-conservative Powell flux functions that uses the normal direction for the average.
-function flux_nonconservative_powell_reduced(u_ll, u_rr,
-                                             normal_direction_ll::AbstractVector,
-					     normal_direction_average::AbstractVector,
-					     equations::IdealGlmMhdEquations2D)
+# Extend the definition of the non-conservative Powell flux functions.
+import Trixi.flux_nonconservative_powell
+function flux_nonconservative_powell(u_ll, u_rr,
+                                     normal_direction_ll::AbstractVector,
+				     equations::IdealGlmMhdEquations2D)
     flux_nonconservative_powell(u_ll, u_rr, normal_direction_ll, normal_direction_ll,
                                 equations)
 end
-volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell_reduced)
+volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
 solver = DGSEM(polydeg = 3,
-               surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell_reduced),
+               surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell),
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 ###########

From 1e64eca8a25b22bd394fcd03fcfaa55fbccc88d5 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Thu, 7 Mar 2024 17:07:40 +0000
Subject: [PATCH 21/58] Added explicit periodicity to the meshes.

---
 examples/structured_2d_dgsem/elixir_mhd_coupled.jl | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
index ec111c5fe64..112d9ec9188 100644
--- a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
+++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
@@ -35,7 +35,7 @@ cells_per_dimension = (32, 64)
 import Trixi.flux_nonconservative_powell
 function flux_nonconservative_powell(u_ll, u_rr,
                                      normal_direction_ll::AbstractVector,
-				     equations::IdealGlmMhdEquations2D)
+                                     equations::IdealGlmMhdEquations2D)
     flux_nonconservative_powell(u_ll, u_rr, normal_direction_ll, normal_direction_ll,
                                 equations)
 end
@@ -53,7 +53,8 @@ coordinates_min1 = (-1 / sin(pi / 4), -1 / sin(pi / 4))
 coordinates_max1 = (0.0, 1 / sin(pi / 4))
 mesh1 = StructuredMesh(cells_per_dimension,
                        coordinates_min1,
-                       coordinates_max1)
+                       coordinates_max1,
+                       periodicity = (false, true))
 
 coupling_function1 = (x, u, equations_other, equations_own) -> u
 boundary_conditions1 = (x_neg = BoundaryConditionCoupled(2, (:end, :i_forward), Float64,
@@ -75,7 +76,8 @@ coordinates_min2 = (0.0, -1 / sin(pi / 4))
 coordinates_max2 = (1 / sin(pi / 4), 1 / sin(pi / 4))
 mesh2 = StructuredMesh(cells_per_dimension,
                        coordinates_min2,
-                       coordinates_max2)
+                       coordinates_max2,
+                       periodicity = (false, true))
 
 coupling_function2 = (x, u, equations_other, equations_own) -> u
 boundary_conditions2 = (x_neg = BoundaryConditionCoupled(1, (:end, :i_forward), Float64,

From fd0055f66d038bf168413d02e4f0680d28960325 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Mon, 11 Mar 2024 13:59:31 +0000
Subject: [PATCH 22/58] Applied autoformatter.

---
 src/semidiscretization/semidiscretization_coupled.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl
index 2da840be1dc..186bd92e954 100644
--- a/src/semidiscretization/semidiscretization_coupled.jl
+++ b/src/semidiscretization/semidiscretization_coupled.jl
@@ -436,7 +436,7 @@ function (boundary_condition::BoundaryConditionCoupled)(u_inner, orientation, di
 
     # Calculate boundary flux
     if surface_flux_function isa Tuple
-	# In case of conservative and non-conservative fluxes add the non-conservative one with a factor of 1/2.
+        # In case of conservative and non-conservative fluxes add the non-conservative one with a factor of 1/2.
         if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
             flux = surface_flux_function[1](u_inner, u_boundary, orientation,
                                             equations) +

From 0ddaaaefe80c32ae29af429d7725a3fb9bb40386 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Tue, 12 Mar 2024 11:25:25 +0000
Subject: [PATCH 23/58] Applied updated autoformatter.

---
 src/semidiscretization/semidiscretization_coupled.jl | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl
index 186bd92e954..148c83840d4 100644
--- a/src/semidiscretization/semidiscretization_coupled.jl
+++ b/src/semidiscretization/semidiscretization_coupled.jl
@@ -480,9 +480,7 @@ function allocate_coupled_boundary_condition(boundary_condition, direction, mesh
 end
 
 # In 2D
-function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{
-                                                                                          2
-                                                                                          },
+function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2},
                                              direction, mesh, equations, dg::DGSEM)
     if direction in (1, 2)
         cell_size = size(mesh, 2)

From 8d8cfca716ed973a491cbbec16edc2823acf07a0 Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Wed, 13 Mar 2024 10:21:19 +0000
Subject: [PATCH 24/58] Fixed issue with the GlmSpeedCallback printing format.

---
 src/callbacks_step/glm_speed.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index f266e96e333..3d5133cffe7 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -44,7 +44,7 @@ function Base.show(io::IO, ::MIME"text/plain",
         setup = [
             "GLM wave speed scaling" => glm_speed_callback.glm_scale,
             "Expected CFL number" => glm_speed_callback.cfl,
-            #            "Coupled semidiscretization indices" => glm_speed_callback.semi_indices,
+	    "Coupled semis indices" => glm_speed_callback.semi_indices,
         ]
         summary_box(io, "GlmSpeedCallback", setup)
     end

From 7c3248c9061b09feeab2d012683b36b44391454f Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Wed, 13 Mar 2024 10:28:38 +0000
Subject: [PATCH 25/58] Applied autoformatter.

---
 src/callbacks_step/glm_speed.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 3d5133cffe7..b7f0323ba02 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -44,7 +44,7 @@ function Base.show(io::IO, ::MIME"text/plain",
         setup = [
             "GLM wave speed scaling" => glm_speed_callback.glm_scale,
             "Expected CFL number" => glm_speed_callback.cfl,
-	    "Coupled semis indices" => glm_speed_callback.semi_indices,
+            "Coupled semis indices" => glm_speed_callback.semi_indices,
         ]
         summary_box(io, "GlmSpeedCallback", setup)
     end

From 2114a81268f42f60326e3c0bbb583fcea42dbbfc Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Thu, 14 Mar 2024 11:03:38 +0000
Subject: [PATCH 26/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>
---
 src/callbacks_step/glm_speed.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index b7f0323ba02..508db686171 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -14,7 +14,7 @@ The `cfl` number should be set to the same value as for the time step size calcu
 `glm_scale` ensures that the GLM wave speed is lower than the fastest physical waves in the MHD
 solution and should thus be set to a value within the interval [0,1]. Note that `glm_scale = 0`
 deactivates the divergence cleaning.
-In case of coupling semidiscretizations, specify for which semi_index the divergence
+In case of coupled semidiscretizations, specify for which semi_index the divergence
 cleaning should be applied.
 """
 struct GlmSpeedCallback{RealT <: Real}

From 630c91e71f4c133393ec689244851b91359c17bd Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Thu, 14 Mar 2024 11:05:51 +0000
Subject: [PATCH 27/58] Clarified what we mean by semi_index.

---
 src/callbacks_step/glm_speed.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 508db686171..437170f0640 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -14,8 +14,8 @@ The `cfl` number should be set to the same value as for the time step size calcu
 `glm_scale` ensures that the GLM wave speed is lower than the fastest physical waves in the MHD
 solution and should thus be set to a value within the interval [0,1]. Note that `glm_scale = 0`
 deactivates the divergence cleaning.
-In case of coupled semidiscretizations, specify for which semi_index the divergence
-cleaning should be applied.
+In case of coupled semidiscretizations, specify for which semi_index, i.e. index of the
+semidiscretization, the divergence cleaning should be applied.
 """
 struct GlmSpeedCallback{RealT <: Real}
     glm_scale::RealT

From 55cf04358762040975d61c6c5ba8fa03dfa87bab Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Thu, 14 Mar 2024 11:06:24 +0000
Subject: [PATCH 28/58] Update
 src/semidiscretization/semidiscretization_coupled.jl

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>
---
 src/semidiscretization/semidiscretization_coupled.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl
index 148c83840d4..e830f5bf3d6 100644
--- a/src/semidiscretization/semidiscretization_coupled.jl
+++ b/src/semidiscretization/semidiscretization_coupled.jl
@@ -438,15 +438,15 @@ function (boundary_condition::BoundaryConditionCoupled)(u_inner, orientation, di
     if surface_flux_function isa Tuple
         # In case of conservative and non-conservative fluxes add the non-conservative one with a factor of 1/2.
         if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
-            flux = surface_flux_function[1](u_inner, u_boundary, orientation,
+            flux = (surface_flux_function[1](u_inner, u_boundary, orientation,
                                             equations) +
                    0.5 *
-                   surface_flux_function[2](u_inner, u_boundary, orientation, equations)
+                   surface_flux_function[2](u_inner, u_boundary, orientation, equations))
         else # u_boundary is "left" of boundary, u_inner is "right" of boundary
-            flux = surface_flux_function[1](u_boundary, u_inner, orientation,
+            flux = (surface_flux_function[1](u_boundary, u_inner, orientation,
                                             equations) +
                    0.5 *
-                   surface_flux_function[2](u_boundary, u_inner, orientation, equations)
+                   surface_flux_function[2](u_boundary, u_inner, orientation, equations))
         end
     else
         if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary

From e8e2abab82d0ca0ceb0867b283c0681a478fa9ea Mon Sep 17 00:00:00 2001
From: iomsn <iomsn1@googlemail.com>
Date: Thu, 14 Mar 2024 12:46:06 +0000
Subject: [PATCH 29/58] Applied autoformatter.

---
 .../semidiscretization_coupled.jl                  | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl
index e830f5bf3d6..4653d103731 100644
--- a/src/semidiscretization/semidiscretization_coupled.jl
+++ b/src/semidiscretization/semidiscretization_coupled.jl
@@ -439,14 +439,16 @@ function (boundary_condition::BoundaryConditionCoupled)(u_inner, orientation, di
         # In case of conservative and non-conservative fluxes add the non-conservative one with a factor of 1/2.
         if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
             flux = (surface_flux_function[1](u_inner, u_boundary, orientation,
-                                            equations) +
-                   0.5 *
-                   surface_flux_function[2](u_inner, u_boundary, orientation, equations))
+                                             equations) +
+                    0.5 *
+                    surface_flux_function[2](u_inner, u_boundary, orientation,
+                                             equations))
         else # u_boundary is "left" of boundary, u_inner is "right" of boundary
             flux = (surface_flux_function[1](u_boundary, u_inner, orientation,
-                                            equations) +
-                   0.5 *
-                   surface_flux_function[2](u_boundary, u_inner, orientation, equations))
+                                             equations) +
+                    0.5 *
+                    surface_flux_function[2](u_boundary, u_inner, orientation,
+                                             equations))
         end
     else
         if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary

From 56792308160ba72395d3d4779efe297bdcccf1d0 Mon Sep 17 00:00:00 2001
From: SimonCan <simon.candelaresi@gmail.com>
Date: Tue, 2 Apr 2024 12:43:19 +0100
Subject: [PATCH 30/58] Added how to to GLM speed callbacks in coupled systems.

---
 docs/src/multi-physics_coupling.md | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/docs/src/multi-physics_coupling.md b/docs/src/multi-physics_coupling.md
index eec92bc21de..1b2ecb57ed5 100644
--- a/docs/src/multi-physics_coupling.md
+++ b/docs/src/multi-physics_coupling.md
@@ -36,6 +36,17 @@ By passing the equations we can make use of their parameters, if they are requir
 Examples can be seen in `examples/structured_2d_dgsem/elixir_advection_coupled.jl`.
 
 
+## GlmSpeedCallback for coupled MHD simulations
+
+When simulating an MHD system and the [`GlmSpeedCallback`](@ref) is required
+we need to specify for which semidiscretization we need the GLM speed updated.
+This can be done with an additional parameter called `semi_indices` which
+is a tuple containing the semidiscretization indices for all systems
+that require the GLM speed updated.
+
+An example elixir can be found at [`examples/structured_2d_dgsem/elixir_mhd_coupled.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/structured_2d_dgsem/elixir_mhd_coupled.jl).
+
+
 ## Warning about binary compatibility
 Currently the coordinate values on the nodes can differ by machine precision when
 simulating the mesh and when splitting the mesh in multiple domains.

From b186b8d9028dfef12bd70117c4aafe4c0fb54df2 Mon Sep 17 00:00:00 2001
From: SimonCan <simon.candelaresi@gmail.com>
Date: Tue, 2 Apr 2024 12:47:21 +0100
Subject: [PATCH 31/58] Added news item about the GLM speed calculation in
 coupled systems.

---
 NEWS.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/NEWS.md b/NEWS.md
index 5516fc9add5..075354102a3 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -11,6 +11,7 @@ for human readability.
 - Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension
   to 1D and 3D on `TreeMesh`.
 - New analysis callback for 2D `P4estMesh` to compute integrated quantities along a boundary surface, e.g., pressure lift and drag coefficients.
+- Optional tuple parameter for `GlmSpeedCallback` called `semi_indices` to specify for which semidiscretization we need to update the GLM speed.
 
 ## Changes when updating to v0.7 from v0.6.x
 

From 7086f764d8abb90e741da9d0c95eff8634ccbfbd Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Mon, 15 Apr 2024 12:48:58 +0100
Subject: [PATCH 32/58] Update NEWS.md

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 NEWS.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/NEWS.md b/NEWS.md
index 9d67564eb60..88afec79987 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -11,7 +11,7 @@ for human readability.
   to 1D and 3D on `TreeMesh`.
 - Implementation of 1D Linearized Euler Equations.
 - New analysis callback for 2D `P4estMesh` to compute integrated quantities along a boundary surface, e.g., pressure lift and drag coefficients.
-- Optional tuple parameter for `GlmSpeedCallback` called `semi_indices` to specify for which semidiscretization we need to update the GLM speed.
+- Optional tuple parameter for `GlmSpeedCallback` called `semi_indices` to specify for which semidiscretization of a `SemidiscretizationCoupled` we need to update the GLM speed.
 
 ## Changes when updating to v0.7 from v0.6.x
 

From 15634bc0377dc2a8e29af52c6dca188205f4fbb0 Mon Sep 17 00:00:00 2001
From: SimonCan <simon.candelaresi@gmail.com>
Date: Mon, 15 Apr 2024 17:11:41 +0100
Subject: [PATCH 33/58] Throw a descriptive error message if the equations of
 subtype AbstractIdealGlmMhdEquations are not included in the list of
 semidiscretizations of AbstractIdealGlmMhdEquations.

---
 src/callbacks_step/glm_speed.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 437170f0640..d9604a7bd68 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -99,6 +99,13 @@ end
         u_modified!(integrator, false)
     end
 
+    # Check that all MHD semidiscretizations received a GLM clening speed update.
+    for (semi_index, semi) in enumerate(semis)
+        if (typeof(semi.equations) <: AbstractIdealGlmMhdEquations && !(semi_index in semi_indices))
+            throw("Equation of semidiscretization $semi_index needs to be included in 'semi_indices' of 'GlmSpeedCallback'.")
+        end
+    end
+
     return nothing
 end
 

From 28a3c06df714edc7d83be2e33e2ff0961a834a14 Mon Sep 17 00:00:00 2001
From: SimonCan <simon.candelaresi@gmail.com>
Date: Tue, 16 Apr 2024 14:23:01 +0100
Subject: [PATCH 34/58] Added checking for semi_indices in case of
 SemidiscretizationCoupled.

---
 src/callbacks_step/glm_speed.jl | 19 ++++++++++++-------
 1 file changed, 12 insertions(+), 7 deletions(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index d9604a7bd68..d8c26ec9387 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -77,6 +77,10 @@ end
 
     @unpack glm_scale, cfl, semi_indices = glm_speed_callback
 
+    if (typeof(semi) <: SemidiscretizationCoupled) && (length(semi_indices) == 0)
+        throw("Since you have more than one semidiscretization you need to specify the 'semi_indices' for which the GLM speed needs to be calculated.")
+    end
+
     # Make sure we can handle multiple semidiscretizationis in coupled simulations.
     if length(semi_indices) == 0
         semis = tuple(semi)
@@ -85,6 +89,14 @@ end
         semis = semi.semis
     end
 
+    # Check that all MHD semidiscretizations received a GLM clening speed update.
+    for (semi_index, semi) in enumerate(semis)
+        if (typeof(semi.equations) <: AbstractIdealGlmMhdEquations &&
+            !(semi_index in semi_indices))
+            throw("Equation of semidiscretization $semi_index needs to be included in 'semi_indices' of 'GlmSpeedCallback'.")
+        end
+    end
+
     for semi_index in semi_indices
         semi = semis[semi_index]
         mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
@@ -99,13 +111,6 @@ end
         u_modified!(integrator, false)
     end
 
-    # Check that all MHD semidiscretizations received a GLM clening speed update.
-    for (semi_index, semi) in enumerate(semis)
-        if (typeof(semi.equations) <: AbstractIdealGlmMhdEquations && !(semi_index in semi_indices))
-            throw("Equation of semidiscretization $semi_index needs to be included in 'semi_indices' of 'GlmSpeedCallback'.")
-        end
-    end
-
     return nothing
 end
 

From 07e6f38684bc28bb43a49847c3684e965c4f55ee Mon Sep 17 00:00:00 2001
From: SimonCan <simon.candelaresi@gmail.com>
Date: Tue, 16 Apr 2024 14:23:54 +0100
Subject: [PATCH 35/58] Added index of semidiscretization for GLM speed
 calculation.

---
 examples/structured_2d_dgsem/elixir_mhd_coupled.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
index 112d9ec9188..bf405c98694 100644
--- a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
+++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
@@ -118,7 +118,7 @@ cfl = 1.0
 
 stepsize_callback = StepsizeCallback(cfl = cfl)
 
-glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl, semi_indices = tuple(1))
+glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl, semi_indices = tuple(1, 2))
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,

From f3b18c2ff94a92dbb913cc18c9c6d54ffde78c79 Mon Sep 17 00:00:00 2001
From: SimonCan <simon.candelaresi@gmail.com>
Date: Wed, 17 Apr 2024 09:15:15 +0100
Subject: [PATCH 36/58] Aplplied autoformatter.

---
 examples/structured_2d_dgsem/elixir_mhd_coupled.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
index bf405c98694..ae16b2c38c2 100644
--- a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
+++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
@@ -118,7 +118,8 @@ cfl = 1.0
 
 stepsize_callback = StepsizeCallback(cfl = cfl)
 
-glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl, semi_indices = tuple(1, 2))
+glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl,
+                                      semi_indices = tuple(1, 2))
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,

From c952f084fec7c285c1db458b6d13588d0b5ab74d Mon Sep 17 00:00:00 2001
From: SimonCan <simon.candelaresi@gmail.com>
Date: Wed, 17 Apr 2024 09:15:31 +0100
Subject: [PATCH 37/58] Applied autoformatter.

---
 src/callbacks_step/glm_speed.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index d8c26ec9387..82972021b06 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -89,7 +89,7 @@ end
         semis = semi.semis
     end
 
-    # Check that all MHD semidiscretizations received a GLM clening speed update.
+    # Check that all MHD semidiscretizations received a GLM cleaning speed update.
     for (semi_index, semi) in enumerate(semis)
         if (typeof(semi.equations) <: AbstractIdealGlmMhdEquations &&
             !(semi_index in semi_indices))

From aa79c2882b96299bbe083fc5df14f78c5323a3f6 Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Wed, 24 Apr 2024 16:07:25 +0100
Subject: [PATCH 38/58] Update docs/src/multi-physics_coupling.md

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 docs/src/multi-physics_coupling.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/multi-physics_coupling.md b/docs/src/multi-physics_coupling.md
index 1b2ecb57ed5..10de8e04303 100644
--- a/docs/src/multi-physics_coupling.md
+++ b/docs/src/multi-physics_coupling.md
@@ -38,7 +38,7 @@ Examples can be seen in `examples/structured_2d_dgsem/elixir_advection_coupled.j
 
 ## GlmSpeedCallback for coupled MHD simulations
 
-When simulating an MHD system and the [`GlmSpeedCallback`](@ref) is required
+When simulating an MHD system and the [`GlmSpeedCallback`](@ref) is required,
 we need to specify for which semidiscretization we need the GLM speed updated.
 This can be done with an additional parameter called `semi_indices` which
 is a tuple containing the semidiscretization indices for all systems

From 4acb302bcf53e57955650c7ed11b7d99b6dbddcd Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Wed, 24 Apr 2024 16:07:38 +0100
Subject: [PATCH 39/58] Update docs/src/multi-physics_coupling.md

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 docs/src/multi-physics_coupling.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/multi-physics_coupling.md b/docs/src/multi-physics_coupling.md
index 10de8e04303..2c82c1b19e5 100644
--- a/docs/src/multi-physics_coupling.md
+++ b/docs/src/multi-physics_coupling.md
@@ -40,7 +40,7 @@ Examples can be seen in `examples/structured_2d_dgsem/elixir_advection_coupled.j
 
 When simulating an MHD system and the [`GlmSpeedCallback`](@ref) is required,
 we need to specify for which semidiscretization we need the GLM speed updated.
-This can be done with an additional parameter called `semi_indices` which
+This can be done with an additional parameter called `semi_indices`, which
 is a tuple containing the semidiscretization indices for all systems
 that require the GLM speed updated.
 

From c3c793698c0a5465ef9b4d2ec9d6ff215900c3e0 Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Wed, 24 Apr 2024 16:07:57 +0100
Subject: [PATCH 40/58] Update
 src/semidiscretization/semidiscretization_coupled.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/semidiscretization/semidiscretization_coupled.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl
index 4653d103731..f6c4ba237e1 100644
--- a/src/semidiscretization/semidiscretization_coupled.jl
+++ b/src/semidiscretization/semidiscretization_coupled.jl
@@ -436,7 +436,8 @@ function (boundary_condition::BoundaryConditionCoupled)(u_inner, orientation, di
 
     # Calculate boundary flux
     if surface_flux_function isa Tuple
-        # In case of conservative and non-conservative fluxes add the non-conservative one with a factor of 1/2.
+        # In case of conservative (index 1) and non-conservative (index 2) fluxes,
+        # add the non-conservative one with a factor of 1/2.
         if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
             flux = (surface_flux_function[1](u_inner, u_boundary, orientation,
                                              equations) +

From 4715e5240770e9ee54c512ccc87dd1904b389be8 Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Wed, 24 Apr 2024 16:08:18 +0100
Subject: [PATCH 41/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/callbacks_step/glm_speed.jl | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 82972021b06..2006a854957 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -14,8 +14,12 @@ The `cfl` number should be set to the same value as for the time step size calcu
 `glm_scale` ensures that the GLM wave speed is lower than the fastest physical waves in the MHD
 solution and should thus be set to a value within the interval [0,1]. Note that `glm_scale = 0`
 deactivates the divergence cleaning.
-In case of coupled semidiscretizations, specify for which semi_index, i.e. index of the
-semidiscretization, the divergence cleaning should be applied.
+
+In case of coupled semidiscretizations, specify for which `semi_index`, i.e. index of the
+semidiscretization, the divergence cleaning should be applied. See also
+[`SemidiscretizationCoupled`](@ref).
+Note: `SemidiscretizationCoupled` and all related features are considered experimental and
+may change at any time.
 """
 struct GlmSpeedCallback{RealT <: Real}
     glm_scale::RealT

From 6cd1baf1b458cd8f0d456c6bf5623e1593dd703a Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Wed, 24 Apr 2024 16:09:10 +0100
Subject: [PATCH 42/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/callbacks_step/glm_speed.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 2006a854957..3ea857f0821 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -48,7 +48,7 @@ function Base.show(io::IO, ::MIME"text/plain",
         setup = [
             "GLM wave speed scaling" => glm_speed_callback.glm_scale,
             "Expected CFL number" => glm_speed_callback.cfl,
-            "Coupled semis indices" => glm_speed_callback.semi_indices,
+            "Selected semidiscretizations" => glm_speed_callback.semi_indices,
         ]
         summary_box(io, "GlmSpeedCallback", setup)
     end

From 8f6c8840a9dd52cf13bccd81863b5b681fa18945 Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Wed, 24 Apr 2024 16:10:34 +0100
Subject: [PATCH 43/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/callbacks_step/glm_speed.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 3ea857f0821..f849f8de9ae 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -110,11 +110,11 @@ end
 
         # c_h is proportional to its own time step divided by the complete MHD time step
         equations.c_h = glm_scale * c_h_deltat / dt
-
-        # avoid re-evaluating possible FSAL stages
-        u_modified!(integrator, false)
     end
 
+    # avoid re-evaluating possible FSAL stages
+    u_modified!(integrator, false)
+
     return nothing
 end
 

From 9653badd7ef69b080d006770a6d4338dc73a7145 Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Wed, 24 Apr 2024 16:11:02 +0100
Subject: [PATCH 44/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/callbacks_step/glm_speed.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index f849f8de9ae..1036b898081 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -24,7 +24,7 @@ may change at any time.
 struct GlmSpeedCallback{RealT <: Real}
     glm_scale::RealT
     cfl::RealT
-    semi_indices::Tuple
+    semi_indices::Vector{Int}
 end
 
 function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:GlmSpeedCallback})

From ba3a92af3e547bb0fb7fbbd26fa08837ddecb891 Mon Sep 17 00:00:00 2001
From: SimonCan <simon.candelaresi@gmail.com>
Date: Wed, 24 Apr 2024 17:21:49 +0100
Subject: [PATCH 45/58] Updated semi_indices for glm speed calculations to
 Vector.

---
 examples/structured_2d_dgsem/elixir_mhd_coupled.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
index ae16b2c38c2..777e56d6fe2 100644
--- a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
+++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
@@ -119,7 +119,7 @@ cfl = 1.0
 stepsize_callback = StepsizeCallback(cfl = cfl)
 
 glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl,
-                                      semi_indices = tuple(1, 2))
+				      semi_indices = Vector([1, 2]))
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,

From b1e0bf50ce085d2ad3da619da1037aff92657560 Mon Sep 17 00:00:00 2001
From: SimonCan <simon.candelaresi@gmail.com>
Date: Wed, 24 Apr 2024 17:22:11 +0100
Subject: [PATCH 46/58] Some refactoring of the GLM speed calculations for
 coupled systems.

---
 src/callbacks_step/glm_speed.jl | 41 ++++++++++++++++++++-------------
 1 file changed, 25 insertions(+), 16 deletions(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 1036b898081..b634fadbdc1 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -54,7 +54,7 @@ function Base.show(io::IO, ::MIME"text/plain",
     end
 end
 
-function GlmSpeedCallback(; glm_scale = 0.5, cfl, semi_indices = ())
+function GlmSpeedCallback(; glm_scale = 0.5, cfl, semi_indices = Vector([1]))
     @assert 0<=glm_scale<=1 "glm_scale must be between 0 and 1"
 
     glm_speed_callback = GlmSpeedCallback(glm_scale, cfl, semi_indices)
@@ -74,27 +74,28 @@ function (glm_speed_callback::GlmSpeedCallback)(u, t, integrator)
     return true
 end
 
-# This method is called as callback after the StepsizeCallback during the time integration.
-@inline function (glm_speed_callback::GlmSpeedCallback)(integrator)
-    dt = get_proposed_dt(integrator)
-    semi = integrator.p
-
+function update_cleaning_speed!(semi, glm_speed_callback, dt)
     @unpack glm_scale, cfl, semi_indices = glm_speed_callback
 
-    if (typeof(semi) <: SemidiscretizationCoupled) && (length(semi_indices) == 0)
-        throw("Since you have more than one semidiscretization you need to specify the 'semi_indices' for which the GLM speed needs to be calculated.")
-    end
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+
+    # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR)
+    c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache)
+
+    # c_h is proportional to its own time step divided by the complete MHD time step
+    equations.c_h = glm_scale * c_h_deltat / dt
+end
+
+function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled, glm_speed_callback, dt)
+    @unpack glm_scale, cfl, semi_indices = glm_speed_callback
 
-    # Make sure we can handle multiple semidiscretizationis in coupled simulations.
     if length(semi_indices) == 0
-        semis = tuple(semi)
-        semi_indices = (1)
-    else
-        semis = semi.semis
+        throw("Since you have more than one semidiscretization you need to specify the 'semi_indices' for which the GLM speed needs to be calculated.")
     end
 
+#     @autoinfiltrate
     # Check that all MHD semidiscretizations received a GLM cleaning speed update.
-    for (semi_index, semi) in enumerate(semis)
+    for (semi_index, semi) in enumerate(semi_coupled.semis)
         if (typeof(semi.equations) <: AbstractIdealGlmMhdEquations &&
             !(semi_index in semi_indices))
             throw("Equation of semidiscretization $semi_index needs to be included in 'semi_indices' of 'GlmSpeedCallback'.")
@@ -102,7 +103,7 @@ end
     end
 
     for semi_index in semi_indices
-        semi = semis[semi_index]
+        semi = semi_coupled.semis[semi_index]
         mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
 
         # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR)
@@ -111,6 +112,14 @@ end
         # c_h is proportional to its own time step divided by the complete MHD time step
         equations.c_h = glm_scale * c_h_deltat / dt
     end
+end
+
+# This method is called as callback after the StepsizeCallback during the time integration.
+@inline function (glm_speed_callback::GlmSpeedCallback)(integrator)
+    dt = get_proposed_dt(integrator)
+    semi = integrator.p
+
+    update_cleaning_speed!(semi, glm_speed_callback, dt)
 
     # avoid re-evaluating possible FSAL stages
     u_modified!(integrator, false)

From 4a856720acf2e86a99f1963d3b52ccf68d6719b7 Mon Sep 17 00:00:00 2001
From: SimonCan <simon.candelaresi@gmail.com>
Date: Wed, 24 Apr 2024 18:09:20 +0100
Subject: [PATCH 47/58] Applied autoformatter.

---
 examples/structured_2d_dgsem/elixir_mhd_coupled.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
index 777e56d6fe2..dcb193bb17a 100644
--- a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
+++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
@@ -119,7 +119,7 @@ cfl = 1.0
 stepsize_callback = StepsizeCallback(cfl = cfl)
 
 glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl,
-				      semi_indices = Vector([1, 2]))
+                                      semi_indices = Vector([1, 2]))
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,

From 6bd2f8710c056a946fe8aa799c02368639cc225f Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Thu, 25 Apr 2024 18:09:34 +0100
Subject: [PATCH 48/58] Update
 examples/structured_2d_dgsem/elixir_mhd_coupled.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 examples/structured_2d_dgsem/elixir_mhd_coupled.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
index dcb193bb17a..d3aa4ecf582 100644
--- a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
+++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
@@ -119,7 +119,7 @@ cfl = 1.0
 stepsize_callback = StepsizeCallback(cfl = cfl)
 
 glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl,
-                                      semi_indices = Vector([1, 2]))
+                                      semi_indices = [1, 2])
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,

From 9aa89a514ddaefde6ba81e22f70378534105288f Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Thu, 25 Apr 2024 18:11:16 +0100
Subject: [PATCH 49/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/callbacks_step/glm_speed.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index b634fadbdc1..0fc56122bc4 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -54,7 +54,7 @@ function Base.show(io::IO, ::MIME"text/plain",
     end
 end
 
-function GlmSpeedCallback(; glm_scale = 0.5, cfl, semi_indices = Vector([1]))
+function GlmSpeedCallback(; glm_scale = 0.5, cfl, semi_indices = Int[])
     @assert 0<=glm_scale<=1 "glm_scale must be between 0 and 1"
 
     glm_speed_callback = GlmSpeedCallback(glm_scale, cfl, semi_indices)

From 7e004c15eaaa4df41cf74d9e75ec246b20346cb8 Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Thu, 25 Apr 2024 18:12:05 +0100
Subject: [PATCH 50/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/callbacks_step/glm_speed.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 0fc56122bc4..85b90b0899a 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -75,7 +75,7 @@ function (glm_speed_callback::GlmSpeedCallback)(u, t, integrator)
 end
 
 function update_cleaning_speed!(semi, glm_speed_callback, dt)
-    @unpack glm_scale, cfl, semi_indices = glm_speed_callback
+    @unpack glm_scale, cfl = glm_speed_callback
 
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
 

From cfd2e4762450e79889cf15a08bf9bb50c337c2f0 Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Thu, 25 Apr 2024 18:12:33 +0100
Subject: [PATCH 51/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/callbacks_step/glm_speed.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 85b90b0899a..6ffd32b79c0 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -84,6 +84,8 @@ function update_cleaning_speed!(semi, glm_speed_callback, dt)
 
     # c_h is proportional to its own time step divided by the complete MHD time step
     equations.c_h = glm_scale * c_h_deltat / dt
+    
+    return semi
 end
 
 function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled, glm_speed_callback, dt)

From 603f9fe917c9d4a61bcb08cc93606d6373ce8673 Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Thu, 25 Apr 2024 18:13:15 +0100
Subject: [PATCH 52/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/callbacks_step/glm_speed.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 6ffd32b79c0..87a3c03c57c 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -95,7 +95,6 @@ function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled, glm_spe
         throw("Since you have more than one semidiscretization you need to specify the 'semi_indices' for which the GLM speed needs to be calculated.")
     end
 
-#     @autoinfiltrate
     # Check that all MHD semidiscretizations received a GLM cleaning speed update.
     for (semi_index, semi) in enumerate(semi_coupled.semis)
         if (typeof(semi.equations) <: AbstractIdealGlmMhdEquations &&

From ca011e3edbff004ffb7159bb6b3a1cbf691a058a Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Thu, 25 Apr 2024 18:13:30 +0100
Subject: [PATCH 53/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/callbacks_step/glm_speed.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 87a3c03c57c..486633b4b56 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -99,7 +99,7 @@ function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled, glm_spe
     for (semi_index, semi) in enumerate(semi_coupled.semis)
         if (typeof(semi.equations) <: AbstractIdealGlmMhdEquations &&
             !(semi_index in semi_indices))
-            throw("Equation of semidiscretization $semi_index needs to be included in 'semi_indices' of 'GlmSpeedCallback'.")
+            error("Equation of semidiscretization $semi_index needs to be included in 'semi_indices' of 'GlmSpeedCallback'.")
         end
     end
 

From a432fb46335a529d4021bf828069a95ee8f37190 Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Thu, 25 Apr 2024 18:14:01 +0100
Subject: [PATCH 54/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/callbacks_step/glm_speed.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 486633b4b56..7c530509fd6 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -120,6 +120,8 @@ end
     dt = get_proposed_dt(integrator)
     semi = integrator.p
 
+    # Call the appropriate update function (this indirection allows to specialize on,
+    # e.g., the semidiscretization type)
     update_cleaning_speed!(semi, glm_speed_callback, dt)
 
     # avoid re-evaluating possible FSAL stages

From 01b374f81cb9bf3b60f24a9098f3c4fc3222c71d Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Thu, 25 Apr 2024 18:14:15 +0100
Subject: [PATCH 55/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/callbacks_step/glm_speed.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 7c530509fd6..b04695ae531 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -113,6 +113,8 @@ function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled, glm_spe
         # c_h is proportional to its own time step divided by the complete MHD time step
         equations.c_h = glm_scale * c_h_deltat / dt
     end
+    
+    return semi_coupled
 end
 
 # This method is called as callback after the StepsizeCallback during the time integration.

From c6d196dca46d156b656f9625d9eec33be198d8e3 Mon Sep 17 00:00:00 2001
From: SimonCan <simon.candelaresi@gmail.com>
Date: Thu, 25 Apr 2024 18:32:37 +0100
Subject: [PATCH 56/58] Moved the update_cleaning_speed function for
 semi_coupled::SemidiscretizationCoupled to
 src/semidiscretization/semidiscretization_coupled.jl.

---
 .../semidiscretization_coupled.jl             | 29 +++++++++++++++++++
 1 file changed, 29 insertions(+)

diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl
index f6c4ba237e1..ed8e461b440 100644
--- a/src/semidiscretization/semidiscretization_coupled.jl
+++ b/src/semidiscretization/semidiscretization_coupled.jl
@@ -348,6 +348,35 @@ function calculate_dt(u_ode, t, cfl_number, semi::SemidiscretizationCoupled)
     return dt
 end
 
+function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled, glm_speed_callback, dt)
+    @unpack glm_scale, cfl, semi_indices = glm_speed_callback
+
+    if length(semi_indices) == 0
+        throw("Since you have more than one semidiscretization you need to specify the 'semi_indices' for which the GLM speed needs to be calculated.")
+    end
+
+    # Check that all MHD semidiscretizations received a GLM cleaning speed update.
+    for (semi_index, semi) in enumerate(semi_coupled.semis)
+        if (typeof(semi.equations) <: AbstractIdealGlmMhdEquations &&
+            !(semi_index in semi_indices))
+            error("Equation of semidiscretization $semi_index needs to be included in 'semi_indices' of 'GlmSpeedCallback'.")
+        end
+    end
+
+    for semi_index in semi_indices
+        semi = semi_coupled.semis[semi_index]
+        mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+
+        # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR)
+        c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache)
+
+        # c_h is proportional to its own time step divided by the complete MHD time step
+        equations.c_h = glm_scale * c_h_deltat / dt
+    end
+    
+    return semi_coupled
+end
+
 ################################################################################
 ### Equations
 ################################################################################

From ff3c0dfac44f2a7e86efa8a9bb3e3c8588328250 Mon Sep 17 00:00:00 2001
From: SimonCan <simon.candelaresi@gmail.com>
Date: Mon, 29 Apr 2024 09:15:51 +0100
Subject: [PATCH 57/58] Applied autoformatter.

---
 src/callbacks_step/glm_speed.jl                      | 7 ++++---
 src/semidiscretization/semidiscretization_coupled.jl | 5 +++--
 2 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index b04695ae531..02a9619e11f 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -84,11 +84,12 @@ function update_cleaning_speed!(semi, glm_speed_callback, dt)
 
     # c_h is proportional to its own time step divided by the complete MHD time step
     equations.c_h = glm_scale * c_h_deltat / dt
-    
+
     return semi
 end
 
-function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled, glm_speed_callback, dt)
+function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled,
+                                glm_speed_callback, dt)
     @unpack glm_scale, cfl, semi_indices = glm_speed_callback
 
     if length(semi_indices) == 0
@@ -113,7 +114,7 @@ function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled, glm_spe
         # c_h is proportional to its own time step divided by the complete MHD time step
         equations.c_h = glm_scale * c_h_deltat / dt
     end
-    
+
     return semi_coupled
 end
 
diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl
index ed8e461b440..33803752c2e 100644
--- a/src/semidiscretization/semidiscretization_coupled.jl
+++ b/src/semidiscretization/semidiscretization_coupled.jl
@@ -348,7 +348,8 @@ function calculate_dt(u_ode, t, cfl_number, semi::SemidiscretizationCoupled)
     return dt
 end
 
-function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled, glm_speed_callback, dt)
+function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled,
+                                glm_speed_callback, dt)
     @unpack glm_scale, cfl, semi_indices = glm_speed_callback
 
     if length(semi_indices) == 0
@@ -373,7 +374,7 @@ function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled, glm_spe
         # c_h is proportional to its own time step divided by the complete MHD time step
         equations.c_h = glm_scale * c_h_deltat / dt
     end
-    
+
     return semi_coupled
 end
 

From 14696002aa253be395f27d289db925b28ab836b1 Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Fri, 3 May 2024 11:32:21 +0100
Subject: [PATCH 58/58] Update src/callbacks_step/glm_speed.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/callbacks_step/glm_speed.jl | 30 ------------------------------
 1 file changed, 30 deletions(-)

diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl
index 02a9619e11f..8ee406af5f9 100644
--- a/src/callbacks_step/glm_speed.jl
+++ b/src/callbacks_step/glm_speed.jl
@@ -88,36 +88,6 @@ function update_cleaning_speed!(semi, glm_speed_callback, dt)
     return semi
 end
 
-function update_cleaning_speed!(semi_coupled::SemidiscretizationCoupled,
-                                glm_speed_callback, dt)
-    @unpack glm_scale, cfl, semi_indices = glm_speed_callback
-
-    if length(semi_indices) == 0
-        throw("Since you have more than one semidiscretization you need to specify the 'semi_indices' for which the GLM speed needs to be calculated.")
-    end
-
-    # Check that all MHD semidiscretizations received a GLM cleaning speed update.
-    for (semi_index, semi) in enumerate(semi_coupled.semis)
-        if (typeof(semi.equations) <: AbstractIdealGlmMhdEquations &&
-            !(semi_index in semi_indices))
-            error("Equation of semidiscretization $semi_index needs to be included in 'semi_indices' of 'GlmSpeedCallback'.")
-        end
-    end
-
-    for semi_index in semi_indices
-        semi = semi_coupled.semis[semi_index]
-        mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
-
-        # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR)
-        c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache)
-
-        # c_h is proportional to its own time step divided by the complete MHD time step
-        equations.c_h = glm_scale * c_h_deltat / dt
-    end
-
-    return semi_coupled
-end
-
 # This method is called as callback after the StepsizeCallback during the time integration.
 @inline function (glm_speed_callback::GlmSpeedCallback)(integrator)
     dt = get_proposed_dt(integrator)