Sc/converters coupling

docs/make.jl

@@ -108,6 +108,7 @@ makedocs(
             ],
             "Time integration" => "time_integration.md",
             "Callbacks" => "callbacks.md",
+            "Coupling" => "multi-physics_coupling.md"
         ],
         "Advanced topics & developers" => [
             "Conventions" =>"conventions.md",

docs/src/multi-physics_coupling.md

@@ -0,0 +1,46 @@
+# [Multi-physics coupling](@id multi-physics-coupling)
+A complex simulation can consist of different spatial domains in which
+different equations are being solved, different numerical methods being used
+or the grid structure is different.
+One example would be a fluid in a tank and an extended hot plate attached to it.
+We would then like to solve the Navier-Stokes equations in the fluid domain
+and the heat conduction equations in the plate.
+The coupling would happen at the interface through the exchange of thermal energy.
+
+
+## Converter coupling
+It may happen that the two systems to be coupled do not share any variables, but
+share some of the physics.
+In such a situation, the same physics is just represented in a different form and with
+a different set of variables.
+This is the case, for instance assuming two domains, if there is a fluid system in one domain
+and a Vlasov system in the other domain.
+In that case we would have variables representing distribution functions of
+the Vlasov system on one side and variables representing the mechanical quantities, like density,
+of the fluid system.
+To translate the fields from one description to the other one needs to use
+converter functions.
+These functions need to be hand tailored by the user in the elixir file where each
+pair of coupled systems requires two coupling functions, one for each direction.
+
+In the general case, we have a system $A$ with $m$ variables
+$u_{A,i}, \: i = 1, \dots, m$ and another
+system $B$ with $n$ variables $u_{B,j}, \: j = 1, \dots, n$.
+We then define two coupling functions, one that transforms $u_A$ into $u_B$
+and one that goes the other way.
+
+In their minimal form they take the position vector $x$, state vector $u$
+and the equations of the two coupled systems
+and return the transformed variables.
+By passing the equations we can make use of their parameters, if they are required.
+Examples can be seen in `examples/structured_2d_dgsem/elixir_advection_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.
+This is an issue coming from the coordinate interpolation on the nodes.
+As a result, running a simulation in a single system and in two coupled domains
+may result in a difference of the order of the machine precision.
+While this is not an issue for most practical problems, it is best to keep this in mind when comparing test runs.
+

examples/structured_2d_dgsem/elixir_advection_coupled.jl

@@ -2,31 +2,38 @@ using OrdinaryDiffEq
 using Trixi
 
 ###############################################################################
-# Coupled semidiscretization of two linear advection systems, which are connected periodically
+# Coupled semidiscretization of four linear advection 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 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:
+# In this elixir, we have a square domain that is divided into a upper-left, lower-left,
+# upper-right and lower-right quarter. On each quarter
+# of the domain, a completely independent SemidiscretizationHyperbolic is created for the
+# linear advection equations. The four systems are coupled in the x and y-direction.
+# For a high-level overview, see also the figure below:
 #
 # (-1,  1)                                   ( 1,  1)
 #     ┌────────────────────┬────────────────────┐
-#     │    ↑ periodic ↑    │    ↑ periodic ↑    │
-#     │                    │                    │
+#     │    ↑ coupled ↑     │    ↑ coupled ↑     │
 #     │                    │                    │
 #     │     =========      │     =========      │
 #     │     system #1      │     system #2      │
 #     │     =========      │     =========      │
 #     │                    │                    │
+#     │<-- coupled         │<-- coupled         │
+#     │         coupled -->│         coupled -->│
 #     │                    │                    │
+#     │    ↓ coupled ↓     │    ↓ coupled ↓     │
+#     ├────────────────────┼────────────────────┤
+#     │    ↑ coupled ↑     │    ↑ coupled ↑     │
 #     │                    │                    │
+#     │     =========      │     =========      │
+#     │     system #3      │     system #4      │
+#     │     =========      │     =========      │
 #     │                    │                    │
-#     │         coupled -->│<-- coupled         │
-#     │                    │                    │
-#     │<-- coupled         │         coupled -->│
-#     │                    │                    │
+#     │<-- coupled         │<-- coupled         │
+#     │         coupled -->│         coupled -->│
 #     │                    │                    │
-#     │    ↓ periodic ↓    │    ↓ periodic ↓    │
+#     │    ↓ coupled  ↓    │    ↓ coupled  ↓    │
 #     └────────────────────┴────────────────────┘
 # (-1, -1)                                   ( 1, -1)
 
@@ -36,60 +43,135 @@ 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))
+# This will be the number of elements for each quarter/semidiscretization.
+cells_per_dimension = (8, 8)
+
+###########
+# system #1
+###########
+
+coordinates_min1 = (-1.0, 0.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)
+mesh1 = StructuredMesh(cells_per_dimension, coordinates_min1, coordinates_max1)
 
-cells_per_dimension1 = cells_per_dimension
+# Define the coupling functions
+coupling_function12 = (x, u, equations_other, equations_own) -> u
+coupling_function13 = (x, u, equations_other, equations_own) -> u
 
-mesh1 = StructuredMesh(cells_per_dimension1, coordinates_min1, coordinates_max1)
+# Define the coupling boundary conditions and the system it is coupled to.
+boundary_conditions_x_neg1 = BoundaryConditionCoupled(2, (:end, :i_forward), Float64,
+                                                      coupling_function12)
+boundary_conditions_x_pos1 = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64,
+                                                      coupling_function12)
+boundary_conditions_y_neg1 = BoundaryConditionCoupled(3, (:i_forward, :end), Float64,
+                                                      coupling_function13)
+boundary_conditions_y_pos1 = BoundaryConditionCoupled(3, (:i_forward, :begin), Float64,
+                                                      coupling_function13)
 
 # 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))
+                                     boundary_conditions = (x_neg = boundary_conditions_x_neg1,
+                                                            x_pos = boundary_conditions_x_pos1,
+                                                            y_neg = boundary_conditions_y_neg1,
+                                                            y_pos = boundary_conditions_y_pos1))
+
+###########
+# system #2
+###########
+
+coordinates_min2 = (0.0, 0.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_dimension, coordinates_min2, coordinates_max2)
 
-mesh2 = StructuredMesh(cells_per_dimension2, coordinates_min2, coordinates_max2)
+# Define the coupling functions
+coupling_function21 = (x, u, equations_other, equations_own) -> u
+coupling_function24 = (x, u, equations_other, equations_own) -> u
 
+# Define the coupling boundary conditions and the system it is coupled to.
+boundary_conditions_x_neg2 = BoundaryConditionCoupled(1, (:end, :i_forward), Float64,
+                                                      coupling_function21)
+boundary_conditions_x_pos2 = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64,
+                                                      coupling_function21)
+boundary_conditions_y_neg2 = BoundaryConditionCoupled(4, (:i_forward, :end), Float64,
+                                                      coupling_function24)
+boundary_conditions_y_pos2 = BoundaryConditionCoupled(4, (:i_forward, :begin), Float64,
+                                                      coupling_function24)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
 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)
+                                     boundary_conditions = (x_neg = boundary_conditions_x_neg2,
+                                                            x_pos = boundary_conditions_x_pos2,
+                                                            y_neg = boundary_conditions_y_neg2,
+                                                            y_pos = boundary_conditions_y_pos2))
+
+###########
+# system #3
+###########
+
+coordinates_min3 = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max3 = (0.0, 0.0) # maximum coordinates (max(x), max(y))
+
+mesh3 = StructuredMesh(cells_per_dimension, coordinates_min3, coordinates_max3)
+
+# Define the coupling functions
+coupling_function34 = (x, u, equations_other, equations_own) -> u
+coupling_function31 = (x, u, equations_other, equations_own) -> u
+
+# Define the coupling boundary conditions and the system it is coupled to.
+boundary_conditions_x_neg3 = BoundaryConditionCoupled(4, (:end, :i_forward), Float64,
+                                                      coupling_function34)
+boundary_conditions_x_pos3 = BoundaryConditionCoupled(4, (:begin, :i_forward), Float64,
+                                                      coupling_function34)
+boundary_conditions_y_neg3 = BoundaryConditionCoupled(1, (:i_forward, :end), Float64,
+                                                      coupling_function31)
+boundary_conditions_y_pos3 = BoundaryConditionCoupled(1, (:i_forward, :begin), Float64,
+                                                      coupling_function31)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi3 = SemidiscretizationHyperbolic(mesh3, equations, initial_condition_convergence_test,
+                                     solver,
+                                     boundary_conditions = (x_neg = boundary_conditions_x_neg3,
+                                                            x_pos = boundary_conditions_x_pos3,
+                                                            y_neg = boundary_conditions_y_neg3,
+                                                            y_pos = boundary_conditions_y_pos3))
+
+###########
+# system #4
+###########
+
+coordinates_min4 = (0.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max4 = (1.0, 0.0) # maximum coordinates (max(x), max(y))
+
+mesh4 = StructuredMesh(cells_per_dimension, coordinates_min4, coordinates_max4)
+
+# Define the coupling functions
+coupling_function43 = (x, u, equations_other, equations_own) -> u
+coupling_function42 = (x, u, equations_other, equations_own) -> u
+
+# Define the coupling boundary conditions and the system it is coupled to.
+boundary_conditions_x_neg4 = BoundaryConditionCoupled(3, (:end, :i_forward), Float64,
+                                                      coupling_function43)
+boundary_conditions_x_pos4 = BoundaryConditionCoupled(3, (:begin, :i_forward), Float64,
+                                                      coupling_function43)
+boundary_conditions_y_neg4 = BoundaryConditionCoupled(2, (:i_forward, :end), Float64,
+                                                      coupling_function42)
+boundary_conditions_y_pos4 = BoundaryConditionCoupled(2, (:i_forward, :begin), Float64,
+                                                      coupling_function42)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi4 = SemidiscretizationHyperbolic(mesh4, equations, initial_condition_convergence_test,
+                                     solver,
+                                     boundary_conditions = (x_neg = boundary_conditions_x_neg4,
+                                                            x_pos = boundary_conditions_x_pos4,
+                                                            y_neg = boundary_conditions_y_neg4,
+                                                            y_pos = boundary_conditions_y_pos4))
+
+# Create a semidiscretization that bundles all the semidiscretizations.
+semi = SemidiscretizationCoupled(semi1, semi2, semi3, semi4)
 
 ###############################################################################
 # ODE solvers, callbacks etc.
@@ -104,7 +186,10 @@ 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)
+analysis_callback3 = AnalysisCallback(semi3, interval = 100)
+analysis_callback4 = AnalysisCallback(semi4, interval = 100)
+analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2,
+                                            analysis_callback3, analysis_callback4)
 
 # The SaveSolutionCallback allows to save the solution to a file in regular intervals
 save_solution = SaveSolutionCallback(interval = 100,