From 95518c5670774dfccfd66db2fc0df15ec91b251f Mon Sep 17 00:00:00 2001
From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com>
Date: Fri, 16 Jun 2023 20:35:24 +0100
Subject: [PATCH 1/2] Initial support for surface coupling of two systems
 (#1452)

* Corrected bugs from coupled to main merger.

* Added polytropic equation.

* Added further polytropic equations and examples.

* Added coupling equations.

* Added coupling equation between Euler.

* Commented debugging bits, like infiltrator.

* Add missing `using`

* Fix `destats` deprecation warning

* Added coupled elixir.

* Removed commented testing code.

* Added other_list to the coupled semi discretisation elixir.

* Removed flux coupling equation.

* Removed surface coupling equation for polytropic Euler.

* Removed polytropic Euler equation.

* Removed any code related to BoundaryConditionCoupledAB.

* Removed flux coupling code.

* Removed numerical fluxes for BoundaryConditionCoupledAB.

* Removed surface fluxes for BoundaryConditionCoupledAB.

* Removed coupled elixir.

* Remove Coupled StructuredMesh from visualization test.

* Remove duplicate function definitions

* make advection elixir go further

* Removed initial_condition_peak.

* Removed src/equations/hyperbolic_diffusion_2d.jl.

* Removed 3d coupling.

* Remopved 3d capability.

* Removed 3d capability.

* Removed 3d plotting of coupled data.

* Remove extra dependencies

* Remove whitespace changes

* Fix type instability

* Some temporary fixes.

* Fix type in semidiscretization

* Removed analysis_callback for simple coupled elixir.

* Removed analysis callbacks for the coupled case.

* Removed AnalysysCallback for coupled elixir.

* Removed polytropic coupling elixir.

* Update summary output

* Update src/solvers/dgsem_structured/dg_2d.jl

* Format summary

* Fix save solution callback

* Remove unused code

* Move timeit call before dispatch on semi

* Avoid copy on calculcate_dt

* Avoid copy on save_solution_file

* Remove unnnecessary override of wrap_array

* Undo changes to analysis callback

* Remove equations_list

* Remove unused functions

* nmeshes -> nsystems

* Further cleanup

* Move BoundaryConditionCoupled to the correct location

* Visualize setup

* Formatting improvmenets

* Change 1:nsystems(semi) to eachsystem(semi)

* Remove redundant ndofs(...) function

* copy_to_coupled_boundary --> copy_to_coupled_boundary!

* Move all SemidiscretizationCoupled-specific code to semi/semi_coupled.jl

* Use uEltype for BCCoupled

* Add comment

* I --> Indices

* Add comment

* Remove Base.summary for SemiCoupled since it appears to be unused

* Add parens

* Int64 -> Int

* Add xref for Documenter.jl

* Fixup comment

* Remove unused `total_volume`

* Remove obsolete comment

* summary_semi --> print_summary_semi for clarity

* Make SemiCoupled ctor more convenient

* Fix docstring

* Add description to elixir

* Rename elixir

* Remove unused kwarg

* Fix argument order and simplify interface for IO functions

* Explicitly return nothing in functions that should do - nothing

* Update comment

* Add AnalysisCallback to coupled semidiscretization (#1505)

* Add AnalysisCallback to coupled semidiscretization

* First non-crashing version of the AnalysisCallbackCoupled

* Added comment to offending line in the analysis callback.

* Fix stupid bug

* Rename variable for easier testing

* Clean up code

* Remove type instability

* Prettify output

* Add test

* Enable `convergence_test` for SemidiscretizationCoupled

* Increased the frequency of the solution write out for a more usable animation.

* Reverted analysis intervals.

* Updated the l2 and linf errors for the elixir_advection_basic_coupled.jl test
to reflect the increased simulation time.

* Corrected bracket typo in structured_2d test.

* Renamed plural variable names to lists.

* Further renaiming plural variable names.

* Added convergence_test for elixir_advection_basic_coupled.

* Fix coverage for convergence_test

* Add test for analysis_callback(sol)

* Split timers between systems

* fully specialize on rhs and parameters when constructing an ODEProblem

* switch example back to OrdinaryDiffEq.jl

* Reverted coupled example to use Trixi.solve instead of OrdinaryDiffEq solve.
This should fix issues with LoadError: Failed to precompile OrdinaryDiffEq
in the thread_legacy test.

* Changed Julia version in project toml to 1.9
to fix OrdinaryDiffEq issues on the github test.

* Change 1:nsystems(semi) to eachsystem(semi)

* Use `get_system_u_ode`

* Move all SemidiscretizationCoupled-specific code to semi/semi_coupled.jl

* Changed file name name of elixir_advection_basic_coupled.jl to elixir_advection_coupled.jl.

* Reverted Julia version to 1.8 in Project toml file.

* Apply suggestions from code review

* -use type SciMLBase.FullSpecialize instead of instance

* Use get_system_u_ode instead of manual view

* Reorder elixir ingredients

* Make comment reflect code again

* Use solve from OrdinaryDiffEq

* Use more precise type for array

* Test EOCs for each system separately

* Allow test to run for the full duration

---------

Co-authored-by: SimonCan <simon.candelaresi@gmail.com>
Co-authored-by: Hendrik Ranocha <mail@ranocha.de>

* Remove unused `total_volume(...)`

* Make `save_solution_file` work for SemiEulerGravity again (and make it multi-system aware)

* Update src/semidiscretization/semidiscretization_euler_gravity.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

* Apply formatting

---------

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
Co-authored-by: Hendrik Ranocha <mail@ranocha.de>
Co-authored-by: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com>
Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Co-authored-by: Lucas Gemein <74359570+NichtLucas@users.noreply.github.com>
---
 .../elixir_advection_coupled.jl               | 117 ++++
 src/Trixi.jl                                  |   8 +-
 src/auxiliary/special_elixirs.jl              |  16 +-
 src/callbacks_step/analysis.jl                |  30 +-
 src/callbacks_step/save_solution.jl           |  91 +--
 src/callbacks_step/stepsize.jl                |  19 +-
 src/callbacks_step/summary.jl                 |  22 +-
 src/meshes/mesh_io.jl                         |   8 +-
 src/semidiscretization/semidiscretization.jl  |   6 +-
 .../semidiscretization_coupled.jl             | 610 ++++++++++++++++++
 .../semidiscretization_euler_gravity.jl       |  17 +-
 test/test_special_elixirs.jl                  |   7 +
 test/test_structured_2d.jl                    |  13 +
 13 files changed, 890 insertions(+), 74 deletions(-)
 create mode 100644 examples/structured_2d_dgsem/elixir_advection_coupled.jl
 create mode 100644 src/semidiscretization/semidiscretization_coupled.jl

diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl
new file mode 100644
index 00000000000..1e54e411db6
--- /dev/null
+++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl
@@ -0,0 +1,117 @@
+using OrdinaryDiffEq
+using Trixi
+
+
+###############################################################################
+# Coupled semidiscretization of two linear advection systems, which are connected periodically
+#
+# In this elixir, we have a square domain that is divided into a left half and a right half. On each
+# half of the domain, a completely independent SemidiscretizationHyperbolic is created for the
+# linear advection equations. The two systems are coupled in the x-direction and have periodic
+# boundaries in the y-direction. For a high-level overview, see also the figure below:
+#
+# (-1,  1)                                   ( 1,  1)
+#     ┌────────────────────┬────────────────────┐
+#     │    ↑ periodic ↑    │    ↑ periodic ↑    │
+#     │                    │                    │
+#     │                    │                    │
+#     │     =========      │     =========      │
+#     │     system #1      │     system #2      │
+#     │     =========      │     =========      │
+#     │                    │                    │
+#     │                    │                    │
+#     │                    │                    │
+#     │                    │                    │
+#     │         coupled -->│<-- coupled         │
+#     │                    │                    │
+#     │<-- coupled         │         coupled -->│
+#     │                    │                    │
+#     │                    │                    │
+#     │    ↓ periodic ↓    │    ↓ periodic ↓    │
+#     └────────────────────┴────────────────────┘
+# (-1, -1)                                   ( 1, -1)
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+# First mesh is the left half of a [-1,1]^2 square
+coordinates_min1 = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max1 = ( 0.0,  1.0) # maximum coordinates (max(x), max(y))
+
+# Define identical resolution as a variable such that it is easier to change from `trixi_include`
+cells_per_dimension = (8, 16)
+
+cells_per_dimension1 = cells_per_dimension
+
+mesh1 = StructuredMesh(cells_per_dimension1, coordinates_min1, coordinates_max1)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, solver,
+                                     boundary_conditions=(
+                                       # Connect left boundary with right boundary of right mesh
+                                       x_neg=BoundaryConditionCoupled(2, (:end, :i_forward), Float64),
+                                       # Connect right boundary with left boundary of right mesh
+                                       x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward),  Float64),
+                                       y_neg=boundary_condition_periodic,
+                                       y_pos=boundary_condition_periodic))
+
+
+# Second mesh is the right half of a [-1,1]^2 square
+coordinates_min2 = (0.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max2 = (1.0,  1.0) # maximum coordinates (max(x), max(y))
+
+cells_per_dimension2 = cells_per_dimension
+
+mesh2 = StructuredMesh(cells_per_dimension2, coordinates_min2, coordinates_max2)
+
+semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, solver,
+                                     boundary_conditions=(
+                                       # Connect left boundary with right boundary of left mesh
+                                       x_neg=BoundaryConditionCoupled(1, (:end, :i_forward), Float64),
+                                       # Connect right boundary with left boundary of left mesh
+                                       x_pos=BoundaryConditionCoupled(1, (:begin, :i_forward),  Float64),
+                                       y_neg=boundary_condition_periodic,
+                                       y_pos=boundary_condition_periodic))
+
+# Create a semidiscretization that bundles semi1 and semi2
+semi = SemidiscretizationCoupled(semi1, semi2)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 2.0
+ode = semidiscretize(semi, (0.0, 2.0));
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback1 = AnalysisCallback(semi1, interval=100)
+analysis_callback2 = AnalysisCallback(semi2, interval=100)
+analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     solution_variables=cons2prim)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl=1.6)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+
+# Print the timer summary
+summary_callback()
diff --git a/src/Trixi.jl b/src/Trixi.jl
index d5579aeea33..86e349c7dad 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -117,6 +117,7 @@ include("semidiscretization/semidiscretization.jl")
 include("semidiscretization/semidiscretization_hyperbolic.jl")
 include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl")
 include("semidiscretization/semidiscretization_euler_acoustics.jl")
