diff --git a/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl index 45ad75c8d6a..a7f157698bb 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl @@ -55,7 +55,7 @@ end dg::DG, indices...) x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...) - return initial_condition(x, t, equations) + return initial_condition_mach3_flow(x, t, equations) end # Supersonic outflow boundary condition. diff --git a/src/solvers/dgsem_p4est/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_p4est/dg_2d_subcell_limiters.jl index 7bf1e1a0438..ad36ba1bba2 100644 --- a/src/solvers/dgsem_p4est/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_p4est/dg_2d_subcell_limiters.jl @@ -86,34 +86,40 @@ flux_primary = flux(u_primary, normal_direction, equations) flux_secondary = flux(u_secondary, normal_direction, equations) - # TODO: Correct order? flux_primary - flux_secondary - bar_state = 0.5 * (u_primary + u_secondary) - - 0.5 * (flux_primary - flux_secondary) / lambda + bar_state_primary = 0.5 * (u_primary + u_secondary) - + 0.5 * (flux_secondary - flux_primary) / lambda + if isodd(secondary_direction) + bar_state_secondary = 0.5 * (u_primary + u_secondary) - + 0.5 * (flux_secondary - flux_primary) / lambda + else + bar_state_secondary = 0.5 * (u_primary + u_secondary) - + 0.5 * (flux_primary - flux_secondary) / lambda + end if primary_direction == 1 - set_node_vars!(bar_states1, bar_state, equations, dg, i_primary, - j_primary, primary_element) + set_node_vars!(bar_states1, bar_state_primary, equations, dg, + i_primary, j_primary, primary_element) elseif primary_direction == 2 - set_node_vars!(bar_states1, bar_state, equations, dg, i_primary + 1, - j_primary, primary_element) + set_node_vars!(bar_states1, bar_state_primary, equations, dg, + i_primary + 1, j_primary, primary_element) elseif primary_direction == 3 - set_node_vars!(bar_states2, bar_state, equations, dg, i_primary, - j_primary, primary_element) + set_node_vars!(bar_states2, bar_state_primary, equations, dg, + i_primary, j_primary, primary_element) else # primary_direction == 4 - set_node_vars!(bar_states2, bar_state, equations, dg, i_primary, - j_primary + 1, primary_element) + set_node_vars!(bar_states2, bar_state_primary, equations, dg, + i_primary, j_primary + 1, primary_element) end if secondary_direction == 1 - set_node_vars!(bar_states1, bar_state, equations, dg, i_secondary, - j_secondary, secondary_element) + set_node_vars!(bar_states1, bar_state_secondary, equations, dg, + i_secondary, j_secondary, secondary_element) elseif secondary_direction == 2 - set_node_vars!(bar_states1, bar_state, equations, dg, i_secondary + 1, - j_secondary, secondary_element) + set_node_vars!(bar_states1, bar_state_secondary, equations, dg, + i_secondary + 1, j_secondary, secondary_element) elseif secondary_direction == 3 - set_node_vars!(bar_states2, bar_state, equations, dg, i_secondary, - j_secondary, secondary_element) + set_node_vars!(bar_states2, bar_state_secondary, equations, dg, + i_secondary, j_secondary, secondary_element) else # secondary_direction == 4 - set_node_vars!(bar_states2, bar_state, equations, dg, i_secondary, - j_secondary + 1, secondary_element) + set_node_vars!(bar_states2, bar_state_secondary, equations, dg, + i_secondary, j_secondary + 1, secondary_element) end end end @@ -127,14 +133,14 @@ end @inline function calc_lambdas_bar_states_boundary!(u, t, limiter, boundary_conditions::BoundaryConditionPeriodic, equations, dg, cache; - calc_bar_states = calc_bar_states) + 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, equations, dg, cache; - calc_bar_states = calc_bar_states) + calc_bar_states = true) (; boundary_condition_types, boundary_indices) = boundary_conditions (; contravariant_vectors) = cache.elements @@ -189,7 +195,7 @@ end # TODO: Correct order? bar_state = 0.5 * (u_inner + u_outer) - - 0.5 * (flux_inner - flux_outer) / lambda + 0.5 * (flux_outer - flux_inner) / lambda if direction == 1 set_node_vars!(bar_states1, bar_state, equations, dg, i_node, j_node, element) diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index a5166e89251..582504598b4 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -1842,11 +1842,12 @@ end boundary_condition::typeof(boundary_condition_slip_wall), normal_direction::AbstractVector, direction, equations, dg, indices...) - u_rotate = rotate_to_x(u_inner, normal_direction, equations) + factor = (normal_direction[1] * u_inner[2] + normal_direction[2] * u_inner[3]) + u_normal = factor * SVector(normal_direction / norm(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 @@ -1864,13 +1865,13 @@ end @inline function get_boundary_outer_state(u_inner, cache, t, boundary_condition::BoundaryConditionCharacteristic, - orientation_or_normal, direction, equations, + orientation::Integer, direction, equations, dg::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_or_normal, + u_inner, orientation, direction, x, t, equations) return u_outer @@ -1878,14 +1879,16 @@ end @inline function get_boundary_outer_state(u_inner, cache, t, boundary_condition::BoundaryConditionCharacteristic, - normal_direction::AbstractVector, equations, - dg::DG, indices...) + normal_direction::AbstractVector, direction, + equations, dg::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, + u_inner, + normal_direction / + norm(normal_direction), direction, x, t, equations) return u_outer