GLM speed calculations for coupled semidiscretizations

NEWS.md

@@ -11,6 +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 of a `SemidiscretizationCoupled` we need to update the GLM speed.
 
 ## Changes when updating to v0.7 from v0.6.x
 

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.

examples/structured_2d_dgsem/elixir_mhd_coupled.jl

@@ -0,0 +1,136 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# 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, an independent SemidiscretizationHyperbolic is created for
+# 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)
+#     ┌────────────────────┬────────────────────┐
+#     │    ↑ periodic ↑    │    ↑ periodic ↑    │
+#     │                    │                    │
+#     │     =========      │     =========      │
+#     │     system #1      │     system #2      │
+#     │     =========      │     =========      │
+#     │                    │                    │
+#     │<-- coupled         │<-- coupled         │
+#     │         coupled -->│         coupled -->│
+#     │                    │                    │
+#     │    ↓ periodic ↓    │    ↓ periodic ↓    │
+#     └────────────────────┴────────────────────┘
+# (-2, -2)                                   ( 2, -2)
+
+gamma = 5 / 3
+equations = IdealGlmMhdEquations2D(gamma)
+
+cells_per_dimension = (32, 64)
+
+# 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)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###########
+# system #1
+###########
+
+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,
+                       periodicity = (false, true))
+
+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_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,
+                       periodicity = (false, true))
+
+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 = [1, 2])
+
+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

src/callbacks_step/glm_speed.jl

@@ -6,26 +6,34 @@
 #! 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.
 The `cfl` number should be set to the same value as for the time step size calculation. The
 `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. 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
     cfl::RealT
+    semi_indices::Vector{Int}
 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_indices = glm_speed_callback
+    print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ", cfl=", cfl, "semi_indices=",
+          semi_indices, ")")
 end
 
 function Base.show(io::IO, ::MIME"text/plain",
@@ -40,15 +48,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,
+            "Selected semidiscretizations" => 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 = Int[])
     @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),
@@ -65,19 +74,29 @@ 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
-    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+function update_cleaning_speed!(semi, glm_speed_callback, dt)
     @unpack glm_scale, cfl = glm_speed_callback
 
+    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
 
+    return semi
+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
+
+    # 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
     u_modified!(integrator, false)
 

src/semidiscretization/semidiscretization_coupled.jl

@@ -348,6 +348,36 @@ 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
 ################################################################################
@@ -435,10 +465,28 @@ 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
+        # 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) +
+                    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

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