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] 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