diff --git a/NEWS.md b/NEWS.md index 3159cdcfaa5..948cc53531c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -13,6 +13,9 @@ for human readability. - Support for second-order finite volume subcell volume integral (`VolumeIntegralPureLGLFiniteVolumeO2`) and stabilized DG-FV blending volume integral (`VolumeIntegralShockCapturingRRG`) on 3D meshes for conservative systems ([#2734], [#2755]). +- Extended 3D support for subcell limiting with `P4estMesh` was added ([#2733]). + In the new version, local (minimum or maximum) limiting for nonlinear variables (using + the keyword `local_onesided_variables_nonlinear` in `SubcellLimiterIDP()`) is supported. #### Changed @@ -26,17 +29,16 @@ for human readability. ## Changes in the v0.13 lifecycle #### Added -- Initial 3D support for subcell limiting with `P4estMesh` was added ([#2582] and [#2647]). +- Initial 3D support for subcell limiting with `P4estMesh` was added ([#2582], [#2647], [#2688], [#2722]). In the new version, IDP positivity limiting for conservative variables (using the keyword `positivity_variables_cons` in `SubcellLimiterIDP()`) and nonlinear variables (using `positivity_variables_nonlinear`) is supported. - `BoundsCheckCallback` is not supported in 3D yet. - Optimized 2D and 3D kernels for nonconservative fluxes with `P4estMesh` were added ([#2653], [#2663]). The optimized kernel can be enabled via the trait `Trixi.combine_conservative_and_nonconservative_fluxes(flux, equations)`. When the trait is set to `Trixi.True()`, a single method has to be defined, that computes and returns the tuple `flux_cons(u_ll, u_rr) + 0.5f0 * flux_noncons(u_ll, u_rr)` and `flux_cons(u_ll, u_rr) + 0.5f0 * flux_noncons(u_rr, u_ll)`. -- Added support for `source_terms_parabolic`, which allows users to specify gradient-dependent source terms when solving parabolic equations ([#2721]). +- Added support for `source_terms_parabolic`, which allows users to specify gradient-dependent source terms when solving parabolic equations ([#2721]). - Support for second-order finite volume subcell volume integral (`VolumeIntegralPureLGLFiniteVolumeO2`) and stabilized DG-FV blending volume integral (`VolumeIntegralShockCapturingRRG`) on 1D & 2D meshes for conservative systems ([#2022], [#2695], [#2718]). diff --git a/examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl b/examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl index 7ffaea44ad4..721688fc6d2 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl @@ -44,7 +44,9 @@ polydeg = 3 basis = LobattoLegendreBasis(polydeg) limiter_idp = SubcellLimiterIDP(equations, basis; positivity_variables_cons = ["rho"], - positivity_variables_nonlinear = [pressure]) + positivity_variables_nonlinear = [pressure], + local_onesided_variables_nonlinear = [], + max_iterations_newton = 25) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) diff --git a/examples/p4est_3d_dgsem/elixir_mhd_shockcapturing_subcell.jl b/examples/p4est_3d_dgsem/elixir_mhd_shockcapturing_subcell.jl index 40f45c0e49b..8c15733ab8b 100644 --- a/examples/p4est_3d_dgsem/elixir_mhd_shockcapturing_subcell.jl +++ b/examples/p4est_3d_dgsem/elixir_mhd_shockcapturing_subcell.jl @@ -117,7 +117,7 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -stage_callbacks = (SubcellLimiterIDPCorrection(),) +stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback()) 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 diff --git a/src/callbacks_stage/subcell_bounds_check_3d.jl b/src/callbacks_stage/subcell_bounds_check_3d.jl index e001ad142e5..5267f061d85 100644 --- a/src/callbacks_stage/subcell_bounds_check_3d.jl +++ b/src/callbacks_stage/subcell_bounds_check_3d.jl @@ -19,6 +19,26 @@ # `@batch` here to allow a possible redefinition of `@threaded` without creating errors here. # See also https://github.com/trixi-framework/Trixi.jl/pull/1888#discussion_r1537785293. + if local_onesided + for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear + key = Symbol(string(variable), "_", string(min_or_max)) + deviation = idp_bounds_delta_local[key] + sign_ = min_or_max(1.0, -1.0) + @batch reduction=(max, deviation) for element in eachelement(solver, cache) + for k in eachnode(solver), j in eachnode(solver), i in eachnode(solver) + v = variable(get_node_vars(u, equations, solver, i, j, k, element), + equations) + # Note: We always save the absolute deviations >= 0 and therefore use the + # `max` operator for lower and upper bounds. The different directions of + # upper and lower bounds are considered with `sign_`. + deviation = max(deviation, + sign_ * + (v - variable_bounds[key][i, j, k, element])) + end + end + idp_bounds_delta_local[key] = deviation + end + end if positivity for v in limiter.positivity_variables_cons key = Symbol(string(v), "_min") diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index 0b877d466be..205c43fe3e6 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -1791,6 +1791,27 @@ end return SVector(w1, w2, w3, w4, w5) end +# Transformation from conservative variables u to entropy vector ds_0/du, +# using the modified specific entropy of Guermond et al. (2019): s_0 = p * rho^(-gamma) / (gamma-1). +# Note: This is *not* the "conventional" specific entropy s = ln(p / rho^(gamma)). +@inline function cons2entropy_guermond_etal(u, equations::CompressibleEulerEquations3D) + rho, rho_v1, rho_v2, rho_v3, rho_e = u + + inv_rho = 1 / rho + v_square_rho = (rho_v1^2 + rho_v2^2 + rho_v3^2) * inv_rho + inv_rho_gammap1 = inv_rho^(equations.gamma + 1) + + # The derivative vector for the modified specific entropy of Guermond et al. + w1 = inv_rho_gammap1 * + (0.5f0 * (equations.gamma + 1) * v_square_rho - equations.gamma * rho_e) + w2 = -rho_v1 * inv_rho_gammap1 + w3 = -rho_v2 * inv_rho_gammap1 + w4 = -rho_v3 * inv_rho_gammap1 + w5 = inv_rho^equations.gamma + + return SVector(w1, w2, w3, w4, w5) +end + @inline function entropy2cons(w, equations::CompressibleEulerEquations3D) # See Hughes, Franca, Mallet (1986) A new finite element formulation for CFD # [DOI: 10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) @@ -1880,6 +1901,34 @@ end return S end +@doc raw""" + entropy_guermond_etal(u, equations::CompressibleEulerEquations3D) + +Calculate the modified specific entropy of Guermond et al. (2019): +```math +s_0 = p * \rho^{-\gamma} / (\gamma-1). +``` +Note: This is *not* the "conventional" specific entropy ``s = ln(p / \rho^\gamma)``. +- Guermond at al. (2019) + Invariant domain preserving discretization-independent schemes and convex limiting for hyperbolic systems. + [DOI: 10.1016/j.cma.2018.11.036](https://doi.org/10.1016/j.cma.2018.11.036) +""" +@inline function entropy_guermond_etal(u, equations::CompressibleEulerEquations3D) + rho, rho_v1, rho_v2, rho_v3, rho_e = u + + # Modified specific entropy from Guermond et al. (2019) + s = (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho) * + (1 / rho)^equations.gamma + + return s +end + +# Transformation from conservative variables u to d(s)/d(u) +@inline function gradient_conservative(::typeof(entropy_guermond_etal), + u, equations::CompressibleEulerEquations3D) + return cons2entropy_guermond_etal(u, equations) +end + # Default entropy is the mathematical entropy @inline function entropy(cons, equations::CompressibleEulerEquations3D) return entropy_math(cons, equations) diff --git a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl index b5f6b729f26..a15a0173366 100644 --- a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl +++ b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl @@ -9,6 +9,187 @@ # IDP Limiting ############################################################################### +############################################################################### +# Calculation of local bounds using low-order FV solution + +@inline function calc_bounds_onesided!(var_minmax, min_or_max, variable, + u::AbstractArray{<:Any, 5}, t, semi) + mesh, equations, dg, cache = mesh_equations_solver_cache(semi) + # Calc bounds inside elements + @threaded for element in eachelement(dg, cache) + # Reset bounds + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + if min_or_max === max + var_minmax[i, j, k, element] = typemin(eltype(var_minmax)) + else + var_minmax[i, j, k, element] = typemax(eltype(var_minmax)) + end + end + + # Calculate bounds at Gauss-Lobatto nodes + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + var = variable(get_node_vars(u, equations, dg, i, j, k, element), equations) + var_minmax[i, j, k, element] = min_or_max(var_minmax[i, j, k, element], var) + + if i > 1 + var_minmax[i - 1, j, k, element] = min_or_max(var_minmax[i - 1, j, k, + element], var) + end + if i < nnodes(dg) + var_minmax[i + 1, j, k, element] = min_or_max(var_minmax[i + 1, j, k, + element], var) + end + if j > 1 + var_minmax[i, j - 1, k, element] = min_or_max(var_minmax[i, j - 1, k, + element], var) + end + if j < nnodes(dg) + var_minmax[i, j + 1, k, element] = min_or_max(var_minmax[i, j + 1, k, + element], var) + end + if k > 1 + var_minmax[i, j, k - 1, element] = min_or_max(var_minmax[i, j, k - 1, + element], var) + end + if k < nnodes(dg) + var_minmax[i, j, k + 1, element] = min_or_max(var_minmax[i, j, k + 1, + element], var) + end + end + end + + # Values at element boundary + calc_bounds_onesided_interface!(var_minmax, min_or_max, variable, u, t, semi, mesh) + + return nothing +end + +function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, semi, + mesh::P4estMesh{3}) + _, equations, dg, cache = mesh_equations_solver_cache(semi) + (; boundary_conditions) = semi + + (; neighbor_ids, node_indices) = cache.interfaces + index_range = eachnode(dg) + + # Calc bounds 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] + + # 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,k indexing + i_primary_start, i_primary_step_i, i_primary_step_j = index_to_start_step_3d(primary_indices[1], + index_range) + j_primary_start, j_primary_step_i, j_primary_step_j = index_to_start_step_3d(primary_indices[2], + index_range) + k_primary_start, k_primary_step_i, k_primary_step_j = index_to_start_step_3d(primary_indices[3], + index_range) + + i_primary = i_primary_start + j_primary = j_primary_start + k_primary = k_primary_start + + i_secondary_start, i_secondary_step_i, i_secondary_step_j = index_to_start_step_3d(secondary_indices[1], + index_range) + j_secondary_start, j_secondary_step_i, j_secondary_step_j = index_to_start_step_3d(secondary_indices[2], + index_range) + k_secondary_start, k_secondary_step_i, k_secondary_step_j = index_to_start_step_3d(secondary_indices[3], + index_range) + + i_secondary = i_secondary_start + j_secondary = j_secondary_start + k_secondary = k_secondary_start + + for j in eachnode(dg) + for i in eachnode(dg) + var_primary = variable(get_node_vars(u, equations, dg, i_primary, + j_primary, k_primary, + primary_element), equations) + var_secondary = variable(get_node_vars(u, equations, dg, i_secondary, + j_secondary, k_secondary, + secondary_element), + equations) + + var_minmax[i_primary, j_primary, k_primary, primary_element] = minmax(var_minmax[i_primary, + j_primary, + k_primary, + primary_element], + var_secondary) + var_minmax[i_secondary, j_secondary, k_secondary, secondary_element] = minmax(var_minmax[i_secondary, + j_secondary, + k_secondary, + secondary_element], + var_primary) + + # Increment the primary element indices + i_primary += i_primary_step_i + j_primary += j_primary_step_i + k_primary += k_primary_step_i + # Increment the secondary element surface indices + i_secondary += i_secondary_step_i + j_secondary += j_secondary_step_i + k_secondary += k_secondary_step_i + end + # Increment the primary element indices + i_primary += i_primary_step_j + j_primary += j_primary_step_j + k_primary += k_primary_step_j + # Increment the secondary element surface indices + i_secondary += i_secondary_step_j + j_secondary += j_secondary_step_j + k_secondary += k_secondary_step_j + end + end + + # Calc bounds at physical boundaries + calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t, + boundary_conditions, + mesh, equations, dg, cache) + + return nothing +end + +@inline function calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t, + boundary_conditions::BoundaryConditionPeriodic, + mesh::P4estMesh{3}, + equations, dg, cache) + return nothing +end + +############################################################################## +# Local one-sided limiting of nonlinear variables + +@inline function idp_local_onesided!(alpha, limiter, u::AbstractArray{<:Real, 5}, + t, dt, semi, + variable, min_or_max) + mesh, equations, dg, cache = mesh_equations_solver_cache(semi) + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + var_minmax = variable_bounds[Symbol(string(variable), "_", string(min_or_max))] + calc_bounds_onesided!(var_minmax, min_or_max, variable, u, t, semi) + + # Perform Newton's bisection method to find new alpha + @threaded for element in eachelement(dg, cache) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + inverse_jacobian = get_inverse_jacobian(cache.elements.inverse_jacobian, + mesh, i, j, k, element) + u_local = get_node_vars(u, equations, dg, i, j, k, element) + newton_loops_alpha!(alpha, var_minmax[i, j, k, element], + u_local, i, j, k, element, + variable, min_or_max, + initial_check_local_onesided_newton_idp, + final_check_local_onesided_newton_idp, + inverse_jacobian, dt, equations, dg, cache, limiter) + end + end + + return nothing +end + ############################################################################### # Global positivity limiting of conservative variables diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 12de969bf2f..f29256f056b 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -371,21 +371,50 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end -@trixi_testset "elixir_euler_sedov_sc_subcell.jl" begin +@trixi_testset "elixir_euler_sedov_sc_subcell.jl (positivity bounds)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_sc_subcell.jl"), l2=[ - 0.1942700455652903, - 0.07557644365785855, - 0.07557644365785836, - 0.07557644365785698, - 0.3713893635249306 + 0.19427117014566905, + 0.07557679851688888, + 0.07557679851688896, + 0.07557679851688917, + 0.3713893373335463 ], linf=[ - 2.7542157588958798, - 1.8885700263691245, - 1.888570026369125, - 1.8885700263691252, - 4.9712792944452096 + 2.754292081408987, + 1.888627709320322, + 1.8886277093203232, + 1.888627709320322, + 4.971280431903264 + ], + tspan=(0.0, 0.3),) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from + # OrdinaryDiffEq.jl + # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 + @test_allocations(Trixi.rhs!, semi, sol, 15_000) +end + +@trixi_testset "elixir_euler_sedov_sc_subcell.jl (local bounds)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_sc_subcell.jl"), + local_onesided_variables_nonlinear=[(entropy_guermond_etal, + min)], + max_iterations_newton=40, + l2=[ + 0.19153085678321066, + 0.07411109384422779, + 0.07411109384410808, + 0.07411109384406232, + 0.36714268468314665 + ], + linf=[ + 1.4037775549639524, + 1.339590863739464, + 1.3395908637591605, + 1.3395908637371077, + 4.824252924073932 ], tspan=(0.0, 0.3)) # Ensure that we do not have excessive memory allocations