+include("semidiscretization/semidiscretization_coupled.jl")
 include("callbacks_step/callbacks_step.jl")
 include("callbacks_stage/callbacks_stage.jl")
 include("semidiscretization/semidiscretization_euler_gravity.jl")
@@ -184,7 +185,8 @@ export boundary_condition_do_nothing,
        boundary_condition_noslip_wall,
        boundary_condition_slip_wall,
        boundary_condition_wall,
-       BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal
+       BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal,
+       BoundaryConditionCoupled
 
 export initial_condition_convergence_test, source_terms_convergence_test
 export source_terms_harmonic
@@ -229,12 +231,14 @@ export SemidiscretizationEulerAcoustics
 export SemidiscretizationEulerGravity, ParametersEulerGravity,
        timestep_gravity_erk52_3Sstar!, timestep_gravity_carpenter_kennedy_erk54_2N!
 
+export SemidiscretizationCoupled
+
 export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        SaveRestartCallback, SaveSolutionCallback, TimeSeriesCallback, VisualizationCallback,
        AveragingCallback,
        AMRCallback, StepsizeCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
-       TrivialCallback
+       TrivialCallback, AnalysisCallbackCoupled
 
 export load_mesh, load_time
 
diff --git a/src/auxiliary/special_elixirs.jl b/src/auxiliary/special_elixirs.jl
index da73b42e572..25bca8939ce 100644
--- a/src/auxiliary/special_elixirs.jl
+++ b/src/auxiliary/special_elixirs.jl
@@ -85,9 +85,21 @@ function convergence_test(mod::Module, elixir::AbstractString, iterations; kwarg
         println("#"^100)
     end
 
-    # number of variables
-    _, equations, _, _ = mesh_equations_solver_cache(mod.semi)
+    # Use raw error values to compute EOC
+    analyze_convergence(errors, iterations, mod.semi)
+end
+
+# Analyze convergence for any semidiscretization
+# Note: this intermediate method is to allow dispatching on the semidiscretization
+function analyze_convergence(errors, iterations, semi::AbstractSemidiscretization)
+    _, equations, _, _ = mesh_equations_solver_cache(semi)
     variablenames = varnames(cons2cons, equations)
+    analyze_convergence(errors, iterations, variablenames)
+end
+
+# This method is called with the collected error values to actually compute and print the EOC
+function analyze_convergence(errors, iterations,
+                             variablenames::Union{Tuple, AbstractArray})
     nvariables = length(variablenames)
 
     # Reshape errors to get a matrix where the i-th row represents the i-th iteration
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 2e038401df7..7fa2e21a244 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -84,11 +84,13 @@ function Base.show(io::IO, ::MIME"text/plain",
     end
 end
 
+# This is the convenience constructor that gets called from the elixirs
 function AnalysisCallback(semi::AbstractSemidiscretization; kwargs...)
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     AnalysisCallback(mesh, equations, solver, cache; kwargs...)
 end
 
+# This is the actual constructor
 function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache;
                           interval = 0,
                           save_analysis = false,
@@ -132,9 +134,18 @@ function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache;
                      initialize = initialize!)
 end
 
+# This method gets called from OrdinaryDiffEq's `solve(...)`
 function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, t,
                      integrator) where {Condition, Affect! <: AnalysisCallback}
     semi = integrator.p
+    du_ode = first(get_tmp_cache(integrator))
+    initialize!(cb, u_ode, du_ode, t, integrator, semi)
+end
+
+# This is the actual initialization method
+# Note: we have this indirection to allow initializing a callback from the AnalysisCallbackCoupled
+function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, du_ode, t,
+                     integrator, semi) where {Condition, Affect! <: AnalysisCallback}
     initial_state_integrals = integrate(u_ode, semi)
     _, equations, _, _ = mesh_equations_solver_cache(semi)
 
@@ -202,13 +213,21 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, t,
     # Note: For details see the actual callback function below
     analysis_callback.start_gc_time = Base.gc_time_ns()
 
-    analysis_callback(integrator)
+    analysis_callback(u_ode, du_ode, integrator, semi)
     return nothing
 end
 
