Sc/converters coupling

docs/make.jl

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

docs/src/coupling.md

@@ -0,0 +1,33 @@
+# [Coupling](@id coupling-id)
+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.
+
+Another type of coupling is bulk or volume coupling.
+There we have at least two systems that share all or parts of the domain.
+We could, for instance, have a Maxwell system and a fluid system.
+The coupling would then occur through the Lorentz force.
+
+
+## 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 one system with `m` variables `u_i` and another
+system with `n` variables `v_j`.
+We then define two coupling functions, one that transforms `u_i` into `v_i`
+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_converter.jl`
+and in `src/coupling_converters/coupling_converters_2d.jl`.

examples/structured_2d_dgsem/elixir_advection_coupled.jl

@@ -48,13 +48,15 @@ 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
+boundary_conditions_x_neg1 = BoundaryConditionCoupled(2, (:end, :i_forward), Float64, coupling_converter_identity(equations))
+boundary_conditions_x_pos1 = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64, coupling_converter_identity(equations))
+
 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),
+                                       x_neg=boundary_conditions_x_neg1,
                                        # Connect right boundary with left boundary of right mesh
-                                       x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward),  Float64),
+                                       x_pos=boundary_conditions_x_pos1,
                                        y_neg=boundary_condition_periodic,
                                        y_pos=boundary_condition_periodic))
 
@@ -67,12 +69,15 @@ cells_per_dimension2 = cells_per_dimension
 
 mesh2 = StructuredMesh(cells_per_dimension2, coordinates_min2, coordinates_max2)
 
+boundary_conditions_x_neg2 = BoundaryConditionCoupled(1, (:end, :i_forward), Float64, coupling_converter_identity(equations))
+boundary_conditions_x_pos2 = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64, coupling_converter_identity(equations))
+
 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),
+                                       x_neg=boundary_conditions_x_neg2,
                                        # Connect right boundary with left boundary of left mesh
-                                       x_pos=BoundaryConditionCoupled(1, (:begin, :i_forward),  Float64),
+                                       x_pos=boundary_conditions_x_pos2,
                                        y_neg=boundary_condition_periodic,
                                        y_pos=boundary_condition_periodic))
 

examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl

@@ -0,0 +1,134 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# 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
+# 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)
+
+# The user can define their own coupling functions.
+coupling_function1 = (x, u) -> (sign(x[2] - 0.0)*0.1 + 1.0)/1.1 * 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))
+
+cells_per_dimension2 = cells_per_dimension
+
+mesh2 = StructuredMesh(cells_per_dimension2, coordinates_min2, coordinates_max2)
+
+coupling_function2 = (x, u) -> (sign(x[2] - 0.0)*0.1 + 1.0)/1.1 * 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
+ode = semidiscretize(semi, (0.0, 20.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 = 1,
+                                     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()

src/Trixi.jl

@@ -123,6 +123,7 @@ include("callbacks_step/callbacks_step.jl")
 include("callbacks_stage/callbacks_stage.jl")
 include("semidiscretization/semidiscretization_euler_gravity.jl")
 include("time_integration/time_integration.jl")
+include("coupling_converters/coupling_converters.jl")
 
 # `trixi_include` and special elixirs such as `convergence_test`
 include("auxiliary/special_elixirs.jl")
@@ -189,6 +190,9 @@ export boundary_condition_do_nothing,
        BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal,
        BoundaryConditionCoupled
 
+export coupling_converter_identity,
+       coupling_converter_heaviside_2d, coupling_converter_linear_2d
+
 export initial_condition_convergence_test, source_terms_convergence_test
 export source_terms_harmonic
 export initial_condition_poisson_nonperiodic, source_terms_poisson_nonperiodic,

src/coupling_converters/coupling_converters.jl

@@ -0,0 +1,25 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+@doc raw"""
+    coupling_converter_identity(semi::AbstractSemidiscretization, tspan)
+    Identity coupling converter function.
+
+The coupling is given as a linear function.
+```math
+c(x) = u(x)
+```
+"""
+function coupling_converter_identity(equations::AbstractEquations)
+    return (x, u) -> u
+end
+
+####################################################################################################
+# Include files with actual implementations for different systems of equations.
+
+include("coupling_converters_2d.jl")
+end # @muladd

src/coupling_converters/coupling_converters_2d.jl

@@ -0,0 +1,35 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+@doc raw"""
+    Coupling converter function for a system of two LinearScalarAdvectionEquation2D.
+
+The coupling is given as a Heaviside step.
+```math
+c(x) = {c_0, for x \ge x_0 \times s
+        0, for x < x_0 \times s}
+```
+Here, `s` is the sign of the step function, x_0 the save_position
+of the step and c_0 the amplitude.
+"""
+function coupling_converter_heaviside_2d(x_0, c_0, s,
+                                         equations::LinearScalarAdvectionEquation2D)
+    return (x, u) -> c_0 * (s * sign(x[2] - x_0) + 1.0) / 2.0 * u
+end
+
+@doc raw"""
+    Coupling converter function for a system of two LinearScalarAdvectionEquation2D.
+
+The coupling is given as a linear function.
+```math
+c(x) = x * u(x)
+```
+"""
+function coupling_converter_linear_2d(equations::LinearScalarAdvectionEquation2D)
+    return (x, u) -> x[2] * u
+end
+end # @muladd

src/semidiscretization/semidiscretization_coupled.jl

@@ -364,15 +364,17 @@ 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, CouplingFunction}
     # 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::CouplingFunction
+
+    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 +386,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
 
@@ -495,9 +500,13 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{
 
         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]]
+                x = cache.elements.node_coordinates[:, i_node, j_node,
+                                                    linear_indices[i_cell, j_cell]]
+                converted_u = boundary_condition.coupling_converter(x,
+                                                                    u[:, i_node, j_node,
+                                                                      linear_indices[i_cell,
+                                                                                     j_cell]])
+                boundary_condition.u_boundary[v, i, cell] = converted_u[v]
             end
             i_node += i_node_step
             j_node += j_node_step

test/test_structured_2d.jl

@@ -22,14 +22,13 @@ 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],
-      coverage_override = (maxiters=10^5,))
+      linf = [6.314906965543265e-5, 6.314906965410039e-5])
+  end
 
-    @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
+  @trixi_testset "elixir_advection_coupled_converter.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_coupled_converter.jl"),    
+      l2   = [0.3495477674652473, 0.3472339065154432],
+      linf = [0.5569080960939969, 0.537610538307045])
   end
 
   @trixi_testset "elixir_advection_extended.jl" begin