diff --git a/src/Trixi.jl b/src/Trixi.jl index d9b9245918e..39532b4c167 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -272,7 +272,7 @@ export load_mesh, load_time, load_timestep, load_timestep!, load_dt, load_adaptive_time_integrator! export ControllerThreeLevel, ControllerThreeLevelCombined, - IndicatorLöhner, IndicatorLoehner, IndicatorMax + IndicatorLöhner, IndicatorLoehner, IndicatorMax, IndicatorClamp export PositivityPreservingLimiterZhangShu diff --git a/src/solvers/dgsem_tree/indicators.jl b/src/solvers/dgsem_tree/indicators.jl index 04b3bc98aec..5d11e3484a4 100644 --- a/src/solvers/dgsem_tree/indicators.jl +++ b/src/solvers/dgsem_tree/indicators.jl @@ -259,4 +259,55 @@ function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorMax) summary_box(io, "IndicatorMax", setup) end end + +""" + IndicatorClamp(equations::AbstractEquations, basis; min=0.0, max=1.0, variable) + IndicatorClamp(semi::AbstractSemidiscretization; min=0.0, max=1.0, variable) + +A simple indicator returning 1.0 when the element average of `variable` lies within [min,max]. +Returns -1.0 otherwise. +""" +struct IndicatorClamp{RealT <: Real, Variable, Cache <: NamedTuple} <: AbstractIndicator + min::RealT + max::RealT + variable::Variable + cache::Cache +end + +function IndicatorClamp(equations::AbstractEquations, basis; min = 0.0, max = 1.0, + variable) + cache = create_cache(IndicatorClamp, equations, basis) + IndicatorClamp{typeof(min), typeof(variable), typeof(cache)}(min, max, variable, + cache) +end + +function IndicatorClamp(semi::AbstractSemidiscretization; min = 0.0, max = 1.0, + variable) + cache = create_cache(IndicatorClamp, semi) + return IndicatorClamp{typeof(min), typeof(variable), typeof(cache)}(min, max, + variable, cache) +end + +function Base.show(io::IO, indicator::IndicatorClamp) + @nospecialize indicator # reduce precompilation time + + print(io, "IndicatorClamp(") + print(io, "min=", indicator.min, ", max=", indicator.max, ", variable=", + indicator.variable, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorClamp) + @nospecialize indicator # reduce precompilation time + + if get(io, :compact, false) + show(io, indicator) + else + setup = [ + "indicator variable" => indicator.variable, + "min" => indicator.min, + "max" => indicator.max + ] + summary_box(io, "IndicatorClamp", setup) + end +end end # @muladd diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index 9d11d05edcf..a3d0ed3a5d1 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -196,4 +196,40 @@ function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 3}, return alpha end + +function create_cache(::Type{IndicatorClamp}, equations::AbstractEquations{1}, + basis::LobattoLegendreBasis) + alpha = Vector{real(basis)}() + return (; alpha, basis.weights) +end + +function create_cache(type::Type{IndicatorClamp}, mesh, equations::AbstractEquations{1}, + dg::DGSEM, cache) + cache = create_cache(type, equations, dg.basis) +end + +function (indicator_clamp::IndicatorClamp)(u::AbstractArray{<:Any, 3}, + mesh, equations, dg::DGSEM, cache; + kwargs...) + @unpack alpha, weights = indicator_clamp.cache + resize!(alpha, nelements(dg, cache)) + + @threaded for element in eachelement(dg, cache) + mean = zero(real(dg.basis)) + + for i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, element) + mean += indicator_clamp.variable(u_local, equations) * weights[i] + end + mean *= 0.5 # Divide by reference element length + + if indicator_clamp.min <= mean <= indicator_clamp.max + alpha[element] = 1.0 + else + alpha[element] = -1.0 + end + end + + return alpha +end end # @muladd diff --git a/src/solvers/dgsem_tree/indicators_2d.jl b/src/solvers/dgsem_tree/indicators_2d.jl index 6ef65c4e9da..4f9cfd5633f 100644 --- a/src/solvers/dgsem_tree/indicators_2d.jl +++ b/src/solvers/dgsem_tree/indicators_2d.jl @@ -229,4 +229,41 @@ function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 4}, return alpha end + +function create_cache(::Type{IndicatorClamp}, equations::AbstractEquations{2}, + basis::LobattoLegendreBasis) + alpha = Vector{real(basis)}() + return (; alpha, basis.weights) +end + +function create_cache(type::Type{IndicatorClamp}, mesh, equations::AbstractEquations{2}, + dg::DGSEM, cache) + cache = create_cache(type, equations, dg.basis) +end + +function (indicator_clamp::IndicatorClamp)(u::AbstractArray{<:Any, 4}, + mesh, equations, dg::DGSEM, cache; + kwargs...) + @unpack alpha, weights = indicator_clamp.cache + resize!(alpha, nelements(dg, cache)) + + @threaded for element in eachelement(dg, cache) + mean = zero(real(dg.basis)) + + for j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, element) + mean += indicator_clamp.variable(u_local, equations) * weights[i] * + weights[j] + end + mean *= 0.25 # Divide by reference element area + + if indicator_clamp.min <= mean <= indicator_clamp.max + alpha[element] = 1.0 + else + alpha[element] = -1.0 + end + end + + return alpha +end end # @muladd diff --git a/src/solvers/dgsem_tree/indicators_3d.jl b/src/solvers/dgsem_tree/indicators_3d.jl index 6d4ee63618a..49123345714 100644 --- a/src/solvers/dgsem_tree/indicators_3d.jl +++ b/src/solvers/dgsem_tree/indicators_3d.jl @@ -255,4 +255,41 @@ function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 5}, return alpha end + +function create_cache(::Type{IndicatorClamp}, equations::AbstractEquations{3}, + basis::LobattoLegendreBasis) + alpha = Vector{real(basis)}() + return (; alpha, basis.weights) +end + +function create_cache(type::Type{IndicatorClamp}, mesh, equations::AbstractEquations{3}, + dg::DGSEM, cache) + cache = create_cache(type, equations, dg.basis) +end + +function (indicator_clamp::IndicatorClamp)(u::AbstractArray{<:Any, 5}, + mesh, equations, dg::DGSEM, cache; + kwargs...) + @unpack alpha, weights = indicator_clamp.cache + resize!(alpha, nelements(dg, cache)) + + @threaded for element in eachelement(dg, cache) + mean = zero(real(dg.basis)) + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, k, element) + mean += indicator_clamp.variable(u_local, equations) * weights[i] * + weights[j] * weights[k] + end + mean *= 0.125 # Divide by reference element volume + + if indicator_clamp.min <= mean <= indicator_clamp.max + alpha[element] = 1.0 + else + alpha[element] = -1.0 + end + end + + return alpha +end end # @muladd diff --git a/test/test_unit.jl b/test/test_unit.jl index 5831122ffe2..7fd17d1fad6 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -433,6 +433,9 @@ end indicator_max = IndicatorMax("variable", (; cache = nothing)) @test_nowarn show(stdout, indicator_max) + + indicator_clamp = IndicatorClamp(0.0, 1.0, "variable", (; cache = nothing)) + @test_nowarn show(stdout, indicator_clamp) end @timed_testset "LBM 2D constructor" begin