diff --git a/examples/structured_2d_dgsem/elixir_euler_double_mach.jl b/examples/structured_2d_dgsem/elixir_euler_double_mach.jl index 2df2489d290..8153d2634b8 100644 --- a/examples/structured_2d_dgsem/elixir_euler_double_mach.jl +++ b/examples/structured_2d_dgsem/elixir_euler_double_mach.jl @@ -63,7 +63,12 @@ boundary_condition_inflow_outflow = BoundaryConditionCharacteristic(initial_cond direction, x, t, equations) # Calculate boundary flux - flux = surface_flux_function(u_boundary, u_inner, normal_direction, equations) + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, orientation_or_normal, + equations) + else + flux = surface_flux_function(u_boundary, u_inner, normal_direction, equations) + end else # x[1] >= 1 / 6 # Use the free slip wall BC otherwise flux = boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t, @@ -77,14 +82,14 @@ 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...) if x[1] < 1 / 6 # BoundaryConditionCharacteristic u_outer = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection, u_inner, - normal_direction / - norm(normal_direction), + normal_direction, direction, x, t, equations) else # if x[1] >= 1 / 6 # boundary_condition_slip_wall @@ -120,7 +125,7 @@ volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_fv = surface_flux) solver = DGSEM(basis, surface_flux, volume_integral) -initial_refinement_level = 6 +initial_refinement_level = 4 cells_per_dimension = (4 * 2^initial_refinement_level, 2^initial_refinement_level) coordinates_min = (0.0, 0.0) coordinates_max = (4.0, 1.0) 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 21f5f2dc3be..e85c263429f 100644 --- a/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl +++ b/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl @@ -63,7 +63,12 @@ boundary_condition_inflow_outflow = BoundaryConditionCharacteristic(initial_cond direction, x, t, equations) # Calculate boundary flux - flux = surface_flux_function(u_boundary, u_inner, normal_direction, equations) + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, orientation_or_normal, + equations) + else + flux = surface_flux_function(u_boundary, u_inner, normal_direction, equations) + end else # x[1] >= 1 / 6 # Use the free slip wall BC otherwise flux = boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t, @@ -77,14 +82,14 @@ 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...) if x[1] < 1 / 6 # BoundaryConditionCharacteristic u_outer = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection, u_inner, - normal_direction / - norm(normal_direction), + normal_direction, direction, x, t, equations) else # if x[1] >= 1 / 6 # boundary_condition_slip_wall @@ -123,7 +128,7 @@ volume_integral = VolumeIntegralSubcellLimiting(limiter_mcl; volume_flux_fv = surface_flux) solver = DGSEM(basis, surface_flux, volume_integral) -initial_refinement_level = 6 +initial_refinement_level = 4 cells_per_dimension = (4 * 2^initial_refinement_level, 2^initial_refinement_level) coordinates_min = (0.0, 0.0) coordinates_max = (4.0, 1.0) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 516f22ff33c..31d1c8a8184 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -226,9 +226,10 @@ end x, t, surface_flux_function, equations) - u_boundary = characteristic_boundary_value_function(boundary_condition.outer_boundary_value_function, - u_inner, orientation_or_normal, - direction, x, t, equations) + u_boundary = boundary_condition.boundary_value_function(boundary_condition.outer_boundary_value_function, + u_inner, + orientation_or_normal, + direction, x, t, equations) # Calculate boundary flux if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary diff --git a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl index d1044db1d44..c6f3dcc1f10 100644 --- a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl @@ -245,7 +245,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) @@ -270,8 +271,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) @@ -301,7 +302,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) @@ -326,8 +328,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 ce4b5440e86..aff68284b40 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -651,7 +651,7 @@ end return nothing end -@inline function calc_lambdas_bar_states!(u, t, mesh::TreeMesh, +@inline function calc_lambdas_bar_states!(u, t, mesh::TreeMesh{2}, nonconservative_terms, equations, limiter, dg, cache, boundary_conditions; calc_bar_states = true) @@ -761,7 +761,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) @@ -781,8 +782,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, @@ -807,7 +808,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) @@ -827,8 +829,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, @@ -1831,7 +1833,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 @@ -1839,7 +1841,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 @@ -1852,9 +1854,9 @@ end @inline function get_boundary_outer_state(u_inner, cache, t, boundary_condition::BoundaryConditionDirichlet, - orientation_or_normal, direction, equations, - dg, indices...) - @unpack node_coordinates = cache.elements + orientation_or_normal, direction, mesh, + equations, dg, indices...) + (; node_coordinates) = cache.elements x = get_node_coords(node_coordinates, equations, dg, indices...) u_outer = boundary_condition.boundary_value_function(x, t, equations) @@ -1864,30 +1866,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, - 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), + u_inner, orientation_or_normal, direction, x, t, equations) return u_outer diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index c37276c5802..8168ce6ab5f 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,