diff --git a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl index f2e97826880..55efabcd346 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl @@ -1,6 +1,6 @@ using Trixi -using OrdinaryDiffEq, Plots +using OrdinaryDiffEq using Downloads: download ############################################################################### @@ -9,14 +9,14 @@ using Downloads: download equations = CompressibleEulerEquations2D(1.4) @inline function initial_condition_mach2_flow(x, t, equations::CompressibleEulerEquations2D) - # set the freestream flow parameters - rho_freestream = 1.4 - v1 = 2.0 - v2 = 0.0 - p_freestream = 1.0 - - prim = SVector(rho_freestream, v1, v2, p_freestream) - return prim2cons(prim, equations) + # set the freestream flow parameters + rho_freestream = 1.4 + v1 = 2.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_mach2_flow @@ -28,10 +28,10 @@ initial_condition = initial_condition_mach2_flow normal_direction::AbstractVector, x, t, surface_flux_function, equations::CompressibleEulerEquations2D) - u_boundary = initial_condition_mach2_flow(x, t, equations) - flux = Trixi.flux(u_boundary, normal_direction, equations) + u_boundary = initial_condition_mach2_flow(x, t, equations) + flux = Trixi.flux(u_boundary, normal_direction, equations) - return flux + return flux end # Supersonic outflow boundary condition. @@ -40,15 +40,15 @@ end @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) + flux = Trixi.flux(u_inner, normal_direction, equations) - return flux + return flux end d = 3 surface_flux = flux_lax_friedrichs -volume_flux = flux_ranocha +volume_flux = flux_ranocha basis = LobattoLegendreBasis(d) shock_indicator = IndicatorHennemannGassner(equations, basis, @@ -67,17 +67,16 @@ solver = DGSEM(polydeg = d, surface_flux = surface_flux, mesh_file = joinpath(@__DIR__, "NACA6412.inp") isfile(mesh_file) || download("https://gist.github.com/DanielDoehring/b434005800a1b0c0b4b50a8772283019/raw/1f6649b3370163d75d268176b51775f8685dd81d/NACA6412.inp", - mesh_file) + mesh_file) boundary_symbols = [:PhysicalLine10, :PhysicalLine20, :PhysicalLine30, :PhysicalLine40] -mesh = P4estMesh{2}(mesh_file, polydeg = d, boundary_symbols=boundary_symbols) +mesh = P4estMesh{2}(mesh_file, polydeg = d, boundary_symbols = boundary_symbols) boundary_conditions = Dict(:PhysicalLine10 => boundary_condition_supersonic_inflow, :PhysicalLine20 => boundary_condition_outflow, # Right boundary :PhysicalLine30 => boundary_condition_slip_wall, # Airfoil - :PhysicalLine40 => boundary_condition_outflow # Top/Bottom - ) + :PhysicalLine40 => boundary_condition_outflow) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) @@ -96,25 +95,10 @@ callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback) -############################################################################### -# run the simulation - -# Positivity limiter might be necessary for examples with strong shocks. Very sensitive -# to the order of the limiter variables, pressure must come first. -stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-7, 1.0e-6), - variables = (pressure, Trixi.density)) - # Run the simulation ############################################################################### -sol = solve(ode, #SSPRK104(; thread = OrdinaryDiffEq.True(), stage_limiter!); - SSPRK104(; thread = OrdinaryDiffEq.True()); +sol = solve(ode, SSPRK104(; thread = OrdinaryDiffEq.True()); dt = 42.0, # overwritten callback = callbacks); summary_callback() # print the timer summary - -plot(sol) - -pd = PlotData2D(sol) -plot(pd["rho"]) -plot!(getmesh(pd)) \ No newline at end of file diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 16bcc38bd27..8138ca6c1af 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -482,7 +482,7 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, node_set_dict = parse_node_sets(meshfile, boundary_symbols) # Read in all elements with associated nodes to specify the boundaries element_node_matrix = parse_elements(meshfile, n_trees, n_dimensions) - + # Initialize boundary information matrix with symbol for no boundary / internal connection boundary_names = fill(Symbol("---"), 2 * n_dimensions, n_trees) @@ -494,22 +494,26 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, # and search for "Node ordering and face numbering on elements" for boundary in keys(node_set_dict) # Check bottom edge - if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] + if tree_nodes[1] in node_set_dict[boundary] && + tree_nodes[2] in node_set_dict[boundary] # Bottom boundary is position 3 in p4est indexing boundary_names[3, tree] = boundary end # Check right edge - if tree_nodes[2] in node_set_dict[boundary] && tree_nodes[3] in node_set_dict[boundary] + if tree_nodes[2] in node_set_dict[boundary] && + tree_nodes[3] in node_set_dict[boundary] # Right boundary is position 2 in p4est indexing boundary_names[2, tree] = boundary end # Check top edge - if tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] + if tree_nodes[3] in node_set_dict[boundary] && + tree_nodes[4] in node_set_dict[boundary] # Top boundary is position 4 in p4est indexing boundary_names[4, tree] = boundary end # Check left edge - if tree_nodes[4] in node_set_dict[boundary] && tree_nodes[1] in node_set_dict[boundary] + if tree_nodes[4] in node_set_dict[boundary] && + tree_nodes[1] in node_set_dict[boundary] # Left boundary is position 1 in p4est indexing boundary_names[1, tree] = boundary end @@ -523,38 +527,50 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, # https://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node26.html for boundary in keys(node_set_dict) # Check "front face" (y_min) - if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] && - tree_nodes[5] in node_set_dict[boundary] && tree_nodes[6] in node_set_dict[boundary] + if tree_nodes[1] in node_set_dict[boundary] && + tree_nodes[2] in node_set_dict[boundary] && + tree_nodes[5] in node_set_dict[boundary] && + tree_nodes[6] in node_set_dict[boundary] # Front face is position 3 in p4est indexing boundary_names[3, tree] = boundary end # Check "back face" (y_max) - if tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] && - tree_nodes[7] in node_set_dict[boundary] && tree_nodes[8] in node_set_dict[boundary] + if tree_nodes[3] in node_set_dict[boundary] && + tree_nodes[4] in node_set_dict[boundary] && + tree_nodes[7] in node_set_dict[boundary] && + tree_nodes[8] in node_set_dict[boundary] # Front face is position 4 in p4est indexing boundary_names[4, tree] = boundary end # Check "left face" (x_min) - if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] && - tree_nodes[5] in node_set_dict[boundary] && tree_nodes[8] in node_set_dict[boundary] + if tree_nodes[1] in node_set_dict[boundary] && + tree_nodes[4] in node_set_dict[boundary] && + tree_nodes[5] in node_set_dict[boundary] && + tree_nodes[8] in node_set_dict[boundary] # Left face is position 1 in p4est indexing boundary_names[1, tree] = boundary end # Check "right face" (x_max) - if tree_nodes[2] in node_set_dict[boundary] && tree_nodes[3] in node_set_dict[boundary] && - tree_nodes[6] in node_set_dict[boundary] && tree_nodes[7] in node_set_dict[boundary] + if tree_nodes[2] in node_set_dict[boundary] && + tree_nodes[3] in node_set_dict[boundary] && + tree_nodes[6] in node_set_dict[boundary] && + tree_nodes[7] in node_set_dict[boundary] # Right face is position 2 in p4est indexing boundary_names[2, tree] = boundary end # Check "bottom face" (z_min) - if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] && - tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] + if tree_nodes[1] in node_set_dict[boundary] && + tree_nodes[2] in node_set_dict[boundary] && + tree_nodes[3] in node_set_dict[boundary] && + tree_nodes[4] in node_set_dict[boundary] # Bottom face is position 5 in p4est indexing boundary_names[5, tree] = boundary end # Check "top face" (z_max) - if tree_nodes[5] in node_set_dict[boundary] && tree_nodes[6] in node_set_dict[boundary] && - tree_nodes[7] in node_set_dict[boundary] && tree_nodes[8] in node_set_dict[boundary] + if tree_nodes[5] in node_set_dict[boundary] && + tree_nodes[6] in node_set_dict[boundary] && + tree_nodes[7] in node_set_dict[boundary] && + tree_nodes[8] in node_set_dict[boundary] # Top face is position 6 in p4est indexing boundary_names[6, tree] = boundary end @@ -568,8 +584,13 @@ end function parse_elements(meshfile, n_trees, n_dims) @assert n_dims ∈ [2, 3] "Only 2D and 3D meshes are supported" - element_types = n_dims == 2 ? ["*ELEMENT, type=CPS4", "*ELEMENT, type=C2D4", "*ELEMENT, type=S4"] : - ["*ELEMENT, type=C3D8"] + element_types = n_dims == 2 ? + [ + "*ELEMENT, type=CPS4", + "*ELEMENT, type=C2D4", + "*ELEMENT, type=S4", + ] : + ["*ELEMENT, type=C3D8"] # 2D quads: 4 nodes + element index, 3D hexes: 8 nodes + element index expected_content_length = n_dims == 2 ? 5 : 9 @@ -621,7 +642,7 @@ function parse_node_sets(meshfile, boundary_symbols) end elseif current_symbol !== nothing # Read only if there was already a nodeset specified # There is always a trailing comma, remove the corresponding empty string - append!(current_nodes, parse.(Int64, split(line, ",")[1:end-1])) + append!(current_nodes, parse.(Int64, split(line, ",")[1:(end - 1)])) end end if current_symbol !== nothing diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index cebc2917d52..d921090be77 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -391,6 +391,27 @@ end end end +@trixi_testset "elixir_euler_airfoil_mach2.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_airfoil_mach2.jl"), + l2=[ + 1.682112049215034e-10, 3.439606229463632e-10, + 1.973815936404324e-10, 7.666245265417042e-10, + ], + linf=[ + 2.636251394960709e-8, 4.4272578936244145e-8, + 2.8049958258193942e-8, 1.0125492977408612e-7, + ], + tspan=(0.0, 0.1)) + # 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)) < 1000 + end +end + @trixi_testset "elixir_eulergravity_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_convergence.jl"), l2=[ diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 4a2d2112c99..6596ecd911f 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -234,6 +234,29 @@ end end end +@trixi_testset "elixir_euler_free_stream_boundaries.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_free_stream_boundaries.jl"), + l2=[ + 3.3824205900551404e-16, 5.761375062854377e-16, + 1.0712466724440564e-15, 1.2802903179864802e-15, + 2.310819174832704e-15, + ], + linf=[ + 2.1094237467877974e-15, 3.3584246494910985e-15, + 9.2148511043888e-15, 1.2101430968414206e-14, + 1.9539925233402755e-14, + ]) + # 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)) < 1000 + end +end + @trixi_testset "elixir_euler_free_stream_extruded.jl with HLLC FLux" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_extruded.jl"), l2=[