Skip to content

Commit

Permalink
Clean up routine get_boundary_outer_state
Browse files Browse the repository at this point in the history
  • Loading branch information
@bennibolm
bennibolm committed Nov 22, 2023
1 parent b72ed16 commit 9c0a288
Show file tree
Hide file tree
Showing 4 changed files with 197 additions and 104 deletions.
90 changes: 88 additions & 2 deletions examples/structured_2d_dgsem/elixir_euler_double_mach.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,11 +7,97 @@ using Trixi
gamma = 1.4
equations = CompressibleEulerEquations2D(gamma)

initial_condition = Trixi.initial_condition_double_mach_reflection
"""
initial_condition_double_mach_reflection(x, t, equations::CompressibleEulerEquations2D)
Compressible Euler setup for a double Mach reflection problem.
Involves strong shock interactions as well as steady / unsteady flow structures.
Also exercises special boundary conditions along the bottom of the domain that is a mixture of
Dirichlet and slip wall.
See Section IV c on the paper below for details.
- Paul Woodward and Phillip Colella (1984)
The Numerical Simulation of Two-Dimensional Fluid Flows with Strong Shocks.
[DOI: 10.1016/0021-9991(84)90142-6](https://doi.org/10.1016/0021-9991(84)90142-6)
"""
@inline function initial_condition_double_mach_reflection(x, t,
equations::CompressibleEulerEquations2D)
if x[1] < 1 / 6 + (x[2] + 20 * t) / sqrt(3)
phi = pi / 6
sin_phi, cos_phi = sincos(phi)

rho = 8.0
v1 = 8.25 * cos_phi
v2 = -8.25 * sin_phi
p = 116.5
else
rho = 1.4
v1 = 0.0
v2 = 0.0
p = 1.0
end

prim = SVector(rho, v1, v2, p)
return prim2cons(prim, equations)
end

initial_condition = initial_condition_double_mach_reflection

boundary_condition_inflow_outflow = BoundaryConditionCharacteristic(initial_condition)

boundary_conditions = (y_neg = Trixi.boundary_condition_mixed_dirichlet_wall,
# Special mixed boundary condition type for the :y_neg of the domain.
# It is charachteristic-based when x < 1/6 and a slip wall when x >= 1/6
# Note: Only for StructuredMesh
@inline function boundary_condition_mixed_characteristic_wall(u_inner,
normal_direction::AbstractVector,
direction,
x, t, surface_flux_function,
equations::CompressibleEulerEquations2D)
if x[1] < 1 / 6
# From the BoundaryConditionCharacteristic
# get the external state of the solution
u_boundary = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection,
u_inner,
normal_direction,
direction, x, t,
equations)
# Calculate boundary flux
flux = surface_flux_function(u_boundary, u_inner, normal_direction, equations)
else # x[1] >= 1 / 6
# Use the free slip wall BC otherwise
flux = boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t,
surface_flux_function, equations)
end

return flux
end

# Note: Only for StructuredMesh
@inline function Trixi.get_boundary_outer_state(u_inner, cache, t,
boundary_condition::typeof(boundary_condition_mixed_characteristic_wall),
normal_direction::AbstractVector, direction,
equations,
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)

else # if x[1] >= 1 / 6 # boundary_condition_slip_wall
u_rotate = Trixi.rotate_to_x(u_inner, normal_direction, equations)

u_outer = SVector(u_inner[1],
u_inner[2] - 2.0 * u_rotate[2],
u_inner[3] - 2.0 * u_rotate[3],
u_inner[4])
end
return u_outer
end

boundary_conditions = (y_neg = boundary_condition_mixed_characteristic_wall,
y_pos = boundary_condition_inflow_outflow,
x_pos = boundary_condition_inflow_outflow,
x_neg = boundary_condition_inflow_outflow)
Expand Down
90 changes: 88 additions & 2 deletions examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,11 +7,97 @@ using Trixi
gamma = 1.4
equations = CompressibleEulerEquations2D(gamma)

initial_condition = Trixi.initial_condition_double_mach_reflection
"""
initial_condition_double_mach_reflection(x, t, equations::CompressibleEulerEquations2D)
Compressible Euler setup for a double Mach reflection problem.
Involves strong shock interactions as well as steady / unsteady flow structures.
Also exercises special boundary conditions along the bottom of the domain that is a mixture of
Dirichlet and slip wall.
See Section IV c on the paper below for details.
- Paul Woodward and Phillip Colella (1984)
The Numerical Simulation of Two-Dimensional Fluid Flows with Strong Shocks.
[DOI: 10.1016/0021-9991(84)90142-6](https://doi.org/10.1016/0021-9991(84)90142-6)
"""
@inline function initial_condition_double_mach_reflection(x, t,
equations::CompressibleEulerEquations2D)
if x[1] < 1 / 6 + (x[2] + 20 * t) / sqrt(3)
phi = pi / 6
sin_phi, cos_phi = sincos(phi)

rho = 8.0
v1 = 8.25 * cos_phi
v2 = -8.25 * sin_phi
p = 116.5
else
rho = 1.4
v1 = 0.0
v2 = 0.0
p = 1.0
end

prim = SVector(rho, v1, v2, p)
return prim2cons(prim, equations)
end

initial_condition = initial_condition_double_mach_reflection

boundary_condition_inflow_outflow = BoundaryConditionCharacteristic(initial_condition)

