From 96ca8e25d5e1697107042b37a94c1dc525be5f38 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 22 Nov 2023 13:06:27 +0100 Subject: [PATCH] Add elixir with supersonic flow --- ...ir_euler_supersonic_cylinder_sc_subcell.jl | 138 ++++++++++++++++++ .../dgsem_p4est/subcell_limiters_2d.jl | 13 +- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 6 +- 3 files changed, 150 insertions(+), 7 deletions(-) create mode 100644 examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl diff --git a/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl new file mode 100644 index 00000000000..b97862ace17 --- /dev/null +++ b/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 diff --git a/src/solvers/dgsem_p4est/subcell_limiters_2d.jl b/src/solvers/dgsem_p4est/subcell_limiters_2d.jl index e7a8a704478..1acdb05aa7a 100644 --- a/src/solvers/dgsem_p4est/subcell_limiters_2d.jl +++ b/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 diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 007144e0bfd..1e9d63d4602 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/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]