-# TODO: Taal refactor, allow passing an IO object (which could be devnull to avoid cluttering the console)
+# This method gets called from OrdinaryDiffEq's `solve(...)`
 function (analysis_callback::AnalysisCallback)(integrator)
     semi = integrator.p
+    du_ode = first(get_tmp_cache(integrator))
+    u_ode = integrator.u
+    analysis_callback(u_ode, du_ode, integrator, semi)
+end
+
+# This method gets called internally as the main entry point to the AnalysiCallback
+# TODO: Taal refactor, allow passing an IO object (which could be devnull to avoid cluttering the console)
+function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi)
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     @unpack dt, t = integrator
     iter = integrator.stats.naccept
@@ -300,15 +319,14 @@ function (analysis_callback::AnalysisCallback)(integrator)
         end
 
         # Calculate current time derivative (needed for semidiscrete entropy time derivative, residual, etc.)
-        du_ode = first(get_tmp_cache(integrator))
         # `integrator.f` is usually just a call to `rhs!`
         # However, we want to allow users to modify the ODE RHS outside of Trixi.jl
         # and allow us to pass a combined ODE RHS to OrdinaryDiffEq, e.g., for
         # hyperbolic-parabolic systems.
-        @notimeit timer() integrator.f(du_ode, integrator.u, semi, t)
-        u = wrap_array(integrator.u, mesh, equations, solver, cache)
+        @notimeit timer() integrator.f(du_ode, u_ode, semi, t)
+        u = wrap_array(u_ode, mesh, equations, solver, cache)
         du = wrap_array(du_ode, mesh, equations, solver, cache)
-        l2_error, linf_error = analysis_callback(io, du, u, integrator.u, t, semi)
+        l2_error, linf_error = analysis_callback(io, du, u, u_ode, t, semi)
 
         mpi_println("─"^100)
         mpi_println()
diff --git a/src/callbacks_step/save_solution.jl b/src/callbacks_step/save_solution.jl
index 55f17bbc1c7..1fe0d6b1e15 100644
--- a/src/callbacks_step/save_solution.jl
+++ b/src/callbacks_step/save_solution.jl
@@ -141,14 +141,7 @@ function initialize_save_cb!(solution_callback::SaveSolutionCallback, u, t, inte
     mpi_isroot() && mkpath(solution_callback.output_directory)
 
     semi = integrator.p
-    mesh, _, _, _ = mesh_equations_solver_cache(semi)
-    @trixi_timeit timer() "I/O" begin
-        if mesh.unsaved_changes
-            mesh.current_filename = save_mesh_file(mesh,
-                                                   solution_callback.output_directory)
-            mesh.unsaved_changes = false
-        end
-    end
+    @trixi_timeit timer() "I/O" save_mesh(semi, solution_callback.output_directory)
 
     if solution_callback.save_initial_solution
         solution_callback(integrator)
@@ -157,6 +150,16 @@ function initialize_save_cb!(solution_callback::SaveSolutionCallback, u, t, inte
     return nothing
 end
 
+# Save mesh for a general semidiscretization (default)
+function save_mesh(semi::AbstractSemidiscretization, output_directory, timestep = 0)
+    mesh, _, _, _ = mesh_equations_solver_cache(semi)
+
+    if mesh.unsaved_changes
+        mesh.current_filename = save_mesh_file(mesh, output_directory)
+        mesh.unsaved_changes = false
+    end
+end
+
 # this method is called to determine whether the callback should be activated
 function (solution_callback::SaveSolutionCallback)(u, t, integrator)
     @unpack interval_or_dt, save_final_solution = solution_callback
@@ -174,41 +177,15 @@ end
 # this method is called when the callback is activated
 function (solution_callback::SaveSolutionCallback)(integrator)
     u_ode = integrator.u
-    @unpack t, dt = integrator
-    iter = integrator.stats.naccept
     semi = integrator.p
-    mesh, _, _, _ = mesh_equations_solver_cache(semi)
+    iter = integrator.stats.naccept
 
     @trixi_timeit timer() "I/O" begin
-        @trixi_timeit timer() "save mesh" if mesh.unsaved_changes
-            mesh.current_filename = save_mesh_file(mesh,
-                                                   solution_callback.output_directory,
-                                                   iter)
-            mesh.unsaved_changes = false
-        end
-
-        element_variables = Dict{Symbol, Any}()
-        @trixi_timeit timer() "get element variables" begin
-            get_element_variables!(element_variables, u_ode, semi)
-            callbacks = integrator.opts.callback
-            if callbacks isa CallbackSet
-                for cb in callbacks.continuous_callbacks
-                    get_element_variables!(element_variables, u_ode, semi, cb;
-                                           t = integrator.t,
-                                           iter = integrator.stats.naccept)
-                end
-                for cb in callbacks.discrete_callbacks
-                    get_element_variables!(element_variables, u_ode, semi, cb;
-                                           t = integrator.t,
-                                           iter = integrator.stats.naccept)
-                end
-            end
-        end
-
-        @trixi_timeit timer() "save solution" save_solution_file(u_ode, t, dt, iter,
-                                                                 semi,
-                                                                 solution_callback,
-                                                                 element_variables)
+        # Call high-level functions that dispatch on semidiscretization type
+        @trixi_timeit timer() "save mesh" save_mesh(semi,
+                                                    solution_callback.output_directory,
+                                                    iter)
+        save_solution_file(semi, u_ode, solution_callback, integrator)
     end
 
     # avoid re-evaluating possible FSAL stages
@@ -216,13 +193,43 @@ function (solution_callback::SaveSolutionCallback)(integrator)
     return nothing
 end
 
+@inline function save_solution_file(semi::AbstractSemidiscretization, u_ode,
+                                    solution_callback,
+                                    integrator; system = "")
+    @unpack t, dt = integrator
+    iter = integrator.stats.naccept
+
+    element_variables = Dict{Symbol, Any}()
+    @trixi_timeit timer() "get element variables" begin
+        get_element_variables!(element_variables, u_ode, semi)
+        callbacks = integrator.opts.callback
+        if callbacks isa CallbackSet
+            for cb in callbacks.continuous_callbacks
+                get_element_variables!(element_variables, u_ode, semi, cb;
+                                       t = integrator.t, iter = iter)
+            end
+            for cb in callbacks.discrete_callbacks
+                get_element_variables!(element_variables, u_ode, semi, cb;
+                                       t = integrator.t, iter = iter)
+            end
+        end
+    end
+
+    @trixi_timeit timer() "save solution" save_solution_file(u_ode, t, dt, iter, semi,
+                                                             solution_callback,
+                                                             element_variables,
+                                                             system = system)
+end
+
 @inline function save_solution_file(u_ode, t, dt, iter,
                                     semi::AbstractSemidiscretization, solution_callback,
-                                    element_variables = Dict{Symbol, Any}())
+                                    element_variables = Dict{Symbol, Any}();
+                                    system = "")
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     u = wrap_array_native(u_ode, mesh, equations, solver, cache)
     save_solution_file(u, t, dt, iter, mesh, equations, solver, cache,
-                       solution_callback, element_variables)
+                       solution_callback,
+                       element_variables; system = system)
 end
 
 # TODO: Taal refactor, move save_mesh_file?
diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl
index 9e9f2d4885b..8b5cb958318 100644
--- a/src/callbacks_step/stepsize.jl
+++ b/src/callbacks_step/stepsize.jl
@@ -64,14 +64,11 @@ end
         t = integrator.t
         u_ode = integrator.u
         semi = integrator.p
-        mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
         @unpack cfl_number = stepsize_callback
-        u = wrap_array(u_ode, mesh, equations, solver, cache)
 
-        dt = @trixi_timeit timer() "calculate dt" begin
-            cfl_number * max_dt(u, t, mesh, have_constant_speed(equations), equations,
-                   solver, cache)
-        end
+        # Dispatch based on semidiscretization
+        dt = @trixi_timeit timer() "calculate dt" calculate_dt(u_ode, t, cfl_number,
+                                                               semi)
 
         set_proposed_dt!(integrator, dt)
         integrator.opts.dtmax = dt
@@ -83,6 +80,16 @@ end
     return nothing
 end
 
+# General case for a single semidiscretization
+function calculate_dt(u_ode, t, cfl_number, semi::AbstractSemidiscretization)
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+    u = wrap_array(u_ode, mesh, equations, solver, cache)
+
+    dt = cfl_number * max_dt(u, t, mesh,
+                have_constant_speed(equations), equations,
+                solver, cache)
+end
+
 # Time integration methods from the DiffEq ecosystem without adaptive time stepping on their own
 # such as `CarpenterKennedy2N54` require passing `dt=...` in `solve(ode, ...)`. Since we don't have
 # an integrator at this stage but only the ODE, this method will be used there. It's called in
diff --git a/src/callbacks_step/summary.jl b/src/callbacks_step/summary.jl
index a73b2a1913b..08e13d0b98d 100644
--- a/src/callbacks_step/summary.jl
+++ b/src/callbacks_step/summary.jl
@@ -152,15 +152,7 @@ function initialize_summary_callback(cb::DiscreteCallback, u, t, integrator)
                            :indentation_level => 0)
 
     semi = integrator.p
-    show(io_context, MIME"text/plain"(), semi)
-    println(io, "\n")
-    mesh, equations, solver, _ = mesh_equations_solver_cache(semi)
-    show(io_context, MIME"text/plain"(), mesh)
-    println(io, "\n")
-    show(io_context, MIME"text/plain"(), equations)
-    println(io, "\n")
-    show(io_context, MIME"text/plain"(), solver)
-    println(io, "\n")
+    print_summary_semidiscretization(io_context, semi)
 
     callbacks = integrator.opts.callback
     if callbacks isa CallbackSet
@@ -208,6 +200,18 @@ function initialize_summary_callback(cb::DiscreteCallback, u, t, integrator)
     return nothing
 end
 
+function print_summary_semidiscretization(io::IO, semi::AbstractSemidiscretization)
+    show(io, MIME"text/plain"(), semi)
+    println(io, "\n")
+    mesh, equations, solver, _ = mesh_equations_solver_cache(semi)
+    show(io, MIME"text/plain"(), mesh)
+    println(io, "\n")
+    show(io, MIME"text/plain"(), equations)
+    println(io, "\n")
+    show(io, MIME"text/plain"(), solver)
+    println(io, "\n")
+end
+
 function (cb::DiscreteCallback{Condition, Affect!})(io::IO = stdout) where {Condition,
                                                                             Affect! <:
                                                                             typeof(summary_callback)
diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl
index b9c462fa15a..ede85d80106 100644
--- a/src/meshes/mesh_io.jl
+++ b/src/meshes/mesh_io.jl
@@ -95,11 +95,15 @@ end
 # of the mesh, like its size and the type of boundary mapping function.
 # Then, within Trixi2Vtk, the StructuredMesh and its node coordinates are reconstructured from
 # these attributes for plotting purposes
-function save_mesh_file(mesh::StructuredMesh, output_directory)
+function save_mesh_file(mesh::StructuredMesh, output_directory; system = "")
     # Create output directory (if it does not exist)
     mkpath(output_directory)
 
-    filename = joinpath(output_directory, "mesh.h5")
+    if isempty(system)
+        filename = joinpath(output_directory, "mesh.h5")
+    else
+        filename = joinpath(output_directory, @sprintf("mesh_%s.h5", system))
+    end
 
     # Open file (clobber existing content)
     h5open(filename, "w") do file
diff --git a/src/semidiscretization/semidiscretization.jl b/src/semidiscretization/semidiscretization.jl
index 8fef66d261e..ac312c57c89 100644
--- a/src/semidiscretization/semidiscretization.jl
+++ b/src/semidiscretization/semidiscretization.jl
@@ -76,7 +76,8 @@ function semidiscretize(semi::AbstractSemidiscretization, tspan)
     #       mpi_isparallel() && MPI.Barrier(mpi_comm())
     #       See https://github.com/trixi-framework/Trixi.jl/issues/328
     iip = true # is-inplace, i.e., we modify a vector when calling rhs!
-    return ODEProblem{iip}(rhs!, u0_ode, tspan, semi)
+    specialize = SciMLBase.FullSpecialize # specialize on rhs! and parameters (semi)
+    return ODEProblem{iip, specialize}(rhs!, u0_ode, tspan, semi)
 end
 
 """
@@ -93,7 +94,8 @@ function semidiscretize(semi::AbstractSemidiscretization, tspan,
     #       mpi_isparallel() && MPI.Barrier(mpi_comm())
     #       See https://github.com/trixi-framework/Trixi.jl/issues/328
     iip = true # is-inplace, i.e., we modify a vector when calling rhs!
-    return ODEProblem{iip}(rhs!, u0_ode, tspan, semi)
+    specialize = SciMLBase.FullSpecialize # specialize on rhs! and parameters (semi)
+    return ODEProblem{iip, specialize}(rhs!, u0_ode, tspan, semi)
 end
 
 """
diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl
new file mode 100644
index 00000000000..b7adff78425
--- /dev/null
+++ b/src/semidiscretization/semidiscretization_coupled.jl
@@ -0,0 +1,610 @@
+"""
+    SemidiscretizationCoupled
+
+A struct used to bundle multiple semidiscretizations.
+[`semidiscretize`](@ref) will return an `ODEProblem` that synchronizes time steps between the semidiscretizations.
+Each call of `rhs!` will call `rhs!` for each semidiscretization individually.
+The semidiscretizations can be coupled by gluing meshes together using [`BoundaryConditionCoupled`](@ref).
+
+!!! warning "Experimental code"
+    This is an experimental feature and can change any time.
+"""
+struct SemidiscretizationCoupled{S, Indices, EquationList} <: AbstractSemidiscretization
+    semis::S
+    u_indices::Indices # u_ode[u_indices[i]] is the part of u_ode corresponding to semis[i]
+    performance_counter::PerformanceCounter
+end
+
+"""
+    SemidiscretizationCoupled(semis...)
+
+Create a coupled semidiscretization that consists of the semidiscretizations passed as arguments.
+"""
+function SemidiscretizationCoupled(semis...)
+    @assert all(semi -> ndims(semi) == ndims(semis[1]), semis) "All semidiscretizations must have the same dimension!"
+
+    # Number of coefficients for each semidiscretization
+    n_coefficients = zeros(Int, length(semis))
+    for i in 1:length(semis)
+        _, equations, _, _ = mesh_equations_solver_cache(semis[i])
+        n_coefficients[i] = ndofs(semis[i]) * nvariables(equations)
+    end
+
+    # Compute range of coefficients associated with each semidiscretization and allocate coupled BCs
+    u_indices = Vector{UnitRange{Int}}(undef, length(semis))
+    for i in 1:length(semis)
+        offset = sum(n_coefficients[1:(i - 1)]) + 1
+        u_indices[i] = range(offset, length = n_coefficients[i])
+
+        allocate_coupled_boundary_conditions(semis[i])
+    end
+
+    performance_counter = PerformanceCounter()
+
+    SemidiscretizationCoupled{typeof(semis), typeof(u_indices), typeof(performance_counter)
+                              }(semis, u_indices, performance_counter)
+end
+
+function Base.show(io::IO, semi::SemidiscretizationCoupled)
+    @nospecialize semi # reduce precompilation time
+
+    print(io, "SemidiscretizationCoupled($(semi.semis))")
+end
+
+function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationCoupled)
+    @nospecialize semi # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, semi)
+    else
+        summary_header(io, "SemidiscretizationCoupled")
+        summary_line(io, "#spatial dimensions", ndims(semi.semis[1]))
+        summary_line(io, "#systems", nsystems(semi))
+        for i in eachsystem(semi)
+            summary_line(io, "system", i)
+            mesh, equations, solver, _ = mesh_equations_solver_cache(semi.semis[i])
+            summary_line(increment_indent(io), "mesh", mesh |> typeof |> nameof)
+            summary_line(increment_indent(io), "equations", equations |> typeof |> nameof)
+            summary_line(increment_indent(io), "initial condition",
+                         semi.semis[i].initial_condition)
+            # no boundary conditions since that could be too much
+            summary_line(increment_indent(io), "source terms", semi.semis[i].source_terms)
+            summary_line(increment_indent(io), "solver", solver |> typeof |> nameof)
+        end
+        summary_line(io, "total #DOFs", ndofs(semi))
+        summary_footer(io)
+    end
+end
+
+function print_summary_semidiscretization(io::IO, semi::SemidiscretizationCoupled)
+    show(io, MIME"text/plain"(), semi)
+    println(io, "\n")
+    for i in eachsystem(semi)
+        mesh, equations, solver, _ = mesh_equations_solver_cache(semi.semis[i])
+        summary_header(io, "System #$i")
+
+        summary_line(io, "mesh", mesh |> typeof |> nameof)
+        show(increment_indent(io), MIME"text/plain"(), mesh)
+
+        summary_line(io, "equations", equations |> typeof |> nameof)
+        show(increment_indent(io), MIME"text/plain"(), equations)
+
+        summary_line(io, "solver", solver |> typeof |> nameof)
+        show(increment_indent(io), MIME"text/plain"(), solver)
+
+        summary_footer(io)
+        println(io, "\n")
+    end
+end
+
+@inline Base.ndims(semi::SemidiscretizationCoupled) = ndims(semi.semis[1])
+
+@inline nsystems(semi::SemidiscretizationCoupled) = length(semi.semis)
+
+@inline eachsystem(semi::SemidiscretizationCoupled) = Base.OneTo(nsystems(semi))
+
+@inline Base.real(semi::SemidiscretizationCoupled) = promote_type(real.(semi.semis)...)
+
+@inline Base.eltype(semi::SemidiscretizationCoupled) = promote_type(eltype.(semi.semis)...)
+
+@inline function ndofs(semi::SemidiscretizationCoupled)
+    sum(ndofs, semi.semis)
+end
+
+@inline function nelements(semi::SemidiscretizationCoupled)
+    return sum(semi.semis) do semi_
+        mesh, equations, solver, cache = mesh_equations_solver_cache(semi_)
+
+        nelements(mesh, solver, cache)
+    end
+end
+
+function compute_coefficients(t, semi::SemidiscretizationCoupled)
+    @unpack u_indices = semi
+
+    u_ode = Vector{real(semi)}(undef, u_indices[end][end])
+
+    for i in eachsystem(semi)
+        # Call `compute_coefficients` in `src/semidiscretization/semidiscretization.jl`
+        u_ode[u_indices[i]] .= compute_coefficients(t, semi.semis[i])
+    end
+
+    return u_ode
+end
+
+@inline function get_system_u_ode(u_ode, index, semi::SemidiscretizationCoupled)
+    @view u_ode[semi.u_indices[index]]
+end
+
+function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t)
+    @unpack u_indices = semi
+
+    time_start = time_ns()
+
+    @trixi_timeit timer() "copy to coupled boundaries" begin
+        for semi_ in semi.semis
+            copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi)
+        end
+    end
+
+    # Call rhs! for each semidiscretization
+    for i in eachsystem(semi)
+        u_loc = get_system_u_ode(u_ode, i, semi)
+        du_loc = get_system_u_ode(du_ode, i, semi)
+
+        @trixi_timeit timer() "system #$i" rhs!(du_loc, u_loc, semi.semis[i], t)
+    end
+
+    runtime = time_ns() - time_start
+    put!(semi.performance_counter, runtime)
+
+    return nothing
+end
+
+################################################################################
+### AnalysisCallback
+################################################################################
+
+"""
+    AnalysisCallbackCoupled(semi, callbacks...)
+
+Combine multiple analysis callbacks for coupled simulations with a
+[`SemidiscretizationCoupled`](@ref). For each coupled system, an indididual
+[`AnalysisCallback`](@ref) **must** be created and passed to the `AnalysisCallbackCoupled` **in
+order**, i.e., in the same sequence as the indidvidual semidiscretizations are stored in the
+`SemidiscretizationCoupled`.
+
+!!! warning "Experimental code"
+    This is an experimental feature and can change any time.
+"""
+struct AnalysisCallbackCoupled{CB}
+    callbacks::CB
+end
+
+function Base.show(io::IO, ::MIME"text/plain",
+                   cb_coupled::DiscreteCallback{<:Any, <:AnalysisCallbackCoupled})
+    @nospecialize cb_coupled # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, cb_coupled)
+    else
+        analysis_callback_coupled = cb_coupled.affect!
+
+        summary_header(io, "AnalysisCallbackCoupled")
+        for (i, cb) in enumerate(analysis_callback_coupled.callbacks)
+            summary_line(io, "Callback #$i", "")
+            show(increment_indent(io), MIME"text/plain"(), cb)
+        end
+        summary_footer(io)
+    end
+end
+
+# Convenience constructor for the coupled callback that gets called directly from the elixirs
+function AnalysisCallbackCoupled(semi_coupled, callbacks...)
+    if length(callbacks) != nsystems(semi_coupled)
+        error("an AnalysisCallbackCoupled requires one AnalysisCallback for each semidiscretization")
+    end
+
+    analysis_callback_coupled = AnalysisCallbackCoupled{typeof(callbacks)}(callbacks)
+
+    # This callback is triggered if any of its subsidiary callbacks' condition is triggered
+    condition = (u, t, integrator) -> any(callbacks) do callback
+        callback.condition(u, t, integrator)
+    end
+
+    DiscreteCallback(condition, analysis_callback_coupled,
+                     save_positions = (false, false),
+                     initialize = initialize!)
+end
+
+# This method gets called during initialization from OrdinaryDiffEq's `solve(...)`
+function initialize!(cb_coupled::DiscreteCallback{Condition, Affect!}, u_ode_coupled, t,
+                     integrator) where {Condition, Affect! <: AnalysisCallbackCoupled}
+    analysis_callback_coupled = cb_coupled.affect!
+    semi_coupled = integrator.p
+    du_ode_coupled = first(get_tmp_cache(integrator))
+
+    # Loop over coupled systems' callbacks and initialize them individually
+    for i in eachsystem(semi_coupled)
+        cb = analysis_callback_coupled.callbacks[i]
+        semi = semi_coupled.semis[i]
+        u_ode = get_system_u_ode(u_ode_coupled, i, semi_coupled)
+        du_ode = get_system_u_ode(du_ode_coupled, i, semi_coupled)
+        initialize!(cb, u_ode, du_ode, t, integrator, semi)
+    end
+end
+
+# This method gets called from OrdinaryDiffEq's `solve(...)`
+function (analysis_callback_coupled::AnalysisCallbackCoupled)(integrator)
+    semi_coupled = integrator.p
+    u_ode_coupled = integrator.u
+    du_ode_coupled = first(get_tmp_cache(integrator))
+
+    # Loop over coupled systems' callbacks and call them individually
+    for i in eachsystem(semi_coupled)
+        @unpack condition = analysis_callback_coupled.callbacks[i]
+        analysis_callback = analysis_callback_coupled.callbacks[i].affect!
+        u_ode = get_system_u_ode(u_ode_coupled, i, semi_coupled)
+
+        # Check condition and skip callback if it is not yet its turn
+        if !condition(u_ode, integrator.t, integrator)
+            continue
+        end
+
+        semi = semi_coupled.semis[i]
+        du_ode = get_system_u_ode(du_ode_coupled, i, semi_coupled)
+        analysis_callback(u_ode, du_ode, integrator, semi)
+    end
+end
+
+# used for error checks and EOC analysis
+function (cb::DiscreteCallback{Condition, Affect!})(sol) where {Condition,
+                                                                Affect! <:
+                                                                AnalysisCallbackCoupled}
+    semi_coupled = sol.prob.p
+    u_ode_coupled = sol.u[end]
+    @unpack callbacks = cb.affect!
+
+    uEltype = real(semi_coupled)
+    l2_error_collection = uEltype[]
+    linf_error_collection = uEltype[]
+    for i in eachsystem(semi_coupled)
+        analysis_callback = callbacks[i].affect!
+        @unpack analyzer = analysis_callback
+        cache_analysis = analysis_callback.cache
+
+        semi = semi_coupled.semis[i]
+        u_ode = get_system_u_ode(u_ode_coupled, i, semi_coupled)
+
+        l2_error, linf_error = calc_error_norms(u_ode, sol.t[end], analyzer, semi,
+                                                cache_analysis)
+        append!(l2_error_collection, l2_error)
+        append!(linf_error_collection, linf_error)
+    end
+
+    (; l2 = l2_error_collection, linf = linf_error_collection)
+end
+
+################################################################################
+### SaveSolutionCallback
+################################################################################
+
+# Save mesh for a coupled semidiscretization, which contains multiple meshes internally
+function save_mesh(semi::SemidiscretizationCoupled, output_directory, timestep = 0)
+    for i in eachsystem(semi)
+        mesh, _, _, _ = mesh_equations_solver_cache(semi.semis[i])
+
+        if mesh.unsaved_changes
+            mesh.current_filename = save_mesh_file(mesh, output_directory, system = i)
+            mesh.unsaved_changes = false
+        end
+    end
+end
+
+@inline function save_solution_file(semi::SemidiscretizationCoupled, u_ode,
+                                    solution_callback,
+                                    integrator)
+    @unpack semis = semi
+
+    for i in eachsystem(semi)
+        u_ode_slice = get_system_u_ode(u_ode, i, semi)
+        save_solution_file(semis[i], u_ode_slice, solution_callback, integrator, system = i)
+    end
+end
+
+################################################################################
+### StepsizeCallback
+################################################################################
+
+# In case of coupled system, use minimum timestep over all systems
+function calculate_dt(u_ode, t, cfl_number, semi::SemidiscretizationCoupled)
+    dt = minimum(eachsystem(semi)) do i
+        u_ode_slice = get_system_u_ode(u_ode, i, semi)
+        calculate_dt(u_ode_slice, t, cfl_number, semi.semis[i])
+    end
+
+    return dt
+end
+
+################################################################################
+### Equations
+################################################################################
+
+"""
+    BoundaryConditionCoupled(other_semi_index, indices, uEltype)
+
+Boundary condition to glue two meshes together. Solution values at the boundary
+of another mesh will be used as boundary values. This requires the use
+of [`SemidiscretizationCoupled`](@ref). The other mesh is specified by `other_semi_index`,
+which is the index of the mesh in the tuple of semidiscretizations.
+
+Note that the elements and nodes of the two meshes at the coupled boundary must coincide.
+This is currently only implemented for [`StructuredMesh`](@ref).
+
+# Arguments
+- `other_semi_index`: the index in `SemidiscretizationCoupled` of the semidiscretization
+                      from which the values are copied
+- `indices::Tuple`: node/cell indices at the boundary of the mesh in the other
+                    semidiscretization. See examples below.
+- `uEltype::Type`: element type of solution
+
+# Examples
+```julia
+# Connect the left boundary of mesh 2 to our boundary such that our positive
+# boundary direction will match the positive y direction of the other boundary
+BoundaryConditionCoupled(2, (:begin, :i), Float64)
+
+# Connect the same two boundaries oppositely oriented
+BoundaryConditionCoupled(2, (:begin, :i_backwards), Float64)
+
+# Using this as y_neg boundary will connect `our_cells[i, 1, j]` to `other_cells[j, end-i, end]`
+BoundaryConditionCoupled(2, (:j, :i_backwards, :end), Float64)
+```
+
+!!! warning "Experimental code"
+    This is an experimental feature and can change any time.
+"""
+mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indices}
+    # NDIMST2M1 == NDIMS * 2 - 1
+    # Buffer for boundary values: [variable, nodes_i, nodes_j, cell_i, cell_j]
+    u_boundary        :: Array{uEltype, NDIMST2M1} # NDIMS * 2 - 1
+    other_semi_index  :: Int
+    other_orientation :: Int
+    indices           :: Indices
+
+    function BoundaryConditionCoupled(other_semi_index, indices, uEltype)
+        NDIMS = length(indices)
+        u_boundary = Array{uEltype, NDIMS * 2 - 1}(undef, ntuple(_ -> 0, NDIMS * 2 - 1))
+
+        if indices[1] in (:begin, :end)
+            other_orientation = 1
+        elseif indices[2] in (:begin, :end)
+            other_orientation = 2
+        else # indices[3] in (:begin, :end)
+            other_orientation = 3
+        end
+
+        new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices)}(u_boundary, other_semi_index,
+                                                            other_orientation, indices)
+    end
+end
+
+function Base.eltype(boundary_condition::BoundaryConditionCoupled)
+    eltype(boundary_condition.u_boundary)
+end
+
+function (boundary_condition::BoundaryConditionCoupled)(u_inner, orientation, direction,
+                                                        cell_indices, surface_node_indices,
+                                                        surface_flux_function, equations)
+    # get_node_vars(boundary_condition.u_boundary, equations, solver, surface_node_indices..., cell_indices...),
+    # but we don't have a solver here
+    u_boundary = SVector(ntuple(v -> boundary_condition.u_boundary[v,
+                                                                   surface_node_indices...,
+                                                                   cell_indices...],
+                                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)
+    end
+
+    return flux
+end
+
+function allocate_coupled_boundary_conditions(semi::AbstractSemidiscretization)
+    n_boundaries = 2 * ndims(semi)
+    mesh, equations, solver, _ = mesh_equations_solver_cache(semi)
+
+    for direction in 1:n_boundaries
+        boundary_condition = semi.boundary_conditions[direction]
+
+        allocate_coupled_boundary_condition(boundary_condition, direction, mesh, equations,
+                                            solver)
+    end
+end
+
+# Don't do anything for other BCs than BoundaryConditionCoupled
+function allocate_coupled_boundary_condition(boundary_condition, direction, mesh, equations,
+                                             solver)
+    return nothing
+end
+
+# In 2D
+function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2
+                                                                                          },
+                                             direction, mesh, equations, dg::DGSEM)
+    if direction in (1, 2)
+        cell_size = size(mesh, 2)
+    else
+        cell_size = size(mesh, 1)
+    end
+
+    uEltype = eltype(boundary_condition)
+    boundary_condition.u_boundary = Array{uEltype, 3}(undef, nvariables(equations),
+                                                      nnodes(dg),
+                                                      cell_size)
+end
+
+# Don't do anything for other BCs than BoundaryConditionCoupled
+function copy_to_coupled_boundary!(boundary_condition, u_ode, semi)
+    return nothing
+end
+
+function copy_to_coupled_boundary!(boundary_conditions::Union{Tuple, NamedTuple}, u_ode,
+                                   semi)
+    for boundary_condition in boundary_conditions
+        copy_to_coupled_boundary!(boundary_condition, u_ode, semi)
+    end
+end
+
+# In 2D
+function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{2}, u_ode,
+                                   semi)
+    @unpack u_indices = semi
+    @unpack other_semi_index, other_orientation, indices = boundary_condition
+
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi.semis[other_semi_index])
+    u = wrap_array(get_system_u_ode(u_ode, other_semi_index, semi), mesh, equations, solver,
+                   cache)
+
+    linear_indices = LinearIndices(size(mesh))
+
+    if other_orientation == 1
+        cells = axes(mesh, 2)
+    else # other_orientation == 2
+        cells = axes(mesh, 1)
+    end
+
+    # Copy solution data to the coupled boundary using "delayed indexing" with
+    # a start value and a step size to get the correct face and orientation.
+    node_index_range = eachnode(solver)
+    i_node_start, i_node_step = index_to_start_step_2d(indices[1], node_index_range)
+    j_node_start, j_node_step = index_to_start_step_2d(indices[2], node_index_range)
+
+    i_cell_start, i_cell_step = index_to_start_step_2d(indices[1], axes(mesh, 1))
+    j_cell_start, j_cell_step = index_to_start_step_2d(indices[2], axes(mesh, 2))
+
+    i_cell = i_cell_start
+    j_cell = j_cell_start
+
+    for cell in cells
+        i_node = i_node_start
+        j_node = j_node_start
+
+        for i in eachnode(solver)
+            for v in 1:size(u, 1)
+                boundary_condition.u_boundary[v, i, cell] = u[v, i_node, j_node,
+                                                              linear_indices[i_cell,
+                                                                             j_cell]]
+            end
+            i_node += i_node_step
+            j_node += j_node_step
+        end
+        i_cell += i_cell_step
+        j_cell += j_cell_step
+    end
+end
+
+################################################################################
+### DGSEM/structured
+################################################################################
+
+@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation,
+                                                  boundary_condition::BoundaryConditionCoupled,
+                                                  mesh::StructuredMesh, equations,
+                                                  surface_integral, dg::DG, cache,
+                                                  direction, node_indices,
+                                                  surface_node_indices, element)
+    @unpack node_coordinates, contravariant_vectors, inverse_jacobian = cache.elements
+    @unpack surface_flux = surface_integral
+
+    cell_indices = get_boundary_indices(element, orientation, mesh)
+
+    u_inner = get_node_vars(u, equations, dg, node_indices..., element)
+
+    # If the mapping is orientation-reversing, the contravariant vectors' orientation
+    # is reversed as well. The normal vector must be oriented in the direction
+    # from `left_element` to `right_element`, or the numerical flux will be computed
+    # incorrectly (downwind direction).
+    sign_jacobian = sign(inverse_jacobian[node_indices..., element])
+
+    # Contravariant vector Ja^i is the normal vector
+    normal = sign_jacobian * get_contravariant_vector(orientation, contravariant_vectors,
+                                      node_indices..., element)
+
+    # If the mapping is orientation-reversing, the normal vector will be reversed (see above).
+    # However, the flux now has the wrong sign, since we need the physical flux in normal direction.
+    flux = sign_jacobian * boundary_condition(u_inner, normal, direction, cell_indices,
+                              surface_node_indices, surface_flux, equations)
+
+    for v in eachvariable(equations)
+        surface_flux_values[v, surface_node_indices..., direction, element] = flux[v]
+    end
+end
+
+function get_boundary_indices(element, orientation, mesh::StructuredMesh{2})
+    cartesian_indices = CartesianIndices(size(mesh))
+    if orientation == 1
+        # Get index of element in y-direction
+        cell_indices = (cartesian_indices[element][2],)
+    else # orientation == 2
+        # Get index of element in x-direction
+        cell_indices = (cartesian_indices[element][1],)
+    end
+
+    return cell_indices
+end
+
+################################################################################
+### Special elixirs
+################################################################################
+
+# Analyze convergence for SemidiscretizationCoupled
+function analyze_convergence(errors_coupled, iterations,
+                             semi_coupled::SemidiscretizationCoupled)
+    # Extract errors: the errors are currently stored as
+    # | iter 1 sys 1 var 1...n | iter 1 sys 2 var 1...n | ... | iter 2 sys 1 var 1...n | ...
+    # but for calling `analyze_convergence` below, we need the following layout
+    # sys n: | iter 1 var 1...n | iter 1 var 1...n | ... | iter 2 var 1...n | ...
+    # That is, we need to extract and join the data for a single system
+    errors = Dict{Symbol, Vector{Float64}}[]
+    for i in eachsystem(semi_coupled)
+        push!(errors, Dict(:l2 => Float64[], :linf => Float64[]))
+    end
+    offset = 0
+    for iter in 1:iterations, i in eachsystem(semi_coupled)
+        # Extract information on current semi
+        semi = semi_coupled.semis[i]
+        _, equations, _, _ = mesh_equations_solver_cache(semi)
+        variablenames = varnames(cons2cons, equations)
+
+        # Compute offset
+        first = offset + 1
+        last = offset + length(variablenames)
+        offset += length(variablenames)
+
+        # Append errors to appropriate storage
+        append!(errors[i][:l2], errors_coupled[:l2][first:last])
+        append!(errors[i][:linf], errors_coupled[:linf][first:last])
+    end
+
+    eoc_mean_values = Vector{Dict{Symbol, Any}}(undef, nsystems(semi_coupled))
+    for i in eachsystem(semi_coupled)
+        # Use visual cues to separate output from multiple systems
+        println()
+        println("="^100)
+        println("# System $i")
+        println("="^100)
+
+        # Extract information on current semi
+        semi = semi_coupled.semis[i]
+        _, equations, _, _ = mesh_equations_solver_cache(semi)
+        variablenames = varnames(cons2cons, equations)
+
+        eoc_mean_values[i] = analyze_convergence(errors[i], iterations, variablenames)
+    end
+
+    return eoc_mean_values
+end
diff --git a/src/semidiscretization/semidiscretization_euler_gravity.jl b/src/semidiscretization/semidiscretization_euler_gravity.jl
index 665f2be9bfa..8fe9de1d2b2 100644
--- a/src/semidiscretization/semidiscretization_euler_gravity.jl
+++ b/src/semidiscretization/semidiscretization_euler_gravity.jl
@@ -477,18 +477,29 @@ end
 @inline function save_solution_file(u_ode, t, dt, iter,
                                     semi::SemidiscretizationEulerGravity,
                                     solution_callback,
-                                    element_variables = Dict{Symbol, Any}())
+                                    element_variables = Dict{Symbol, Any}();
+                                    system = "")
+    # If this is called already as part of a multi-system setup (i.e., system is non-empty),
+    # we build a combined system name
+    if !isempty(system)
+        system_euler = system * "_euler"
+        system_gravity = system * "_gravity"
+    else
+        system_euler = "euler"
+        system_gravity = "gravity"
+    end
+
     u_euler = wrap_array_native(u_ode, semi.semi_euler)
     filename_euler = save_solution_file(u_euler, t, dt, iter,
                                         mesh_equations_solver_cache(semi.semi_euler)...,
                                         solution_callback, element_variables,
-                                        system = "euler")
+                                        system = system_euler)
 
     u_gravity = wrap_array_native(semi.cache.u_ode, semi.semi_gravity)
     filename_gravity = save_solution_file(u_gravity, t, dt, iter,
                                           mesh_equations_solver_cache(semi.semi_gravity)...,
                                           solution_callback, element_variables,
-                                          system = "gravity")
+                                          system = system_gravity)
 
     return filename_euler, filename_gravity
 end
