From bcfc1e952c5e408e55c7c394b966afcc234e11a4 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 28 Nov 2023 17:22:43 +0100 Subject: [PATCH] Fix normalization for Charactersitics and slip wall --- .../elixir_euler_double_mach.jl | 15 ++++--- .../elixir_euler_double_mach_MCL.jl | 15 ++++--- .../dg_2d_subcell_limiters.jl | 2 +- .../dgsem_tree/dg_2d_subcell_limiters.jl | 45 ++++++++++++------- src/time_integration/methods_SSP.jl | 4 +- test/test_structured_2d.jl | 32 ++++++------- 6 files changed, 64 insertions(+), 49 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_double_mach.jl b/examples/structured_2d_dgsem/elixir_euler_double_mach.jl index 1b268422e1f..d854ab45f76 100644 --- a/examples/structured_2d_dgsem/elixir_euler_double_mach.jl +++ b/examples/structured_2d_dgsem/elixir_euler_double_mach.jl @@ -76,22 +76,23 @@ end @inline function Trixi.get_boundary_outer_state(u_inner, cache, t, boundary_condition::typeof(boundary_condition_mixed_characteristic_wall), normal_direction::AbstractVector, direction, - equations, + equations::CompressibleEulerEquations2D, dg, indices...) x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...) if x[1] < 1 / 6 # BoundaryConditionCharacteristic u_outer = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection, u_inner, - normal_direction, - direction, x, t, - equations) + normal_direction / + norm(normal_direction), + direction, x, t, equations) else # if x[1] >= 1 / 6 # boundary_condition_slip_wall - u_rotate = Trixi.rotate_to_x(u_inner, normal_direction, equations) + factor = (normal_direction[1] * u_inner[2] + normal_direction[2] * u_inner[3]) + u_normal = (factor / sum(normal_direction .^ 2)) * normal_direction u_outer = SVector(u_inner[1], - u_inner[2] - 2.0 * u_rotate[2], - u_inner[3] - 2.0 * u_rotate[3], + u_inner[2] - 2.0 * u_normal[1], + u_inner[3] - 2.0 * u_normal[2], u_inner[4]) end return u_outer diff --git a/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl b/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl index ad537b49684..f5abd10ca87 100644 --- a/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl +++ b/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl @@ -76,22 +76,23 @@ end @inline function Trixi.get_boundary_outer_state(u_inner, cache, t, boundary_condition::typeof(boundary_condition_mixed_characteristic_wall), normal_direction::AbstractVector, direction, - equations, + equations::CompressibleEulerEquations2D, dg, indices...) x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...) if x[1] < 1 / 6 # BoundaryConditionCharacteristic u_outer = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection, u_inner, - normal_direction, - direction, x, t, - equations) + normal_direction / + norm(normal_direction), + direction, x, t, equations) else # if x[1] >= 1 / 6 # boundary_condition_slip_wall - u_rotate = Trixi.rotate_to_x(u_inner, normal_direction, equations) + factor = (normal_direction[1] * u_inner[2] + normal_direction[2] * u_inner[3]) + u_normal = (factor / sum(normal_direction .^ 2)) * normal_direction u_outer = SVector(u_inner[1], - u_inner[2] - 2.0 * u_rotate[2], - u_inner[3] - 2.0 * u_rotate[3], + u_inner[2] - 2.0 * u_normal[1], + u_inner[3] - 2.0 * u_normal[2], u_inner[4]) end return u_outer diff --git a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl index 4db5d5a3bf9..d1044db1d44 100644 --- a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl @@ -109,7 +109,7 @@ return nothing end -@inline function calc_lambdas_bar_states!(u, t, mesh::StructuredMesh, +@inline function calc_lambdas_bar_states!(u, t, mesh::StructuredMesh{2}, nonconservative_terms, equations, limiter, dg, cache, boundary_conditions; calc_bar_states = true) diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index fe751a26881..ce4b5440e86 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -1831,29 +1831,25 @@ end @inline function get_boundary_outer_state(u_inner, cache, t, boundary_condition::typeof(boundary_condition_slip_wall), orientation::Integer, direction, - equations, dg, indices...) + equations::CompressibleEulerEquations2D, + dg, indices...) return SVector(u_inner[1], -u_inner[2], -u_inner[3], u_inner[4]) end @inline function get_boundary_outer_state(u_inner, cache, t, boundary_condition::typeof(boundary_condition_slip_wall), - normal_direction::AbstractVector, - direction, equations, dg, indices...) - u_rotate = rotate_to_x(u_inner, normal_direction, equations) + normal_direction::AbstractVector, direction, + equations::CompressibleEulerEquations2D, + dg, indices...) + factor = (normal_direction[1] * u_inner[2] + normal_direction[2] * u_inner[3]) + u_normal = (factor / sum(normal_direction .^ 2)) * normal_direction return SVector(u_inner[1], - u_inner[2] - 2.0 * u_rotate[2], - u_inner[3] - 2.0 * u_rotate[3], + u_inner[2] - 2.0 * u_normal[1], + u_inner[3] - 2.0 * u_normal[2], u_inner[4]) end -# Default implementation of `get_boundary_outer_state` returns inner value. -@inline function get_boundary_outer_state(u_inner, cache, t, boundary_condition, - orientation_or_normal, direction, equations, - dg, indices...) - return u_inner -end - @inline function get_boundary_outer_state(u_inner, cache, t, boundary_condition::BoundaryConditionDirichlet, orientation_or_normal, direction, equations, @@ -1868,13 +1864,30 @@ end @inline function get_boundary_outer_state(u_inner, cache, t, boundary_condition::BoundaryConditionCharacteristic, - orientation_or_normal, direction, equations, + orientation::Integer, direction, equations, dg, indices...) - @unpack node_coordinates = cache.elements + (; node_coordinates) = cache.elements x = get_node_coords(node_coordinates, equations, dg, indices...) u_outer = boundary_condition.boundary_value_function(boundary_condition.outer_boundary_value_function, - u_inner, orientation_or_normal, + u_inner, orientation, + direction, x, t, equations) + + return u_outer +end + +@inline function get_boundary_outer_state(u_inner, cache, t, + boundary_condition::BoundaryConditionCharacteristic, + normal_direction::AbstractVector, direction, + equations, dg, indices...) + (; node_coordinates) = cache.elements + + x = get_node_coords(node_coordinates, equations, dg, indices...) + + u_outer = boundary_condition.boundary_value_function(boundary_condition.outer_boundary_value_function, + u_inner, + normal_direction / + norm(normal_direction), direction, x, t, equations) return u_outer diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index 7443e6ed262..bc5e46e1931 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -269,8 +269,8 @@ function calc_normal_directions!(container_bar_states, mesh::TreeMesh, equations nothing end -function calc_normal_directions!(container_bar_states, mesh::StructuredMesh, equations, - dg, cache) +function calc_normal_directions!(container_bar_states, mesh::StructuredMesh{2}, + equations, dg, cache) (; weights, derivative_matrix) = dg.basis (; contravariant_vectors) = cache.elements diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 50f032e9e0a..8061d62705a 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -695,16 +695,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shock_upstream_sc_subcell.jl"), l2=[ - 1.2351650580585378, - 1.1268555757641623, - 1.723421391322672, - 11.715156856828612, + 1.2360964282793594, + 1.1280398030799386, + 1.7267648591561005, + 11.7190241252948, ], linf=[ - 5.378363125399197, - 6.5624165241817245, - 10.008377369949162, - 51.24533429408666, + 5.393863681425582, + 6.476830142076269, + 10.13110957325047, + 49.922573759404, ], tspan=(0.0, 0.5)) # Ensure that we do not have excessive memory allocations @@ -720,16 +720,16 @@ end @trixi_testset "elixir_euler_shock_upstream_MCL.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shock_upstream_MCL.jl"), l2=[ - 1.2607430289878367, - 1.1565837325291748, - 1.779179030245946, - 11.891223800389836, + 1.2655219565562812, + 1.1981291403309833, + 1.7964926709729876, + 11.968876601483794, ], linf=[ - 5.68876088478026, - 8.16555442595031, - 10.859100194836678, - 50.258224089906975, + 5.675285252773072, + 10.063582862371597, + 12.3837111905078, + 51.34700027665198, ], tspan=(0.0, 0.5)) # Ensure that we do not have excessive memory allocations