boundary_conditions = (y_neg = Trixi.boundary_condition_mixed_dirichlet_wall,
# Special mixed boundary condition type for the :y_neg of the domain.
# It is charachteristic-based when x < 1/6 and a slip wall when x >= 1/6
# Note: Only for StructuredMesh
@inline function boundary_condition_mixed_characteristic_wall(u_inner,
normal_direction::AbstractVector,
direction,
x, t, surface_flux_function,
equations::CompressibleEulerEquations2D)
if x[1] < 1 / 6
# From the BoundaryConditionCharacteristic
# get the external state of the solution
u_boundary = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection,
u_inner,
normal_direction,
direction, x, t,
equations)
# Calculate boundary flux
flux = surface_flux_function(u_boundary, u_inner, normal_direction, equations)
else # x[1] >= 1 / 6
# Use the free slip wall BC otherwise
flux = boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t,
surface_flux_function, equations)
end

return flux
end

# Note: Only for StructuredMesh
@inline function Trixi.get_boundary_outer_state(u_inner, cache, t,
boundary_condition::typeof(boundary_condition_mixed_characteristic_wall),
normal_direction::AbstractVector, direction,
equations,
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)

else # if x[1] >= 1 / 6 # boundary_condition_slip_wall
u_rotate = Trixi.rotate_to_x(u_inner, normal_direction, equations)

u_outer = SVector(u_inner[1],
u_inner[2] - 2.0 * u_rotate[2],
u_inner[3] - 2.0 * u_rotate[3],
u_inner[4])
end
return u_outer
end

boundary_conditions = (y_neg = boundary_condition_mixed_characteristic_wall,
y_pos = boundary_condition_inflow_outflow,
x_pos = boundary_condition_inflow_outflow,
x_neg = boundary_condition_inflow_outflow)
Expand Down
65 changes: 0 additions & 65 deletions src/equations/compressible_euler_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -480,71 +480,6 @@ end
return cons
end

"""
initial_condition_double_mach_reflection(x, t, equations::CompressibleEulerEquations2D)
Compressible Euler setup for a double Mach reflection problem.
Involves strong shock interactions as well as steady / unsteady flow structures.
Also exercises special boundary conditions along the bottom of the domain that is a mixture of
Dirichlet and slip wall.
See Section IV c on the paper below for details.
- Paul Woodward and Phillip Colella (1984)
The Numerical Simulation of Two-Dimensional Fluid Flows with Strong Shocks.
[DOI: 10.1016/0021-9991(84)90142-6](https://doi.org/10.1016/0021-9991(84)90142-6)
"""
@inline function initial_condition_double_mach_reflection(x, t,
equations::CompressibleEulerEquations2D)
if x[1] < 1 / 6 + (x[2] + 20 * t) / sqrt(3)
phi = pi / 6
sin_phi, cos_phi = sincos(phi)

rho = 8.0
v1 = 8.25 * cos_phi
v2 = -8.25 * sin_phi
p = 116.5
else
rho = 1.4
v1 = 0.0
v2 = 0.0
p = 1.0
end

prim = SVector(rho, v1, v2, p)
return prim2cons(prim, equations)
end

# Special mixed boundary condition type for the :Bottom of the domain.
# It is charachteristic when x < 1/6 and a slip wall when x >= 1/6
@inline function boundary_condition_mixed_dirichlet_wall(u_inner,
normal_direction::AbstractVector,
direction,
x, t, surface_flux_function,
equations::CompressibleEulerEquations2D)
# Note: Only for StructuredMesh
if x[1] < 1 / 6
# # From the BoundaryConditionDirichlet
# # get the external value of the solution
# u_boundary = initial_condition_double_mach_reflection(x, t, equations)

# From the BoundaryConditionCharacteristic
# get the external state of the solution
u_boundary = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection,
u_inner,
normal_direction,
direction, x, t,
equations)
# Calculate boundary flux
flux = surface_flux_function(u_boundary, u_inner, normal_direction, equations)
else # x[1] >= 1 / 6
# Use the free slip wall BC otherwise
flux = boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t,
surface_flux_function, equations)
end

return flux
end

# Calculate 2D flux for a single point
@inline function flux(u, orientation::Integer, equations::CompressibleEulerEquations2D)
rho, rho_v1, rho_v2, rho_e = u
Expand Down
56 changes: 21 additions & 35 deletions src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1813,44 +1813,30 @@ end
return nothing
end

@inline function get_boundary_outer_state(u_inner, cache, t,
boundary_condition::typeof(boundary_condition_slip_wall),
orientation::Integer, direction,
equations, 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)

return SVector(u_inner[1],
u_inner[2] - 2.0 * u_rotate[2],
u_inner[3] - 2.0 * u_rotate[3],
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...)
if boundary_condition == boundary_condition_slip_wall #boundary_condition_reflecting_euler_wall
if orientation_or_normal isa AbstractArray
u_rotate = rotate_to_x(u_inner, orientation_or_normal, equations)

return SVector(u_inner[1],
u_inner[2] - 2.0 * u_rotate[2],
u_inner[3] - 2.0 * u_rotate[3],
u_inner[4])
else # orientation_or_normal isa Integer
return SVector(u_inner[1], -u_inner[2], -u_inner[3], u_inner[4])
end
elseif boundary_condition == boundary_condition_mixed_dirichlet_wall
x = 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,
orientation_or_normal,
direction, x, t,
equations)

return u_outer
else # x[1] >= 1 / 6 # boundary_condition_slip_wall
if orientation_or_normal isa AbstractArray
u_rotate = rotate_to_x(u_inner, orientation_or_normal, equations)

return SVector(u_inner[1],
u_inner[2] - 2.0 * u_rotate[2],
u_inner[3] - 2.0 * u_rotate[3],
u_inner[4])
else # orientation_or_normal isa Integer
return SVector(u_inner[1], -u_inner[2], -u_inner[3], u_inner[4])
end
end
end

error("rr")
return u_inner
end

Expand Down

0 comments on commit 9c0a288

Please sign in to comment.