Add elixir with supersonic flow

examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl

@@ -0,0 +1,138 @@
+# Channel flow around a cylinder at Mach 3
+#
+# Boundary conditions are supersonic Mach 3 inflow at the left portion of the domain
+# and supersonic outflow at the right portion of the domain. The top and bottom of the
+# channel as well as the cylinder are treated as Euler slip wall boundaries.
+# This flow results in strong shock reflections / interactions as well as Kelvin-Helmholtz
+# instabilities at later times as two Mach stems form above and below the cylinder.
+#
+# For complete details on the problem setup see Section 5.7 of the paper:
+# - Jean-Luc Guermond, Murtazo Nazarov, Bojan Popov, and Ignacio Tomas (2018)
+#   Second-Order Invariant Domain Preserving Approximation of the Euler Equations using Convex Limiting.
+#   [DOI: 10.1137/17M1149961](https://doi.org/10.1137/17M1149961)
+#
+# Keywords: supersonic flow, shock capturing, AMR, unstructured curved mesh, positivity preservation, compressible Euler, 2D
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+@inline function initial_condition_mach3_flow(x, t, equations::CompressibleEulerEquations2D)
+    # set the freestream flow parameters
+    rho_freestream = 1.4
+    v1 = 3.0
+    v2 = 0.0
+    p_freestream = 1.0
+
+    prim = SVector(rho_freestream, v1, v2, p_freestream)
+    return prim2cons(prim, equations)
+end
+
+initial_condition = initial_condition_mach3_flow
+
+# Supersonic inflow boundary condition.
+# Calculate the boundary flux entirely from the external solution state, i.e., set
+# external solution state values for everything entering the domain.
+@inline function boundary_condition_supersonic_inflow(u_inner,
+                                                      normal_direction::AbstractVector,
+                                                      x, t, surface_flux_function,
+                                                      equations::CompressibleEulerEquations2D)
+    u_boundary = initial_condition_mach3_flow(x, t, equations)
+    flux = Trixi.flux(u_boundary, normal_direction, equations)
+
+    return flux
+end
+
+# Supersonic outflow boundary condition.
+# Calculate the boundary flux entirely from the internal solution state. Analogous to supersonic inflow
+# except all the solution state values are set from the internal solution as everything leaves the domain
+@inline function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t,
+                                            surface_flux_function,
+                                            equations::CompressibleEulerEquations2D)
+    flux = Trixi.flux(u_inner, normal_direction, equations)
+
+    return flux
+end
+
+# boundary_condition_inflow_outflow = BoundaryConditionCharacteristic(initial_condition)
+
+# boundary_conditions = Dict(:Bottom => boundary_condition_slip_wall,
+#                            :Circle => boundary_condition_slip_wall,
+#                            :Top => boundary_condition_slip_wall,
+#                            :Right => boundary_condition_inflow_outflow,
+#                            :Left => boundary_condition_inflow_outflow)
+
+boundary_conditions = Dict(:Bottom => boundary_condition_slip_wall,
+                           :Circle => boundary_condition_slip_wall,
+                           :Top => boundary_condition_slip_wall,
+                           :Right => boundary_condition_outflow,
+                           :Left => boundary_condition_supersonic_inflow)
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                local_minmax_variables_cons = ["rho"],
+                                positivity_variables_cons = ["rho"],
+                                positivity_variables_nonlinear = [pressure],
+                                positivity_correction_factor = 0.5,
+                                spec_entropy = true,
+                                bar_states = false,
+                                max_iterations_newton = 1000,
+                                newton_tolerances = (1.0e-14, 1.0e-15))
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+# Get the unstructured quad mesh from a file (downloads the file if not available locally)
+default_mesh_file = joinpath(@__DIR__, "abaqus_cylinder_in_channel.inp")
+isfile(default_mesh_file) ||
+    download("https://gist.githubusercontent.com/andrewwinters5000/a08f78f6b185b63c3baeff911a63f628/raw/addac716ea0541f588b9d2bd3f92f643eb27b88f/abaqus_cylinder_in_channel.inp",
+             default_mesh_file)
+mesh_file = default_mesh_file
+
+mesh = P4estMesh{2}(mesh_file, initial_refinement_level = 1)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.4)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback,
+                        save_solution)
+
+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
+                  save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

src/solvers/dgsem_p4est/subcell_limiters_2d.jl

@@ -167,17 +167,20 @@ function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, sem
         j_secondary = j_secondary_start
 
         for node in eachnode(dg)
-            var_primary = variable(get_node_vars(u, equations, dg, i_primary, j_primary, primary_element), equations)
-            var_secondary = variable(get_node_vars(u, equations, dg, i_secondary, j_secondary, secondary_element), equations)
+            var_primary = variable(get_node_vars(u, equations, dg, i_primary, j_primary,
+                                                 primary_element), equations)
+            var_secondary = variable(get_node_vars(u, equations, dg, i_secondary,
+                                                   j_secondary, secondary_element),
+                                     equations)
 
             var_minmax[i_primary, j_primary, primary_element] = minmax(var_minmax[i_primary,
                                                                                   j_primary,
                                                                                   primary_element],
                                                                        var_secondary)
             var_minmax[i_secondary, j_secondary, secondary_element] = minmax(var_minmax[i_secondary,
-                                                                                           j_secondary,
-                                                                                           secondary_element],
-                                                                                var_primary)
+                                                                                        j_secondary,
+                                                                                        secondary_element],
+                                                                             var_primary)
 
             # Increment primary element indices
             i_primary += i_primary_step

src/solvers/dgsem_tree/subcell_limiters_2d.jl

@@ -475,7 +475,8 @@ end
 end
 
 @inline function idp_spec_entropy!(alpha, limiter, u, t, dt, semi,
-                                   mesh::Union{StructuredMesh{2}, P4estMesh{2}}, elements)
+                                   mesh::Union{StructuredMesh{2}, P4estMesh{2}},
+                                   elements)
     _, equations, dg, cache = mesh_equations_solver_cache(semi)
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
     s_min = variable_bounds[:spec_entropy_min]
@@ -526,7 +527,8 @@ end
 end
 
 @inline function idp_math_entropy!(alpha, limiter, u, t, dt, semi,
-                                   mesh::Union{StructuredMesh{2}, P4estMesh{2}}, elements)
+                                   mesh::Union{StructuredMesh{2}, P4estMesh{2}},
+                                   elements)
     _, equations, dg, cache = mesh_equations_solver_cache(semi)
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
     s_max = variable_bounds[:math_entropy_max]