From bb074193b9f8b1f66be706ef2a5db8cad7fa6b6b Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 26 Jun 2023 16:31:21 +0100 Subject: [PATCH 001/117] Added coupling converters. --- src/Trixi.jl | 2 ++ .../coupling_converters.jl | 13 ++++++++++ .../coupling_converters_2d.jl | 25 +++++++++++++++++++ 3 files changed, 40 insertions(+) create mode 100644 src/coupling_converters/coupling_converters.jl create mode 100644 src/coupling_converters/coupling_converters_2d.jl diff --git a/src/Trixi.jl b/src/Trixi.jl index 66878f4b459..85b536352bd 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -189,6 +189,8 @@ export boundary_condition_do_nothing, BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal, BoundaryConditionCoupled +export coupling_converter_heaviside_2d + export initial_condition_convergence_test, source_terms_convergence_test export source_terms_harmonic export initial_condition_poisson_nonperiodic, source_terms_poisson_nonperiodic, diff --git a/src/coupling_converters/coupling_converters.jl b/src/coupling_converters/coupling_converters.jl new file mode 100644 index 00000000000..a82b849f284 --- /dev/null +++ b/src/coupling_converters/coupling_converters.jl @@ -0,0 +1,13 @@ +# 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 + +#################################################################################################### +# Include files with actual implementations for different systems of equations. + +include("coupling_converters_2d.jl") + +end # @muladd diff --git a/src/coupling_converters/coupling_converters_2d.jl b/src/coupling_converters/coupling_converters_2d.jl new file mode 100644 index 00000000000..f18058632c0 --- /dev/null +++ b/src/coupling_converters/coupling_converters_2d.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 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, x_0, c_0, s, + equations_left::LinearScalarAdvectionEquation2D, + equation_right::LinearScalarAdvectionEquation2D) + return c_0 * (s*sign(x[2] - x_0) + 1.0)/2.0 +end + +end # @muladd From 95892bf8348689fa19e303f0d37c71cb563c3a05 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 29 Jun 2023 14:23:06 +0100 Subject: [PATCH 002/117] Added generic converter_function for structured 2d meshes. --- src/Trixi.jl | 1 + .../coupling_converters_2d.jl | 20 +++++++++++++++---- .../semidiscretization_coupled.jl | 14 ++++++++----- 3 files changed, 26 insertions(+), 9 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 85b536352bd..0759ba129d2 100644 --- a/src/Trixi.jl +++ b/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") diff --git a/src/coupling_converters/coupling_converters_2d.jl b/src/coupling_converters/coupling_converters_2d.jl index f18058632c0..e360ffc4eac 100644 --- a/src/coupling_converters/coupling_converters_2d.jl +++ b/src/coupling_converters/coupling_converters_2d.jl @@ -5,6 +5,7 @@ @muladd begin #! format: noindent + @doc raw""" Coupling converter function for a system of two LinearScalarAdvectionEquation2D. @@ -16,10 +17,21 @@ c(x) = {c_0, for x \ge 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, x_0, c_0, s, - equations_left::LinearScalarAdvectionEquation2D, - equation_right::LinearScalarAdvectionEquation2D) - return c_0 * (s*sign(x[2] - x_0) + 1.0)/2.0 +function coupling_converter_heaviside_2d(x_0, c_0, s) + 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() + return (x, u) -> x[2]*u end end # @muladd diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index b7adff78425..babf4527a57 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -371,8 +371,9 @@ mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indic other_semi_index :: Int other_orientation :: Int indices :: Indices + coupling_converter :: Function - function BoundaryConditionCoupled(other_semi_index, indices, uEltype) + 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 +386,7 @@ 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 +496,12 @@ 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 From 4e82276d4717eca3a78bc4d083d5dc47ca516240 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 3 Jul 2023 16:33:47 +0100 Subject: [PATCH 003/117] Added example elixir for coupling converters. --- .../elixir_advection_coupled_converter.jl | 127 ++++++++++++++++++ 1 file changed, 127 insertions(+) create mode 100644 examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl b/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl new file mode 100644 index 00000000000..5546209aa19 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl @@ -0,0 +1,127 @@ +using OrdinaryDiffEq +using Trixi +using Trixi2Vtk + + +############################################################################### +# Coupled semidiscretization of two linear advection systems, which are connected periodically +# +# 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) +# coupling_function1 = (x, u) -> (sign(x[2] - 0.0)*0.1 + 1.0)/1.1 * u + +# 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, coupling_function1), + # Connect right boundary with left boundary of right mesh + x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward), Float64, coupling_function1), + 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 + +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, coupling_function2), + # Connect right boundary with left boundary of left mesh + x_pos=BoundaryConditionCoupled(1, (:begin, :i_forward), Float64, coupling_function2), + 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() + +# Convert the snapshot to vtk data. +trixi2vtk("out/solution_*.h5") From 7355ed02d0425a7735154d839f6880c73f7dfcc1 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 3 Jul 2023 17:24:58 +0100 Subject: [PATCH 004/117] Cleaned up converter coupling elixir. --- .../elixir_advection_coupled_converter.jl | 23 +++++++++++-------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl b/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl index 5546209aa19..8aada6bdb2b 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl @@ -1,10 +1,10 @@ using OrdinaryDiffEq using Trixi -using Trixi2Vtk - ############################################################################### -# 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 @@ -50,15 +50,19 @@ 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=BoundaryConditionCoupled(2, (:end, :i_forward), Float64, coupling_function1), + x_neg=boundary_conditions_x_neg1, # Connect right boundary with left boundary of right mesh - x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward), Float64, coupling_function1), + x_pos=boundary_conditions_x_pos1, y_neg=boundary_condition_periodic, y_pos=boundary_condition_periodic)) @@ -74,12 +78,15 @@ 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=BoundaryConditionCoupled(1, (:end, :i_forward), Float64, coupling_function2), + x_neg=boundary_conditions_x_neg2, # Connect right boundary with left boundary of left mesh - x_pos=BoundaryConditionCoupled(1, (:begin, :i_forward), Float64, coupling_function2), + x_pos=boundary_conditions_x_pos2, y_neg=boundary_condition_periodic, y_pos=boundary_condition_periodic)) @@ -123,5 +130,3 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), # Print the timer summary summary_callback() -# Convert the snapshot to vtk data. -trixi2vtk("out/solution_*.h5") From b8cabf41796ed89fb8a233732884b2b94946174a Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 3 Jul 2023 17:25:17 +0100 Subject: [PATCH 005/117] Added equations in coupling converters. --- src/coupling_converters/coupling_converters_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coupling_converters/coupling_converters_2d.jl b/src/coupling_converters/coupling_converters_2d.jl index e360ffc4eac..82ee854d826 100644 --- a/src/coupling_converters/coupling_converters_2d.jl +++ b/src/coupling_converters/coupling_converters_2d.jl @@ -17,7 +17,7 @@ c(x) = {c_0, for x \ge 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) +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 @@ -30,7 +30,7 @@ The coupling is given as a linear function. c(x) = x * u(x) ``` """ -function coupling_converter_linear_2d() +function coupling_converter_linear_2d(equations::LinearScalarAdvectionEquation2D) return (x, u) -> x[2]*u end From 6d521a7f292c8d6a78cfc616bf473c7f8128a06e Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 4 Jul 2023 14:48:28 +0100 Subject: [PATCH 006/117] Added converter functions. --- .../elixir_advection_coupled.jl | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index 1e54e411db6..1e3146906b9 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/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)) From ab8e7507fb103a2859d4a4234151297428ad37e1 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 4 Jul 2023 14:49:24 +0100 Subject: [PATCH 007/117] Added identity converter function. --- src/Trixi.jl | 3 ++- .../coupling_converters.jl | 25 +++++++++++++++++++ .../coupling_converters_2d.jl | 1 - test/test_structured_2d.jl | 7 ++++++ 4 files changed, 34 insertions(+), 2 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 0759ba129d2..613b45e5e1f 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -190,7 +190,8 @@ export boundary_condition_do_nothing, BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal, BoundaryConditionCoupled -export coupling_converter_heaviside_2d +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 diff --git a/src/coupling_converters/coupling_converters.jl b/src/coupling_converters/coupling_converters.jl index a82b849f284..53ebf900e41 100644 --- a/src/coupling_converters/coupling_converters.jl +++ b/src/coupling_converters/coupling_converters.jl @@ -5,6 +5,31 @@ @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. diff --git a/src/coupling_converters/coupling_converters_2d.jl b/src/coupling_converters/coupling_converters_2d.jl index 82ee854d826..34011f4e1f7 100644 --- a/src/coupling_converters/coupling_converters_2d.jl +++ b/src/coupling_converters/coupling_converters_2d.jl @@ -5,7 +5,6 @@ @muladd begin #! format: noindent - @doc raw""" Coupling converter function for a system of two LinearScalarAdvectionEquation2D. diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 16fc72f0a46..d593f6e5a21 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -32,6 +32,13 @@ isdir(outdir) && rm(outdir, recursive=true) end end + @trixi_testset "elixir_advection_coupled_converters.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_coupled_converters.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], From 24773c7830f87b1b8775ec41d04c7809959de911 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 5 Jul 2023 15:59:58 +0100 Subject: [PATCH 008/117] Autoformat for converter coupling implementation. --- .../coupling_converters.jl | 2 -- .../coupling_converters_2d.jl | 9 +++---- .../semidiscretization_coupled.jl | 26 ++++++++++++------- 3 files changed, 20 insertions(+), 17 deletions(-) diff --git a/src/coupling_converters/coupling_converters.jl b/src/coupling_converters/coupling_converters.jl index 53ebf900e41..5856cb373ea 100644 --- a/src/coupling_converters/coupling_converters.jl +++ b/src/coupling_converters/coupling_converters.jl @@ -29,10 +29,8 @@ 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 diff --git a/src/coupling_converters/coupling_converters_2d.jl b/src/coupling_converters/coupling_converters_2d.jl index 34011f4e1f7..d527fe96894 100644 --- a/src/coupling_converters/coupling_converters_2d.jl +++ b/src/coupling_converters/coupling_converters_2d.jl @@ -16,11 +16,11 @@ c(x) = {c_0, for x \ge 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 +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. @@ -30,7 +30,6 @@ c(x) = x * u(x) ``` """ function coupling_converter_linear_2d(equations::LinearScalarAdvectionEquation2D) - return (x, u) -> x[2]*u + return (x, u) -> x[2] * u end - end # @muladd diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index babf4527a57..b24546394a5 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -367,13 +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 - coupling_converter :: Function - - function BoundaryConditionCoupled(other_semi_index, indices, uEltype, coupling_converter) + 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)) @@ -386,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, coupling_converter) + other_orientation, indices, + coupling_converter) end end @@ -496,8 +498,12 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ for i in eachnode(solver) for v in 1:size(u, 1) - 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]]) + 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, From c7da2cc809d507e17e39a27c781fcf6e0f78031a Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 5 Jul 2023 16:04:50 +0100 Subject: [PATCH 009/117] Added coupled converter elixir. --- .../elixir_advection_coupled_converter.jl | 76 ++++++++++--------- 1 file changed, 40 insertions(+), 36 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl b/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl index 8aada6bdb2b..18009ef4f92 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl @@ -36,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) @@ -53,23 +53,25 @@ 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) +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)) - +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 @@ -78,17 +80,20 @@ 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)) +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) @@ -104,29 +109,28 @@ ode = semidiscretize(semi, (0.0, 20.0)); 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=1, - solution_variables=cons2prim) +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) +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() - From 98bf2cf462eab663e3de55aca25ce9eaff93ed7a Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 5 Jul 2023 16:06:47 +0100 Subject: [PATCH 010/117] Corrected file name of coupled converters test. --- test/test_structured_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index d593f6e5a21..2182c197c8b 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -32,8 +32,8 @@ isdir(outdir) && rm(outdir, recursive=true) end end - @trixi_testset "elixir_advection_coupled_converters.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_coupled_converters.jl"), + @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,)) From d97302dc48df0561584094d42952a6b216550d6d Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 18 Jul 2023 16:36:05 +0100 Subject: [PATCH 011/117] Removed redundant doc string. --- src/coupling_converters/coupling_converters.jl | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/coupling_converters/coupling_converters.jl b/src/coupling_converters/coupling_converters.jl index 5856cb373ea..3c2e08e34ac 100644 --- a/src/coupling_converters/coupling_converters.jl +++ b/src/coupling_converters/coupling_converters.jl @@ -5,18 +5,6 @@ @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. From 4a22bfea3274641cdacceea286fd768d4d876ee8 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 18 Jul 2023 16:40:06 +0100 Subject: [PATCH 012/117] Added function signature in doc string. --- src/coupling_converters/coupling_converters.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/coupling_converters/coupling_converters.jl b/src/coupling_converters/coupling_converters.jl index 3c2e08e34ac..b4ecb96bf62 100644 --- a/src/coupling_converters/coupling_converters.jl +++ b/src/coupling_converters/coupling_converters.jl @@ -6,6 +6,7 @@ #! format: noindent @doc raw""" + coupling_converter_identity(semi::AbstractSemidiscretization, tspan) Identity coupling converter function. The coupling is given as a linear function. From f1f6ee8970b13e2b09a09f31ea10507c2cf501b4 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 18 Jul 2023 16:44:57 +0100 Subject: [PATCH 013/117] Removed coverage_override in coupled tests. --- test/test_structured_2d.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 2182c197c8b..5875ef829ac 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -22,8 +22,7 @@ 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]) @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin errors = analysis_callback(sol) @@ -35,8 +34,7 @@ isdir(outdir) && rm(outdir, recursive=true) @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,)) + linf = [1.5703274355039958, 1.6235401582169442]) end @trixi_testset "elixir_advection_extended.jl" begin From 7a7def09b11262b1f0da17dead4084f4081be10f Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 18 Jul 2023 16:46:25 +0100 Subject: [PATCH 014/117] Removed old commented code. --- src/semidiscretization/semidiscretization_coupled.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index b24546394a5..940a1125983 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -505,9 +505,6 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ 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 From 9045c5d95674a41fd68e162ce41304bbed015d47 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Wed, 19 Jul 2023 12:21:56 +0100 Subject: [PATCH 015/117] Update make.jl Added interface coupling docs to the main menu. --- docs/make.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/make.jl b/docs/make.jl index 5069e4dc49a..59827f1e99b 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -113,6 +113,7 @@ makedocs( "Testing" => "testing.md", "Performance" => "performance.md", "Parallelization" => "parallelization.md", + "Coupling" => "coupling.md" ], "Troubleshooting and FAQ" => "troubleshooting.md", "Reference" => [ From 31cf07f42a7b60f976d53e8777a1c8794f22ef58 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Wed, 19 Jul 2023 12:27:20 +0100 Subject: [PATCH 016/117] Update make.jl Moved converter coupling section. --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index 59827f1e99b..a95581f0350 100644 --- a/docs/make.jl +++ b/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", @@ -113,7 +114,6 @@ makedocs( "Testing" => "testing.md", "Performance" => "performance.md", "Parallelization" => "parallelization.md", - "Coupling" => "coupling.md" ], "Troubleshooting and FAQ" => "troubleshooting.md", "Reference" => [ From 2f19315ea8361a17c7ba2d37e783f574b239c0e4 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Wed, 19 Jul 2023 12:36:22 +0100 Subject: [PATCH 017/117] Create coupling.md --- docs/src/coupling.md | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/src/coupling.md diff --git a/docs/src/coupling.md b/docs/src/coupling.md new file mode 100644 index 00000000000..e0a9d916e88 --- /dev/null +++ b/docs/src/coupling.md @@ -0,0 +1 @@ +# [Coupling](@id coupling-id) From 26d8dd3871a163ace1e97168f453a8a6e85c7d57 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Wed, 19 Jul 2023 13:05:08 +0100 Subject: [PATCH 018/117] Update coupling.md Added some documentation on coupling converters. --- docs/src/coupling.md | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/docs/src/coupling.md b/docs/src/coupling.md index e0a9d916e88..cce3594e140 100644 --- a/docs/src/coupling.md +++ b/docs/src/coupling.md @@ -1 +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`. From 05579d12722f440728f4e2fcba2ec5d43c33c76b Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 20 Jul 2023 11:24:19 +0100 Subject: [PATCH 019/117] Removed troublesome AnalysisCallbackCoupled from test. --- test/test_structured_2d.jl | 6 ------ 1 file changed, 6 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 5875ef829ac..5ffe73cbd34 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -23,12 +23,6 @@ isdir(outdir) && rm(outdir, recursive=true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_coupled.jl"), l2 = [7.816742843181738e-6, 7.816742843196112e-6], linf = [6.314906965543265e-5, 6.314906965410039e-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 - end end @trixi_testset "elixir_advection_coupled_converter.jl" begin From 84c872e1fcdd0612b1a162875d14b1d616387d71 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 20 Jul 2023 12:37:41 +0100 Subject: [PATCH 020/117] Chenged coupling converter function. --- .../structured_2d_dgsem/elixir_advection_coupled_converter.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl b/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl index 18009ef4f92..766fbff8a14 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl @@ -49,9 +49,8 @@ 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 +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) From 78f994874b28b70baf52f2b3a3e0942d69baa0ca Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 20 Jul 2023 12:46:26 +0100 Subject: [PATCH 021/117] Changed coupling converter function and updated tests. --- .../structured_2d_dgsem/elixir_advection_coupled_converter.jl | 3 +-- test/test_structured_2d.jl | 4 ++-- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl b/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl index 766fbff8a14..0648952faaf 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl @@ -76,8 +76,7 @@ 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 +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) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 5ffe73cbd34..3cc55aa8da7 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -27,8 +27,8 @@ isdir(outdir) && rm(outdir, recursive=true) @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]) + l2 = [0.3495477674652473, 0.3472339065154432], + linf = [0.5569080960939969, 0.537610538307045]) end @trixi_testset "elixir_advection_extended.jl" begin From 2e350469ba58eaccb203cf3327354475929f17fc Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 24 Jul 2023 14:45:54 +0100 Subject: [PATCH 022/117] Sepcialized coupling function call. --- src/semidiscretization/semidiscretization_coupled.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 940a1125983..00818e40b69 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -364,14 +364,14 @@ 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 - coupling_converter::Function + coupling_converter::CouplingFunction function BoundaryConditionCoupled(other_semi_index, indices, uEltype, coupling_converter) @@ -386,9 +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, - coupling_converter) + new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices), typeof(coupling_converter)}(u_boundary, + other_semi_index, + other_orientation, + indices, + coupling_converter) end end From 6bc7c3be6fb72d9b8ac9695ee4b76014eecde5f7 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 25 Jul 2023 14:56:15 +0100 Subject: [PATCH 023/117] Removed volume coupling from documentation to avoit confusion. --- docs/src/coupling.md | 5 ----- 1 file changed, 5 deletions(-) diff --git a/docs/src/coupling.md b/docs/src/coupling.md index cce3594e140..a10b90f7ee8 100644 --- a/docs/src/coupling.md +++ b/docs/src/coupling.md @@ -7,11 +7,6 @@ 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 From b4decb7a070209590a28c3704133fd9e7930ef83 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 25 Jul 2023 14:57:51 +0100 Subject: [PATCH 024/117] Update src/coupling_converters/coupling_converters.jl Co-authored-by: Hendrik Ranocha --- src/coupling_converters/coupling_converters.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/coupling_converters/coupling_converters.jl b/src/coupling_converters/coupling_converters.jl index b4ecb96bf62..7e7d815d9a7 100644 --- a/src/coupling_converters/coupling_converters.jl +++ b/src/coupling_converters/coupling_converters.jl @@ -7,7 +7,8 @@ @doc raw""" coupling_converter_identity(semi::AbstractSemidiscretization, tspan) - Identity coupling converter function. + +Identity coupling converter function. The coupling is given as a linear function. ```math From 54d8180173953b9e4b0525d504c2da67488ac210 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 25 Jul 2023 14:59:34 +0100 Subject: [PATCH 025/117] Removed redundant converter function for coupling. --- .../coupling_converters_2d.jl | 35 ------------------- 1 file changed, 35 deletions(-) delete mode 100644 src/coupling_converters/coupling_converters_2d.jl diff --git a/src/coupling_converters/coupling_converters_2d.jl b/src/coupling_converters/coupling_converters_2d.jl deleted file mode 100644 index d527fe96894..00000000000 --- a/src/coupling_converters/coupling_converters_2d.jl +++ /dev/null @@ -1,35 +0,0 @@ -# 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 From 160c7bf895d4639affedee70049e05309a63b427 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 25 Jul 2023 15:00:33 +0100 Subject: [PATCH 026/117] Removed redundant coupling converter file mentioned in some files. --- docs/src/coupling.md | 3 +-- src/coupling_converters/coupling_converters.jl | 1 - 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/docs/src/coupling.md b/docs/src/coupling.md index a10b90f7ee8..10194801cab 100644 --- a/docs/src/coupling.md +++ b/docs/src/coupling.md @@ -24,5 +24,4 @@ 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 can be seen in `examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl`. diff --git a/src/coupling_converters/coupling_converters.jl b/src/coupling_converters/coupling_converters.jl index b4ecb96bf62..b97801883d6 100644 --- a/src/coupling_converters/coupling_converters.jl +++ b/src/coupling_converters/coupling_converters.jl @@ -21,5 +21,4 @@ end #################################################################################################### # Include files with actual implementations for different systems of equations. -include("coupling_converters_2d.jl") end # @muladd From cab33c274eec963344faecef0b8c46f6d7e395ab Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 25 Jul 2023 15:25:40 +0100 Subject: [PATCH 027/117] Autoreformatted. --- src/semidiscretization/semidiscretization_coupled.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 00818e40b69..e16750fa32a 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -364,7 +364,8 @@ 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, CouplingFunction} +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 From 024106c3d2f6c44d40c38f078a52e628c1626153 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 25 Jul 2023 17:05:45 +0100 Subject: [PATCH 028/117] Removed old coupled elixir and replaced it with one using converter functions. --- docs/src/coupling.md | 2 +- .../elixir_advection_coupled.jl | 88 +++++++----- .../elixir_advection_coupled_converter.jl | 134 ------------------ test/test_structured_2d.jl | 6 - 4 files changed, 51 insertions(+), 179 deletions(-) delete mode 100644 examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl diff --git a/docs/src/coupling.md b/docs/src/coupling.md index 10194801cab..f715bfbdc0e 100644 --- a/docs/src/coupling.md +++ b/docs/src/coupling.md @@ -24,4 +24,4 @@ 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`. +Examples can be seen in `examples/structured_2d_dgsem/elixir_advection_coupled.jl`. diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index 1e3146906b9..0648952faaf 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/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,38 +49,49 @@ cells_per_dimension1 = cells_per_dimension mesh1 = StructuredMesh(cells_per_dimension1, coordinates_min1, coordinates_max1) -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=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)) - +# 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)) +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) -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)) +coupling_function2 = (x, u) -> (sign(x[2] - 0.0)*0.1 + 1.0)/1.1 * u -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)) +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) @@ -88,35 +100,35 @@ 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, 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_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 = 1, + 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() diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl b/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl deleted file mode 100644 index 0648952faaf..00000000000 --- a/examples/structured_2d_dgsem/elixir_advection_coupled_converter.jl +++ /dev/null @@ -1,134 +0,0 @@ -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() diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 3cc55aa8da7..634d243cb8b 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -20,12 +20,6 @@ isdir(outdir) && rm(outdir, recursive=true) 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]) - 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]) From c6d6da673db628870f4008209d040f4c4554fcf9 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 1 Aug 2023 10:04:16 +0100 Subject: [PATCH 029/117] Updated errors for coupled tests. --- test/test_structured_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 75937ba82ad..6adce677ee4 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -23,8 +23,8 @@ 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 = [0.3495477674652473, 0.3472339065154432], + linf = [0.5569080960939969, 0.537610538307045], coverage_override = (maxiters=10^5,)) @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin From 62d1408797a8553cb9b1c392f568e9cd50b0763b Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 1 Aug 2023 15:07:03 +0100 Subject: [PATCH 030/117] Corrected test results for coupled equations. --- test/test_structured_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 6adce677ee4..0bf0a3923f2 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -29,8 +29,8 @@ isdir(outdir) && rm(outdir, recursive=true) @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 ≈ [0.3495477674652473, 0.3472339065154432] rtol=1.0e-4 + @test errors.linf ≈ [0.5569080960939969, 0.537610538307045] rtol=1.0e-4 end end From 09fd6919c3a24643e61614af497927db85e410d8 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 2 Aug 2023 12:10:15 +0100 Subject: [PATCH 031/117] Corrected comment. --- examples/structured_2d_dgsem/elixir_advection_coupled.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index 0648952faaf..b5eeb9543de 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl @@ -99,7 +99,7 @@ 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, 20.0)); # At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup From df67f00d78bc389c6fc5c55422ffd271a2367498 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 2 Aug 2023 12:10:48 +0100 Subject: [PATCH 032/117] Removed coupled test from special tests. --- test/test_special_elixirs.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_special_elixirs.jl b/test/test_special_elixirs.jl index 23017059eaa..6cea1497ba1 100644 --- a/test/test_special_elixirs.jl +++ b/test/test_special_elixirs.jl @@ -32,6 +32,7 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none", @timed_testset "structured_2d_dgsem coupled" begin mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_coupled.jl"), 3) + println("mean_convergence = ", mean_convergence) @test isapprox(mean_convergence[1][:l2], [4.0], rtol=0.05) @test isapprox(mean_convergence[2][:l2], [4.0], rtol=0.05) end From 35980a1bcb3d87d485f529af525537ebfce9a477 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 2 Aug 2023 12:12:11 +0100 Subject: [PATCH 033/117] Removed coupled test from specials. --- test/test_special_elixirs.jl | 7 ------- 1 file changed, 7 deletions(-) diff --git a/test/test_special_elixirs.jl b/test/test_special_elixirs.jl index 6cea1497ba1..d97b99bf75a 100644 --- a/test/test_special_elixirs.jl +++ b/test/test_special_elixirs.jl @@ -30,13 +30,6 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none", @test isapprox(mean_convergence[:l2], [4.0], rtol=0.05) end - @timed_testset "structured_2d_dgsem coupled" begin - mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_coupled.jl"), 3) - println("mean_convergence = ", mean_convergence) - @test isapprox(mean_convergence[1][:l2], [4.0], rtol=0.05) - @test isapprox(mean_convergence[2][:l2], [4.0], rtol=0.05) - end - @timed_testset "p4est_2d_dgsem" begin # Run convergence test on unrefined mesh no_refine = @cfunction((p4est, which_tree, quadrant) -> Cint(0), Cint, (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p4est_quadrant_t})) From aac84302e390612b0162408ed4a809ab716b2be1 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 2 Aug 2023 17:01:27 +0100 Subject: [PATCH 034/117] Chaned the coupling function to the identity. --- examples/structured_2d_dgsem/elixir_advection_coupled.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index b5eeb9543de..75bac2e759b 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl @@ -50,7 +50,7 @@ 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 +coupling_function1 = coupling_converter_identity(equations) boundary_conditions_x_neg1 = BoundaryConditionCoupled(2, (:end, :i_forward), Float64, coupling_function1) @@ -76,7 +76,7 @@ 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 +coupling_function2 = coupling_converter_identity(equations) boundary_conditions_x_neg2 = BoundaryConditionCoupled(1, (:end, :i_forward), Float64, coupling_function2) From eaeec5b27347bc369c96e19bcaad3a939793f36d Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 2 Aug 2023 17:01:57 +0100 Subject: [PATCH 035/117] Updated coupling tests. --- test/test_structured_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 0bf0a3923f2..ccb881d20c0 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -23,8 +23,8 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "elixir_advection_coupled.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_coupled.jl"), - l2 = [0.3495477674652473, 0.3472339065154432], - linf = [0.5569080960939969, 0.537610538307045], + 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 From f1619383fb519a6b1692af194e3695a835f9eb65 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 3 Aug 2023 08:41:16 +0100 Subject: [PATCH 036/117] Updated errors for coupled test. --- test/test_structured_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index ccb881d20c0..ba452d5ff4d 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -29,8 +29,8 @@ isdir(outdir) && rm(outdir, recursive=true) @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin errors = analysis_callback(sol) - @test errors.l2 ≈ [0.3495477674652473, 0.3472339065154432] rtol=1.0e-4 - @test errors.linf ≈ [0.5569080960939969, 0.537610538307045] 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 end end From 569da04831b7541e2cdf3043c7b2dfc133a51688 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 7 Sep 2023 19:24:03 +0100 Subject: [PATCH 037/117] Added advice about binary compatability for coupled equations in the documentation. --- docs/src/coupling.md | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/docs/src/coupling.md b/docs/src/coupling.md index f715bfbdc0e..070abac6ca0 100644 --- a/docs/src/coupling.md +++ b/docs/src/coupling.md @@ -25,3 +25,13 @@ 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 Compatability +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. + From 93fc57187230686daf62d205deea1f179f0c4af7 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 19 Sep 2023 13:01:42 +0100 Subject: [PATCH 038/117] Typo. --- docs/src/coupling.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/coupling.md b/docs/src/coupling.md index 070abac6ca0..ca9a26ce360 100644 --- a/docs/src/coupling.md +++ b/docs/src/coupling.md @@ -27,7 +27,7 @@ and return the transformed variables. Examples can be seen in `examples/structured_2d_dgsem/elixir_advection_coupled.jl`. -## Warning about Binary Compatability +## 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. From ed1fec950a12f24c54b241cd858faf233075a62f Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 5 Oct 2023 16:17:15 +0100 Subject: [PATCH 039/117] Added numerical fluxes. --- .../compressible_euler_multicomponent_2d.jl | 173 ++++++++++++++++++ 1 file changed, 173 insertions(+) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 7b437f4a1b4..47eaa5924ee 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -270,6 +270,29 @@ end return vcat(f_other, f_rho) end +# Calculate 1D flux for a single point +@inline function flux(u, normal_direction::AbstractVector, + equations::CompressibleEulerMulticomponentEquations2D) + rho_v1, rho_v2, rho_e = u + + rho = density(u, equations) + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + gamma = totalgamma(u, equations) + p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) + + f_rho = densities(u, v_normal, equations) + f1 = rho_v1 * v_normal + p * normal_direction[1] + f2 = rho_v2 * v_normal + p * normal_direction[2] + f3 = (rho_e + p) * v_normal + + f_other = SVector{3, real(equations)}(f1, f2, f3) + + return vcat(f_other, f_rho) +end + """ flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations2D) @@ -446,6 +469,76 @@ See also return vcat(f_other, f_rho) end +@inline function flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerMulticomponentEquations2D) + # Unpack left and right state + @unpack gammas, gas_constants, cv = equations + rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll + rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr + rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], + u_rr[i + 3]) + for i in eachcomponent(equations)) + rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + + u_rr[i + 3]) + for i in eachcomponent(equations)) + + # Iterating over all partial densities + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + # Calculating gamma + gamma = totalgamma(0.5 * (u_ll + u_rr), equations) + inv_gamma_minus_one = 1 / (gamma - 1) + + # extract velocities + v1_ll = rho_v1_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_ll = rho_v2_ll / rho_ll + v2_rr = rho_v2_rr / rho_rr + v2_avg = 0.5 * (v2_ll + v2_rr) + velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # helpful variables + help1_ll = zero(v1_ll) + help1_rr = zero(v1_rr) + enth_ll = zero(v1_ll) + enth_rr = zero(v1_rr) + for i in eachcomponent(equations) + enth_ll += u_ll[i + 3] * gas_constants[i] + enth_rr += u_rr[i + 3] * gas_constants[i] + help1_ll += u_ll[i + 3] * cv[i] + help1_rr += u_rr[i + 3] * cv[i] + end + + # temperature and pressure + T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll + T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr + p_ll = T_ll * enth_ll + p_rr = T_rr * enth_rr + p_avg = 0.5 * (p_ll + p_rr) + inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) + + f_rho_sum = zero(T_rr) + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5 * (v_dot_n_ll + v_dot_n_rr) + for i in eachcomponent(equations)) + for i in eachcomponent(equations) + f_rho_sum += f_rho[i] + end + f1 = f_rho_sum * v1_avg + p_avg * normal_direction[1] + f2 = f_rho_sum * v2_avg + p_avg * normal_direction[2] + f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) + + # momentum and energy flux + f_other = SVector{3, real(equations)}(f1, f2, f3) + + return vcat(f_other, f_rho) +end + + # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D) @@ -476,6 +569,30 @@ end λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) end +@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerMulticomponentEquations2D) + rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll + rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr + + # Get the density and gas gamma + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + gamma_ll = totalgamma(u_ll, equations) + gamma_rr = totalgamma(u_rr, equations) + + # Get the velocities based on direction + v_ll = rho_v1_ll / rho_ll * normal_direction[1] + rho_v1_ll / rho_ll * normal_direction[2] + v_rr = rho_v1_rr / rho_rr * normal_direction[1] + rho_v1_rr / rho_rr * normal_direction[2] + + # Compute the sound speeds on the left and right + p_ll = (gamma_ll - 1) * (rho_e_ll - 1 / 2 * (rho_v1_ll^2 + rho_v2_ll^2) / rho_ll) + c_ll = sqrt(gamma_ll * p_ll / rho_ll) + p_rr = (gamma_rr - 1) * (rho_e_rr - 1 / 2 * (rho_v1_rr^2 + rho_v2_rr^2) / rho_rr) + c_rr = sqrt(gamma_rr * p_rr / rho_rr) + + λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) +end + @inline function max_abs_speeds(u, equations::CompressibleEulerMulticomponentEquations2D) rho_v1, rho_v2, rho_e = u @@ -491,6 +608,62 @@ end return (abs(v1) + c, abs(v2) + c) end +@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerMulticomponentEquations2D) + rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll + rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr + + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + + gamma = totalgamma(u_ll, equations) + p_ll = (gamma - 1) * (rho_e_ll - 1 / 2 * rho_ll * (v1_ll^2 + v2_ll^2)) + p_rr = (gamma - 1) * (rho_e_rr - 1 / 2 * rho_rr * (v1_rr^2 + v2_rr^2)) + + if orientation == 1 # x-direction + λ_min = v1_ll - sqrt(gamma * p_ll / rho_ll) + λ_max = v1_rr + sqrt(gamma * p_rr / rho_rr) + else # y-direction + λ_min = v2_ll - sqrt(gamma * p_ll / rho_ll) + λ_max = v2_rr + sqrt(gamma * p_rr / rho_rr) + end + + return λ_min, λ_max +end + +@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerMulticomponentEquations2D) + rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll + rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr + + rho_ll = density(u_ll, equations) + rho_rr = density(u_rr, equations) + + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + + gamma = totalgamma(u_ll, equations) + p_ll = (gamma - 1) * (rho_e_ll - 1 / 2 * rho_ll * (v1_ll^2 + v2_ll^2)) + p_rr = (gamma - 1) * (rho_e_rr - 1 / 2 * rho_rr * (v1_rr^2 + v2_rr^2)) + + v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_normal_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + norm_ = norm(normal_direction) + # The v_normals are already scaled by the norm + λ_min = v_normal_ll - sqrt(gamma * p_ll / rho_ll) * norm_ + λ_max = v_normal_rr + sqrt(gamma * p_rr / rho_rr) * norm_ + + return λ_min, λ_max +end + # Convert conservative variables to primitive @inline function cons2prim(u, equations::CompressibleEulerMulticomponentEquations2D) rho_v1, rho_v2, rho_e = u From 3ea484d11db92359c771cc53a6aeee586e961c35 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 6 Oct 2023 16:51:10 +0100 Subject: [PATCH 040/117] Corrected rs copy routine. Now loop over this semi's components. --- src/semidiscretization/semidiscretization_coupled.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index e16750fa32a..f2d8c05b350 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -500,7 +500,7 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ j_node = j_node_start for i in eachnode(solver) - for v in 1:size(u, 1) + for v in 1:size(boundary_condition.u_boundary, 1) x = cache.elements.node_coordinates[:, i_node, j_node, linear_indices[i_cell, j_cell]] converted_u = boundary_condition.coupling_converter(x, From bfd4ec14c810a8c56b75efb43511ec84109fa428 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 11 Oct 2023 17:22:43 +0100 Subject: [PATCH 041/117] Reformatted equations source file. --- .../compressible_euler_multicomponent_2d.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 47eaa5924ee..a13acf7b1fc 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -522,7 +522,8 @@ end inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) f_rho_sum = zero(T_rr) - f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5 * + (v_dot_n_ll + v_dot_n_rr) for i in eachcomponent(equations)) for i in eachcomponent(equations) f_rho_sum += f_rho[i] @@ -538,7 +539,6 @@ end return vcat(f_other, f_rho) end - # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D) @@ -581,8 +581,10 @@ end gamma_rr = totalgamma(u_rr, equations) # Get the velocities based on direction - v_ll = rho_v1_ll / rho_ll * normal_direction[1] + rho_v1_ll / rho_ll * normal_direction[2] - v_rr = rho_v1_rr / rho_rr * normal_direction[1] + rho_v1_rr / rho_rr * normal_direction[2] + v_ll = rho_v1_ll / rho_ll * normal_direction[1] + + rho_v1_ll / rho_ll * normal_direction[2] + v_rr = rho_v1_rr / rho_rr * normal_direction[1] + + rho_v1_rr / rho_rr * normal_direction[2] # Compute the sound speeds on the left and right p_ll = (gamma_ll - 1) * (rho_e_ll - 1 / 2 * (rho_v1_ll^2 + rho_v2_ll^2) / rho_ll) @@ -643,7 +645,7 @@ end rho_ll = density(u_ll, equations) rho_rr = density(u_rr, equations) - + v1_ll = rho_v1_ll / rho_ll v2_ll = rho_v2_ll / rho_ll v1_rr = rho_v1_rr / rho_rr @@ -652,7 +654,7 @@ end gamma = totalgamma(u_ll, equations) p_ll = (gamma - 1) * (rho_e_ll - 1 / 2 * rho_ll * (v1_ll^2 + v2_ll^2)) p_rr = (gamma - 1) * (rho_e_rr - 1 / 2 * rho_rr * (v1_rr^2 + v2_rr^2)) - + v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] v_normal_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] From c4bc8b3edd0201839de99ea9634ae4da662674d6 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 16 Oct 2023 10:41:01 +0100 Subject: [PATCH 042/117] Removed problemating include of time_integration.jl. --- src/Trixi.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 3142adcef20..0bbc83819e8 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -125,7 +125,6 @@ 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("time_integration/time_integration.jl") include("coupling_converters/coupling_converters.jl") # `trixi_include` and special elixirs such as `convergence_test` From cdd984a3b02ea378097c9a9dbddc1d65655d8cfa Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 16 Oct 2023 17:06:29 +0100 Subject: [PATCH 043/117] Removed export of deleted methods. --- src/Trixi.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 0bbc83819e8..483e1ea5a28 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -199,8 +199,7 @@ export boundary_condition_do_nothing, BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal, BoundaryConditionCoupled -export coupling_converter_identity, - coupling_converter_heaviside_2d, coupling_converter_linear_2d +export coupling_converter_identity export initial_condition_convergence_test, source_terms_convergence_test export source_terms_harmonic From 6eb92b6664f9d6b846bdd9cc5f8cbd6141a9ea5f Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 20 Oct 2023 14:52:47 +0100 Subject: [PATCH 044/117] Reverted to old version of compressible Euler multicomponent with no support for structured grid. --- .../compressible_euler_multicomponent_2d.jl | 175 ------------------ 1 file changed, 175 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index a13acf7b1fc..7b437f4a1b4 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -270,29 +270,6 @@ end return vcat(f_other, f_rho) end -# Calculate 1D flux for a single point -@inline function flux(u, normal_direction::AbstractVector, - equations::CompressibleEulerMulticomponentEquations2D) - rho_v1, rho_v2, rho_e = u - - rho = density(u, equations) - - v1 = rho_v1 / rho - v2 = rho_v2 / rho - v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] - gamma = totalgamma(u, equations) - p = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2)) - - f_rho = densities(u, v_normal, equations) - f1 = rho_v1 * v_normal + p * normal_direction[1] - f2 = rho_v2 * v_normal + p * normal_direction[2] - f3 = (rho_e + p) * v_normal - - f_other = SVector{3, real(equations)}(f1, f2, f3) - - return vcat(f_other, f_rho) -end - """ flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations2D) @@ -469,76 +446,6 @@ See also return vcat(f_other, f_rho) end -@inline function flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressibleEulerMulticomponentEquations2D) - # Unpack left and right state - @unpack gammas, gas_constants, cv = equations - rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll - rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr - rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 3], - u_rr[i + 3]) - for i in eachcomponent(equations)) - rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 3] + - u_rr[i + 3]) - for i in eachcomponent(equations)) - - # Iterating over all partial densities - rho_ll = density(u_ll, equations) - rho_rr = density(u_rr, equations) - - # Calculating gamma - gamma = totalgamma(0.5 * (u_ll + u_rr), equations) - inv_gamma_minus_one = 1 / (gamma - 1) - - # extract velocities - v1_ll = rho_v1_ll / rho_ll - v1_rr = rho_v1_rr / rho_rr - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_ll = rho_v2_ll / rho_ll - v2_rr = rho_v2_rr / rho_rr - v2_avg = 0.5 * (v2_ll + v2_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr) - v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] - v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] - - # helpful variables - help1_ll = zero(v1_ll) - help1_rr = zero(v1_rr) - enth_ll = zero(v1_ll) - enth_rr = zero(v1_rr) - for i in eachcomponent(equations) - enth_ll += u_ll[i + 3] * gas_constants[i] - enth_rr += u_rr[i + 3] * gas_constants[i] - help1_ll += u_ll[i + 3] * cv[i] - help1_rr += u_rr[i + 3] * cv[i] - end - - # temperature and pressure - T_ll = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll - T_rr = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr - p_ll = T_ll * enth_ll - p_rr = T_rr * enth_rr - p_avg = 0.5 * (p_ll + p_rr) - inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - - f_rho_sum = zero(T_rr) - f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * 0.5 * - (v_dot_n_ll + v_dot_n_rr) - for i in eachcomponent(equations)) - for i in eachcomponent(equations) - f_rho_sum += f_rho[i] - end - f1 = f_rho_sum * v1_avg + p_avg * normal_direction[1] - f2 = f_rho_sum * v2_avg + p_avg * normal_direction[2] - f3 = f_rho_sum * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one) + - 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll) - - # momentum and energy flux - f_other = SVector{3, real(equations)}(f1, f2, f3) - - return vcat(f_other, f_rho) -end - # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D) @@ -569,32 +476,6 @@ end λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) end -@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressibleEulerMulticomponentEquations2D) - rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll - rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr - - # Get the density and gas gamma - rho_ll = density(u_ll, equations) - rho_rr = density(u_rr, equations) - gamma_ll = totalgamma(u_ll, equations) - gamma_rr = totalgamma(u_rr, equations) - - # Get the velocities based on direction - v_ll = rho_v1_ll / rho_ll * normal_direction[1] + - rho_v1_ll / rho_ll * normal_direction[2] - v_rr = rho_v1_rr / rho_rr * normal_direction[1] + - rho_v1_rr / rho_rr * normal_direction[2] - - # Compute the sound speeds on the left and right - p_ll = (gamma_ll - 1) * (rho_e_ll - 1 / 2 * (rho_v1_ll^2 + rho_v2_ll^2) / rho_ll) - c_ll = sqrt(gamma_ll * p_ll / rho_ll) - p_rr = (gamma_rr - 1) * (rho_e_rr - 1 / 2 * (rho_v1_rr^2 + rho_v2_rr^2) / rho_rr) - c_rr = sqrt(gamma_rr * p_rr / rho_rr) - - λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) -end - @inline function max_abs_speeds(u, equations::CompressibleEulerMulticomponentEquations2D) rho_v1, rho_v2, rho_e = u @@ -610,62 +491,6 @@ end return (abs(v1) + c, abs(v2) + c) end -@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerMulticomponentEquations2D) - rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll - rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr - - rho_ll = density(u_ll, equations) - rho_rr = density(u_rr, equations) - - v1_ll = rho_v1_ll / rho_ll - v2_ll = rho_v2_ll / rho_ll - v1_rr = rho_v1_rr / rho_rr - v2_rr = rho_v2_rr / rho_rr - - gamma = totalgamma(u_ll, equations) - p_ll = (gamma - 1) * (rho_e_ll - 1 / 2 * rho_ll * (v1_ll^2 + v2_ll^2)) - p_rr = (gamma - 1) * (rho_e_rr - 1 / 2 * rho_rr * (v1_rr^2 + v2_rr^2)) - - if orientation == 1 # x-direction - λ_min = v1_ll - sqrt(gamma * p_ll / rho_ll) - λ_max = v1_rr + sqrt(gamma * p_rr / rho_rr) - else # y-direction - λ_min = v2_ll - sqrt(gamma * p_ll / rho_ll) - λ_max = v2_rr + sqrt(gamma * p_rr / rho_rr) - end - - return λ_min, λ_max -end - -@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, - equations::CompressibleEulerMulticomponentEquations2D) - rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll - rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr - - rho_ll = density(u_ll, equations) - rho_rr = density(u_rr, equations) - - v1_ll = rho_v1_ll / rho_ll - v2_ll = rho_v2_ll / rho_ll - v1_rr = rho_v1_rr / rho_rr - v2_rr = rho_v2_rr / rho_rr - - gamma = totalgamma(u_ll, equations) - p_ll = (gamma - 1) * (rho_e_ll - 1 / 2 * rho_ll * (v1_ll^2 + v2_ll^2)) - p_rr = (gamma - 1) * (rho_e_rr - 1 / 2 * rho_rr * (v1_rr^2 + v2_rr^2)) - - v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] - v_normal_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] - - norm_ = norm(normal_direction) - # The v_normals are already scaled by the norm - λ_min = v_normal_ll - sqrt(gamma * p_ll / rho_ll) * norm_ - λ_max = v_normal_rr + sqrt(gamma * p_rr / rho_rr) * norm_ - - return λ_min, λ_max -end - # Convert conservative variables to primitive @inline function cons2prim(u, equations::CompressibleEulerMulticomponentEquations2D) rho_v1, rho_v2, rho_e = u From 8d26189fdd657f1ba06df7563a7f22e898efe390 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 20 Oct 2023 15:39:02 +0100 Subject: [PATCH 045/117] Renamed documentation file for multi-physics coupling. --- docs/src/{coupling.md => multi-physics_coupling.md} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename docs/src/{coupling.md => multi-physics_coupling.md} (100%) diff --git a/docs/src/coupling.md b/docs/src/multi-physics_coupling.md similarity index 100% rename from docs/src/coupling.md rename to docs/src/multi-physics_coupling.md From afbcc9da0f0944b6d17b38ac3af97704cb4464eb Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 20 Oct 2023 15:40:34 +0100 Subject: [PATCH 046/117] Renamed doc reference. --- docs/src/multi-physics_coupling.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/multi-physics_coupling.md b/docs/src/multi-physics_coupling.md index ca9a26ce360..da1a66d344f 100644 --- a/docs/src/multi-physics_coupling.md +++ b/docs/src/multi-physics_coupling.md @@ -1,4 +1,4 @@ -# [Coupling](@id coupling-id) +# [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. From e730c9aa7b07c3b0b68bb31ded0556b4367e1e9a Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Fri, 20 Oct 2023 15:41:18 +0100 Subject: [PATCH 047/117] Update src/semidiscretization/semidiscretization_coupled.jl Co-authored-by: Michael Schlottke-Lakemper --- src/semidiscretization/semidiscretization_coupled.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index d3a9003595c..342a8a8302d 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -365,7 +365,7 @@ BoundaryConditionCoupled(2, (:j, :i_backwards, :end), Float64) This is an experimental feature and can change any time. """ mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indices, - CouplingFunction} + 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 From 1cbed944d9541bf96e02e1029b4c8bc8ba49bb7e Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 24 Oct 2023 14:20:01 +0100 Subject: [PATCH 048/117] Update docs/src/multi-physics_coupling.md Co-authored-by: Michael Schlottke-Lakemper --- docs/src/multi-physics_coupling.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/multi-physics_coupling.md b/docs/src/multi-physics_coupling.md index da1a66d344f..9cf16176304 100644 --- a/docs/src/multi-physics_coupling.md +++ b/docs/src/multi-physics_coupling.md @@ -17,9 +17,9 @@ 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` +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` From fefbe5972587cfac6acdeb5d74ef2a848795e140 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 24 Oct 2023 14:20:13 +0100 Subject: [PATCH 049/117] Update docs/src/multi-physics_coupling.md Co-authored-by: Michael Schlottke-Lakemper --- docs/src/multi-physics_coupling.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/multi-physics_coupling.md b/docs/src/multi-physics_coupling.md index 9cf16176304..6d7fc254d62 100644 --- a/docs/src/multi-physics_coupling.md +++ b/docs/src/multi-physics_coupling.md @@ -22,7 +22,7 @@ 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` +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`. From f9ab6427e0a90bfc02f7b2842f3aedaff1bacd2e Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 24 Oct 2023 14:20:35 +0100 Subject: [PATCH 050/117] Update docs/src/multi-physics_coupling.md Co-authored-by: Michael Schlottke-Lakemper --- docs/src/multi-physics_coupling.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/multi-physics_coupling.md b/docs/src/multi-physics_coupling.md index 6d7fc254d62..4ca33e8da28 100644 --- a/docs/src/multi-physics_coupling.md +++ b/docs/src/multi-physics_coupling.md @@ -27,7 +27,7 @@ and return the transformed variables. Examples can be seen in `examples/structured_2d_dgsem/elixir_advection_coupled.jl`. -## Warning about Binary Compatibility +## 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. From b1f63dcd784dea3f7bbda265d91aaa0c9d8dcb6a Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 24 Oct 2023 14:22:59 +0100 Subject: [PATCH 051/117] Reinstated structured_2d_dgsem coupled in special tests. --- test/test_special_elixirs.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/test/test_special_elixirs.jl b/test/test_special_elixirs.jl index 4932b63a80a..c05dfbdfca1 100644 --- a/test/test_special_elixirs.jl +++ b/test/test_special_elixirs.jl @@ -30,6 +30,12 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none", @test isapprox(mean_convergence[:l2], [4.0], rtol=0.05) end + @timed_testset "structured_2d_dgsem coupled" begin + mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_coupled.jl"), 3) + @test isapprox(mean_convergence[1][:l2], [4.0], rtol=0.05) + @test isapprox(mean_convergence[2][:l2], [4.0], rtol=0.05) + end + @timed_testset "p4est_2d_dgsem" begin # Run convergence test on unrefined mesh no_refine = @cfunction((p4est, which_tree, quadrant) -> Cint(0), Cint, (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p4est_quadrant_t})) From 4c3254b74990db8eb6503ee639f2194362969dca Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 24 Oct 2023 14:26:11 +0100 Subject: [PATCH 052/117] Update examples/structured_2d_dgsem/elixir_advection_coupled.jl Co-authored-by: Michael Schlottke-Lakemper --- examples/structured_2d_dgsem/elixir_advection_coupled.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index 75bac2e759b..440506eddb3 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl @@ -112,7 +112,7 @@ 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, +save_solution = SaveSolutionCallback(interval = 100, solution_variables = cons2prim) # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step From 64e8dda071c672ca0b7d7ae61f6fdf0e19077a76 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 24 Oct 2023 14:47:55 +0100 Subject: [PATCH 053/117] Renamed CouplingFunction to CouplingConverter. --- src/semidiscretization/semidiscretization_coupled.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 342a8a8302d..cbabcd39728 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -372,7 +372,7 @@ mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indic other_semi_index::Int other_orientation::Int indices::Indices - coupling_converter::CouplingFunction + coupling_converter::CouplingConverter function BoundaryConditionCoupled(other_semi_index, indices, uEltype, coupling_converter) From bd59f724009bf38eb682e8e003278b02dfc007ca Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 24 Oct 2023 15:41:09 +0100 Subject: [PATCH 054/117] Update src/semidiscretization/semidiscretization_coupled.jl Co-authored-by: Michael Schlottke-Lakemper --- src/semidiscretization/semidiscretization_coupled.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 342a8a8302d..6e06610373c 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -498,6 +498,7 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ for cell in cells i_node = i_node_start j_node = j_node_start + element_id = linear_indices[i_cell, j_cell] for i in eachnode(solver) for v in 1:size(boundary_condition.u_boundary, 1) From 2331cc103e1cade1360a9e031f65190093c5f610 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 24 Oct 2023 16:55:32 +0100 Subject: [PATCH 055/117] Cleaned the copy of coupled boundary values. --- .../semidiscretization_coupled.jl | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index cbabcd39728..e997564f881 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -470,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) @@ -499,15 +501,12 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ i_node = i_node_start j_node = j_node_start - for i in eachnode(solver) + for element_id in eachnode(solver) for v in 1:size(boundary_condition.u_boundary, 1) - 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] + 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 end i_node += i_node_step j_node += j_node_step From 35bffba8b6985f06bc124c018d83afec8ad85bd3 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 24 Oct 2023 16:55:52 +0100 Subject: [PATCH 056/117] Reduced time span for example coupling elixir. --- examples/structured_2d_dgsem/elixir_advection_coupled.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index 75bac2e759b..5424d3a8e3a 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl @@ -100,7 +100,7 @@ semi = SemidiscretizationCoupled(semi1, semi2) # ODE solvers, callbacks etc. # Create ODE problem with time span from 0.0 to 20.0 -ode = semidiscretize(semi, (0.0, 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 From 488cef86ff304042887520b4bfe5c43ddcb40eba Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 24 Oct 2023 16:58:03 +0100 Subject: [PATCH 057/117] Removed redundant loop. --- src/semidiscretization/semidiscretization_coupled.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index e997564f881..ee92992fc89 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -502,12 +502,10 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ j_node = j_node_start for element_id in eachnode(solver) - for v in 1:size(boundary_condition.u_boundary, 1) - 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 - end + 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 From 81a10d01c20a2eb7d762342c80b83c1b2f9f9f6c Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 24 Oct 2023 16:58:34 +0100 Subject: [PATCH 058/117] Applied formatter. --- src/semidiscretization/semidiscretization_coupled.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index ee92992fc89..2975db19a6a 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -502,9 +502,10 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ j_node = j_node_start 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) + 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 From 0579dc127c1a54f84a00e7da290d372baff37c41 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 24 Oct 2023 17:04:49 +0100 Subject: [PATCH 059/117] Removed default coupling covnerter function. --- src/Trixi.jl | 2 -- .../coupling_converters.jl | 25 ------------------- 2 files changed, 27 deletions(-) delete mode 100644 src/coupling_converters/coupling_converters.jl diff --git a/src/Trixi.jl b/src/Trixi.jl index b0317438b99..3eb4ee3229b 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -200,8 +200,6 @@ export boundary_condition_do_nothing, BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal, BoundaryConditionCoupled -export coupling_converter_identity - export initial_condition_convergence_test, source_terms_convergence_test export source_terms_harmonic export initial_condition_poisson_nonperiodic, source_terms_poisson_nonperiodic, diff --git a/src/coupling_converters/coupling_converters.jl b/src/coupling_converters/coupling_converters.jl deleted file mode 100644 index 774986a503e..00000000000 --- a/src/coupling_converters/coupling_converters.jl +++ /dev/null @@ -1,25 +0,0 @@ -# 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. - -end # @muladd From baf28fdec6dd29fdef7c3826014672d8caccd471 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 24 Oct 2023 17:06:25 +0100 Subject: [PATCH 060/117] Moved coupling converter function into elixir. --- examples/structured_2d_dgsem/elixir_advection_coupled.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index 28fb5705ff5..5beedf58da8 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl @@ -50,7 +50,7 @@ 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 = coupling_converter_identity(equations) +coupling_function1 = (x, u) -> u boundary_conditions_x_neg1 = BoundaryConditionCoupled(2, (:end, :i_forward), Float64, coupling_function1) @@ -76,7 +76,7 @@ cells_per_dimension2 = cells_per_dimension mesh2 = StructuredMesh(cells_per_dimension2, coordinates_min2, coordinates_max2) -coupling_function2 = coupling_converter_identity(equations) +coupling_function2 = (x, u) -> u boundary_conditions_x_neg2 = BoundaryConditionCoupled(1, (:end, :i_forward), Float64, coupling_function2) From a4971a5aa8a4decb6f069c1115030ab0227ed3a1 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Wed, 25 Oct 2023 10:06:09 +0200 Subject: [PATCH 061/117] Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- docs/src/multi-physics_coupling.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/multi-physics_coupling.md b/docs/src/multi-physics_coupling.md index 4ca33e8da28..98a05422332 100644 --- a/docs/src/multi-physics_coupling.md +++ b/docs/src/multi-physics_coupling.md @@ -17,12 +17,12 @@ 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$ +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$ +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`. From 3491b1a4a81aef6e2349eeb45309a290698c8811 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Wed, 25 Oct 2023 10:07:13 +0200 Subject: [PATCH 062/117] Update docs/make.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index ac7013a5a82..259df9c1d0e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -107,7 +107,7 @@ makedocs( ], "Time integration" => "time_integration.md", "Callbacks" => "callbacks.md", - "Coupling" => "coupling.md" + "Coupling" => "multi-physics_coupling.md" ], "Advanced topics & developers" => [ "Conventions" =>"conventions.md", From 0ff77e2abd37aca6304677fc3f64469517fbe209 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 25 Oct 2023 09:44:05 +0100 Subject: [PATCH 063/117] Removed coupling_converters.jl from the include. --- src/Trixi.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 3eb4ee3229b..457d9dc336d 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -125,7 +125,6 @@ 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") From 19935605cdbf916e044ab1d90e6e8811a25021bc Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 25 Oct 2023 12:32:38 +0100 Subject: [PATCH 064/117] Corrected introduced issue with coupling boundary copy. The latest change to clean up the boundary copying introduced a bug related to the determination of the wrong node indices. This is now corrected. --- src/semidiscretization/semidiscretization_coupled.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 50c2f4af1ad..fefecd4c8d2 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -505,7 +505,7 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ 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) + u_node = u[:, i_node, j_node, linear_indices[i_cell, j_cell]] converted_u = coupling_converter(x, u_node) boundary_condition.u_boundary[:, element_id, cell] = converted_u i_node += i_node_step From d5129e249abbd9373e2c94b6fc15490fb3fd0c5d Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 25 Oct 2023 12:36:53 +0100 Subject: [PATCH 065/117] Corrected comment on final simulation time. --- examples/structured_2d_dgsem/elixir_advection_coupled.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index 5beedf58da8..bc2469392bf 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl @@ -99,7 +99,7 @@ semi = SemidiscretizationCoupled(semi1, semi2) ############################################################################### # ODE solvers, callbacks etc. -# Create ODE problem with time span from 0.0 to 20.0 +# Create ODE problem with time span from 0.0 to 2.0 ode = semidiscretize(semi, (0.0, 2.0)); # At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup From 21f203f3493d33eb38215103f794125e40ef8402 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 25 Oct 2023 12:43:25 +0100 Subject: [PATCH 066/117] Updated errors for coupled test to reflect changed final simulation time. --- test/test_structured_2d.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index ec6b23470af..c3192f52b2c 100644 --- a/test/test_structured_2d.jl +++ b/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 = [8.336466488044185e-6, 8.336466488028504e-6], - linf = [6.406463687191888e-5, 6.406463686947639e-5], + l2 = [7.816742843181738e-6, 7.816742843196112e-6], + linf = [6.314906965543265e-5, 6.314906965410039e-5], coverage_override = (maxiters=10^5,)) @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin errors = analysis_callback(sol) - @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 + @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 # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 4f92bc678f23bad738f342760d6ba010b6a76a9f Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 30 Oct 2023 11:39:48 +0000 Subject: [PATCH 067/117] Added miladd. --- src/semidiscretization/semidiscretization_coupled.jl | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index fefecd4c8d2..e5b7bf96c9f 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 @@ -150,8 +157,7 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) # 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) - + du_loc = get_system_u_ode(du_ode, i, semi) @trixi_timeit timer() "system #$i" rhs!(du_loc, u_loc, semi.semis[i], t) end @@ -617,3 +623,5 @@ function analyze_convergence(errors_coupled, iterations, return eoc_mean_values end + +end # @muladd \ No newline at end of file From 8789c4bc53ed31c82b934be83519deb0260ede69 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 3 Nov 2023 13:51:58 +0000 Subject: [PATCH 068/117] Corrected coordinate finding in semidiscretization_coupled. --- .../semidiscretization_coupled.jl | 41 +++++++++++++++---- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index e5b7bf96c9f..a986274f7a2 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -147,20 +147,36 @@ 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 - @trixi_timeit timer() "copy to coupled boundaries" begin - for semi_ in semi.semis - copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi) - end + for semi_ in semi.semis + copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi) 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) + du_loc = get_system_u_ode(du_ode, i, semi) + rhs!(du_loc, u_loc, semi.semis[i], t) @trixi_timeit timer() "system #$i" rhs!(du_loc, u_loc, semi.semis[i], t) end + # semi_ = semi.semis[1] + # copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi) + # semi_ = semi.semis[2] + # copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi) + + # Call rhs! for each semidiscretization + # i = 1 + # 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.semis[i], t) + # i = 2 + # 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.semis[i], t) + runtime = time_ns() - time_start put!(semi.performance_counter, runtime) @@ -469,6 +485,14 @@ function copy_to_coupled_boundary!(boundary_conditions::Union{Tuple, NamedTuple} for boundary_condition in boundary_conditions copy_to_coupled_boundary!(boundary_condition, u_ode, semi) end + # boundary_condition = boundary_conditions[1] + # copy_to_coupled_boundary!(boundary_condition, u_ode, semi) + # boundary_condition = boundary_conditions[2] + # copy_to_coupled_boundary!(boundary_condition, u_ode, semi) + # boundary_condition = boundary_conditions[3] + # copy_to_coupled_boundary!(boundary_condition, u_ode, semi) + # boundary_condition = boundary_conditions[4] + # copy_to_coupled_boundary!(boundary_condition, u_ode, semi) end # In 2D @@ -479,6 +503,7 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ @unpack coupling_converter, u_boundary = boundary_condition mesh, equations, solver, cache = mesh_equations_solver_cache(semi.semis[other_semi_index]) + # @unpack 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) @@ -509,14 +534,16 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ 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) + x = cache.elements.node_coordinates[:, i_node, j_node, + linear_indices[i_cell, j_cell]] u_node = u[:, i_node, j_node, linear_indices[i_cell, j_cell]] converted_u = coupling_converter(x, u_node) + # u_boundary[:, element_id, cell] = converted_u boundary_condition.u_boundary[:, element_id, cell] = converted_u i_node += i_node_step j_node += j_node_step end + i_cell += i_cell_step j_cell += j_cell_step end From 87609bcd1a50791da29fa423ec02725059488b9f Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 3 Nov 2023 16:06:39 +0000 Subject: [PATCH 069/117] Fixed issued related to memory allocation. --- .../semidiscretization_coupled.jl | 49 +++++-------------- 1 file changed, 13 insertions(+), 36 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index a986274f7a2..4a759199fec 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -143,6 +143,13 @@ end @view u_ode[semi.u_indices[index]] end +@inline function call_rhs!(u_ode, du_ode, i, t, semi::SemidiscretizationCoupled) + 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.semis[i], t) + # @trixi_timeit timer() "system #$i" rhs!(du_loc, u_loc, semi.semis[i], t) +end + function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) @unpack u_indices = semi @@ -150,32 +157,10 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) # @trixi_timeit timer() "copy to coupled boundaries" begin - for semi_ in semi.semis - copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi) - end + foreach(semi_ -> copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi), semi.semis) # 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) - rhs!(du_loc, u_loc, semi.semis[i], t) - @trixi_timeit timer() "system #$i" rhs!(du_loc, u_loc, semi.semis[i], t) - end - - # semi_ = semi.semis[1] - # copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi) - # semi_ = semi.semis[2] - # copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi) - - # Call rhs! for each semidiscretization - # i = 1 - # 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.semis[i], t) - # i = 2 - # 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.semis[i], t) + foreach(index -> call_rhs!(u_ode, du_ode, index, t, semi), nsystems(semi)) runtime = time_ns() - time_start put!(semi.performance_counter, runtime) @@ -482,17 +467,7 @@ 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) - end - # boundary_condition = boundary_conditions[1] - # copy_to_coupled_boundary!(boundary_condition, u_ode, semi) - # boundary_condition = boundary_conditions[2] - # copy_to_coupled_boundary!(boundary_condition, u_ode, semi) - # boundary_condition = boundary_conditions[3] - # copy_to_coupled_boundary!(boundary_condition, u_ode, semi) - # boundary_condition = boundary_conditions[4] - # copy_to_coupled_boundary!(boundary_condition, u_ode, semi) + foreach(boundary_condition -> copy_to_coupled_boundary!(boundary_condition, u_ode, semi), boundary_conditions) end # In 2D @@ -534,7 +509,9 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ element_id = linear_indices[i_cell, j_cell] for element_id in eachnode(solver) - x = cache.elements.node_coordinates[:, i_node, j_node, + # x = get_node_vars(node_coordinates, equations, solver, i_node, j_node, + # element_id) + x = @view cache.elements.node_coordinates[:, i_node, j_node, linear_indices[i_cell, j_cell]] u_node = u[:, i_node, j_node, linear_indices[i_cell, j_cell]] converted_u = coupling_converter(x, u_node) From 2906208eaf611614d8ab27ec21d3b6b8d25667ae Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 3 Nov 2023 17:49:19 +0000 Subject: [PATCH 070/117] Corrected loop over semidiscretization. --- src/semidiscretization/semidiscretization_coupled.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 4a759199fec..c64e3e00dfb 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -160,7 +160,7 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) foreach(semi_ -> copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi), semi.semis) # Call rhs! for each semidiscretization - foreach(index -> call_rhs!(u_ode, du_ode, index, t, semi), nsystems(semi)) + foreach(index -> call_rhs!(u_ode, du_ode, index, t, semi), eachsystem(semi)) runtime = time_ns() - time_start put!(semi.performance_counter, runtime) @@ -512,11 +512,11 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ # x = get_node_vars(node_coordinates, equations, solver, i_node, j_node, # element_id) x = @view cache.elements.node_coordinates[:, i_node, j_node, - linear_indices[i_cell, j_cell]] + linear_indices[i_cell, j_cell]] u_node = u[:, i_node, j_node, linear_indices[i_cell, j_cell]] converted_u = coupling_converter(x, u_node) - # u_boundary[:, element_id, cell] = converted_u - boundary_condition.u_boundary[:, element_id, cell] = converted_u + u_boundary[:, element_id, cell] = converted_u + # boundary_condition.u_boundary[:, element_id, cell] = converted_u i_node += i_node_step j_node += j_node_step end From 9ee224f16de0384460be1b8004ca3f588dd297d9 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 3 Nov 2023 21:03:51 +0000 Subject: [PATCH 071/117] Removed commented out code. --- src/semidiscretization/semidiscretization_coupled.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index c64e3e00dfb..6e0bd7681ec 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -478,7 +478,6 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ @unpack coupling_converter, u_boundary = boundary_condition mesh, equations, solver, cache = mesh_equations_solver_cache(semi.semis[other_semi_index]) - # @unpack 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) @@ -509,14 +508,11 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ 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) x = @view cache.elements.node_coordinates[:, i_node, j_node, linear_indices[i_cell, j_cell]] u_node = u[:, i_node, j_node, linear_indices[i_cell, j_cell]] converted_u = coupling_converter(x, u_node) u_boundary[:, element_id, cell] = converted_u - # boundary_condition.u_boundary[:, element_id, cell] = converted_u i_node += i_node_step j_node += j_node_step end From dd847e4370a9807651cda0d474af89cc5dfa5522 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 7 Nov 2023 23:39:50 +0000 Subject: [PATCH 072/117] Fixed type instability with loops over semidiscretizations using lispy tuple programming. --- .../semidiscretization_coupled.jl | 21 ++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 6e0bd7681ec..635b25dada7 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -143,11 +143,14 @@ end @view u_ode[semi.u_indices[index]] end -@inline function call_rhs!(u_ode, du_ode, i, t, semi::SemidiscretizationCoupled) +# Fix for the type instability. +@inline function call_rhs!(u_ode, du_ode, t, semi, i, semi_, semi_tuple...) 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.semis[i], t) - # @trixi_timeit timer() "system #$i" rhs!(du_loc, u_loc, semi.semis[i], t) + if length(semi_tuple) > 0 + call_rhs!(u_ode, du_ode, t, semi, i+1, semi.semis) + end end function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) @@ -160,7 +163,13 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) foreach(semi_ -> copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi), semi.semis) # Call rhs! for each semidiscretization - foreach(index -> call_rhs!(u_ode, du_ode, index, t, semi), eachsystem(semi)) + call_rhs!(u_ode, du_ode, t, semi, 1, semi.semis...) + + # 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) + # rhs!(du_loc, u_loc, semi.semis[i], t) + # end runtime = time_ns() - time_start put!(semi.performance_counter, runtime) @@ -470,6 +479,10 @@ function copy_to_coupled_boundary!(boundary_conditions::Union{Tuple, NamedTuple} foreach(boundary_condition -> copy_to_coupled_boundary!(boundary_condition, u_ode, semi), boundary_conditions) end +function mesh_equations_solver_cache_coupled!(semi_, other_semi_index, mesh, equations, solver, cache) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi.semis[other_semi_index]) +end + # In 2D function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{2}, u_ode, semi) @@ -477,6 +490,8 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ @unpack other_semi_index, other_orientation, indices = boundary_condition @unpack coupling_converter, u_boundary = boundary_condition + # foreach(index -> call_rhs!(u_ode, du_ode, index, t, semi), eachsystem(semi)) + # foreach(index -> mesh_equations_solver_cache_coupled!(index, _index, mesh, equations, solver, cache), eachsystem(semi)) 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, From a54a9e8ec9cad67448dbe2e74362e6b384b22ecd Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 8 Nov 2023 00:22:06 +0000 Subject: [PATCH 073/117] Removed obsolete code. --- src/semidiscretization/semidiscretization_coupled.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 635b25dada7..8337bf733ea 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -490,8 +490,6 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ @unpack other_semi_index, other_orientation, indices = boundary_condition @unpack coupling_converter, u_boundary = boundary_condition - # foreach(index -> call_rhs!(u_ode, du_ode, index, t, semi), eachsystem(semi)) - # foreach(index -> mesh_equations_solver_cache_coupled!(index, _index, mesh, equations, solver, cache), eachsystem(semi)) 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, From eb14606fd9a7e562534b2856080bd67350948827 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 8 Nov 2023 14:26:32 +0000 Subject: [PATCH 074/117] Fixed another typa instability in coupled semidiscretization. --- .../semidiscretization_coupled.jl | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 8337bf733ea..0da12b6430f 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -143,13 +143,14 @@ end @view u_ode[semi.u_indices[index]] end -# Fix for the type instability. -@inline function call_rhs!(u_ode, du_ode, t, semi, i, semi_, semi_tuple...) +# Rcursive call of the RHS for the semidiscretizations using Lispy tuple programming +# to avoid type instability. +@inline function rhs!(u_ode, du_ode, t, semi, i, semi_, semi_tuple...) 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.semis[i], t) if length(semi_tuple) > 0 - call_rhs!(u_ode, du_ode, t, semi, i+1, semi.semis) + rhs!(u_ode, du_ode, t, semi, i+1, semi.semis) end end @@ -163,7 +164,7 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) foreach(semi_ -> copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi), semi.semis) # Call rhs! for each semidiscretization - call_rhs!(u_ode, du_ode, t, semi, 1, semi.semis...) + rhs!(u_ode, du_ode, t, semi, 1, semi.semis...) # for i in eachsystem(semi) # u_loc = get_system_u_ode(u_ode, i, semi) @@ -479,8 +480,12 @@ function copy_to_coupled_boundary!(boundary_conditions::Union{Tuple, NamedTuple} foreach(boundary_condition -> copy_to_coupled_boundary!(boundary_condition, u_ode, semi), boundary_conditions) end -function mesh_equations_solver_cache_coupled!(semi_, other_semi_index, mesh, equations, solver, cache) - mesh, equations, solver, cache = mesh_equations_solver_cache(semi.semis[other_semi_index]) +function mesh_equations_solver_cache(other_semi_index, i, semi_, semi_tuple...) + if i == other_semi_index + return mesh_equations_solver_cache(semi_) + else + mesh_equations_solver_cache(other_semi_index, i+1, semi_tuple...) + end end # In 2D @@ -490,7 +495,7 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ @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]) + mesh, equations, solver, cache = mesh_equations_solver_cache(other_semi_index, 1, semi.semis...) @unpack node_coordinates = cache.elements u = wrap_array(get_system_u_ode(u_ode, other_semi_index, semi), mesh, equations, solver, cache) From 20a1676b3d76dd25e8a26fc515fe002feb98e6ba Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 8 Nov 2023 16:22:45 +0000 Subject: [PATCH 075/117] Cleaning up of the coupled semidiscretization. --- .../semidiscretization_coupled.jl | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 0da12b6430f..613e30db902 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -148,9 +148,9 @@ end @inline function rhs!(u_ode, du_ode, t, semi, i, semi_, semi_tuple...) 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.semis[i], t) - if length(semi_tuple) > 0 - rhs!(u_ode, du_ode, t, semi, i+1, semi.semis) + rhs!(du_loc, u_loc, semi_, t) + if i < nsystems(semi) + rhs!(u_ode, du_ode, t, semi, i+1, semi_tuple...) end end @@ -159,19 +159,11 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) time_start = time_ns() - # @trixi_timeit timer() "copy to coupled boundaries" begin - foreach(semi_ -> copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi), semi.semis) # Call rhs! for each semidiscretization rhs!(u_ode, du_ode, t, semi, 1, semi.semis...) - # 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) - # rhs!(du_loc, u_loc, semi.semis[i], t) - # end - runtime = time_ns() - time_start put!(semi.performance_counter, runtime) From 58ec8fbc9ac496910144c98c1b8a5b61b7aff3b9 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 8 Nov 2023 16:29:43 +0000 Subject: [PATCH 076/117] Autoformatted coupled semidiscretization. --- .../semidiscretization_coupled.jl | 63 ++++++++++++------- 1 file changed, 40 insertions(+), 23 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 613e30db902..4c77f57e88c 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -4,7 +4,7 @@ # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. @muladd begin #! format: noindent - + """ SemidiscretizationCoupled @@ -48,7 +48,8 @@ function SemidiscretizationCoupled(semis...) performance_counter = PerformanceCounter() - SemidiscretizationCoupled{typeof(semis), typeof(u_indices), typeof(performance_counter) + SemidiscretizationCoupled{typeof(semis), typeof(u_indices), + typeof(performance_counter) }(semis, u_indices, performance_counter) end @@ -71,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)) @@ -112,7 +115,9 @@ 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) @@ -150,7 +155,7 @@ end du_loc = get_system_u_ode(du_ode, i, semi) rhs!(du_loc, u_loc, semi_, t) if i < nsystems(semi) - rhs!(u_ode, du_ode, t, semi, i+1, semi_tuple...) + rhs!(u_ode, du_ode, t, semi, i + 1, semi_tuple...) end end @@ -158,8 +163,9 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) @unpack u_indices = semi time_start = time_ns() - - foreach(semi_ -> copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi), semi.semis) + + foreach(semi_ -> copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi), + semi.semis) # Call rhs! for each semidiscretization rhs!(u_ode, du_ode, t, semi, 1, semi.semis...) @@ -317,7 +323,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 @@ -409,8 +416,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, @@ -435,19 +444,22 @@ 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 # In 2D -function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2 +function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{ + 2 }, direction, mesh, equations, dg::DGSEM) if direction in (1, 2) @@ -469,27 +481,31 @@ end function copy_to_coupled_boundary!(boundary_conditions::Union{Tuple, NamedTuple}, u_ode, semi) - foreach(boundary_condition -> copy_to_coupled_boundary!(boundary_condition, u_ode, semi), boundary_conditions) + foreach(boundary_condition -> copy_to_coupled_boundary!(boundary_condition, u_ode, + semi), 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 - mesh_equations_solver_cache(other_semi_index, i+1, semi_tuple...) + 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, +function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{2}, + u_ode, 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(other_semi_index, 1, semi.semis...) + mesh, equations, solver, cache = mesh_equations_solver_cache(other_semi_index, 1, + semi.semis...) @unpack node_coordinates = cache.elements - u = wrap_array(get_system_u_ode(u_ode, other_semi_index, semi), mesh, equations, solver, + u = wrap_array(get_system_u_ode(u_ode, other_semi_index, semi), mesh, equations, + solver, cache) linear_indices = LinearIndices(size(mesh)) @@ -536,7 +552,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, @@ -556,7 +573,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). @@ -633,5 +651,4 @@ function analyze_convergence(errors_coupled, iterations, return eoc_mean_values end - -end # @muladd \ No newline at end of file +end # @muladd From 49b121e7cdd9979f42ba7457ed8ad0cbb0cdc42c Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 14 Nov 2023 11:26:25 +0000 Subject: [PATCH 077/117] Fixed last type instability in coupling. --- src/semidiscretization/semidiscretization_coupled.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 4c77f57e88c..66d831f0f73 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -479,10 +479,16 @@ function copy_to_coupled_boundary!(boundary_condition, u_ode, semi) return nothing end +function copy_to_coupled_boundary!(u_ode, semi, i, boundary_condition, boundary_conditions...) + copy_to_coupled_boundary!(boundary_condition, u_ode, semi) + if i < length(boundary_conditions) + copy_to_coupled_boundary!(u_ode, semi, i + 1, boundary_conditions...) + end +end + function copy_to_coupled_boundary!(boundary_conditions::Union{Tuple, NamedTuple}, u_ode, semi) - foreach(boundary_condition -> copy_to_coupled_boundary!(boundary_condition, u_ode, - semi), boundary_conditions) + copy_to_coupled_boundary!(u_ode, semi, 1, boundary_conditions...) end function mesh_equations_solver_cache(other_semi_index, i, semi_, semi_tuple...) From 3b5d49201ee82f3a1633733012f40ca724ceeb5a Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 14 Nov 2023 11:28:34 +0000 Subject: [PATCH 078/117] Autoformatter on semidiscretization. --- src/semidiscretization/semidiscretization_coupled.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 66d831f0f73..0d80e8195fb 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -479,7 +479,8 @@ function copy_to_coupled_boundary!(boundary_condition, u_ode, semi) return nothing end -function copy_to_coupled_boundary!(u_ode, semi, i, boundary_condition, boundary_conditions...) +function copy_to_coupled_boundary!(u_ode, semi, i, boundary_condition, + boundary_conditions...) copy_to_coupled_boundary!(boundary_condition, u_ode, semi) if i < length(boundary_conditions) copy_to_coupled_boundary!(u_ode, semi, i + 1, boundary_conditions...) From 538e40194473d4613b1f5d00d81692a75a389cc2 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 15 Nov 2023 17:08:52 +0000 Subject: [PATCH 079/117] Fixed bug in boundary values copy that arose when coupling multiple systems. --- src/semidiscretization/semidiscretization_coupled.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 0d80e8195fb..3455994224d 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -479,17 +479,17 @@ function copy_to_coupled_boundary!(boundary_condition, u_ode, semi) return nothing end -function copy_to_coupled_boundary!(u_ode, semi, i, boundary_condition, +function copy_to_coupled_boundary!(u_ode, semi, i, n_boundaries, boundary_condition, boundary_conditions...) copy_to_coupled_boundary!(boundary_condition, u_ode, semi) - if i < length(boundary_conditions) - copy_to_coupled_boundary!(u_ode, semi, i + 1, boundary_conditions...) + if i < n_boundaries + copy_to_coupled_boundary!(u_ode, semi, i + 1, n_boundaries, boundary_conditions...) end end function copy_to_coupled_boundary!(boundary_conditions::Union{Tuple, NamedTuple}, u_ode, semi) - copy_to_coupled_boundary!(u_ode, semi, 1, boundary_conditions...) + copy_to_coupled_boundary!(u_ode, semi, 1, length(boundary_conditions), boundary_conditions...) end function mesh_equations_solver_cache(other_semi_index, i, semi_, semi_tuple...) From e6b621a13bdfef1ecfb8d09366df52b9023aa046 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 15 Nov 2023 18:11:12 +0000 Subject: [PATCH 080/117] aplpied autoformatter on coupled semidiscretization. --- src/semidiscretization/semidiscretization_coupled.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 3455994224d..775b316bebc 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -483,13 +483,15 @@ function copy_to_coupled_boundary!(u_ode, semi, i, n_boundaries, boundary_condit boundary_conditions...) copy_to_coupled_boundary!(boundary_condition, u_ode, semi) if i < n_boundaries - copy_to_coupled_boundary!(u_ode, semi, i + 1, n_boundaries, boundary_conditions...) + copy_to_coupled_boundary!(u_ode, semi, i + 1, n_boundaries, + boundary_conditions...) end end function copy_to_coupled_boundary!(boundary_conditions::Union{Tuple, NamedTuple}, u_ode, semi) - copy_to_coupled_boundary!(u_ode, semi, 1, length(boundary_conditions), boundary_conditions...) + copy_to_coupled_boundary!(u_ode, semi, 1, length(boundary_conditions), + boundary_conditions...) end function mesh_equations_solver_cache(other_semi_index, i, semi_, semi_tuple...) From 2a7fea389a2f64d36183ff83fe76ce03d5d5f675 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 23 Nov 2023 12:28:15 +0000 Subject: [PATCH 081/117] Extended the structured 2d example elixir for the coupled advection to 4 semidiscretizations. This hase two purpuses: 1. Users are given an example fro 2d coupling avoiding common pitfalls. 2. This increases the code coverege for the test. --- .../elixir_advection_coupled.jl | 178 +++++++++++++----- 1 file changed, 130 insertions(+), 48 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index cfc2b0b0c63..b7e8edd9f32 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl @@ -2,33 +2,38 @@ 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. +# 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) @@ -38,62 +43,136 @@ 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)) +# These will be the number of elements for each quarter/semidiscretization. +cells_per_dimension = (8, 8) -# Define identical resolution as a variable such that it is easier to change from `trixi_include` -cells_per_dimension = (8, 16) +########### +# system #1 +########### -cells_per_dimension1 = cells_per_dimension +coordinates_min1 = (-1.0, 0.0) # minimum coordinates (min(x), min(y)) +coordinates_max1 = (0.0, 1.0) # maximum coordinates (max(x), max(y)) -mesh1 = StructuredMesh(cells_per_dimension1, coordinates_min1, coordinates_max1) -# The user can define their own coupling functions. -coupling_function1 = (x, u) -> u +mesh1 = StructuredMesh(cells_per_dimension, coordinates_min1, coordinates_max1) +# Define the coupling functions +coupling_function12 = (x, u) -> u +coupling_function13 = (x, u) -> u + +# Define the coupling boundary conditions and the system it is coupled to. boundary_conditions_x_neg1 = BoundaryConditionCoupled(2, (:end, :i_forward), Float64, - coupling_function1) + coupling_function12) boundary_conditions_x_pos1 = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64, - coupling_function1) + 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 = boundary_conditions_x_neg1, - # Connect right boundary with left boundary of right mesh + boundary_conditions = (x_neg = boundary_conditions_x_neg1, x_pos = boundary_conditions_x_pos1, - y_neg = boundary_condition_periodic, - y_pos = boundary_condition_periodic)) + y_neg = boundary_conditions_y_neg1, + y_pos = boundary_conditions_y_pos1)) -# 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)) +########### +# system #2 +########### -cells_per_dimension2 = cells_per_dimension +coordinates_min2 = (0.0, 0.0) # minimum coordinates (min(x), min(y)) +coordinates_max2 = (1.0, 1.0) # maximum coordinates (max(x), max(y)) -mesh2 = StructuredMesh(cells_per_dimension2, coordinates_min2, coordinates_max2) +mesh2 = StructuredMesh(cells_per_dimension, coordinates_min2, coordinates_max2) -coupling_function2 = (x, u) -> u +# Define the coupling functions +coupling_function21 = (x, u) -> u +coupling_function24 = (x, u) -> u +# Define the coupling boundary conditions and the system it is coupled to. boundary_conditions_x_neg2 = BoundaryConditionCoupled(1, (:end, :i_forward), Float64, - coupling_function2) + coupling_function21) boundary_conditions_x_pos2 = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64, - coupling_function2) + 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 = boundary_conditions_x_neg2, - # Connect right boundary with left boundary of left mesh + boundary_conditions = (x_neg = boundary_conditions_x_neg2, x_pos = boundary_conditions_x_pos2, - y_neg = boundary_condition_periodic, - y_pos = boundary_condition_periodic)) + 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) -> u +coupling_function31 = (x, u) -> 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 +########### + +# First mesh is the left half of a [-1,1]^2 square +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) -> u +coupling_function42 = (x, u) -> 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 semi1 and semi2 -semi = SemidiscretizationCoupled(semi1, semi2) +# Create a semidiscretization that bundles all the semidiscretizations. +semi = SemidiscretizationCoupled(semi1, semi2, semi3, semi4) ############################################################################### # ODE solvers, callbacks etc. @@ -108,7 +187,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, From e8c2d5277224efa36c152d99581b6081a47ddc72 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 23 Nov 2023 12:29:43 +0000 Subject: [PATCH 082/117] Updated test results for coupled advection in 2d to reflect the 4 semidiscretizations that are now used. --- test/test_structured_2d.jl | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 1addc29e3e6..f7fca31ba2c 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -33,8 +33,18 @@ 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 From 4f76579106d9d81b4f05e70bbc2a7ed159da1cae Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 24 Nov 2023 14:52:02 +0000 Subject: [PATCH 083/117] Added correct errors for tests for the coupled adveciton equations in structured 2d. --- test/test_structured_2d.jl | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index f7fca31ba2c..199e5b278f9 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -49,8 +49,18 @@ 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 + @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 From 418796a95dfdc5d8a0eea525ecb617674c93a7d7 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 11 Dec 2023 18:00:57 +0000 Subject: [PATCH 084/117] Update examples/structured_2d_dgsem/elixir_advection_coupled.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- examples/structured_2d_dgsem/elixir_advection_coupled.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index b7e8edd9f32..e471e52c659 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl @@ -43,7 +43,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) -# These will be the number of elements for each quarter/semidiscretization. +# This will be the number of elements for each quarter/semidiscretization. cells_per_dimension = (8, 8) ########### From fa4852bb939038abb04e8a9c5d9888839e8f3982 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 11 Dec 2023 18:01:34 +0000 Subject: [PATCH 085/117] Update examples/structured_2d_dgsem/elixir_advection_coupled.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- examples/structured_2d_dgsem/elixir_advection_coupled.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index e471e52c659..bde4b72b07f 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl @@ -143,7 +143,6 @@ semi3 = SemidiscretizationHyperbolic(mesh3, equations, initial_condition_converg # system #4 ########### -# First mesh is the left half of a [-1,1]^2 square coordinates_min4 = (0.0, -1.0) # minimum coordinates (min(x), min(y)) coordinates_max4 = (1.0, 0.0) # maximum coordinates (max(x), max(y)) From 9b55294b9ce144e0e9bd4fe555da8413733d1c47 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 11 Dec 2023 18:01:48 +0000 Subject: [PATCH 086/117] Update src/semidiscretization/semidiscretization_coupled.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- src/semidiscretization/semidiscretization_coupled.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index de5faeeefcf..958690a42c7 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -148,7 +148,7 @@ end @view u_ode[semi.u_indices[index]] end -# Rcursive call of the RHS for the semidiscretizations using Lispy tuple programming +# Recursive call of the RHS for the semidiscretizations using Lispy tuple programming # to avoid type instability. @inline function rhs!(u_ode, du_ode, t, semi, i, semi_, semi_tuple...) u_loc = get_system_u_ode(u_ode, i, semi) From c830e22facdb7cb29a6e11ca3a5b1a08fefe9932 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 11 Dec 2023 18:03:02 +0000 Subject: [PATCH 087/117] Update src/semidiscretization/semidiscretization_coupled.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- src/semidiscretization/semidiscretization_coupled.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 958690a42c7..f58cfbf4aec 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -168,7 +168,11 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) semi.semis) # Call rhs! for each semidiscretization - rhs!(u_ode, du_ode, t, semi, 1, semi.semis...) + 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 runtime = time_ns() - time_start put!(semi.performance_counter, runtime) From e064a56f198eb447d6c43fe2e3c76eb1acab6dd6 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 11 Dec 2023 18:04:05 +0000 Subject: [PATCH 088/117] Update src/semidiscretization/semidiscretization_coupled.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- src/semidiscretization/semidiscretization_coupled.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index f58cfbf4aec..c1cd09b65a8 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -389,10 +389,10 @@ mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indic # 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 - coupling_converter::CouplingConverter + other_semi_index :: Int + other_orientation :: Int + indices :: Indices + coupling_converter :: CouplingConverter function BoundaryConditionCoupled(other_semi_index, indices, uEltype, coupling_converter) From 9d8f4e54a816949dba698ee258a55c758de4d250 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 11 Dec 2023 18:04:28 +0000 Subject: [PATCH 089/117] Update src/semidiscretization/semidiscretization_coupled.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- src/semidiscretization/semidiscretization_coupled.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index c1cd09b65a8..38c1f9f34cb 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -407,11 +407,10 @@ mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indic other_orientation = 3 end - new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices), typeof(coupling_converter)}(u_boundary, - other_semi_index, - other_orientation, - indices, - coupling_converter) + new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices), + typeof(coupling_converter)}(u_boundary, + other_semi_index, other_orientation, + indices, coupling_converter) end end From 6c60c3ec4a4737fcabddc876e0043fa3e8c82b98 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 11 Dec 2023 18:07:35 +0000 Subject: [PATCH 090/117] Update src/semidiscretization/semidiscretization_coupled.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- src/semidiscretization/semidiscretization_coupled.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 38c1f9f34cb..2681e8ceb36 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -499,6 +499,7 @@ 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 From 388a30990ff50d7872e6440ce95b7438e66c997c Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 11 Dec 2023 18:09:05 +0000 Subject: [PATCH 091/117] Update src/semidiscretization/semidiscretization_coupled.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- src/semidiscretization/semidiscretization_coupled.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 2681e8ceb36..8c97f155169 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -545,9 +545,10 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ element_id = linear_indices[i_cell, j_cell] for element_id in eachnode(solver) - x = @view cache.elements.node_coordinates[:, i_node, j_node, - linear_indices[i_cell, j_cell]] - u_node = u[:, i_node, j_node, linear_indices[i_cell, j_cell]] + x = get_node_coords(node_coordinates, equations, solver, i_node, j_node, + linear_indices[i_cell, j_cell]) + u_node = get_node_vars(u, equations, solver, i_node, j_node, + linear_indices[i_cell, j_cell]] converted_u = coupling_converter(x, u_node) u_boundary[:, element_id, cell] = converted_u i_node += i_node_step From 0e7fd1960e227ec8c634908f2c3fd541ecaa6544 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 11 Dec 2023 18:09:30 +0000 Subject: [PATCH 092/117] Update src/semidiscretization/semidiscretization_coupled.jl Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- src/semidiscretization/semidiscretization_coupled.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 8c97f155169..f03fae7b574 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -550,7 +550,9 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ u_node = get_node_vars(u, equations, solver, i_node, j_node, linear_indices[i_cell, j_cell]] converted_u = coupling_converter(x, u_node) - u_boundary[:, element_id, cell] = converted_u + @inbounds for i in eachindex(converted_u) + u_boundary[i, element_id, cell] = converted_u[i] + end i_node += i_node_step j_node += j_node_step end From b4493a229c997181e26e9d0becf07be71102e055 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 11 Dec 2023 18:53:59 +0000 Subject: [PATCH 093/117] Corrected foreach_enumerate implementation. --- .../semidiscretization_coupled.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 38c1f9f34cb..956a9d9b62c 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -148,6 +148,22 @@ 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 + # Recursive call of the RHS for the semidiscretizations using Lispy tuple programming # to avoid type instability. @inline function rhs!(u_ode, du_ode, t, semi, i, semi_, semi_tuple...) From 34e6fcd0f933821d4095a22ac9a38c2268e282aa Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Dec 2023 20:41:51 +0100 Subject: [PATCH 094/117] Fix closing parens --- src/semidiscretization/semidiscretization_coupled.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index b4bef70746f..add68c71ee7 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -564,7 +564,7 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ x = get_node_coords(node_coordinates, equations, solver, i_node, j_node, linear_indices[i_cell, j_cell]) u_node = get_node_vars(u, equations, solver, i_node, j_node, - linear_indices[i_cell, j_cell]] + linear_indices[i_cell, j_cell]) converted_u = coupling_converter(x, u_node) @inbounds for i in eachindex(converted_u) u_boundary[i, element_id, cell] = converted_u[i] From 4c2b90b2a1b6afbd201674779136d9667fdfc152 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Dec 2023 20:43:45 +0100 Subject: [PATCH 095/117] Remove unused recursive rhs! --- src/semidiscretization/semidiscretization_coupled.jl | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index add68c71ee7..7b15531ebec 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -164,17 +164,6 @@ end foreach_enumerate(func, remaining_collection, index + 1) end -# Recursive call of the RHS for the semidiscretizations using Lispy tuple programming -# to avoid type instability. -@inline function rhs!(u_ode, du_ode, t, semi, i, semi_, semi_tuple...) - 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) - if i < nsystems(semi) - rhs!(u_ode, du_ode, t, semi, i + 1, semi_tuple...) - end -end - function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) @unpack u_indices = semi From f0911c05848945741ab999be1873e3c4bf6c57e6 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Dec 2023 22:22:25 +0100 Subject: [PATCH 096/117] Pass equations to converter function --- .../elixir_advection_coupled.jl | 16 ++-- .../semidiscretization_coupled.jl | 87 +++++++++---------- 2 files changed, 51 insertions(+), 52 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index bde4b72b07f..43b68f21b03 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl @@ -56,8 +56,8 @@ coordinates_max1 = (0.0, 1.0) # maximum coordinates (max(x), max(y)) mesh1 = StructuredMesh(cells_per_dimension, coordinates_min1, coordinates_max1) # Define the coupling functions -coupling_function12 = (x, u) -> u -coupling_function13 = (x, u) -> u +coupling_function12 = (x, u, equations_other, equations_own) -> u +coupling_function13 = (x, u, equations_other, equations_own) -> u # Define the coupling boundary conditions and the system it is coupled to. boundary_conditions_x_neg1 = BoundaryConditionCoupled(2, (:end, :i_forward), Float64, @@ -87,8 +87,8 @@ coordinates_max2 = (1.0, 1.0) # maximum coordinates (max(x), max(y)) mesh2 = StructuredMesh(cells_per_dimension, coordinates_min2, coordinates_max2) # Define the coupling functions -coupling_function21 = (x, u) -> u -coupling_function24 = (x, u) -> u +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, @@ -118,8 +118,8 @@ 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) -> u -coupling_function31 = (x, u) -> u +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, @@ -149,8 +149,8 @@ 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) -> u -coupling_function42 = (x, u) -> u +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, diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 7b15531ebec..6956a19c445 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -169,8 +169,9 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) time_start = time_ns() - foreach(semi_ -> copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi), - semi.semis) + foreach(semi.semis) do semi_ + copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi, semi_) + end # Call rhs! for each semidiscretization foreach_enumerate(semi.semis) do (i, semi_) @@ -356,7 +357,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 @@ -372,18 +373,20 @@ 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" @@ -481,65 +484,58 @@ 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, i, n_boundaries, boundary_condition, - boundary_conditions...) - copy_to_coupled_boundary!(boundary_condition, u_ode, semi) +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, i + 1, 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) - copy_to_coupled_boundary!(u_ode, semi, 1, length(boundary_conditions), + 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 + 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, equations, solver, cache = mesh_equations_solver_cache(other_semi_index, 1, - semi.semis...) - @unpack node_coordinates = cache.elements - u = wrap_array(get_system_u_ode(u_ode, other_semi_index, semi), mesh, equations, - solver, - cache) + mesh_own, equations_own, solver_own, cache_own = mesh_equations_solver_cache(semi) + + semi_other = semi_coupled.semis[other_semi_index] + mesh_other, equations_other, solver_other, cache_other = mesh_equations_solver_cache(semi_other) + + 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 @@ -549,15 +545,18 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ j_node = j_node_start element_id = linear_indices[i_cell, j_cell] - for element_id in eachnode(solver) - x = get_node_coords(node_coordinates, equations, solver, i_node, j_node, - linear_indices[i_cell, j_cell]) - u_node = get_node_vars(u, equations, solver, i_node, j_node, - linear_indices[i_cell, j_cell]) - converted_u = coupling_converter(x, u_node) - @inbounds for i in eachindex(converted_u) - u_boundary[i, element_id, cell] = converted_u[i] + 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 From dbf2aea61636be00821b0814b6091f339fc4e92e Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Dec 2023 22:54:26 +0100 Subject: [PATCH 097/117] Apply formatting --- src/semidiscretization/semidiscretization_coupled.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 6956a19c445..9786a0b6e30 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -396,7 +396,7 @@ mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indic 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 + u_boundary :: Array{uEltype, NDIMST2M1} # NDIMS * 2 - 1 other_semi_index :: Int other_orientation :: Int indices :: Indices @@ -518,7 +518,8 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ 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) + u_other = wrap_array(u_ode_other, mesh_other, equations_other, solver_other, + cache_other) linear_indices = LinearIndices(size(mesh_other)) @@ -546,11 +547,13 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ 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, + 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, + u_node_converted = coupling_converter(x_other, u_node_other, + equations_other, equations_own) for i in eachindex(u_node_converted) From 401e835783b4e3fa478a82a4efe230370b762588 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 9 Jan 2024 16:27:23 +0000 Subject: [PATCH 098/117] Reverted copy_to_coupled_boundary to previou version to avoid type instability. --- .../semidiscretization_coupled.jl | 51 ++++++++----------- 1 file changed, 20 insertions(+), 31 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 9786a0b6e30..7ebda74c79d 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -506,37 +506,34 @@ end # In 2D function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{2}, u_ode, - semi_coupled, semi) - @unpack u_indices = semi_coupled + semi) + @unpack u_indices = semi @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) - - semi_other = semi_coupled.semis[other_semi_index] - mesh_other, equations_other, solver_other, cache_other = mesh_equations_solver_cache(semi_other) + mesh, equations, solver, cache = mesh_equations_solver_cache(other_semi_index, 1, + semi.semis...) + @unpack node_coordinates = cache.elements + 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_other)) + linear_indices = LinearIndices(size(mesh)) if other_orientation == 1 - cells = axes(mesh_other, 2) + cells = axes(mesh, 2) else # other_orientation == 2 - cells = axes(mesh_other, 1) + cells = axes(mesh, 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_other) + node_index_range = eachnode(solver) 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_other, 1)) - j_cell_start, j_cell_step = index_to_start_step_2d(indices[2], axes(mesh_other, 2)) + 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 = i_cell_start j_cell = j_cell_start @@ -546,20 +543,12 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ j_node = j_node_start 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 - + for element_id in eachnode(solver) + x = @view cache.elements.node_coordinates[:, i_node, j_node, + linear_indices[i_cell, j_cell]] + u_node = u[:, i_node, j_node, linear_indices[i_cell, j_cell]] + converted_u = coupling_converter(x, u_node) + u_boundary[:, element_id, cell] = converted_u i_node += i_node_step j_node += j_node_step end From 6472d36b3fb72dcd855e87133ca29b5844b21b81 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 11 Jan 2024 15:30:44 +0000 Subject: [PATCH 099/117] Corrected computation of coupled semidiscretizations and fixed memory issue. --- .../semidiscretization_coupled.jl | 65 +++++++++++++------ 1 file changed, 44 insertions(+), 21 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 7ebda74c79d..ad7e2e26270 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -469,7 +469,9 @@ function allocate_coupled_boundary_condition(boundary_condition, direction, mesh end # In 2D -function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2}, +function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{ + 2 + }, direction, mesh, equations, dg::DGSEM) if direction in (1, 2) cell_size = size(mesh, 2) @@ -503,37 +505,50 @@ function copy_to_coupled_boundary!(boundary_conditions::Union{Tuple, NamedTuple} 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 + 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, equations, solver, cache = mesh_equations_solver_cache(other_semi_index, 1, - semi.semis...) - @unpack node_coordinates = cache.elements - u = wrap_array(get_system_u_ode(u_ode, other_semi_index, semi), mesh, equations, - solver, - cache) + 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...) + + 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 @@ -543,12 +558,20 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ j_node = j_node_start element_id = linear_indices[i_cell, j_cell] - for element_id in eachnode(solver) - x = @view cache.elements.node_coordinates[:, i_node, j_node, - linear_indices[i_cell, j_cell]] - u_node = u[:, i_node, j_node, linear_indices[i_cell, j_cell]] - converted_u = coupling_converter(x, u_node) - u_boundary[:, element_id, cell] = converted_u + 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 From 52f7645bcd382f512caa759b6f23716c6daf28a8 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 15 Jan 2024 09:57:18 +0000 Subject: [PATCH 100/117] Removed redundant nelements function, as it is no longer used. --- src/semidiscretization/semidiscretization_coupled.jl | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index ad7e2e26270..7bd16d7701b 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -123,14 +123,6 @@ end 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 From 9bfc7ccd29062828a176ae3ecb6e9d67e337308f Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 15 Jan 2024 15:19:58 +0000 Subject: [PATCH 101/117] Applied autoformatter. --- src/semidiscretization/semidiscretization_coupled.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 7bd16d7701b..db61d8027fe 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -461,9 +461,7 @@ function allocate_coupled_boundary_condition(boundary_condition, direction, mesh end # In 2D -function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{ - 2 - }, +function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2}, direction, mesh, equations, dg::DGSEM) if direction in (1, 2) cell_size = size(mesh, 2) From 6eccc51c5edcc16806ca2da111ce992aa0340ae1 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 29 Jan 2024 16:16:08 +0000 Subject: [PATCH 102/117] Added max_abs_speed_naive( and max_abs_speed_naive for PolytropicEulerEquations2D. --- src/equations/polytropic_euler_2d.jl | 42 ++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index f5d2f7b0bad..6cf29ee7f6c 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -301,6 +301,48 @@ end return abs(v1) + c, abs(v2) + c end +# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the +# maximum velocity magnitude plus the maximum speed of sound +@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::PolytropicEulerEquations2D) + rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations) + + # Get the velocity value in the appropriate direction + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + else # orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + end + # Calculate sound speeds + c_ll = sqrt(equations.gamma * equations.gamma * rho_ll^(equations.gamma - 1)) + c_rr = sqrt(equations.gamma * equations.gamma * rho_rr^(equations.gamma - 1)) + + λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) +end + +@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::PolytropicEulerEquations2D) + rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations) + + # Calculate normal velocities and sound speed + # left + v_ll = (v1_ll * normal_direction[1] + + + v2_ll * normal_direction[2]) + c_ll = sqrt(equations.gamma * equations.gamma * rho_ll^(equations.gamma - 1)) + # right + v_rr = (v1_rr * normal_direction[1] + + + v2_rr * normal_direction[2]) + c_rr = sqrt(equations.gamma * equations.gamma * rho_rr^(equations.gamma - 1)) + + return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction) +end + # Convert conservative variables to primitive @inline function cons2prim(u, equations::PolytropicEulerEquations2D) rho, rho_v1, rho_v2 = u From 673a51afb7bfcc69c798df64159e61713848a0fe Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 30 Jan 2024 11:38:37 +0000 Subject: [PATCH 103/117] Reverted coupling elixir to main branch version. The modified version should be part of a different PR. --- .../elixir_advection_coupled.jl | 191 +++++------------- 1 file changed, 53 insertions(+), 138 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index 43b68f21b03..2a56d23f4c0 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl @@ -2,38 +2,31 @@ using OrdinaryDiffEq using Trixi ############################################################################### -# 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. +# Coupled semidiscretization of two linear advection systems, which are connected periodically # -# 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: +# 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) # ┌────────────────────┬────────────────────┐ -# │ ↑ coupled ↑ │ ↑ coupled ↑ │ +# │ ↑ periodic ↑ │ ↑ periodic ↑ │ +# │ │ │ # │ │ │ # │ ========= │ ========= │ # │ system #1 │ system #2 │ # │ ========= │ ========= │ # │ │ │ -# │<-- coupled │<-- coupled │ -# │ coupled -->│ coupled -->│ # │ │ │ -# │ ↓ coupled ↓ │ ↓ coupled ↓ │ -# ├────────────────────┼────────────────────┤ -# │ ↑ coupled ↑ │ ↑ coupled ↑ │ # │ │ │ -# │ ========= │ ========= │ -# │ system #3 │ system #4 │ -# │ ========= │ ========= │ # │ │ │ -# │<-- coupled │<-- coupled │ -# │ coupled -->│ coupled -->│ +# │ coupled -->│<-- coupled │ +# │ │ │ +# │<-- coupled │ coupled -->│ +# │ │ │ # │ │ │ -# │ ↓ coupled ↓ │ ↓ coupled ↓ │ +# │ ↓ periodic ↓ │ ↓ periodic ↓ │ # └────────────────────┴────────────────────┘ # (-1, -1) ( 1, -1) @@ -43,135 +36,60 @@ 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) -# 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)) +# 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)) -mesh1 = StructuredMesh(cells_per_dimension, coordinates_min1, coordinates_max1) +# Define identical resolution as a variable such that it is easier to change from `trixi_include` +cells_per_dimension = (8, 16) -# Define the coupling functions -coupling_function12 = (x, u, equations_other, equations_own) -> u -coupling_function13 = (x, u, equations_other, equations_own) -> u +cells_per_dimension1 = cells_per_dimension -# 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) +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 = (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)) + 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)) coordinates_max2 = (1.0, 1.0) # maximum coordinates (max(x), max(y)) -mesh2 = StructuredMesh(cells_per_dimension, coordinates_min2, coordinates_max2) +cells_per_dimension2 = cells_per_dimension -# Define the coupling functions -coupling_function21 = (x, u, equations_other, equations_own) -> u -coupling_function24 = (x, u, equations_other, equations_own) -> u +mesh2 = StructuredMesh(cells_per_dimension2, coordinates_min2, coordinates_max2) -# 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 = (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) + 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) ############################################################################### # ODE solvers, callbacks etc. @@ -186,10 +104,7 @@ 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_callback3 = AnalysisCallback(semi3, interval = 100) -analysis_callback4 = AnalysisCallback(semi4, interval = 100) -analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2, - analysis_callback3, analysis_callback4) +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, From ea742858d22dbb9e63ca8f1098d029f35b306b3b Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 30 Jan 2024 11:39:24 +0000 Subject: [PATCH 104/117] Removed some modified coupling code as this should be part of a different PR. --- docs/make.jl | 1 - docs/src/multi-physics_coupling.md | 37 ------------------------------ 2 files changed, 38 deletions(-) delete mode 100644 docs/src/multi-physics_coupling.md diff --git a/docs/make.jl b/docs/make.jl index 259df9c1d0e..df8ac04be12 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -107,7 +107,6 @@ 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 deleted file mode 100644 index 98a05422332..00000000000 --- a/docs/src/multi-physics_coupling.md +++ /dev/null @@ -1,37 +0,0 @@ -# [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. - From 7975e9d9bac122831d026cd7a223af0ef06d891a Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 30 Jan 2024 11:40:38 +0000 Subject: [PATCH 105/117] Reverted changes on ooupling semidiscretization as this should be part of a different PR. --- .../semidiscretization_coupled.jl | 199 ++++++------------ 1 file changed, 68 insertions(+), 131 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index db61d8027fe..4fc77cb6a60 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -1,10 +1,3 @@ -# 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 @@ -48,9 +41,8 @@ function SemidiscretizationCoupled(semis...) performance_counter = PerformanceCounter() - SemidiscretizationCoupled{typeof(semis), typeof(u_indices), - typeof(performance_counter)}(semis, u_indices, - performance_counter) + SemidiscretizationCoupled{typeof(semis), typeof(u_indices), typeof(performance_counter) + }(semis, u_indices, performance_counter) end function Base.show(io::IO, semi::SemidiscretizationCoupled) @@ -72,13 +64,11 @@ 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)) @@ -115,14 +105,20 @@ end @inline Base.real(semi::SemidiscretizationCoupled) = promote_type(real.(semi.semis)...) -@inline function Base.eltype(semi::SemidiscretizationCoupled) - promote_type(eltype.(semi.semis)...) -end +@inline Base.eltype(semi::SemidiscretizationCoupled) = promote_type(eltype.(semi.semis)...) @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 @@ -140,36 +136,23 @@ 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() - foreach(semi.semis) do semi_ - copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi, semi_) + @trixi_timeit timer() "copy to coupled boundaries" begin + for semi_ in semi.semis + copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi) + end end # Call rhs! for each semidiscretization - foreach_enumerate(semi.semis) do (i, semi_) + 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) - rhs!(du_loc, u_loc, semi_, t) + + @trixi_timeit timer() "system #$i" rhs!(du_loc, u_loc, semi.semis[i], t) end runtime = time_ns() - time_start @@ -325,8 +308,7 @@ 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 @@ -349,7 +331,7 @@ end ################################################################################ """ - BoundaryConditionCoupled(other_semi_index, indices, uEltype, coupling_converter) + BoundaryConditionCoupled(other_semi_index, indices, uEltype) 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 @@ -365,37 +347,33 @@ 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, fun) +BoundaryConditionCoupled(2, (:begin, :i), Float64) # Connect the same two boundaries oppositely oriented -BoundaryConditionCoupled(2, (:begin, :i_backwards), Float64, fun) +BoundaryConditionCoupled(2, (:begin, :i_backwards), Float64) # 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, fun) +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, - CouplingConverter} +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 - coupling_converter :: CouplingConverter - - function BoundaryConditionCoupled(other_semi_index, indices, uEltype, - coupling_converter) + 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)) @@ -407,10 +385,8 @@ mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indic other_orientation = 3 end - new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices), - typeof(coupling_converter)}(u_boundary, - other_semi_index, other_orientation, - indices, coupling_converter) + new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices)}(u_boundary, other_semi_index, + other_orientation, indices, coupling_converter) end end @@ -419,10 +395,8 @@ 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, @@ -447,21 +421,20 @@ 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 # In 2D -function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2}, +function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2 + }, direction, mesh, equations, dg::DGSEM) if direction in (1, 2) cell_size = size(mesh, 2) @@ -476,69 +449,43 @@ 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_coupled, semi) +function copy_to_coupled_boundary!(boundary_condition, u_ode, 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_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...) + semi) + for boundary_condition in boundary_conditions + copy_to_coupled_boundary!(boundary_condition, u_ode, semi) end end # In 2D -function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{2}, - u_ode, - semi_coupled, semi) - @unpack u_indices = semi_coupled +function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{2}, u_ode, + semi) + @unpack u_indices = semi @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...) - 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) + 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) - linear_indices = LinearIndices(size(mesh_other)) + linear_indices = LinearIndices(size(mesh)) if other_orientation == 1 - cells = axes(mesh_other, 2) + cells = axes(mesh, 2) else # other_orientation == 2 - cells = axes(mesh_other, 1) + cells = axes(mesh, 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_other) + node_index_range = eachnode(solver) 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_other, 1)) - j_cell_start, j_cell_step = index_to_start_step_2d(indices[2], axes(mesh_other, 2)) + 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 = i_cell_start j_cell = j_cell_start @@ -546,26 +493,19 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ for cell in cells i_node = i_node_start j_node = j_node_start - 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 + for i in eachnode(solver) + for v in 1:size(u, 1) + 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 end - i_cell += i_cell_step j_cell += j_cell_step end @@ -575,8 +515,7 @@ 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, @@ -596,8 +535,7 @@ 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). @@ -674,4 +612,3 @@ function analyze_convergence(errors_coupled, iterations, return eoc_mean_values end -end # @muladd From 35f7c0dbfa277d69a6300f7e7892b4b128ae4995 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 30 Jan 2024 11:41:47 +0000 Subject: [PATCH 106/117] Reverted changes partaining the coupling PR. --- test/test_structured_2d.jl | 55 +++----------------------------------- 1 file changed, 4 insertions(+), 51 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 522510a42e3..1addc29e3e6 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -33,34 +33,14 @@ end @trixi_testset "elixir_advection_coupled.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_coupled.jl"), - 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, - ], + l2=[7.816742843181738e-6, 7.816742843196112e-6], + linf=[6.314906965543265e-5, 6.314906965410039e-5], coverage_override=(maxiters = 10^5,)) @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin errors = analysis_callback(sol) - @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 + @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 # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -631,33 +611,6 @@ end end end -@trixi_testset "elixir_euler_warm_bubble.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_warm_bubble.jl"), - l2=[ - 0.00019387402388722496, - 0.03086514388623955, - 0.04541427917165, - 43.892826583444716, - ], - linf=[ - 0.0015942305974430138, - 0.17449778969139373, - 0.3729704262394843, - 307.6706958565337, - ], - cells_per_dimension=(32, 16), - tspan=(0.0, 10.0)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 - end -end - @trixi_testset "elixir_eulerpolytropic_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"), l2=[ From f437ec95438fa765315bf8dd8fa567b4a4f74b43 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 30 Jan 2024 11:50:36 +0000 Subject: [PATCH 107/117] Removed changes partaining coupling PR. --- .../semidiscretization_coupled.jl | 22 ++++++++----------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 4fc77cb6a60..0941ae6a8ca 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -41,8 +41,9 @@ function SemidiscretizationCoupled(semis...) performance_counter = PerformanceCounter() - SemidiscretizationCoupled{typeof(semis), typeof(u_indices), typeof(performance_counter) - }(semis, u_indices, performance_counter) + SemidiscretizationCoupled{typeof(semis), typeof(u_indices), + typeof(performance_counter)}(semis, u_indices, + performance_counter) end function Base.show(io::IO, semi::SemidiscretizationCoupled) @@ -371,9 +372,8 @@ mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indic other_semi_index :: Int other_orientation :: Int indices :: Indices - coupling_converter :: Function - function BoundaryConditionCoupled(other_semi_index, indices, uEltype, coupling_converter) + function BoundaryConditionCoupled(other_semi_index, indices, uEltype) NDIMS = length(indices) u_boundary = Array{uEltype, NDIMS * 2 - 1}(undef, ntuple(_ -> 0, NDIMS * 2 - 1)) @@ -386,7 +386,7 @@ 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, coupling_converter) + other_orientation, indices) end end @@ -433,8 +433,7 @@ function allocate_coupled_boundary_condition(boundary_condition, direction, mesh end # In 2D -function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2 - }, +function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2}, direction, mesh, equations, dg::DGSEM) if direction in (1, 2) cell_size = size(mesh, 2) @@ -496,12 +495,9 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ for i in eachnode(solver) for v in 1:size(u, 1) - 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]] + 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 From b659ce232fc39a9c22108bb9ede9c845dd644106 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 30 Jan 2024 11:51:25 +0000 Subject: [PATCH 108/117] REverted to version including elixir_euler_warm_bubble.jl tests. --- test/test_structured_2d.jl | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 1addc29e3e6..e5d45ebcc07 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -611,6 +611,33 @@ end end end +@trixi_testset "elixir_euler_warm_bubble.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_warm_bubble.jl"), + l2=[ + 0.00019387402388722496, + 0.03086514388623955, + 0.04541427917165, + 43.892826583444716, + ], + linf=[ + 0.0015942305974430138, + 0.17449778969139373, + 0.3729704262394843, + 307.6706958565337, + ], + cells_per_dimension=(32, 16), + tspan=(0.0, 10.0)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + @trixi_testset "elixir_eulerpolytropic_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"), l2=[ From 8c5eb3bba4f4a337d328b4ae2b25e4d9877dc042 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 30 Jan 2024 11:55:56 +0000 Subject: [PATCH 109/117] Update src/equations/polytropic_euler_2d.jl Co-authored-by: Andrew Winters --- src/equations/polytropic_euler_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index 6cf29ee7f6c..62a09114d0c 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -317,8 +317,8 @@ end v_rr = v2_rr end # Calculate sound speeds - c_ll = sqrt(equations.gamma * equations.gamma * rho_ll^(equations.gamma - 1)) - c_rr = sqrt(equations.gamma * equations.gamma * rho_rr^(equations.gamma - 1)) + c_ll = sqrt(equations.gamma * equations.kappa * rho_ll^(equations.gamma - 1)) + c_rr = sqrt(equations.gamma * equations.kappa * rho_rr^(equations.gamma - 1)) λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) end From 66457c2bb9bec7ca0c80eb019bdc68f1faa34d5a Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 30 Jan 2024 11:56:31 +0000 Subject: [PATCH 110/117] Update src/equations/polytropic_euler_2d.jl Co-authored-by: Andrew Winters --- src/equations/polytropic_euler_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index 62a09114d0c..341963cb327 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -333,7 +333,7 @@ end v_ll = (v1_ll * normal_direction[1] + v2_ll * normal_direction[2]) - c_ll = sqrt(equations.gamma * equations.gamma * rho_ll^(equations.gamma - 1)) + c_ll = sqrt(equations.gamma * equations.kappa * rho_ll^(equations.gamma - 1)) # right v_rr = (v1_rr * normal_direction[1] + From b891b5bfe487739a874b7f790eff8cb82f84612d Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 30 Jan 2024 11:57:01 +0000 Subject: [PATCH 111/117] Update src/equations/polytropic_euler_2d.jl Co-authored-by: Andrew Winters --- src/equations/polytropic_euler_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index 341963cb327..2d050c5ac19 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -338,7 +338,7 @@ end v_rr = (v1_rr * normal_direction[1] + v2_rr * normal_direction[2]) - c_rr = sqrt(equations.gamma * equations.gamma * rho_rr^(equations.gamma - 1)) + c_rr = sqrt(equations.gamma * equations.kappa * rho_rr^(equations.gamma - 1)) return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction) end From 2007993d004934faac7c9aeadc9d1a1cd3fbd1db Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 30 Jan 2024 12:38:38 +0000 Subject: [PATCH 112/117] Update src/equations/polytropic_euler_2d.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Andrés Rueda-Ramírez --- src/equations/polytropic_euler_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index 2d050c5ac19..4d9b6c82b86 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -316,7 +316,7 @@ end v_ll = v2_ll v_rr = v2_rr end - # Calculate sound speeds + # Calculate sound speeds (we have p = kappa * rho^gamma) c_ll = sqrt(equations.gamma * equations.kappa * rho_ll^(equations.gamma - 1)) c_rr = sqrt(equations.gamma * equations.kappa * rho_rr^(equations.gamma - 1)) From a168e888378329256d51249f65f1d0ad555a761d Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 30 Jan 2024 12:39:02 +0000 Subject: [PATCH 113/117] Update src/equations/polytropic_euler_2d.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Andrés Rueda-Ramírez --- src/equations/polytropic_euler_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index 4d9b6c82b86..72e50c1b76c 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -328,7 +328,7 @@ end rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations) rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations) - # Calculate normal velocities and sound speed + # Calculate normal velocities and sound speed (we have p = kappa * rho^gamma) # left v_ll = (v1_ll * normal_direction[1] + From e0fd66e7256592b924004f891eff7a37be002761 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Wed, 31 Jan 2024 11:52:20 +0000 Subject: [PATCH 114/117] Update src/equations/polytropic_euler_2d.jl Co-authored-by: Daniel Doehring --- src/equations/polytropic_euler_2d.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index 72e50c1b76c..b3b888ba95d 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -330,14 +330,12 @@ end # Calculate normal velocities and sound speed (we have p = kappa * rho^gamma) # left - v_ll = (v1_ll * normal_direction[1] - + - v2_ll * normal_direction[2]) + v_ll = v1_ll * normal_direction[1] + + v2_ll * normal_direction[2] c_ll = sqrt(equations.gamma * equations.kappa * rho_ll^(equations.gamma - 1)) # right - v_rr = (v1_rr * normal_direction[1] - + - v2_rr * normal_direction[2]) + v_rr = v1_rr * normal_direction[1] + + v2_rr * normal_direction[2] c_rr = sqrt(equations.gamma * equations.kappa * rho_rr^(equations.gamma - 1)) return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction) From 35c8787d12d18bdbd78c6d0f115f4fd877541979 Mon Sep 17 00:00:00 2001 From: iomsn Date: Thu, 1 Feb 2024 15:37:15 +0000 Subject: [PATCH 115/117] Added consistency and rotation test for LAx-friedrich fluxes for polytropic equations in 2d. --- test/test_unit.jl | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/test/test_unit.jl b/test/test_unit.jl index e8a8effbe29..e68623f61d1 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -857,6 +857,30 @@ end end end +@timed_testset "Consistency check for Lax-Friedrich flux: Polytropic CEE" begin + for gamma in [1.4, 1.0, 5 / 3] + kappa = 0.5 # Scaling factor for the pressure. + equations = PolytropicEulerEquations2D(gamma, kappa) + u = SVector(1.1, -0.5, 2.34) + + orientations = [1, 2] + for orientation in orientations + @test flux_lax_friedrichs(u, u, orientation, equations) ≈ + flux(u, orientation, equations) + end + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions + @test flux_lax_friedrichs(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + end +end + @timed_testset "Consistency check for HLL flux with Davis wave speed estimates: LEE" begin flux_hll = FluxHLL(min_max_speed_davis) From debca0e52b9282404016f0facd48222889279acc Mon Sep 17 00:00:00 2001 From: iomsn Date: Mon, 5 Feb 2024 18:32:59 +0000 Subject: [PATCH 116/117] Applied auto-formatter on polytropic 2d equation. --- src/equations/polytropic_euler_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index b3b888ba95d..233b599cd20 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -334,7 +334,7 @@ end v2_ll * normal_direction[2] c_ll = sqrt(equations.gamma * equations.kappa * rho_ll^(equations.gamma - 1)) # right - v_rr = v1_rr * normal_direction[1] + + v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] c_rr = sqrt(equations.gamma * equations.kappa * rho_rr^(equations.gamma - 1)) From 4f5cc57f80a07457c25a6661f7f96694946ca821 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Thu, 15 Feb 2024 14:02:17 +0000 Subject: [PATCH 117/117] Update src/equations/polytropic_euler_2d.jl Co-authored-by: Andrew Winters --- src/equations/polytropic_euler_2d.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index 233b599cd20..e900fd64073 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -330,12 +330,12 @@ end # Calculate normal velocities and sound speed (we have p = kappa * rho^gamma) # left - v_ll = v1_ll * normal_direction[1] + - v2_ll * normal_direction[2] + v_ll = (v1_ll * normal_direction[1] + + v2_ll * normal_direction[2]) c_ll = sqrt(equations.gamma * equations.kappa * rho_ll^(equations.gamma - 1)) # right - v_rr = v1_rr * normal_direction[1] + - v2_rr * normal_direction[2] + v_rr = (v1_rr * normal_direction[1] + + v2_rr * normal_direction[2]) c_rr = sqrt(equations.gamma * equations.kappa * rho_rr^(equations.gamma - 1)) return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction)