diff --git a/docs/make.jl b/docs/make.jl index dee87371bd1..7fce3b31e24 100644 --- a/docs/make.jl +++ b/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", diff --git a/docs/src/multi-physics_coupling.md b/docs/src/multi-physics_coupling.md new file mode 100644 index 00000000000..eec92bc21de --- /dev/null +++ b/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. + diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index 2a56d23f4c0..43b68f21b03 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/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, diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 0941ae6a8ca..dc21dbe9a1e 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -1,3 +1,10 @@ +# 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 + """ SemidiscretizationCoupled @@ -65,11 +72,13 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationCoupled) 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), "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), "source terms", + semi.semis[i].source_terms) summary_line(increment_indent(io), "solver", solver |> typeof |> nameof) end summary_line(io, "total #DOFs per field", ndofs(semi)) @@ -106,20 +115,14 @@ end @inline Base.real(semi::SemidiscretizationCoupled) = promote_type(real.(semi.semis)...) -@inline Base.eltype(semi::SemidiscretizationCoupled) = promote_type(eltype.(semi.semis)...) +@inline function Base.eltype(semi::SemidiscretizationCoupled) + promote_type(eltype.(semi.semis)...) +end @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 @@ -137,23 +140,40 @@ end @view u_ode[semi.u_indices[index]] end +# Same as `foreach(enumerate(something))`, but without allocations. +# +# Note that compile times may increase if this is used with big tuples. +@inline foreach_enumerate(func, collection) = foreach_enumerate(func, collection, 1) +@inline foreach_enumerate(func, collection::Tuple{}, index) = nothing + +@inline function foreach_enumerate(func, collection, index) + element = first(collection) + remaining_collection = Base.tail(collection) + + func((index, element)) + + # Process remaining collection + foreach_enumerate(func, remaining_collection, index + 1) +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) + foreach(semi.semis) do semi_ + copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi, 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) + @trixi_timeit timer() "copy to coupled boundaries" begin + foreach_enumerate(semi.semis) do (i, semi_) + u_loc = get_system_u_ode(u_ode, i, semi) + du_loc = get_system_u_ode(du_ode, i, semi) + rhs!(du_loc, u_loc, semi_, t) + end end runtime = time_ns() - time_start @@ -309,7 +329,8 @@ end 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) + save_solution_file(semis[i], u_ode_slice, solution_callback, integrator, + system = i) end end @@ -332,7 +353,7 @@ end ################################################################################ """ - BoundaryConditionCoupled(other_semi_index, indices, uEltype) + BoundaryConditionCoupled(other_semi_index, indices, uEltype, coupling_converter) 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 @@ -348,32 +369,37 @@ This is currently only implemented for [`StructuredMesh`](@ref). - `indices::Tuple`: node/cell indices at the boundary of the mesh in the other semidiscretization. See examples below. - `uEltype::Type`: element type of solution +- `coupling_converter::CouplingConverter`: function to call for converting the solution + state of one system to the other system # 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) +BoundaryConditionCoupled(2, (:begin, :i), Float64, fun) # Connect the same two boundaries oppositely oriented -BoundaryConditionCoupled(2, (:begin, :i_backwards), Float64) +BoundaryConditionCoupled(2, (:begin, :i_backwards), Float64, fun) # 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) +BoundaryConditionCoupled(2, (:j, :i_backwards, :end), Float64, fun) ``` !!! 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)) @@ -385,8 +411,10 @@ 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 @@ -395,8 +423,10 @@ function Base.eltype(boundary_condition::BoundaryConditionCoupled) end function (boundary_condition::BoundaryConditionCoupled)(u_inner, orientation, direction, - cell_indices, surface_node_indices, - surface_flux_function, equations) + 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, @@ -421,13 +451,15 @@ function allocate_coupled_boundary_conditions(semi::AbstractSemidiscretization) for direction in 1:n_boundaries boundary_condition = semi.boundary_conditions[direction] - allocate_coupled_boundary_condition(boundary_condition, direction, mesh, equations, + 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, +function allocate_coupled_boundary_condition(boundary_condition, direction, mesh, + equations, solver) return nothing end @@ -448,43 +480,69 @@ function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditi end # Don't do anything for other BCs than BoundaryConditionCoupled -function copy_to_coupled_boundary!(boundary_condition, u_ode, semi) +function copy_to_coupled_boundary!(boundary_condition, u_ode, semi_coupled, semi) return nothing end +function copy_to_coupled_boundary!(u_ode, semi_coupled, semi, i, n_boundaries, + boundary_condition, boundary_conditions...) + copy_to_coupled_boundary!(boundary_condition, u_ode, semi_coupled, semi) + if i < n_boundaries + copy_to_coupled_boundary!(u_ode, semi_coupled, semi, i + 1, n_boundaries, + boundary_conditions...) + end +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) + semi_coupled, semi) + copy_to_coupled_boundary!(u_ode, semi_coupled, semi, 1, length(boundary_conditions), + boundary_conditions...) +end + +function mesh_equations_solver_cache(other_semi_index, i, semi_, semi_tuple...) + if i == other_semi_index + return mesh_equations_solver_cache(semi_) + else + # Walk through semidiscretizations until we find `i` + mesh_equations_solver_cache(other_semi_index, i + 1, semi_tuple...) end end # In 2D -function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{2}, u_ode, - semi) - @unpack u_indices = semi +function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{2}, + u_ode, + semi_coupled, semi) + @unpack u_indices = semi_coupled @unpack other_semi_index, other_orientation, indices = boundary_condition + @unpack coupling_converter, u_boundary = boundary_condition + + mesh_own, equations_own, solver_own, cache_own = mesh_equations_solver_cache(semi) + + mesh_other, equations_other, solver_other, cache_other = mesh_equations_solver_cache(other_semi_index, + 1, + semi_coupled.semis...) - 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) + node_coordinates_other = cache_other.elements.node_coordinates + u_ode_other = get_system_u_ode(u_ode, other_semi_index, semi_coupled) + u_other = wrap_array(u_ode_other, mesh_other, equations_other, solver_other, + cache_other) - linear_indices = LinearIndices(size(mesh)) + linear_indices = LinearIndices(size(mesh_other)) if other_orientation == 1 - cells = axes(mesh, 2) + cells = axes(mesh_other, 2) else # other_orientation == 2 - cells = axes(mesh, 1) + cells = axes(mesh_other, 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) + node_index_range = eachnode(solver_other) 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_start, i_cell_step = index_to_start_step_2d(indices[1], axes(mesh_other, 1)) + j_cell_start, j_cell_step = index_to_start_step_2d(indices[2], axes(mesh_other, 2)) i_cell = i_cell_start j_cell = j_cell_start @@ -492,16 +550,26 @@ 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]] + element_id = linear_indices[i_cell, j_cell] + + for element_id in eachnode(solver_other) + x_other = get_node_coords(node_coordinates_other, equations_other, + solver_other, + i_node, j_node, linear_indices[i_cell, j_cell]) + u_node_other = get_node_vars(u_other, equations_other, solver_other, i_node, + j_node, linear_indices[i_cell, j_cell]) + u_node_converted = coupling_converter(x_other, u_node_other, + equations_other, + equations_own) + + for i in eachindex(u_node_converted) + u_boundary[i, element_id, cell] = u_node_converted[i] end + i_node += i_node_step j_node += j_node_step end + i_cell += i_cell_step j_cell += j_cell_step end @@ -511,7 +579,8 @@ end ### DGSEM/structured ################################################################################ -@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, +@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, + orientation, boundary_condition::BoundaryConditionCoupled, mesh::StructuredMesh, equations, surface_integral, dg::DG, cache, @@ -531,7 +600,8 @@ end 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, + 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). @@ -608,3 +678,4 @@ function analyze_convergence(errors_coupled, iterations, return eoc_mean_values end +end # @muladd diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index e5d45ebcc07..522510a42e3 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -33,14 +33,34 @@ 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], + l2=[ + 7.816742843336293e-6, + 7.816742843340186e-6, + 7.816742843025513e-6, + 7.816742843061526e-6, + ], + linf=[ + 6.314906965276812e-5, + 6.314906965187994e-5, + 6.31490696496595e-5, + 6.314906965032563e-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≈[ + 7.816742843336293e-6, + 7.816742843340186e-6, + 7.816742843025513e-6, + 7.816742843061526e-6, + ] rtol=1.0e-4 + @test errors.linf≈[ + 6.314906965276812e-5, + 6.314906965187994e-5, + 6.31490696496595e-5, + 6.314906965032563e-5, + ] rtol=1.0e-4 # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let