Sc/converters coupling

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,136 @@
+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)
+
+coupling_function1 = coupling_converter_heaviside_2d(-0.5, 1.0, 1.0, equations)
+# 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 = coupling_converter_heaviside_2d(-0.5, 1.0, 1.0, equations)
+# 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,36 @@
+# 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_converters
+
+Define converter functions for two coupled systems.
+These should be used together with SemidiscretizationCoupled.
+Using converter functions we can couple two systems that do not
+share any variables.
+This is done by taking the last inner point of system i, apply
+a converter function on the state vector u_i and obtain a state
+vector u_j for the boundary of system j.
+"""
+
+@doc raw"""
+    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

@@ -367,12 +367,14 @@ BoundaryConditionCoupled(2, (:j, :i_backwards, :end), Float64)
 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)
+    u_boundary::Array{uEltype, NDIMST2M1} # NDIMS * 2 - 1
+    other_semi_index::Int
+    other_orientation::Int
+    indices::Indices
+    coupling_converter::Function
+
+    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))
 
@@ -385,7 +387,8 @@ mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indic
         end
 
         new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices)}(u_boundary, other_semi_index,
-                                                            other_orientation, indices)
+                                                            other_orientation, indices,
+                                                            coupling_converter)
     end
 end
 
@@ -495,9 +498,16 @@ 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]
+                # 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

test/test_structured_2d.jl

@@ -32,6 +32,13 @@ isdir(outdir) && rm(outdir, recursive=true)
     end
   end
 
+  @trixi_testset "elixir_advection_coupled_converter.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_coupled_converter.jl"),    
+      l2   = [0.97560175611287, 0.9707973849860191],
+      linf = [1.5703274355039958, 1.6235401582169442],
+      coverage_override = (maxiters=10^5,))
+  end
+
   @trixi_testset "elixir_advection_extended.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
       l2   = [4.220397559713772e-6],