Feature: Simple clamp AMR indicator.

src/Trixi.jl

@@ -271,7 +271,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
 

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

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(typ::Type{IndicatorClamp}, mesh, equations::AbstractEquations{1},
+                      dg::DGSEM, cache)
+    cache = create_cache(typ, 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

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(typ::Type{IndicatorClamp}, mesh, equations::AbstractEquations{2},
+                      dg::DGSEM, cache)
+    cache = create_cache(typ, 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

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(typ::Type{IndicatorClamp}, mesh, equations::AbstractEquations{3},
+                      dg::DGSEM, cache)
+    cache = create_cache(typ, 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

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