diff --git a/test/test_special_elixirs.jl b/test/test_special_elixirs.jl
index 742a3abc376..23017059eaa 100644
--- a/test/test_special_elixirs.jl
+++ b/test/test_special_elixirs.jl
@@ -30,6 +30,12 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none",
         @test isapprox(mean_convergence[:l2], [4.0], rtol=0.05)
       end
 
+      @timed_testset "structured_2d_dgsem coupled" begin
+        mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_coupled.jl"), 3)
+        @test isapprox(mean_convergence[1][:l2], [4.0], rtol=0.05)
+        @test isapprox(mean_convergence[2][:l2], [4.0], rtol=0.05)
+      end
+
       @timed_testset "p4est_2d_dgsem" begin
         # Run convergence test on unrefined mesh
         no_refine = @cfunction((p4est, which_tree, quadrant) -> Cint(0), Cint, (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p4est_quadrant_t}))
@@ -57,6 +63,7 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none",
       @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_basic.jl"), 2, tspan=(0.0, 0.01))
       @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_extended.jl"), 2, initial_refinement_level=0, tspan=(0.0, 0.1))
       @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_basic.jl"), 2, tspan=(0.0, 0.01))
+      @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_coupled.jl"), 2, tspan=(0.0, 0.01))
       @test_nowarn_mod convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_extended.jl"), 2, cells_per_dimension=(1, 1), tspan=(0.0, 0.1))
     end
   end
diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl
index feaf66c4a7f..16fc72f0a46 100644
--- a/test/test_structured_2d.jl
+++ b/test/test_structured_2d.jl
@@ -19,6 +19,19 @@ isdir(outdir) && rm(outdir, recursive=true)
       linf = [6.627000273229378e-5])
   end
 
+  @trixi_testset "elixir_advection_coupled.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_coupled.jl"),    
+      l2   = [7.816742843181738e-6, 7.816742843196112e-6],
+      linf = [6.314906965543265e-5, 6.314906965410039e-5],
+      coverage_override = (maxiters=10^5,))
+
+    @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin
+      errors = analysis_callback(sol)
+      @test errors.l2   ≈ [7.816742843181738e-6, 7.816742843196112e-6] rtol=1.0e-4
+      @test errors.linf ≈ [6.314906965543265e-5, 6.314906965410039e-5] rtol=1.0e-4
+    end
+  end
+
   @trixi_testset "elixir_advection_extended.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
       l2   = [4.220397559713772e-6],

From deb027adefdd88fccf6cec5ce4ca5c76106a0439 Mon Sep 17 00:00:00 2001
From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com>
Date: Fri, 16 Jun 2023 22:53:10 +0200
Subject: [PATCH 2/2] Bump crate-ci/typos from 1.14.12 to 1.15.0 (#1524)

* Bump crate-ci/typos from 1.14.12 to 1.15.0

Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.14.12 to 1.15.0.
- [Release notes](https://github.com/crate-ci/typos/releases)
- [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md)
- [Commits](https://github.com/crate-ci/typos/compare/v1.14.12...v1.15.0)

---
updated-dependencies:
- dependency-name: crate-ci/typos
  dependency-type: direct:production
  update-type: version-update:semver-minor
...

Signed-off-by: dependabot[bot] <support@github.com>

* Fix typos

---------

Signed-off-by: dependabot[bot] <support@github.com>
Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>
Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 .github/workflows/SpellCheck.yml                            | 2 +-
 docs/src/visualization.md                                   | 2 +-
 examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml
index c4ab3a98557..bc324c689bc 100644
--- a/.github/workflows/SpellCheck.yml
+++ b/.github/workflows/SpellCheck.yml
@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v3
       - name: Check spelling
-        uses: crate-ci/typos@v1.14.12
+        uses: crate-ci/typos@v1.15.0
diff --git a/docs/src/visualization.md b/docs/src/visualization.md
index e29313cc080..8f72bb4b1c6 100644
--- a/docs/src/visualization.md
+++ b/docs/src/visualization.md
@@ -339,7 +339,7 @@ create a [`PlotData1D`](@ref) with the keyword argument `curve` set to your list
 Let's give an example of this with the basic advection equation from above by creating
 a plot along the circle marked in green:
 
-![2d-plot-along-cirlce](https://user-images.githubusercontent.com/72009492/130951042-e1849447-8e55-4798-9361-c8badb9f3a49.png)
+![2d-plot-along-circle](https://user-images.githubusercontent.com/72009492/130951042-e1849447-8e55-4798-9361-c8badb9f3a49.png)
 
 We can write a function like this, that outputs a list of points on a circle:
 ```julia
diff --git a/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl
index 42370e861ce..366be700f9f 100644
--- a/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl
+++ b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl
@@ -3,7 +3,7 @@
 # Boundary conditions are supersonic Mach 3 inflow at the left portion of the domain
 # and supersonic outflow at the right portion of the domain. The top and bottom of the
 # channel as well as the cylinder are treated as Euler slip wall boundaries.
-# This flow results in strong shock refletions / interactions as well as Kelvin-Helmholtz
+# This flow results in strong shock reflections / interactions as well as Kelvin-Helmholtz
 # instabilities at later times as two Mach stems form above and below the cylinder.
 #
 # For complete details on the problem setup see Section 5.7 of the paper: