Sc/converters coupling

docs/make.jl

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

docs/src/multi-physics_coupling.md

@@ -0,0 +1,37 @@
+# [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
+We can have the case where the two systems do not share any variables, but
+share some of the physics.
+Here, the same physics is just represented in a different form and with
+different variables.
+This is the case for a fluid system on one side and a Vlasov system on the other.
+To translate the fields from one description to the other one needs to use
+converter functions.
+
+In the general case, we have a system A with $m$ variables $u_{A,i}$ and another
+system B with $n$ variables $u_{B,j}$.
+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$ and state vector $u$
+and return the transformed variables.
+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

@@ -1,9 +1,10 @@
 using OrdinaryDiffEq
 using Trixi
 
-
 ###############################################################################
-# Coupled semidiscretization of two linear advection systems, which are connected periodically
+# Coupled semidiscretization of two linear advection systems using converter functions such that
+# the upper half of the domain is coupled periodically, while the lower half is not coupled
+# and any incoming wave is completely absorbed.
 #
 # 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
@@ -35,11 +36,11 @@ 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)
+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))
+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)
@@ -48,70 +49,86 @@ 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))
+# The user can define their own coupling functions.
+coupling_function1 = (x, u) -> u
 
+boundary_conditions_x_neg1 = BoundaryConditionCoupled(2, (:end, :i_forward), Float64,
+                                                      coupling_function1)
+boundary_conditions_x_pos1 = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64,
+                                                      coupling_function1)
+
+# 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 = boundary_conditions_x_neg1,
+                                                            # Connect right boundary with left boundary of right mesh
+                                                            x_pos = boundary_conditions_x_pos1,
+                                                            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))
+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))
+coupling_function2 = (x, u) -> u
+
+boundary_conditions_x_neg2 = BoundaryConditionCoupled(1, (:end, :i_forward), Float64,
+                                                      coupling_function2)
+boundary_conditions_x_pos2 = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64,
+                                                      coupling_function2)
+
+semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test,
+                                     solver,
+                                     boundary_conditions = (
+                                                            # Connect left boundary with right boundary of left mesh
+                                                            x_neg = boundary_conditions_x_neg2,
+                                                            # Connect right boundary with left boundary of left mesh
+                                                            x_pos = boundary_conditions_x_pos2,
+                                                            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
+# Create ODE problem with time span from 0.0 to 20.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_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)
+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)
+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)
-
+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);
+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()

src/Trixi.jl

@@ -125,6 +125,7 @@ include("time_integration/time_integration.jl")
 include("callbacks_step/callbacks_step.jl")
 include("callbacks_stage/callbacks_stage.jl")
 include("semidiscretization/semidiscretization_euler_gravity.jl")
+include("coupling_converters/coupling_converters.jl")
 
 # `trixi_include` and special elixirs such as `convergence_test`
 include("auxiliary/special_elixirs.jl")

src/semidiscretization/semidiscretization_coupled.jl

@@ -364,15 +364,18 @@ 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}
+mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indices,
+                                        CouplingConverter}
     # 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)
+    u_boundary::Array{uEltype, NDIMST2M1} # NDIMS * 2 - 1
+    other_semi_index::Int
+    other_orientation::Int
+    indices::Indices
+    coupling_converter::CouplingConverter
+
+    function BoundaryConditionCoupled(other_semi_index, indices, uEltype,
+                                      coupling_converter)
         NDIMS = length(indices)
         u_boundary = Array{uEltype, NDIMS * 2 - 1}(undef, ntuple(_ -> 0, NDIMS * 2 - 1))
 
@@ -384,8 +387,11 @@ mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indic
             other_orientation = 3
         end
 
-        new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices)}(u_boundary, other_semi_index,
-                                                            other_orientation, indices)
+        new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices), typeof(coupling_converter)}(u_boundary,
+                                                                                        other_semi_index,
+                                                                                        other_orientation,
+                                                                                        indices,
+                                                                                        coupling_converter)
     end
 end
 
@@ -464,8 +470,10 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{
                                    semi)
     @unpack u_indices = semi
     @unpack other_semi_index, other_orientation, indices = boundary_condition
+    @unpack coupling_converter, u_boundary = boundary_condition
 
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi.semis[other_semi_index])
+    @unpack node_coordinates = cache.elements
     u = wrap_array(get_system_u_ode(u_ode, other_semi_index, semi), mesh, equations, solver,
                    cache)
 
@@ -492,13 +500,14 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{
     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
+        element_id = linear_indices[i_cell, j_cell]
+
+        for element_id in eachnode(solver)
+            x = get_node_vars(node_coordinates, equations, solver, i_node, j_node,
+                              element_id)
+            u_node = get_node_vars(u, equations, solver, i_node, j_node, element_id)
+            converted_u = coupling_converter(x, u_node)
+            boundary_condition.u_boundary[:, element_id, cell] = converted_u
             i_node += i_node_step
             j_node += j_node_step
         end

test/test_structured_2d.jl

@@ -32,14 +32,14 @@ isdir(outdir) && rm(outdir, recursive=true)
 
   @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],
+      l2   = [8.336466488044185e-6, 8.336466488028504e-6],
+      linf = [6.406463687191888e-5, 6.406463686947639e-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
+      @test errors.l2   ≈ [8.336466488044185e-6, 8.336466488028504e-6] rtol=1.0e-4
+      @test errors.linf ≈ [6.406463687191888e-5, 6.406463686947639e-5] rtol=1.0e-4
       # Ensure that we do not have excessive memory allocations 
       # (e.g., from type instabilities) 
       let