From b45e92e8f75db64f81b3a6dce36b1962756a2f82 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Fri, 1 Dec 2023 19:09:56 +0100 Subject: [PATCH] Add mesh as parameter of `get_boundary_outer_state` - Fixes bug that P4estMesh calls the right `characteristic_boundary_state` routine --- .../elixir_euler_double_mach.jl | 1 + .../elixir_euler_double_mach_MCL.jl | 1 + src/equations/equations.jl | 7 ++-- .../dgsem_p4est/dg_2d_subcell_limiters.jl | 20 ++++++------ .../dgsem_p4est/subcell_limiters_2d.jl | 24 +++++++------- .../dg_2d_subcell_limiters.jl | 14 ++++---- .../dgsem_structured/subcell_limiters_2d.jl | 28 +++++++++------- .../dgsem_tree/dg_2d_subcell_limiters.jl | 32 +++++++++++-------- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 17 ++++++---- 9 files changed, 82 insertions(+), 62 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_double_mach.jl b/examples/structured_2d_dgsem/elixir_euler_double_mach.jl index 5e3117aa6a8..ef18aa6a68f 100644 --- a/examples/structured_2d_dgsem/elixir_euler_double_mach.jl +++ b/examples/structured_2d_dgsem/elixir_euler_double_mach.jl @@ -77,6 +77,7 @@ end @inline function Trixi.get_boundary_outer_state(u_inner, cache, t, boundary_condition::typeof(boundary_condition_mixed_characteristic_wall), normal_direction::AbstractVector, direction, + mesh::StructuredMesh{2}, equations::CompressibleEulerEquations2D, dg, indices...) x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...) 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 dc0af3a45da..afca1b7852e 100644 --- a/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl +++ b/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl @@ -77,6 +77,7 @@ end @inline function Trixi.get_boundary_outer_state(u_inner, cache, t, boundary_condition::typeof(boundary_condition_mixed_characteristic_wall), normal_direction::AbstractVector, direction, + mesh::StructuredMesh{2}, equations::CompressibleEulerEquations2D, dg, indices...) x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index b3a93a7382a..6d6e3e92a98 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -227,7 +227,8 @@ end surface_flux_function, equations) u_boundary = boundary_condition.boundary_value_function(boundary_condition.outer_boundary_value_function, - u_inner, orientation_or_normal, + u_inner, + orientation_or_normal, direction, x, t, equations) # Calculate boundary flux @@ -249,8 +250,8 @@ end surface_flux_function, equations) u_boundary = boundary_condition.boundary_value_function(boundary_condition.outer_boundary_value_function, - u_inner, normal_direction, x, t, - equations) + u_inner, normal_direction, + x, t, equations) # Calculate boundary flux flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations) diff --git a/src/solvers/dgsem_p4est/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_p4est/dg_2d_subcell_limiters.jl index 2ebd26dd845..9ad0f25721e 100644 --- a/src/solvers/dgsem_p4est/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_p4est/dg_2d_subcell_limiters.jl @@ -86,7 +86,8 @@ flux_primary = flux(u_primary, normal_direction, equations) flux_secondary = flux(u_secondary, normal_direction, equations) - bar_state = 0.5 * (u_primary + u_secondary) - 0.5 * (flux_secondary - flux_primary) / lambda + bar_state = 0.5 * (u_primary + u_secondary) - + 0.5 * (flux_secondary - flux_primary) / lambda if primary_direction == 1 set_node_vars!(bar_states1, bar_state, equations, dg, i_primary, j_primary, primary_element) @@ -116,23 +117,24 @@ end end - calc_lambdas_bar_states_boundary!(u, t, limiter, boundary_conditions, equations, dg, - cache; calc_bar_states = calc_bar_states) + calc_lambdas_bar_states_boundary!(u, t, limiter, boundary_conditions, + mesh, equations, dg, cache; + calc_bar_states = calc_bar_states) return nothing end @inline function calc_lambdas_bar_states_boundary!(u, t, limiter, boundary_conditions::BoundaryConditionPeriodic, - mesh::P4estMesh, equations, dg, cache; - calc_bar_states = true) + mesh::P4estMesh, equations, dg, + cache; calc_bar_states = true) return nothing end # Calc lambdas and bar states at physical boundaries @inline function calc_lambdas_bar_states_boundary!(u, t, limiter, boundary_conditions, - mesh::P4estMesh{2}, equations, dg, cache; - calc_bar_states = true) + mesh::P4estMesh{2}, equations, dg, + cache; calc_bar_states = true) (; boundary_condition_types, boundary_indices) = boundary_conditions (; contravariant_vectors) = cache.elements @@ -166,8 +168,8 @@ end u_inner = get_node_vars(u, equations, dg, i_node, j_node, element) u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_condition, normal_direction, - direction, equations, dg, i_node, - j_node, element) + direction, mesh, equations, dg, + i_node, j_node, element) lambda = max_abs_speed_naive(u_inner, u_outer, normal_direction, equations) diff --git a/src/solvers/dgsem_p4est/subcell_limiters_2d.jl b/src/solvers/dgsem_p4est/subcell_limiters_2d.jl index e196dc26a02..d014a7b0380 100644 --- a/src/solvers/dgsem_p4est/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_p4est/subcell_limiters_2d.jl @@ -68,20 +68,21 @@ function calc_bounds_twosided_interface!(var_min, var_max, variable, u, t, semi, end calc_bounds_twosided_interface_inner!(var_min, var_max, variable, u, t, - boundary_conditions, equations, dg, cache) + boundary_conditions, + mesh, equations, dg, cache) return nothing end @inline function calc_bounds_twosided_interface_inner!(var_min, var_max, variable, u, t, boundary_conditions::BoundaryConditionPeriodic, - equations, dg, cache) + mesh, equations, dg, cache) return nothing end @inline function calc_bounds_twosided_interface_inner!(var_min, var_max, variable, u, t, - boundary_conditions, equations, - dg, cache) + boundary_conditions, + mesh, equations, dg, cache) (; boundary_condition_types, boundary_indices) = boundary_conditions (; contravariant_vectors) = cache.elements @@ -110,8 +111,8 @@ end u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_condition, normal_direction, - direction, equations, dg, i_node, - j_node, element) + direction, mesh, equations, dg, + i_node, j_node, element) var_outer = u_outer[variable] var_min[i_node, j_node, element] = min(var_min[i_node, j_node, element], @@ -186,7 +187,8 @@ function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, sem end calc_bounds_onesided_interface_inner!(var_minmax, minmax, variable, u, t, - boundary_conditions, equations, dg, cache) + boundary_conditions, + mesh, equations, dg, cache) return nothing end @@ -194,13 +196,13 @@ end @inline function calc_bounds_onesided_interface_inner!(var_minmax, minmax, variable, u, t, boundary_conditions::BoundaryConditionPeriodic, - equations, dg, cache) + mesh, equations, dg, cache) return nothing end @inline function calc_bounds_onesided_interface_inner!(var_minmax, minmax, variable, u, t, boundary_conditions, - equations, dg, cache) + mesh, equations, dg, cache) (; boundary_condition_types, boundary_indices) = boundary_conditions (; contravariant_vectors) = cache.elements @@ -229,8 +231,8 @@ end u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_condition, normal_direction, - direction, equations, dg, i_node, - j_node, element) + direction, mesh, equations, dg, + i_node, j_node, element) var_outer = variable(u_outer, equations) var_minmax[i_node, j_node, element] = minmax(var_minmax[i_node, j_node, diff --git a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl index a485d79075e..3a1a5602141 100644 --- a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl @@ -246,7 +246,8 @@ end u_inner = get_node_vars(u, equations, dg, 1, j, element) u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[1], Ja1, 1, - equations, dg, 1, j, element) + mesh, equations, dg, + 1, j, element) lambda1[1, j, element] = max_abs_speed_naive(u_inner, u_outer, Ja1, equations) @@ -271,8 +272,8 @@ end u_inner = get_node_vars(u, equations, dg, nnodes(dg), j, element) u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[2], Ja1, 2, - equations, dg, nnodes(dg), j, - element) + mesh, equations, dg, + nnodes(dg), j, element) lambda1[nnodes(dg) + 1, j, element] = max_abs_speed_naive(u_inner, u_outer, Ja1, equations) @@ -302,7 +303,8 @@ end u_inner = get_node_vars(u, equations, dg, i, 1, element) u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[3], Ja2, 3, - equations, dg, i, 1, element) + mesh, equations, dg, + i, 1, element) lambda2[i, 1, element] = max_abs_speed_naive(u_inner, u_outer, Ja2, equations) @@ -327,8 +329,8 @@ end u_inner = get_node_vars(u, equations, dg, i, nnodes(dg), element) u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[4], Ja2, 4, - equations, dg, i, nnodes(dg), - element) + mesh, equations, dg, + i, nnodes(dg), element) lambda2[i, nnodes(dg) + 1, element] = max_abs_speed_naive(u_inner, u_outer, Ja2, equations) diff --git a/src/solvers/dgsem_structured/subcell_limiters_2d.jl b/src/solvers/dgsem_structured/subcell_limiters_2d.jl index c5c21c34ffb..83a979a16a8 100644 --- a/src/solvers/dgsem_structured/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_structured/subcell_limiters_2d.jl @@ -61,7 +61,8 @@ function calc_bounds_twosided_interface!(var_min, var_max, variable, u, t, semi, u_inner = get_node_vars(u, equations, dg, 1, j, element) u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[1], Ja1, 1, - equations, dg, 1, j, element) + mesh, equations, dg, + 1, j, element) var_outer = u_outer[variable] var_min[1, j, element] = min(var_min[1, j, element], var_outer) @@ -77,8 +78,8 @@ function calc_bounds_twosided_interface!(var_min, var_max, variable, u, t, semi, u_inner = get_node_vars(u, equations, dg, nnodes(dg), j, element) u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[2], Ja1, 2, - equations, dg, nnodes(dg), j, - element) + mesh, equations, dg, + nnodes(dg), j, element) var_outer = u_outer[variable] var_min[nnodes(dg), j, element] = min(var_min[nnodes(dg), j, element], @@ -97,7 +98,8 @@ function calc_bounds_twosided_interface!(var_min, var_max, variable, u, t, semi, u_inner = get_node_vars(u, equations, dg, i, 1, element) u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[3], Ja2, 3, - equations, dg, i, 1, element) + mesh, equations, dg, + i, 1, element) var_outer = u_outer[variable] var_min[i, 1, element] = min(var_min[i, 1, element], var_outer) @@ -113,8 +115,8 @@ function calc_bounds_twosided_interface!(var_min, var_max, variable, u, t, semi, u_inner = get_node_vars(u, equations, dg, i, nnodes(dg), element) u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[4], Ja2, 4, - equations, dg, i, nnodes(dg), - element) + mesh, equations, dg, + i, nnodes(dg), element) var_outer = u_outer[variable] var_min[i, nnodes(dg), element] = min(var_min[i, nnodes(dg), element], @@ -180,7 +182,8 @@ function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, sem u_inner = get_node_vars(u, equations, dg, 1, j, element) u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[1], Ja1, 1, - equations, dg, 1, j, element) + mesh, equations, dg, + 1, j, element) var_outer = variable(u_outer, equations) var_minmax[1, j, element] = minmax(var_minmax[1, j, element], var_outer) @@ -195,8 +198,8 @@ function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, sem u_inner = get_node_vars(u, equations, dg, nnodes(dg), j, element) u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[2], Ja1, 2, - equations, dg, nnodes(dg), j, - element) + mesh, equations, dg, + nnodes(dg), j, element) var_outer = variable(u_outer, equations) var_minmax[nnodes(dg), j, element] = minmax(var_minmax[nnodes(dg), j, @@ -214,7 +217,8 @@ function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, sem u_inner = get_node_vars(u, equations, dg, i, 1, element) u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[3], Ja2, 3, - equations, dg, i, 1, element) + mesh, equations, dg, + i, 1, element) var_outer = variable(u_outer, equations) var_minmax[i, 1, element] = minmax(var_minmax[i, 1, element], var_outer) @@ -229,8 +233,8 @@ function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, sem u_inner = get_node_vars(u, equations, dg, i, nnodes(dg), element) u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[4], Ja2, 4, - equations, dg, i, nnodes(dg), - element) + mesh, equations, dg, + i, nnodes(dg), element) var_outer = variable(u_outer, equations) var_minmax[i, nnodes(dg), element] = minmax(var_minmax[i, nnodes(dg), diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 15c5bbc3daf..3b40def2a97 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -764,7 +764,8 @@ end u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[1], orientation, 1, - equations, dg, 1, j, element) + mesh, equations, dg, + 1, j, element) lambda1[1, j, element] = max_abs_speed_naive(u_inner, u_outer, orientation, equations) @@ -784,8 +785,8 @@ end u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[2], orientation, 2, - equations, dg, nnodes(dg), j, - element) + mesh, equations, dg, + nnodes(dg), j, element) lambda1[nnodes(dg) + 1, j, element] = max_abs_speed_naive(u_inner, u_outer, orientation, @@ -810,7 +811,8 @@ end u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[3], orientation, 3, - equations, dg, i, 1, element) + mesh, equations, dg, + i, 1, element) lambda2[i, 1, element] = max_abs_speed_naive(u_inner, u_outer, orientation, equations) @@ -830,8 +832,8 @@ end u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[4], orientation, 4, - equations, dg, i, nnodes(dg), - element) + mesh, equations, dg, + i, nnodes(dg), element) lambda2[i, nnodes(dg) + 1, element] = max_abs_speed_naive(u_inner, u_outer, orientation, @@ -1834,7 +1836,7 @@ end @inline function get_boundary_outer_state(u_inner, cache, t, boundary_condition::typeof(boundary_condition_slip_wall), orientation::Integer, direction, - equations::CompressibleEulerEquations2D, + mesh, equations::CompressibleEulerEquations2D, dg, indices...) return SVector(u_inner[1], -u_inner[2], -u_inner[3], u_inner[4]) end @@ -1842,7 +1844,7 @@ end @inline function get_boundary_outer_state(u_inner, cache, t, boundary_condition::typeof(boundary_condition_slip_wall), normal_direction::AbstractVector, direction, - equations::CompressibleEulerEquations2D, + mesh, 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 @@ -1855,8 +1857,8 @@ end @inline function get_boundary_outer_state(u_inner, cache, t, boundary_condition::BoundaryConditionDirichlet, - orientation_or_normal, direction, equations, - dg, indices...) + orientation_or_normal, direction, mesh, + equations, dg, indices...) @unpack node_coordinates = cache.elements x = get_node_coords(node_coordinates, equations, dg, indices...) @@ -1867,13 +1869,15 @@ end @inline function get_boundary_outer_state(u_inner, cache, t, boundary_condition::BoundaryConditionCharacteristic, - orientation::Integer, direction, equations, + orientation_or_normal, direction, + mesh::Union{TreeMesh, StructuredMesh}, + 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, orientation, + u_inner, orientation_or_normal, direction, x, t, equations) return u_outer @@ -1882,7 +1886,7 @@ end @inline function get_boundary_outer_state(u_inner, cache, t, boundary_condition::BoundaryConditionCharacteristic, normal_direction::AbstractVector, direction, - equations, dg, indices...) + mesh::P4estMesh, equations, dg, indices...) (; node_coordinates) = cache.elements x = get_node_coords(node_coordinates, equations, dg, indices...) @@ -1890,7 +1894,7 @@ end u_outer = boundary_condition.boundary_value_function(boundary_condition.outer_boundary_value_function, u_inner, normal_direction, - direction, x, t, equations) + x, t, equations) return u_outer end diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 1e9d63d4602..52009113ad4 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -182,7 +182,8 @@ end u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[boundary_index], orientation, boundary_index, - equations, dg, index..., element) + mesh, equations, dg, + index..., element) var_outer = u_outer[variable] var_min[index..., element] = min(var_min[index..., element], var_outer) @@ -280,7 +281,8 @@ end u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[1], orientation, 1, - equations, dg, 1, j, element) + mesh, equations, dg, + 1, j, element) var_outer = variable(u_outer, equations) var_minmax[1, j, element] = minmax(var_minmax[1, j, element], @@ -292,8 +294,8 @@ end u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[2], orientation, 2, - equations, dg, nnodes(dg), j, - element) + mesh, equations, dg, + nnodes(dg), j, element) var_outer = variable(u_outer, equations) var_minmax[nnodes(dg), j, element] = minmax(var_minmax[nnodes(dg), @@ -308,7 +310,8 @@ end u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[3], orientation, 3, - equations, dg, i, 1, element) + mesh, equations, dg, + i, 1, element) var_outer = variable(u_outer, equations) var_minmax[i, 1, element] = minmax(var_minmax[i, 1, element], @@ -320,8 +323,8 @@ end u_outer = get_boundary_outer_state(u_inner, cache, t, boundary_conditions[4], orientation, 4, - equations, dg, i, nnodes(dg), - element) + mesh, equations, dg, + i, nnodes(dg), element) var_outer = variable(u_outer, equations) var_minmax[i, nnodes(dg), element] = minmax(var_minmax[i,