From 1a09eaa1a707e6e77fa84591dd44e13eaf1f5b90 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> Date: Tue, 5 Dec 2023 10:47:08 +0100 Subject: [PATCH] Add support for P4estMesh to `subcell-limiting` (#121) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Add support for P4estMesh to `subcell-limiting` * Add elixir with supersonic flow * Extract inner code of `characteristic_boundary_value_function` * Add support for `BoundaryConditionDirichlet` for P4estMesh * Fix calculation of `vn` * Remove default routine for `get_boundary_outer_state` * Revise calculation of bar states * Fix allocations with `foreach(enumerate(...))` * Fix normalization in `get_boundary_outer_state` * Add TODO note for future * Remove unnecessary code; format * Adapt test errors after normalizing normal vectors * Add dispatch for equations * Fix normalization; Is done inside boundary_value_function * Fix bug in calculation of bar state at interfaces for P4estMesh * Add mesh as parameter of `get_boundary_outer_state` - Fixes bug that P4estMesh calls the right `characteristic_boundary_state` routine * Add p4est double mach elixir * Fix error * Remove comment * Add test for supersonic flow * fmt * Add comments * Update src/equations/compressible_euler_2d.jl Co-authored-by: Andrés Rueda-Ramírez * fmt * Fix error --------- Co-authored-by: Andrés Rueda-Ramírez --- .../elixir_euler_double_mach.jl | 167 ++++++++++++ .../elixir_euler_free_stream_sc_subcell.jl | 91 +++++++ ...ir_euler_supersonic_cylinder_sc_subcell.jl | 152 +++++++++++ .../elixir_euler_double_mach.jl | 1 - .../elixir_euler_double_mach_MCL.jl | 1 - src/auxiliary/auxiliary.jl | 17 ++ .../subcell_limiter_idp_correction_2d.jl | 3 +- src/callbacks_step/stepsize_dg2d.jl | 9 +- src/equations/compressible_euler_2d.jl | 41 ++- src/equations/equations.jl | 16 ++ src/solvers/dgsem_p4est/dg.jl | 3 + .../dgsem_p4est/dg_2d_subcell_limiters.jl | 212 +++++++++++++++ .../dgsem_p4est/subcell_limiters_2d.jl | 250 ++++++++++++++++++ .../dg_2d_subcell_limiters.jl | 7 +- .../dgsem_tree/dg_2d_subcell_limiters.jl | 27 +- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 15 +- src/time_integration/methods_SSP.jl | 3 +- test/test_p4est_2d.jl | 81 ++++++ 18 files changed, 1066 insertions(+), 30 deletions(-) create mode 100644 examples/p4est_2d_dgsem/elixir_euler_double_mach.jl create mode 100644 examples/p4est_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl create mode 100644 examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl create mode 100644 src/solvers/dgsem_p4est/dg_2d_subcell_limiters.jl create mode 100644 src/solvers/dgsem_p4est/subcell_limiters_2d.jl 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..949bc40e641 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_double_mach.jl @@ -0,0 +1,167 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# 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/p4est_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl b/examples/p4est_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl new file mode 100644 index 00000000000..9d67ad76ce0 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl @@ -0,0 +1,91 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +initial_condition = initial_condition_constant + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +limiter_idp = SubcellLimiterIDP(equations, basis; + positivity_variables_cons = ["rho"], + positivity_variables_nonlinear = [pressure], + positivity_correction_factor = 0.1, + spec_entropy = false, + 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) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but reduced to 2D. +# This particular mesh is unstructured in the yz-plane, but extruded in x-direction. +# Apply the warping mapping in the yz-plane to get a curved 2D mesh that is extruded +# in x-direction to ensure free stream preservation on a non-conforming mesh. +# See https://doi.org/10.1007/s10915-018-00897-9, Section 6. + +# Mapping as described in https://arxiv.org/abs/2012.12040, but reduced to 2D +function mapping(xi_, eta_) + # Transform input variables between -1 and 1 onto [0,3] + xi = 1.5 * xi_ + 1.5 + eta = 1.5 * eta_ + 1.5 + + y = eta + 3 / 8 * (cos(1.5 * pi * (2 * xi - 3) / 3) * + cos(0.5 * pi * (2 * eta - 3) / 3)) + + x = xi + 3 / 8 * (cos(0.5 * pi * (2 * xi - 3) / 3) * + cos(2 * pi * (2 * y - 3) / 3)) + + return SVector(x, y) +end + +trees_per_dimension = (16, 16) + +# Create P4estMesh with 16 x 16 trees and 16 x 16 elements +mesh = P4estMesh(trees_per_dimension, polydeg = 3, + mapping = mapping, + initial_refinement_level = 0, periodicity = true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 10000, + 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 + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary 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..8b6b93e0ac1 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl @@ -0,0 +1,152 @@ +# 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 + +@inline function Trixi.get_boundary_outer_state(u_inner, cache, t, + boundary_condition::typeof(boundary_condition_supersonic_inflow), + normal_direction::AbstractVector, direction, + mesh::P4estMesh{2}, equations, dg, + indices...) + x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...) + + return initial_condition_mach3_flow(x, t, equations) +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 + +@inline function Trixi.get_boundary_outer_state(u_inner, cache, t, + boundary_condition::typeof(boundary_condition_outflow), + orientation_or_normal, direction, + mesh::P4estMesh{2}, equations, dg, + indices...) + return u_inner +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; + positivity_variables_cons = ["rho"], + positivity_variables_nonlinear = [pressure], + spec_entropy = false, + bar_states = false) +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 = 0) + +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/examples/structured_2d_dgsem/elixir_euler_double_mach.jl b/examples/structured_2d_dgsem/elixir_euler_double_mach.jl index 8153d2634b8..c40d79e3cc5 100644 --- a/examples/structured_2d_dgsem/elixir_euler_double_mach.jl +++ b/examples/structured_2d_dgsem/elixir_euler_double_mach.jl @@ -1,7 +1,6 @@ using OrdinaryDiffEq using Trixi -using LinearAlgebra: norm # for use in get_boundary_outer_state ############################################################################### # semidiscretization of the compressible Euler equations 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 e85c263429f..e56392f9b2d 100644 --- a/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl +++ b/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl @@ -1,7 +1,6 @@ using OrdinaryDiffEq using Trixi -using LinearAlgebra: norm # for use in get_boundary_outer_state ############################################################################### # semidiscretization of the compressible Euler equations diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index 1f7d30d6aa8..ac3ff63ade7 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -345,4 +345,21 @@ function register_error_hints() return nothing end + +# Same as `foreach(enumerate(something))`, but without allocations. +# +# Note that compile times may increase if this is used with big tuples. +# TODO: Add comment in the respective PR (here and where it is used: `dg_p4est/dg_2d_subcell_limiters.jl`) +@inline foreach_enumerate(func, collection) = foreach_enumerate(func, collection, 1) +@inline foreach_enumerate(func, collection::Tuple{}, index) = nothing + +@inline function foreach_enumerate(func, collection, index) + element = first(collection) + remaining_collection = Base.tail(collection) + + func((index, element)) + + # Process remaining collection + foreach_enumerate(func, remaining_collection, index + 1) +end end # @muladd diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl index f7c6fd55828..4ab1357de6f 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl +++ b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl @@ -26,7 +26,8 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations, dg, cache) return nothing end -function perform_idp_correction!(u, dt, mesh::StructuredMesh{2}, equations, dg, cache) +function perform_idp_correction!(u, dt, mesh::Union{StructuredMesh{2}, P4estMesh{2}}, + equations, dg, cache) if dg.volume_integral.limiter.smoothness_indicator elements = cache.element_ids_dgfv else diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index f66686ff363..92c21d65e1d 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -43,7 +43,7 @@ function max_dt(u, t, mesh::TreeMesh{2}, return 2 / (nnodes(dg) * max_scaled_speed) end -@inline function max_dt(u, t, mesh::Union{TreeMesh{2}, StructuredMesh{2}}, +@inline function max_dt(u, t, mesh::Union{TreeMesh{2}, StructuredMesh{2}, P4estMesh{2}}, constant_speed::False, equations, semi, dg::DG, cache, limiter::Union{SubcellLimiterIDP, SubcellLimiterMCL}) @unpack inverse_weights = dg.basis @@ -75,7 +75,7 @@ end J = 1 / cache.elements.inverse_jacobian[element] end for j in eachnode(dg), i in eachnode(dg) - if mesh isa StructuredMesh{2} + if (mesh isa StructuredMesh{2}) || (mesh isa P4estMesh{2}) J = 1 / cache.elements.inverse_jacobian[i, j, element] end denom = inverse_weights[i] * @@ -114,8 +114,9 @@ end return 2 / (nnodes(dg) * max_scaled_speed) end -@inline function max_dt_RK(u, t, mesh::StructuredMesh{2}, constant_speed, equations, - dg::DG, cache, limiter, element_ids_dg) +@inline function max_dt_RK(u, t, mesh::Union{StructuredMesh{2}, P4estMesh{2}}, + constant_speed, equations, dg::DG, cache, limiter, + element_ids_dg) @unpack contravariant_vectors, inverse_jacobian = cache.elements max_scaled_speed = nextfloat(zero(t)) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 8d42a507a56..88e136f996f 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -395,6 +395,7 @@ Should be used together with [`StructuredMesh`](@ref). end # TODO: Add docstring when about to merge. +# Using with TreeMesh{2} @inline function characteristic_boundary_value_function(outer_boundary_value_function, u_inner, orientation::Integer, direction, x, t, @@ -414,11 +415,12 @@ end vn = factor * u_inner[3] * srho end - return characteristic_boundary_value_function_inner(outer_boundary_value_function, - u_inner, srho, vn, x, t, - equations) + return calc_characteristic_boundary_value_function(outer_boundary_value_function, + u_inner, srho, vn, x, t, + equations) end +# Using with StructuredMesh{2} @inline function characteristic_boundary_value_function(outer_boundary_value_function, u_inner, normal_direction::AbstractVector, @@ -437,15 +439,34 @@ end (normal_direction[1] * u_inner[2] + normal_direction[2] * u_inner[3]) / norm(normal_direction) - return characteristic_boundary_value_function_inner(outer_boundary_value_function, - u_inner, srho, vn, x, t, - equations) + return calc_characteristic_boundary_value_function(outer_boundary_value_function, + u_inner, srho, vn, x, t, + equations) end -# Inner function to distinguish between different mesh types. -@inline function characteristic_boundary_value_function_inner(outer_boundary_value_function, - u_inner, srho, vn, x, t, - equations::CompressibleEulerEquations2D) +# Using with P4estMesh{2} +@inline function characteristic_boundary_value_function(outer_boundary_value_function, + u_inner, + normal_direction::AbstractVector, + x, t, + equations::CompressibleEulerEquations2D) + # Get inverse of density + srho = 1 / u_inner[1] + + # Get normal velocity + vn = srho * (normal_direction[1] * u_inner[2] + normal_direction[2] * u_inner[3]) / + norm(normal_direction) + + return calc_characteristic_boundary_value_function(outer_boundary_value_function, + u_inner, srho, vn, x, t, + equations) +end + +# Function to compute the outer state of the characteristics-based boundary condition. +# This function is called by all mesh types. +@inline function calc_characteristic_boundary_value_function(outer_boundary_value_function, + u_inner, srho, vn, x, t, + equations::CompressibleEulerEquations2D) # get pressure and Mach from state p = pressure(u_inner, equations) a = sqrt(equations.gamma * p * srho) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 31d1c8a8184..6d6e3e92a98 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -243,6 +243,22 @@ end return flux end +# Characteristic-based boundary condition for use with P4estMesh +@inline function (boundary_condition::BoundaryConditionCharacteristic)(u_inner, + normal_direction, + x, t, + surface_flux_function, + equations) + u_boundary = boundary_condition.boundary_value_function(boundary_condition.outer_boundary_value_function, + u_inner, normal_direction, + x, t, equations) + + # Calculate boundary flux + flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations) + + return flux +end + # operator types used for dispatch on parabolic boundary fluxes struct Gradient end struct Divergence end diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index ec50627d3ef..ad3b0717498 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -52,4 +52,7 @@ include("dg_2d_parabolic.jl") include("dg_3d.jl") include("dg_3d_parabolic.jl") include("dg_parallel.jl") + +include("subcell_limiters_2d.jl") +include("dg_2d_subcell_limiters.jl") end # @muladd diff --git a/src/solvers/dgsem_p4est/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_p4est/dg_2d_subcell_limiters.jl new file mode 100644 index 00000000000..f8a5ea24ec7 --- /dev/null +++ b/src/solvers/dgsem_p4est/dg_2d_subcell_limiters.jl @@ -0,0 +1,212 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@inline function calc_lambdas_bar_states_interface!(u, t, limiter, boundary_conditions, + mesh::P4estMesh{2}, equations, + dg, cache; calc_bar_states = true) + (; contravariant_vectors) = cache.elements + (; lambda1, lambda2, bar_states1, bar_states2) = limiter.cache.container_bar_states + + (; neighbor_ids, node_indices) = cache.interfaces + index_range = eachnode(dg) + + # Calc lambdas and bar states at interfaces and periodic boundaries + for interface in eachinterface(dg, cache) + # Get element and side index information on the primary element + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + primary_direction = indices2direction(primary_indices) + + # Get element and side index information on the secondary element + secondary_element = neighbor_ids[2, interface] + secondary_indices = node_indices[2, interface] + secondary_direction = indices2direction(secondary_indices) + + # Create the local i,j indexing + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_range) + i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], + index_range) + j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], + index_range) + + i_primary = i_primary_start - i_primary_step + j_primary = j_primary_start - j_primary_step + i_secondary = i_secondary_start - i_secondary_step + j_secondary = j_secondary_start - j_secondary_step + + for node in eachnode(dg) + # Increment primary element indices + i_primary += i_primary_step + j_primary += j_primary_step + i_secondary += i_secondary_step + j_secondary += j_secondary_step + + # Get the normal direction on the primary element. + # Contravariant vectors at interfaces in negative coordinate direction + # are pointing inwards. This is handled by `get_normal_direction`. + normal_direction = get_normal_direction(primary_direction, + contravariant_vectors, i_primary, + j_primary, primary_element) + + u_primary = get_node_vars(u, equations, dg, i_primary, j_primary, + primary_element) + u_secondary = get_node_vars(u, equations, dg, i_secondary, j_secondary, + secondary_element) + + lambda = max_abs_speed_naive(u_primary, u_secondary, normal_direction, + equations) + if primary_direction == 1 + lambda1[i_primary, j_primary, primary_element] = lambda + elseif primary_direction == 2 + lambda1[i_primary + 1, j_primary, primary_element] = lambda + elseif primary_direction == 3 + lambda2[i_primary, j_primary, primary_element] = lambda + else # primary_direction == 4 + lambda2[i_primary, j_primary + 1, primary_element] = lambda + end + if secondary_direction == 1 + lambda1[i_secondary, j_secondary, secondary_element] = lambda + elseif secondary_direction == 2 + lambda1[i_secondary + 1, j_secondary, secondary_element] = lambda + elseif secondary_direction == 3 + lambda2[i_secondary, j_secondary, secondary_element] = lambda + else # secondary_direction == 4 + lambda2[i_secondary, j_secondary + 1, secondary_element] = lambda + end + + !calc_bar_states && continue + + flux_primary = flux(u_primary, normal_direction, equations) + flux_secondary = flux(u_secondary, normal_direction, equations) + + bar_state = 0.5 * (u_primary + u_secondary) - + 0.5 * (flux_secondary - flux_primary) / lambda + if primary_direction == 1 + set_node_vars!(bar_states1, bar_state, equations, dg, + i_primary, j_primary, primary_element) + elseif primary_direction == 2 + set_node_vars!(bar_states1, bar_state, equations, dg, + i_primary + 1, j_primary, primary_element) + elseif primary_direction == 3 + set_node_vars!(bar_states2, bar_state, equations, dg, + i_primary, j_primary, primary_element) + else # primary_direction == 4 + set_node_vars!(bar_states2, bar_state, equations, dg, + i_primary, j_primary + 1, primary_element) + end + if secondary_direction == 1 + set_node_vars!(bar_states1, bar_state, equations, dg, + i_secondary, j_secondary, secondary_element) + elseif secondary_direction == 2 + set_node_vars!(bar_states1, bar_state, equations, dg, + i_secondary + 1, j_secondary, secondary_element) + elseif secondary_direction == 3 + set_node_vars!(bar_states2, bar_state, equations, dg, + i_secondary, j_secondary, secondary_element) + else # secondary_direction == 4 + set_node_vars!(bar_states2, bar_state, equations, dg, + i_secondary, j_secondary + 1, secondary_element) + end + end + end + + calc_lambdas_bar_states_boundary!(u, t, limiter, boundary_conditions, + mesh, equations, dg, cache; + calc_bar_states = calc_bar_states) + + return nothing +end + +@inline function calc_lambdas_bar_states_boundary!(u, t, limiter, + boundary_conditions::BoundaryConditionPeriodic, + mesh::P4estMesh{2}, equations, dg, + cache; calc_bar_states = true) + return nothing +end + +# Calc lambdas and bar states at physical boundaries +@inline function calc_lambdas_bar_states_boundary!(u, t, limiter, boundary_conditions, + mesh::P4estMesh{2}, equations, dg, + cache; calc_bar_states = true) + (; boundary_condition_types, boundary_indices) = boundary_conditions + (; contravariant_vectors) = cache.elements + + (; lambda1, lambda2, bar_states1, bar_states2) = limiter.cache.container_bar_states + + (; boundaries) = cache + index_range = eachnode(dg) + + # Allocation-free version of `foreach(enumerate((...))` + foreach_enumerate(boundary_condition_types) do (i, boundary_condition) + for boundary in boundary_indices[i] + element = boundaries.neighbor_ids[boundary] + node_indices = boundaries.node_indices[boundary] + direction = indices2direction(node_indices) + + i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], + index_range) + j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], + index_range) + + i_node = i_node_start - i_node_step + j_node = j_node_start - j_node_step + for node in eachnode(dg) + i_node += i_node_step + j_node += j_node_step + + normal_direction = get_normal_direction(direction, + contravariant_vectors, + i_node, j_node, element) + + u_inner = get_node_vars(u, equations, dg, i_node, j_node, element) + u_outer = get_boundary_outer_state(u_inner, cache, t, + boundary_condition, normal_direction, + direction, mesh, equations, dg, + i_node, j_node, element) + + lambda = max_abs_speed_naive(u_inner, u_outer, normal_direction, + equations) + if direction == 1 + lambda1[i_node, j_node, element] = lambda + elseif direction == 2 + lambda1[i_node + 1, j_node, element] = lambda + elseif direction == 3 + lambda2[i_node, j_node, element] = lambda + else # direction == 4 + lambda2[i_node, j_node + 1, element] = lambda + end + + !calc_bar_states && continue + + flux_inner = flux(u_inner, normal_direction, equations) + flux_outer = flux(u_outer, normal_direction, equations) + + bar_state = 0.5 * (u_inner + u_outer) - + 0.5 * (flux_outer - flux_inner) / lambda + if direction == 1 + set_node_vars!(bar_states1, bar_state, equations, dg, + i_node, j_node, element) + elseif direction == 2 + set_node_vars!(bar_states1, bar_state, equations, dg, + i_node + 1, j_node, element) + elseif direction == 3 + set_node_vars!(bar_states2, bar_state, equations, dg, + i_node, j_node, element) + else # direction == 4 + set_node_vars!(bar_states2, bar_state, equations, dg, + i_node, j_node + 1, element) + end + end + end + end + + return nothing +end +end # @muladd diff --git a/src/solvers/dgsem_p4est/subcell_limiters_2d.jl b/src/solvers/dgsem_p4est/subcell_limiters_2d.jl new file mode 100644 index 00000000000..d014a7b0380 --- /dev/null +++ b/src/solvers/dgsem_p4est/subcell_limiters_2d.jl @@ -0,0 +1,250 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +function calc_bounds_twosided_interface!(var_min, var_max, variable, u, t, semi, + mesh::P4estMesh{2}) + _, equations, dg, cache = mesh_equations_solver_cache(semi) + (; boundary_conditions) = semi + + (; neighbor_ids, node_indices) = cache.interfaces + index_range = eachnode(dg) + + for interface in eachinterface(dg, cache) + # Get element and side index information on the primary element + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + + # Get element and side index information on the secondary element + secondary_element = neighbor_ids[2, interface] + secondary_indices = node_indices[2, interface] + + # Create the local i,j indexing + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_range) + i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], + index_range) + j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], + index_range) + + i_primary = i_primary_start + j_primary = j_primary_start + i_secondary = i_secondary_start + j_secondary = j_secondary_start + + for node in eachnode(dg) + var_primary = u[variable, i_primary, j_primary, primary_element] + var_secondary = u[variable, i_secondary, j_secondary, secondary_element] + + var_min[i_primary, j_primary, primary_element] = min(var_min[i_primary, + j_primary, + primary_element], + var_secondary) + var_max[i_primary, j_primary, primary_element] = max(var_max[i_primary, + j_primary, + primary_element], + var_secondary) + + var_min[i_secondary, j_secondary, secondary_element] = min(var_min[i_secondary, + j_secondary, + secondary_element], + var_primary) + var_max[i_secondary, j_secondary, secondary_element] = max(var_max[i_secondary, + j_secondary, + secondary_element], + var_primary) + + # Increment primary element indices + i_primary += i_primary_step + j_primary += j_primary_step + i_secondary += i_secondary_step + j_secondary += j_secondary_step + end + end + + calc_bounds_twosided_interface_inner!(var_min, var_max, variable, u, t, + boundary_conditions, + mesh, equations, dg, cache) + + return nothing +end + +@inline function calc_bounds_twosided_interface_inner!(var_min, var_max, variable, u, t, + boundary_conditions::BoundaryConditionPeriodic, + mesh, equations, dg, cache) + return nothing +end + +@inline function calc_bounds_twosided_interface_inner!(var_min, var_max, variable, u, t, + boundary_conditions, + mesh, equations, dg, cache) + (; boundary_condition_types, boundary_indices) = boundary_conditions + (; contravariant_vectors) = cache.elements + + (; boundaries) = cache + index_range = eachnode(dg) + + foreach(enumerate(boundary_condition_types)) do (i, boundary_condition) + for boundary in boundary_indices[i] + element = boundaries.neighbor_ids[boundary] + node_indices = boundaries.node_indices[boundary] + direction = indices2direction(node_indices) + + i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], + index_range) + j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], + index_range) + + i_node = i_node_start + j_node = j_node_start + for i in eachnode(dg) + normal_direction = get_normal_direction(direction, + contravariant_vectors, + i_node, j_node, element) + + u_inner = get_node_vars(u, equations, dg, i_node, j_node, element) + + u_outer = get_boundary_outer_state(u_inner, cache, t, + boundary_condition, normal_direction, + direction, mesh, equations, dg, + i_node, j_node, element) + var_outer = u_outer[variable] + + var_min[i_node, j_node, element] = min(var_min[i_node, j_node, element], + var_outer) + var_max[i_node, j_node, element] = max(var_max[i_node, j_node, element], + var_outer) + + i_node += i_node_step + j_node += j_node_step + end + end + end + + return nothing +end + +function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, semi, + mesh::P4estMesh{2}) + _, equations, dg, cache = mesh_equations_solver_cache(semi) + (; boundary_conditions) = semi + + (; neighbor_ids, node_indices) = cache.interfaces + index_range = eachnode(dg) + index_end = last(index_range) + + for interface in eachinterface(dg, cache) + # Get element and side index information on the primary element + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + + # Get element and side index information on the secondary element + secondary_element = neighbor_ids[2, interface] + secondary_indices = node_indices[2, interface] + + # Create the local i,j indexing + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_range) + i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], + index_range) + j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], + index_range) + + i_primary = i_primary_start + j_primary = j_primary_start + i_secondary = i_secondary_start + 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_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) + + # Increment primary element indices + i_primary += i_primary_step + j_primary += j_primary_step + i_secondary += i_secondary_step + j_secondary += j_secondary_step + end + end + + calc_bounds_onesided_interface_inner!(var_minmax, minmax, variable, u, t, + boundary_conditions, + mesh, equations, dg, cache) + + return nothing +end + +@inline function calc_bounds_onesided_interface_inner!(var_minmax, minmax, variable, u, + t, + boundary_conditions::BoundaryConditionPeriodic, + mesh, equations, dg, cache) + return nothing +end + +@inline function calc_bounds_onesided_interface_inner!(var_minmax, minmax, variable, u, + t, boundary_conditions, + mesh, equations, dg, cache) + (; boundary_condition_types, boundary_indices) = boundary_conditions + (; contravariant_vectors) = cache.elements + + (; boundaries) = cache + index_range = eachnode(dg) + + foreach(enumerate(boundary_condition_types)) do (i, boundary_condition) + for boundary in boundary_indices[i] + element = boundaries.neighbor_ids[boundary] + node_indices = boundaries.node_indices[boundary] + direction = indices2direction(node_indices) + + i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], + index_range) + j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], + index_range) + + i_node = i_node_start + j_node = j_node_start + for node in eachnode(dg) + normal_direction = get_normal_direction(direction, + contravariant_vectors, + i_node, j_node, element) + + u_inner = get_node_vars(u, equations, dg, i_node, j_node, element) + + u_outer = get_boundary_outer_state(u_inner, cache, t, + boundary_condition, normal_direction, + direction, mesh, equations, dg, + i_node, j_node, element) + var_outer = variable(u_outer, equations) + + var_minmax[i_node, j_node, element] = minmax(var_minmax[i_node, j_node, + element], + var_outer) + + i_node += i_node_step + j_node += j_node_step + end + end + end + + return nothing +end +end # @muladd diff --git a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl index c6f3dcc1f10..3a1a5602141 100644 --- a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl @@ -10,8 +10,8 @@ # # See also `flux_differencing_kernel!`. @inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, u, - mesh::StructuredMesh{2}, nonconservative_terms::False, - equations, + mesh::Union{StructuredMesh{2}, P4estMesh{2}}, + nonconservative_terms::False, equations, volume_flux, dg::DGSEM, element, cache) (; contravariant_vectors) = cache.elements (; weights, derivative_split) = dg.basis @@ -109,7 +109,8 @@ return nothing end -@inline function calc_lambdas_bar_states!(u, t, mesh::StructuredMesh{2}, +@inline function calc_lambdas_bar_states!(u, t, + mesh::Union{StructuredMesh{2}, P4estMesh{2}}, nonconservative_terms, equations, limiter, dg, cache, boundary_conditions; calc_bar_states = true) diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index aff68284b40..c1afd67306e 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -5,8 +5,9 @@ @muladd begin #! format: noindent -function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}}, equations, - volume_integral::VolumeIntegralSubcellLimiting, dg::DG, uEltype) +function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, P4estMesh{2}}, + equations, volume_integral::VolumeIntegralSubcellLimiting, + dg::DG, uEltype) cache = create_cache(mesh, equations, VolumeIntegralPureLGLFiniteVolume(volume_integral.volume_flux_fv), dg, uEltype) @@ -61,7 +62,8 @@ function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}}, equations, end function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{2}, StructuredMesh{2}}, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + P4estMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralSubcellLimiting, dg::DGSEM, cache, t, boundary_conditions) @@ -124,7 +126,8 @@ function calc_volume_integral!(du, u, end @inline function subcell_limiting_kernel!(du, u, element, - mesh::Union{TreeMesh{2}, StructuredMesh{2}}, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + P4estMesh{2}}, nonconservative_terms, equations, volume_integral, limiter::SubcellLimiterIDP, dg::DGSEM, cache) @@ -1879,4 +1882,20 @@ end return u_outer end + +@inline function get_boundary_outer_state(u_inner, cache, t, + boundary_condition::BoundaryConditionCharacteristic, + normal_direction::AbstractVector, direction, + mesh::P4estMesh, equations, dg, indices...) + (; node_coordinates) = cache.elements + + x = get_node_coords(node_coordinates, equations, dg, indices...) + + u_outer = boundary_condition.boundary_value_function(boundary_condition.outer_boundary_value_function, + u_inner, + normal_direction, + x, t, equations) + + return u_outer +end end # @muladd diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 431780dab8b..5a37ab77390 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -375,7 +375,8 @@ end end @inline function idp_local_minmax!(alpha, limiter, u, t, dt, semi, - mesh::StructuredMesh{2}, elements, variable) + mesh::Union{StructuredMesh{2}, P4estMesh{2}}, + elements, variable) _, _, dg, cache = mesh_equations_solver_cache(semi) (; variable_bounds) = limiter.cache.subcell_limiter_coefficients @@ -476,7 +477,8 @@ end end @inline function idp_spec_entropy!(alpha, limiter, u, t, dt, semi, - mesh::StructuredMesh{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] @@ -527,7 +529,8 @@ end end @inline function idp_math_entropy!(alpha, limiter, u, t, dt, semi, - mesh::StructuredMesh{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] @@ -599,7 +602,8 @@ end return nothing end -@inline function idp_positivity!(alpha, limiter, u, dt, semi, mesh::StructuredMesh{2}, +@inline function idp_positivity!(alpha, limiter, u, dt, semi, + mesh::Union{StructuredMesh{2}, P4estMesh{2}}, elements, variable) _, _, dg, cache = mesh_equations_solver_cache(semi) @@ -694,7 +698,8 @@ end end @inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, - mesh::StructuredMesh{2}, elements, variable) + mesh::Union{StructuredMesh{2}, P4estMesh{2}}, + elements, variable) _, _, dg, cache = mesh_equations_solver_cache(semi) (; variable_bounds) = limiter.cache.subcell_limiter_coefficients diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index bc5e46e1931..e73e475a73f 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -269,7 +269,8 @@ function calc_normal_directions!(container_bar_states, mesh::TreeMesh, equations nothing end -function calc_normal_directions!(container_bar_states, mesh::StructuredMesh{2}, +function calc_normal_directions!(container_bar_states, + mesh::Union{StructuredMesh{2}, P4estMesh{2}}, equations, dg, cache) (; weights, derivative_matrix) = dg.basis (; contravariant_vectors) = cache.elements diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index db34aecc168..d967b945a2b 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -153,6 +153,32 @@ end end end +@trixi_testset "elixir_euler_free_stream_sc_subcell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_free_stream_sc_subcell.jl"), + l2=[ + 1.4663777294625118e-15, + 2.320054900530864e-14, + 3.487555722563465e-14, + 2.008802099296406e-14, + ], + linf=[ + 2.3092638912203256e-14, + 2.0623780461193064e-13, + 2.6795232699328153e-13, + 2.362554596402333e-13, + ], + atol=2.0e-12,) + # 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)) < 15000 + end +end + @trixi_testset "elixir_euler_shockcapturing_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing_ec.jl"), l2=[ @@ -337,6 +363,32 @@ end end end +@trixi_testset "elixir_euler_double_mach.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_double_mach.jl"), + l2=[ + 0.8741784143331414, + 6.669726935141086, + 3.4980245896042237, + 76.33557073504075, + ], + linf=[ + 11.428353668952052, + 142.73486850872337, + 38.91639544604301, + 1651.7541390872523, + ], + initial_refinement_level=1, + 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)) < 15000 + end +end + @trixi_testset "elixir_euler_supersonic_cylinder.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_supersonic_cylinder.jl"), l2=[ @@ -365,6 +417,35 @@ end end end +@trixi_testset "elixir_euler_supersonic_cylinder_sc_subcell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_supersonic_cylinder_sc_subcell.jl"), + l2=[ + 0.01733051773398731, + 0.038254257166961285, + 0.018157981470786955, + 0.12176639664229769, + ], + linf=[ + 1.3534563960399795, + 2.861333164923601, + 2.248472479406512, + 9.797432332463623, + ], + tspan=(0.0, 0.001), + skip_coverage=true) + if @isdefined sol # Skipped in coverage run + # 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)) < 15000 + end + end +end + @trixi_testset "elixir_eulergravity_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_convergence.jl"), l2=[