From bb074193b9f8b1f66be706ef2a5db8cad7fa6b6b Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 26 Jun 2023 16:31:21 +0100 Subject: [PATCH 1/3] 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 2/3] 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 2ff14f5c2c293f9a75e477ec1c9b6078f055c5e3 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 31 Jan 2024 17:43:13 +0000 Subject: [PATCH 3/3] Added max_abs_speed_naive and max_abs_speed_naive for PolytropicEulerEquations2D. --- src/equations/polytropic_euler_2d.jl | 40 ++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index f5d2f7b0bad..b3b888ba95d 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -301,6 +301,46 @@ 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 (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)) + + λ_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 (we have p = kappa * rho^gamma) + # left + 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] + 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 + # Convert conservative variables to primitive @inline function cons2prim(u, equations::PolytropicEulerEquations2D) rho, rho_v1, rho_v2 = u