Skip to content

Commit

Permalink
Revise calculation of bar states
Browse files Browse the repository at this point in the history
  • Loading branch information
@bennibolm
bennibolm committed Nov 24, 2023
1 parent e2fd3a5 commit ba301ed
Show file tree
Hide file tree
Showing 3 changed files with 40 additions and 31 deletions.
2 changes: 1 addition & 1 deletion examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
50 changes: 28 additions & 22 deletions src/solvers/dgsem_p4est/dg_2d_subcell_limiters.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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)
Expand Down
19 changes: 11 additions & 8 deletions src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -1864,28 +1865,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::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
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
Expand Down

0 comments on commit ba301ed

Please sign in to comment.