diff --git a/NEWS.md b/NEWS.md index 88afec79987..0b4279e9e7a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,6 +12,7 @@ for human readability. - Implementation of 1D Linearized Euler Equations. - New analysis callback for 2D `P4estMesh` to compute integrated quantities along a boundary surface, e.g., pressure lift and drag coefficients. - Optional tuple parameter for `GlmSpeedCallback` called `semi_indices` to specify for which semidiscretization of a `SemidiscretizationCoupled` we need to update the GLM speed. +- Subcell local one-sided limiting support for nonlinear variables in 2D for `TreeMesh`. ## Changes when updating to v0.7 from v0.6.x diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 8b8e49a28a8..8a98fdae283 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -106,7 +106,7 @@ positivity_variables_nonlinear = [pressure] # As for the limiting with global bounds you are passing the quantity names of the conservative # variables you want to limit. So, to limit the density with lower and upper local bounds pass # the following. -local_minmax_variables_cons = ["rho"] +local_twosided_variables_cons = ["rho"] # ## Exemplary simulation # How to set up a simulation using the IDP limiting becomes clearer when looking at an exemplary @@ -154,7 +154,7 @@ volume_flux = flux_ranocha # Here, the simulation should contain local limiting for the density using lower and upper bounds. basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; - local_minmax_variables_cons = ["rho"]) + local_twosided_variables_cons = ["rho"]) # The initialized limiter is passed to `VolumeIntegralSubcellLimiting` in addition to the volume # fluxes of the low-order and high-order scheme. diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl index 209aa2ae352..00d3c69f2e6 100644 --- a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl +++ b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl @@ -41,7 +41,9 @@ surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; - local_minmax_variables_cons = ["rho"]) + local_twosided_variables_cons = ["rho"], + local_onesided_variables_nonlinear = [(Trixi.entropy_math, + max)]) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) diff --git a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl index c1ba3d96962..6cbbe4eb4e6 100644 --- a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl @@ -42,7 +42,9 @@ surface_flux = flux_lax_friedrichs volume_flux = flux_chandrashekar basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; - local_minmax_variables_cons = ["rho"]) + local_twosided_variables_cons = ["rho"], + local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, + min)]) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) diff --git a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl index 4b606502ebe..b2d49ecbd48 100644 --- a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl +++ b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl @@ -86,8 +86,8 @@ volume_flux = flux_ranocha basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; - local_minmax_variables_cons = ["rho" * string(i) - for i in eachcomponent(equations)]) + local_twosided_variables_cons = ["rho" * string(i) + for i in eachcomponent(equations)]) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 4dbf44d29c4..ba193ab2997 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -77,22 +77,27 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi return nothing end - (; local_minmax, positivity) = limiter + (; local_twosided, positivity, local_onesided) = limiter (; output_directory) = callback variables = varnames(cons2cons, semi.equations) mkpath(output_directory) open("$output_directory/deviations.txt", "a") do f print(f, "# iter, simu_time") - if local_minmax - for v in limiter.local_minmax_variables_cons + if local_twosided + for v in limiter.local_twosided_variables_cons variable_string = string(variables[v]) print(f, ", " * variable_string * "_min, " * variable_string * "_max") end end + if local_onesided + for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear + print(f, ", " * string(variable) * "_" * string(min_or_max)) + end + end if positivity for v in limiter.positivity_variables_cons - if v in limiter.local_minmax_variables_cons + if v in limiter.local_twosided_variables_cons continue end print(f, ", " * string(variables[v]) * "_min") @@ -120,15 +125,15 @@ end @inline function finalize_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimiterIDP) - (; local_minmax, positivity) = limiter + (; local_twosided, positivity, local_onesided) = limiter (; idp_bounds_delta_global) = limiter.cache variables = varnames(cons2cons, semi.equations) println("─"^100) println("Maximum deviation from bounds:") println("─"^100) - if local_minmax - for v in limiter.local_minmax_variables_cons + if local_twosided + for v in limiter.local_twosided_variables_cons v_string = string(v) println("$(variables[v]):") println("- lower bound: ", @@ -137,9 +142,19 @@ end idp_bounds_delta_global[Symbol(v_string, "_max")]) end end + if local_onesided + for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear + variable_string = string(variable) + minmax_string = string(min_or_max) + println("$variable_string:") + println("- $minmax_string bound: ", + idp_bounds_delta_global[Symbol(variable_string, "_", + minmax_string)]) + end + end if positivity for v in limiter.positivity_variables_cons - if v in limiter.local_minmax_variables_cons + if v in limiter.local_twosided_variables_cons continue end println(string(variables[v]) * ":\n- positivity: ", diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 3a56ea71f62..0f713a296e2 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -8,7 +8,7 @@ @inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache, limiter::SubcellLimiterIDP, time, iter, output_directory, save_errors) - (; local_minmax, positivity) = solver.volume_integral.limiter + (; local_twosided, positivity, local_onesided) = solver.volume_integral.limiter (; variable_bounds) = limiter.cache.subcell_limiter_coefficients (; idp_bounds_delta_local, idp_bounds_delta_global) = limiter.cache @@ -20,8 +20,8 @@ # `@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_minmax - for v in limiter.local_minmax_variables_cons + 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") @@ -45,9 +45,28 @@ 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)) + 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 j in eachnode(solver), i in eachnode(solver) + v = variable(get_node_vars(u, equations, solver, i, j, 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, element])) + end + end + idp_bounds_delta_local[key] = deviation + end + end if positivity for v in limiter.positivity_variables_cons - if v in limiter.local_minmax_variables_cons + if v in limiter.local_twosided_variables_cons continue end key = Symbol(string(v), "_min") @@ -86,16 +105,23 @@ # Print to output file open("$output_directory/deviations.txt", "a") do f print(f, iter, ", ", time) - if local_minmax - for v in limiter.local_minmax_variables_cons + if local_twosided + for v in limiter.local_twosided_variables_cons v_string = string(v) print(f, ", ", idp_bounds_delta_local[Symbol(v_string, "_min")], ", ", idp_bounds_delta_local[Symbol(v_string, "_max")]) end end + if local_onesided + for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear + print(f, ", ", + idp_bounds_delta_local[Symbol(string(variable), "_", + string(min_or_max))][stride_size]) + end + end if positivity for v in limiter.positivity_variables_cons - if v in limiter.local_minmax_variables_cons + if v in limiter.local_twosided_variables_cons continue end print(f, ", ", idp_bounds_delta_local[Symbol(string(v), "_min")]) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index fa7af285aa4..0614066806c 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1886,6 +1886,27 @@ end return SVector(w1, w2, w3, w4) 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::CompressibleEulerEquations2D) + rho, rho_v1, rho_v2, rho_e = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_square = v1^2 + v2^2 + inv_rho_gammap1 = (1 / rho)^(equations.gamma + 1.0) + + # The derivative vector for the modified specific entropy of Guermond et al. + w1 = inv_rho_gammap1 * + (0.5 * rho * (equations.gamma + 1.0) * v_square - equations.gamma * rho_e) + w2 = -rho_v1 * inv_rho_gammap1 + w3 = -rho_v2 * inv_rho_gammap1 + w4 = (1 / rho)^equations.gamma + + return SVector(w1, w2, w3, w4) +end + @inline function entropy2cons(w, equations::CompressibleEulerEquations2D) # 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) @@ -1991,6 +2012,29 @@ end return S end +# Transformation from conservative variables u to d(s)/d(u) +@inline function gradient_conservative(::typeof(entropy_math), + u, equations::CompressibleEulerEquations2D) + return cons2entropy(u, equations) +end + +# Calculate 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 entropy_guermond_etal(u, equations::CompressibleEulerEquations2D) + rho, rho_v1, rho_v2, rho_e = u + + # Modified specific entropy from Guermond et al. (2019) + s = (rho_e - 0.5 * (rho_v1^2 + rho_v2^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::CompressibleEulerEquations2D) + return cons2entropy_guermond_etal(u, equations) +end + # Default entropy is the mathematical entropy @inline function entropy(cons, equations::CompressibleEulerEquations2D) entropy_math(cons, equations) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 4366cd32f08..622410de855 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -295,8 +295,8 @@ of local and symmetric parts. It is equivalent to the non-conservative flux of B et al. (`flux_nonconservative_powell`) for conforming meshes but it yields different results on non-conforming meshes(!). -The two other flux functions with the same name return either the local -or symmetric portion of the non-conservative flux based on the type of the +The two other flux functions with the same name return either the local +or symmetric portion of the non-conservative flux based on the type of the nonconservative_type argument, employing multiple dispatch. They are used to compute the subcell fluxes in dg_2d_subcell_limiters.jl. diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index e433c953779..17c9488316d 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -14,27 +14,36 @@ end """ SubcellLimiterIDP(equations::AbstractEquations, basis; - local_minmax_variables_cons = String[], + local_twosided_variables_cons = String[], positivity_variables_cons = String[], positivity_variables_nonlinear = [], positivity_correction_factor = 0.1, + local_onesided_variables_nonlinear = [], max_iterations_newton = 10, newton_tolerances = (1.0e-12, 1.0e-14), gamma_constant_newton = 2 * ndims(equations)) Subcell invariant domain preserving (IDP) limiting used with [`VolumeIntegralSubcellLimiting`](@ref) including: -- Local maximum/minimum Zalesak-type limiting for conservative variables (`local_minmax_variables_cons`) +- Local two-sided Zalesak-type limiting for conservative variables (`local_twosided_variables_cons`) - Positivity limiting for conservative variables (`positivity_variables_cons`) and nonlinear variables (`positivity_variables_nonlinear`) +- Local one-sided limiting for nonlinear variables, e.g. `entropy_guermond_etal` and `entropy_math` +with `local_onesided_variables_nonlinear` -Conservative variables to be limited are passed as a vector of strings, e.g. `local_minmax_variables_cons = ["rho"]` -and `positivity_variables_cons = ["rho"]`. For nonlinear variables the specific functions are -passed in a vector, e.g. `positivity_variables_nonlinear = [pressure]`. +To use these three limiting options use the following structure: + +***Conservative variables*** to be limited are passed as a vector of strings, e.g. +`local_twosided_variables_cons = ["rho"]` and `positivity_variables_cons = ["rho"]`. +For ***nonlinear variables***, the wanted variable functions are passed within a vector: To ensure +positivity use a plain vector including the desired variables, e.g. `positivity_variables_nonlinear = [pressure]`. +For local one-sided limiting pass the variable function combined with the requested bound +(`min` or `max`) as a tuple. For instance, to impose a lower local bound on the modified specific +entropy by Guermond et al. use `local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, min)]`. The bounds are calculated using the low-order FV solution. The positivity limiter uses `positivity_correction_factor` such that `u^new >= positivity_correction_factor * u^FV`. -The limiting of nonlinear variables uses a Newton-bisection method with a maximum of +Local and global limiting of nonlinear variables uses a Newton-bisection method with a maximum of `max_iterations_newton` iterations, relative and absolute tolerances of `newton_tolerances` and a provisional update constant `gamma_constant_newton` (`gamma_constant_newton>=2*d`, where `d = #dimensions`). See equation (20) of Pazner (2020) and equation (30) of Rueda-Ramírez et al. (2022). @@ -55,14 +64,17 @@ where `d = #dimensions`). See equation (20) of Pazner (2020) and equation (30) o !!! warning "Experimental implementation" This is an experimental feature and may change in future releases. """ -struct SubcellLimiterIDP{RealT <: Real, LimitingVariablesNonlinear, Cache} <: +struct SubcellLimiterIDP{RealT <: Real, LimitingVariablesNonlinear, + LimitingOnesidedVariablesNonlinear, Cache} <: AbstractSubcellLimiter - local_minmax::Bool - local_minmax_variables_cons::Vector{Int} # Local mininum/maximum principles for conservative variables + local_twosided::Bool + local_twosided_variables_cons::Vector{Int} # Local two-sided limiting for conservative variables positivity::Bool positivity_variables_cons::Vector{Int} # Positivity for conservative variables positivity_variables_nonlinear::LimitingVariablesNonlinear # Positivity for nonlinear variables positivity_correction_factor::RealT + local_onesided::Bool + local_onesided_variables_nonlinear::LimitingOnesidedVariablesNonlinear # Local one-sided limiting for nonlinear variables cache::Cache max_iterations_newton::Int newton_tolerances::Tuple{RealT, RealT} # Relative and absolute tolerances for Newton's method @@ -71,32 +83,56 @@ end # this method is used when the limiter is constructed as for shock-capturing volume integrals function SubcellLimiterIDP(equations::AbstractEquations, basis; - local_minmax_variables_cons = String[], + local_twosided_variables_cons = String[], positivity_variables_cons = String[], positivity_variables_nonlinear = [], positivity_correction_factor = 0.1, + local_onesided_variables_nonlinear = [], max_iterations_newton = 10, newton_tolerances = (1.0e-12, 1.0e-14), gamma_constant_newton = 2 * ndims(equations)) - local_minmax = (length(local_minmax_variables_cons) > 0) + local_twosided = (length(local_twosided_variables_cons) > 0) + local_onesided = (length(local_onesided_variables_nonlinear) > 0) positivity = (length(positivity_variables_cons) + length(positivity_variables_nonlinear) > 0) - local_minmax_variables_cons_ = get_variable_index.(local_minmax_variables_cons, - equations) + # When passing `min` or `max` in the elixir, the specific function of Base is used. + # To speed up the simulation, we replace it with `Trixi.min` and `Trixi.max` respectively. + local_onesided_variables_nonlinear_ = Tuple{Function, Function}[] + for (variable, min_or_max) in local_onesided_variables_nonlinear + if min_or_max === Base.max + push!(local_onesided_variables_nonlinear_, (variable, max)) + elseif min_or_max === Base.min + push!(local_onesided_variables_nonlinear_, (variable, min)) + elseif min_or_max === Trixi.max || min_or_max === Trixi.min + push!(local_onesided_variables_nonlinear_, (variable, min_or_max)) + else + error("Parameter $min_or_max is not a valid input. Use `max` or `min` instead.") + end + end + local_onesided_variables_nonlinear_ = Tuple(local_onesided_variables_nonlinear_) + + local_twosided_variables_cons_ = get_variable_index.(local_twosided_variables_cons, + equations) positivity_variables_cons_ = get_variable_index.(positivity_variables_cons, equations) bound_keys = () - if local_minmax - for v in local_minmax_variables_cons_ + if local_twosided + for v in local_twosided_variables_cons_ v_string = string(v) bound_keys = (bound_keys..., Symbol(v_string, "_min"), Symbol(v_string, "_max")) end end + if local_onesided + for (variable, min_or_max) in local_onesided_variables_nonlinear_ + bound_keys = (bound_keys..., + Symbol(string(variable), "_", string(min_or_max))) + end + end for v in positivity_variables_cons_ - if !(v in local_minmax_variables_cons_) + if !(v in local_twosided_variables_cons_) bound_keys = (bound_keys..., Symbol(string(v), "_min")) end end @@ -108,29 +144,36 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis; SubcellLimiterIDP{typeof(positivity_correction_factor), typeof(positivity_variables_nonlinear), - typeof(cache)}(local_minmax, local_minmax_variables_cons_, + typeof(local_onesided_variables_nonlinear_), + typeof(cache)}(local_twosided, local_twosided_variables_cons_, positivity, positivity_variables_cons_, positivity_variables_nonlinear, - positivity_correction_factor, cache, + positivity_correction_factor, + local_onesided, + local_onesided_variables_nonlinear_, + cache, max_iterations_newton, newton_tolerances, gamma_constant_newton) end function Base.show(io::IO, limiter::SubcellLimiterIDP) @nospecialize limiter # reduce precompilation time - (; local_minmax, positivity) = limiter + (; local_twosided, positivity, local_onesided) = limiter print(io, "SubcellLimiterIDP(") - if !(local_minmax || positivity) + if !(local_twosided || positivity || local_onesided) print(io, "No limiter selected => pure DG method") else features = String[] - if local_minmax + if local_twosided push!(features, "local min/max") end if positivity push!(features, "positivity") end + if local_onesided + push!(features, "local onesided") + end join(io, features, ", ") print(io, "Limiter=($features), ") end @@ -140,19 +183,19 @@ end function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP) @nospecialize limiter # reduce precompilation time - (; local_minmax, positivity) = limiter + (; local_twosided, positivity, local_onesided) = limiter if get(io, :compact, false) show(io, limiter) else - if !(local_minmax || positivity) - setup = ["limiter" => "No limiter selected => pure DG method"] + if !(local_twosided || positivity || local_onesided) + setup = ["Limiter" => "No limiter selected => pure DG method"] else - setup = ["limiter" => ""] - if local_minmax + setup = ["Limiter" => ""] + if local_twosided setup = [ setup..., - "" => "Local maximum/minimum limiting for conservative variables $(limiter.local_minmax_variables_cons)", + "" => "Local two-sided limiting for conservative variables $(limiter.local_twosided_variables_cons)", ] end if positivity @@ -163,6 +206,11 @@ function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP) "" => "- with positivity correction factor = $(limiter.positivity_correction_factor)", ] end + if local_onesided + for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear + setup = [setup..., "" => "Local $min_or_max limiting for $variable"] + end + end setup = [ setup..., "Local bounds" => "FV solution", diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 9343cee4397..33ae0599748 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -37,13 +37,17 @@ function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSE # TODO: Do not abuse `reset_du!` but maybe implement a generic `set_zero!` @trixi_timeit timer() "reset alpha" reset_du!(alpha, dg, semi.cache) - if limiter.local_minmax - @trixi_timeit timer() "local min/max limiting" idp_local_minmax!(alpha, limiter, - u, t, dt, semi) + if limiter.local_twosided + @trixi_timeit timer() "local twosided" idp_local_twosided!(alpha, limiter, + u, t, dt, semi) end if limiter.positivity @trixi_timeit timer() "positivity" idp_positivity!(alpha, limiter, u, dt, semi) end + if limiter.local_onesided + @trixi_timeit timer() "local onesided" idp_local_onesided!(alpha, limiter, + u, t, dt, semi) + end # Calculate alpha1 and alpha2 @unpack alpha1, alpha2 = limiter.cache.subcell_limiter_coefficients @@ -164,18 +168,121 @@ end return nothing end +@inline function calc_bounds_onesided!(var_minmax, min_or_max, variable, u, 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 j in eachnode(dg), i in eachnode(dg) + if min_or_max === max + var_minmax[i, j, element] = typemin(eltype(var_minmax)) + else + var_minmax[i, j, element] = typemax(eltype(var_minmax)) + end + end + + # Calculate bounds at Gauss-Lobatto nodes using u + for j in eachnode(dg), i in eachnode(dg) + var = variable(get_node_vars(u, equations, dg, i, j, element), equations) + var_minmax[i, j, element] = min_or_max(var_minmax[i, j, element], var) + + if i > 1 + var_minmax[i - 1, j, element] = min_or_max(var_minmax[i - 1, j, + element], var) + end + if i < nnodes(dg) + var_minmax[i + 1, j, element] = min_or_max(var_minmax[i + 1, j, + element], var) + end + if j > 1 + var_minmax[i, j - 1, element] = min_or_max(var_minmax[i, j - 1, + element], var) + end + if j < nnodes(dg) + var_minmax[i, j + 1, element] = min_or_max(var_minmax[i, j + 1, + element], var) + end + end + end + + # Values at element boundary + calc_bounds_onesided_interface!(var_minmax, min_or_max, variable, u, t, semi, mesh) +end + +@inline function calc_bounds_onesided_interface!(var_minmax, min_or_max, variable, u, t, + semi, mesh::TreeMesh2D) + _, equations, dg, cache = mesh_equations_solver_cache(semi) + (; boundary_conditions) = semi + # Calc bounds at interfaces and periodic boundaries + for interface in eachinterface(dg, cache) + # Get neighboring element ids + left = cache.interfaces.neighbor_ids[1, interface] + right = cache.interfaces.neighbor_ids[2, interface] + + orientation = cache.interfaces.orientations[interface] + + for i in eachnode(dg) + index_left = (nnodes(dg), i) + index_right = (1, i) + if orientation == 2 + index_left = reverse(index_left) + index_right = reverse(index_right) + end + var_left = variable(get_node_vars(u, equations, dg, index_left..., left), + equations) + var_right = variable(get_node_vars(u, equations, dg, index_right..., right), + equations) + + var_minmax[index_right..., right] = min_or_max(var_minmax[index_right..., + right], var_left) + var_minmax[index_left..., left] = min_or_max(var_minmax[index_left..., + left], var_right) + end + end + + # Calc bounds at physical boundaries + for boundary in eachboundary(dg, cache) + element = cache.boundaries.neighbor_ids[boundary] + orientation = cache.boundaries.orientations[boundary] + neighbor_side = cache.boundaries.neighbor_sides[boundary] + + for i in eachnode(dg) + if neighbor_side == 2 # Element is on the right, boundary on the left + index = (1, i) + boundary_index = 1 + else # Element is on the left, boundary on the right + index = (nnodes(dg), i) + boundary_index = 2 + end + if orientation == 2 + index = reverse(index) + boundary_index += 2 + end + u_outer = get_boundary_outer_state(boundary_conditions[boundary_index], + cache, t, equations, dg, + index..., element) + var_outer = variable(u_outer, equations) + + var_minmax[index..., element] = min_or_max(var_minmax[index..., element], + var_outer) + end + end + + return nothing +end + ############################################################################### -# Local minimum/maximum limiting +# Local two-sided limiting of conservative variables -@inline function idp_local_minmax!(alpha, limiter, u, t, dt, semi) - for variable in limiter.local_minmax_variables_cons - idp_local_minmax!(alpha, limiter, u, t, dt, semi, variable) +@inline function idp_local_twosided!(alpha, limiter, u, t, dt, semi) + for variable in limiter.local_twosided_variables_cons + idp_local_twosided!(alpha, limiter, u, t, dt, semi, variable) end return nothing end -@inline function idp_local_minmax!(alpha, limiter, u, t, dt, semi, variable) +@inline function idp_local_twosided!(alpha, limiter, u, t, dt, semi, variable) _, _, dg, cache = mesh_equations_solver_cache(semi) (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes (; inverse_weights) = dg.basis @@ -236,6 +343,40 @@ end return nothing end +############################################################################## +# Local one-sided limiting of nonlinear variables + +@inline function idp_local_onesided!(alpha, limiter, u, t, dt, semi) + for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear + idp_local_onesided!(alpha, limiter, u, t, dt, semi, variable, min_or_max) + end + + return nothing +end + +@inline function idp_local_onesided!(alpha, limiter, u, t, dt, semi, variable, + min_or_max) + _, 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) + inverse_jacobian = cache.elements.inverse_jacobian[element] + for j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, element) + newton_loops_alpha!(alpha, var_minmax[i, j, element], u_local, + i, j, 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 @@ -283,8 +424,8 @@ end end # Compute bound - if limiter.local_minmax && - variable in limiter.local_minmax_variables_cons && + if limiter.local_twosided && + variable in limiter.local_twosided_variables_cons && var_min[i, j, element] >= positivity_correction_factor * var # Local limiting is more restrictive that positivity limiting # => Skip positivity limiting for this node @@ -346,7 +487,7 @@ end # Perform Newton's bisection method to find new alpha newton_loops_alpha!(alpha, var_min[i, j, element], u_local, i, j, element, - variable, initial_check_nonnegative_newton_idp, + variable, min, initial_check_nonnegative_newton_idp, final_check_nonnegative_newton_idp, inverse_jacobian, dt, equations, dg, cache, limiter) end @@ -356,8 +497,9 @@ end end @inline function newton_loops_alpha!(alpha, bound, u, i, j, element, variable, - initial_check, final_check, inverse_jacobian, dt, - equations, dg, cache, limiter) + min_or_max, initial_check, final_check, + inverse_jacobian, dt, equations, dg, cache, + limiter) (; inverse_weights) = dg.basis (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes @@ -367,37 +509,38 @@ end antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[i] * get_node_vars(antidiffusive_flux1_R, equations, dg, i, j, element) - newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, - equations, dt, limiter, antidiffusive_flux) + newton_loop!(alpha, bound, u, i, j, element, variable, min_or_max, initial_check, + final_check, equations, dt, limiter, antidiffusive_flux) # positive xi direction antidiffusive_flux = -gamma_constant_newton * inverse_jacobian * inverse_weights[i] * get_node_vars(antidiffusive_flux1_L, equations, dg, i + 1, j, element) - newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, - equations, dt, limiter, antidiffusive_flux) + newton_loop!(alpha, bound, u, i, j, element, variable, min_or_max, initial_check, + final_check, equations, dt, limiter, antidiffusive_flux) # negative eta direction antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[j] * get_node_vars(antidiffusive_flux2_R, equations, dg, i, j, element) - newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, - equations, dt, limiter, antidiffusive_flux) + newton_loop!(alpha, bound, u, i, j, element, variable, min_or_max, initial_check, + final_check, equations, dt, limiter, antidiffusive_flux) # positive eta direction antidiffusive_flux = -gamma_constant_newton * inverse_jacobian * inverse_weights[j] * get_node_vars(antidiffusive_flux2_L, equations, dg, i, j + 1, element) - newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, - equations, dt, limiter, antidiffusive_flux) + newton_loop!(alpha, bound, u, i, j, element, variable, min_or_max, initial_check, + final_check, equations, dt, limiter, antidiffusive_flux) return nothing end -@inline function newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, - final_check, equations, dt, limiter, antidiffusive_flux) +@inline function newton_loop!(alpha, bound, u, i, j, element, variable, min_or_max, + initial_check, final_check, equations, dt, limiter, + antidiffusive_flux) newton_reltol, newton_abstol = limiter.newton_tolerances beta = 1 - alpha[i, j, element] @@ -411,7 +554,7 @@ end if isvalid(u_curr, equations) goal = goal_function_newton_idp(variable, bound, u_curr, equations) - initial_check(bound, goal, newton_abstol) && return nothing + initial_check(min_or_max, bound, goal, newton_abstol) && return nothing end # Newton iterations @@ -446,7 +589,7 @@ end # Check new beta for condition and update bounds goal = goal_function_newton_idp(variable, bound, u_curr, equations) - if initial_check(bound, goal, newton_abstol) + if initial_check(min_or_max, bound, goal, newton_abstol) # New beta fulfills condition beta_L = beta else @@ -479,18 +622,25 @@ end end new_alpha = 1 - beta - if alpha[i, j, element] > new_alpha + newton_abstol - error("Alpha is getting smaller. old: $(alpha[i, j, element]), new: $new_alpha") - else - alpha[i, j, element] = new_alpha - end + alpha[i, j, element] = new_alpha return nothing end ### Auxiliary routines for Newton's bisection method ### # Initial checks -@inline initial_check_nonnegative_newton_idp(bound, goal, newton_abstol) = goal <= 0 +@inline function initial_check_local_onesided_newton_idp(::typeof(min), bound, + goal, newton_abstol) + goal <= max(newton_abstol, abs(bound) * newton_abstol) +end + +@inline function initial_check_local_onesided_newton_idp(::typeof(max), bound, + goal, newton_abstol) + goal >= -max(newton_abstol, abs(bound) * newton_abstol) +end + +@inline initial_check_nonnegative_newton_idp(min_or_max, bound, goal, newton_abstol) = goal <= + 0 # Goal and d(Goal)d(u) function @inline goal_function_newton_idp(variable, bound, u, equations) = bound - @@ -501,6 +651,12 @@ end end # Final checks +# final check for one-sided local limiting +@inline function final_check_local_onesided_newton_idp(bound, goal, newton_abstol) + abs(goal) < max(newton_abstol, abs(bound) * newton_abstol) +end + +# final check for nonnegativity limiting @inline function final_check_nonnegative_newton_idp(bound, goal, newton_abstol) (goal <= eps()) && (goal > -max(newton_abstol, abs(bound) * newton_abstol)) end diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index d593dc540f7..efd4058031f 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -349,16 +349,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_sc_subcell_nonperiodic.jl"), l2=[ - 0.3517507570120483, - 0.19252291020146015, - 0.19249751956580294, - 0.618717827188004, + 0.3221177942225801, + 0.1798478357478982, + 0.1798364616438908, + 0.6136884131056267, ], linf=[ - 1.6699566795772216, - 1.3608007992899402, - 1.361864507190922, - 2.44022884092527, + 1.343766644801395, + 1.1749593109683463, + 1.1747613085307178, + 2.4216006041018785, ], tspan=(0.0, 0.5), initial_refinement_level=4, @@ -403,16 +403,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_sc_subcell.jl"), l2=[ - 0.432804941135901, - 0.15009019787510924, - 0.15009019787510922, - 0.6160764058367757, + 0.41444427153173785, + 0.1460669409661223, + 0.14606693069201596, + 0.6168046457461059, ], linf=[ - 1.6122663996643651, - 0.8612394422674909, - 0.8612394422674919, - 6.449588561676761, + 1.5720584643579567, + 0.7946656826861964, + 0.7946656525739751, + 6.455520291414711, ], tspan=(0.0, 1.0), initial_refinement_level=4, diff --git a/test/test_unit.jl b/test/test_unit.jl index c72b51113f5..90ee21030d3 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -416,8 +416,9 @@ end indicator_hg = IndicatorHennemannGassner(1.0, 0.0, true, "variable", "cache") @test_nowarn show(stdout, indicator_hg) - limiter_idp = SubcellLimiterIDP(true, [1], true, [1], ["variable"], 0.1, "cache", 1, - (1.0, 1.0), 1.0) + limiter_idp = SubcellLimiterIDP(true, [1], true, [1], ["variable"], 0.1, + true, [(Trixi.entropy_guermond_etal, min)], "cache", + 1, (1.0, 1.0), 1.0) @test_nowarn show(stdout, limiter_idp) indicator_loehner = IndicatorLöhner(1.0, "variable", (; cache = nothing))