diff --git a/NEWS.md b/NEWS.md index 2570d6fad24..edddc607a42 100644 --- a/NEWS.md +++ b/NEWS.md @@ -14,19 +14,22 @@ This is useful for (locally) diffusion-dominated problems. This enables in particular adaptive mesh refinement for that solver-mesh combination ([#2712]). - Added functionality to `ScalarPlotData2D` allowing visualization a field provided by a user-defined scalar function ([#2796]). - Added `NonIdealCompressibleEuler2D` ([#2768]). -- Generalization of `VolumeIntegralShockCapturingHG` and `VolumeIntegralShockCapturingRRG` to support different volume integrals on the +- Generalization of `VolumeIntegralShockCapturingHG` and `VolumeIntegralShockCapturingRRG` to support different volume integrals on the non-stabilized and stabilized elements/cells. The generalized volume integral is called `VolumeIntegralShockCapturingHGType` and takes the three keyword arguments `volume_integral_default`, `volume_integral_blend_high_order`, and `volume_integral_blend_low_order` besides the usual `indicator` argument. In particular, `volume_integral_default` may be e.g. `VolumeIntegralWeakForm` or `VolumeIntegralAdaptive`, i.e., the non-stabilized elements/cells are no longer restricted to `VolumeIntegralFluxDifferencing` only ([#2802]). -- Added `IndicatorEntropyCorrection`. When combined with `VolumeIntegralAdaptive`, this blends together a stabilized and non-stabilized - volume integral based on the violation of a volume entropy condition. `IndicatorEntropyCorrectionShockCapturingCombined` additionally - guides the blending by taking the maximum of the entropy correction indicator and a shock capturing indicator ([#2764]). +- Added `IndicatorEntropyCorrection`. When combined with `VolumeIntegralAdaptive`, this blends together a stabilized and non-stabilized + volume integral based on the violation of a volume entropy condition. `IndicatorEntropyCorrectionShockCapturingCombined` additionally + guides the blending by taking the maximum of the entropy correction indicator and a shock capturing indicator ([#2764]). - The second-order subcell volume integral is no longer limited to reconstruction in primitive variables. Instead, it is possible to reconstruct in custom variables, if functions `cons2recon` and `recon2cons` are provided to `VolumeIntegralPureLGLFiniteVolumeO2` and `VolumeIntegralShockCapturingRRG`([#2817]). - Add Legendre-Gauss basis for DGSEM and implement solver (`VolumeIntegralWeakForm` and `SurfaceIntegralWeakForm` only) support for conforming 1D & 2D `TreeMesh`es ([#1965]). +- Extended 3D support for subcell limiting with `P4estMesh` was added ([#2763]). + In the new version, local (minimum and maximum) limiting for conservative variables (using the + keyword `local_twosided_variables_cons` in `SubcellLimiterIDP()`) is supported. ## Changes when updating to v0.15 from v0.14.x 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 572b0e0de20..210d0d7a0a0 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl @@ -45,6 +45,7 @@ basis = LobattoLegendreBasis(polydeg) limiter_idp = SubcellLimiterIDP(equations, basis; positivity_variables_cons = ["rho"], positivity_variables_nonlinear = [pressure], + local_twosided_variables_cons = [], local_onesided_variables_nonlinear = [], max_iterations_newton = 25) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index c03f701bda1..4a7461b4d12 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -65,6 +65,8 @@ end if positivity for v in limiter.positivity_variables_cons + # Note: If a variable appears here and in the local min/max limiting, the positivity + # lower bound is taken into account there. Skip these variables here. if v in limiter.local_twosided_variables_cons continue end diff --git a/src/callbacks_stage/subcell_bounds_check_3d.jl b/src/callbacks_stage/subcell_bounds_check_3d.jl index d2abb99cee8..52693102cde 100644 --- a/src/callbacks_stage/subcell_bounds_check_3d.jl +++ b/src/callbacks_stage/subcell_bounds_check_3d.jl @@ -19,6 +19,33 @@ # `@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_twosided + for v in limiter.local_twosided_variables_cons + v_string = string(v) + key_min = Symbol(v_string, "_min") + key_max = Symbol(v_string, "_max") + deviation_min = idp_bounds_delta_local[key_min] + deviation_max = idp_bounds_delta_local[key_max] + @batch reduction=((max, deviation_min), (max, deviation_max)) for element in eachelement(solver, + cache) + for k in eachnode(solver), j in eachnode(solver), i in eachnode(solver) + var = u[v, i, j, k, element] + # Note: We always save the absolute deviations >= 0 and therefore use the + # `max` operator for the lower and upper bound. The different directions of + # upper and lower bound are considered in their calculations with a + # different sign. + deviation_min = max(deviation_min, + variable_bounds[key_min][i, j, k, element] - + var) + deviation_max = max(deviation_max, + var - + variable_bounds[key_max][i, j, k, element]) + end + end + idp_bounds_delta_local[key_min] = deviation_min + idp_bounds_delta_local[key_max] = deviation_max + end + end if local_onesided for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear key = Symbol(string(variable), "_", string(min_or_max)) @@ -41,6 +68,11 @@ end if positivity for v in limiter.positivity_variables_cons + # Note: If a variable appears here and in the local min/max limiting, the positivity + # lower bound is taken into account there. Skip these variables here. + if v in limiter.local_twosided_variables_cons + continue + end key = Symbol(string(v), "_min") deviation = idp_bounds_delta_local[key] @batch reduction=(max, deviation) for element in eachelement(solver, cache) diff --git a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl index a15a0173366..4b1e10da09a 100644 --- a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl +++ b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl @@ -12,10 +12,171 @@ ############################################################################### # Calculation of local bounds using low-order FV solution +@inline function calc_bounds_twosided!(var_min, var_max, variable, + u::AbstractArray{<:Any, 5}, t, semi, equations) + mesh, _, dg, cache = mesh_equations_solver_cache(semi) + # Calc bounds inside elements + @threaded for element in eachelement(dg, cache) + # Calculate bounds at Gauss-Lobatto nodes + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + var = u[variable, i, j, k, element] + var_min[i, j, k, element] = var + var_max[i, j, k, element] = var + end + + # Apply values in x direction + for k in eachnode(dg), j in eachnode(dg), i in 2:nnodes(dg) + var = u[variable, i - 1, j, k, element] + var_min[i, j, k, element] = min(var_min[i, j, k, element], var) + var_max[i, j, k, element] = max(var_max[i, j, k, element], var) + + var = u[variable, i, j, k, element] + var_min[i - 1, j, k, element] = min(var_min[i - 1, j, k, element], var) + var_max[i - 1, j, k, element] = max(var_max[i - 1, j, k, element], var) + end + + # Apply values in y direction + for k in eachnode(dg), j in 2:nnodes(dg), i in eachnode(dg) + var = u[variable, i, j - 1, k, element] + var_min[i, j, k, element] = min(var_min[i, j, k, element], var) + var_max[i, j, k, element] = max(var_max[i, j, k, element], var) + + var = u[variable, i, j, k, element] + var_min[i, j - 1, k, element] = min(var_min[i, j - 1, k, element], var) + var_max[i, j - 1, k, element] = max(var_max[i, j - 1, k, element], var) + end + + # Apply values in z direction + for k in 2:nnodes(dg), j in eachnode(dg), i in eachnode(dg) + var = u[variable, i, j, k - 1, element] + var_min[i, j, k, element] = min(var_min[i, j, k, element], var) + var_max[i, j, k, element] = max(var_max[i, j, k, element], var) + + var = u[variable, i, j, k, element] + var_min[i, j, k - 1, element] = min(var_min[i, j, k - 1, element], var) + var_max[i, j, k - 1, element] = max(var_max[i, j, k - 1, element], var) + end + end + + # Values at element boundary + calc_bounds_twosided_interface!(var_min, var_max, variable, + u, t, semi, mesh, equations) + return nothing +end + +function calc_bounds_twosided_interface!(var_min, var_max, 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 = u[variable, i_primary, j_primary, k_primary, + primary_element] + var_secondary = u[variable, i_secondary, j_secondary, k_secondary, + secondary_element] + + var_min[i_primary, j_primary, k_primary, primary_element] = min(var_min[i_primary, + j_primary, + k_primary, + primary_element], + var_secondary) + var_max[i_primary, j_primary, k_primary, primary_element] = max(var_max[i_primary, + j_primary, + k_primary, + primary_element], + var_secondary) + + var_min[i_secondary, j_secondary, k_secondary, secondary_element] = min(var_min[i_secondary, + j_secondary, + k_secondary, + secondary_element], + var_primary) + var_max[i_secondary, j_secondary, k_secondary, secondary_element] = max(var_max[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_twosided_boundary!(var_min, var_max, variable, u, t, + boundary_conditions, + mesh, equations, dg, cache) + + return nothing +end + +@inline function calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t, + boundary_conditions::BoundaryConditionPeriodic, + mesh::P4estMesh{3}, + equations, dg, cache) + return nothing +end + @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 + + # The approach used in `calc_bounds_twosided!` is not used here because it requires more + # evaluations of the variable and is therefore slower. + @threaded for element in eachelement(dg, cache) # Reset bounds for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) @@ -161,6 +322,75 @@ end return nothing end +############################################################################### +# Local two-sided limiting of conservative variables + +@inline function idp_local_twosided!(alpha, limiter, u::AbstractArray{<:Any, 5}, + t, dt, semi, variable) + mesh, equations, dg, cache = mesh_equations_solver_cache(semi) + (; antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R) = cache.antidiffusive_fluxes + (; inverse_weights) = dg.basis + + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + variable_string = string(variable) + var_min = variable_bounds[Symbol(variable_string, "_min")] + var_max = variable_bounds[Symbol(variable_string, "_max")] + calc_bounds_twosided!(var_min, var_max, variable, u, t, semi, equations) + + @threaded for element in eachelement(dg, semi.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) + var = u[variable, i, j, k, element] + # Real Zalesak type limiter + # * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids" + # * Kuzmin et al. (2010). "Failsafe flux limiting and constrained data projections for equations of gas dynamics" + # Note: The Zalesak limiter has to be computed, even if the state is valid, because the correction is + # for each interface, not each node + + Qp = max(0, (var_max[i, j, k, element] - var) / dt) + Qm = min(0, (var_min[i, j, k, element] - var) / dt) + + # Calculate Pp and Pm + # Note: Boundaries of antidiffusive_flux1/2 are constant 0, so they make no difference here. + val_flux1_local = inverse_weights[i] * + antidiffusive_flux1_R[variable, i, j, k, element] + val_flux1_local_ip1 = -inverse_weights[i] * + antidiffusive_flux1_L[variable, i + 1, j, k, element] + val_flux2_local = inverse_weights[j] * + antidiffusive_flux2_R[variable, i, j, k, element] + val_flux2_local_jp1 = -inverse_weights[j] * + antidiffusive_flux2_L[variable, i, j + 1, k, element] + val_flux3_local = inverse_weights[k] * + antidiffusive_flux3_R[variable, i, j, k, element] + val_flux3_local_jp1 = -inverse_weights[k] * + antidiffusive_flux3_L[variable, i, j, k + 1, element] + + Pp = max(0, val_flux1_local) + max(0, val_flux1_local_ip1) + + max(0, val_flux2_local) + max(0, val_flux2_local_jp1) + + max(0, val_flux3_local) + max(0, val_flux3_local_jp1) + Pm = min(0, val_flux1_local) + min(0, val_flux1_local_ip1) + + min(0, val_flux2_local) + min(0, val_flux2_local_jp1) + + min(0, val_flux3_local) + min(0, val_flux3_local_jp1) + + Pp = inverse_jacobian * Pp + Pm = inverse_jacobian * Pm + + # Compute blending coefficient avoiding division by zero + # (as in paper of [Guermond, Nazarov, Popov, Thomas] (4.8)) + Qp = abs(Qp) / + (abs(Pp) + eps(typeof(Qp)) * 100 * abs(var_max[i, j, k, element])) + Qm = abs(Qm) / + (abs(Pm) + eps(typeof(Qm)) * 100 * abs(var_max[i, j, k, element])) + + # Calculate alpha at nodes + alpha[i, j, k, element] = max(alpha[i, j, k, element], 1 - min(1, Qp, Qm)) + end + end + + return nothing +end + ############################################################################## # Local one-sided limiting of nonlinear variables @@ -214,6 +444,13 @@ end end # Compute bound + if limiter.local_twosided && + (variable in limiter.local_twosided_variables_cons) && + (var_min[i, j, k, element] >= positivity_correction_factor * var) + # Local limiting is more restrictive that positivity limiting + # => Skip positivity limiting for this node + continue + end var_min[i, j, k, element] = positivity_correction_factor * var # Real one-sided Zalesak-type limiter diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index f29256f056b..b5908bc8b03 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -387,7 +387,7 @@ end 1.888627709320322, 4.971280431903264 ], - tspan=(0.0, 0.3),) + 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 @@ -399,22 +399,23 @@ 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_twosided_variables_cons=["rho"], local_onesided_variables_nonlinear=[(entropy_guermond_etal, min)], - max_iterations_newton=40, + max_iterations_newton=30, l2=[ - 0.19153085678321066, - 0.07411109384422779, - 0.07411109384410808, - 0.07411109384406232, - 0.36714268468314665 + 0.16504564013491585, + 0.06461384162458203, + 0.06461384162461223, + 0.06461384162461678, + 0.36193245790622036 ], linf=[ - 1.4037775549639524, - 1.339590863739464, - 1.3395908637591605, - 1.3395908637371077, - 4.824252924073932 + 0.9138327077620716, + 0.5707102472596818, + 0.5707102472739252, + 0.5707102472781822, + 4.777595503303726 ], tspan=(0.0, 0.3)) # Ensure that we do not have excessive memory allocations