diff --git a/examples/p4est_2d_dgsem/elixir_euler_double_mach.jl b/examples/p4est_2d_dgsem/elixir_euler_double_mach.jl new file mode 100644 index 00000000000..d56d1d0f8b5 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_double_mach.jl @@ -0,0 +1,168 @@ + +using OrdinaryDiffEq +using Trixi +using LinearAlgebra: norm # for use in get_boundary_outer_state + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +""" + 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) + +# 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 P4estMesh +@inline function boundary_condition_mixed_characteristic_wall(u_inner, + normal_direction::AbstractVector, + 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, + x, t, + equations) + # Calculate boundary flux + flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations) + else # x[1] >= 1 / 6 + # Use the free slip wall BC otherwise + flux = boundary_condition_slip_wall(u_inner, normal_direction, x, t, + surface_flux_function, equations) + end + + return flux +end + +# Note: Only for P4estMesh +@inline function Trixi.get_boundary_outer_state(u_inner, cache, t, + boundary_condition::typeof(boundary_condition_mixed_characteristic_wall), + normal_direction::AbstractVector, direction, + mesh::P4estMesh{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, + x, t, equations) + + else # if x[1] >= 1 / 6 # boundary_condition_slip_wall + factor = (normal_direction[1] * u_inner[2] + normal_direction[2] * u_inner[3]) + u_normal = (factor / sum(normal_direction .^ 2)) * normal_direction + + u_outer = SVector(u_inner[1], + u_inner[2] - 2.0 * u_normal[1], + u_inner[3] - 2.0 * u_normal[2], + u_inner[4]) + end + return u_outer +end + +boundary_conditions = Dict(: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) + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) + +limiter_idp = SubcellLimiterIDP(equations, basis; + local_minmax_variables_cons = ["rho"], + spec_entropy = true, + positivity_correction_factor = 0.1, + max_iterations_newton = 100, + bar_states = true) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +initial_refinement_level = 4 +trees_per_dimension = (4 * 2^initial_refinement_level, 2^initial_refinement_level) +coordinates_min = (0.0, 0.0) +coordinates_max = (4.0, 1.0) +mesh = P4estMesh(trees_per_dimension, polydeg = polydeg, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + initial_refinement_level = 0, + periodicity = false) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.2) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (entropy,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.9) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback, + save_solution) + +############################################################################### +# run the simulation + +stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false)) + +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/structured_2d_dgsem/elixir_euler_double_mach.jl b/examples/structured_2d_dgsem/elixir_euler_double_mach.jl index ef18aa6a68f..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, @@ -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 afca1b7852e..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, @@ -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/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 59ce2257a3e..4adb1a6fb20 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -363,6 +363,33 @@ end end end +@trixi_testset "elixir_euler_double_mach.jl" begin + # Same results as for StructuredMesh with initial_refinement_level=2 + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_double_mach.jl"), + l2=[ + 0.87417841433288, + 6.669726935171785, + 3.4980245896465387, + 76.33557073534843, + ], + linf=[ + 11.428353671462515, + 142.73486852796972, + 38.91639544578682, + 1651.7541392659086, + ], + initial_refinement_level=2, + tspan=(0.0, 0.05)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 10000 + end +end + @trixi_testset "elixir_euler_supersonic_cylinder.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_supersonic_cylinder.jl"), l2=[