Feature: Simple clamp AMR indicator.

src/Trixi.jl

@@ -222,7 +222,7 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
 export load_mesh, load_time
 
 export ControllerThreeLevel, ControllerThreeLevelCombined,
-       IndicatorLöhner, IndicatorLoehner, IndicatorMax,
+       IndicatorLöhner, IndicatorLoehner, IndicatorMax, IndicatorClamp
        IndicatorNeuralNetwork, NeuralNetworkPerssonPeraire, NeuralNetworkRayHesthaven, NeuralNetworkCNN
 
 export PositivityPreservingLimiterZhangShu

src/solvers/dgsem_tree/indicators.jl

@@ -226,6 +226,52 @@ function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorMax)
   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
+
 """
     IndicatorNeuralNetwork
 

src/solvers/dgsem_tree/indicators_1d.jl

@@ -201,6 +201,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, indicator_threaded, 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,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::Variable = 0.0
+    for i in eachnode(dg)
+      u_local = get_node_vars(u, equations, dg, i, element)
+      mean += indicator_clamp.variable(u_local, equations) * weights[i]*0.5
+    end
+
+    if indicator_clamp.min <= mean <= indicator_clamp.max
+      alpha[element] =  1.0
+    else
+      alpha[element] = -1.0
+    end
+  end
+
+  return alpha
+end
+
 # this method is used when the indicator is constructed as for shock-capturing volume integrals
 # empty cache is default
 function create_cache(::Type{<:IndicatorNeuralNetwork},
@@ -402,4 +436,4 @@ function (indicator_ann::IndicatorNeuralNetwork{NeuralNetworkRayHesthaven})(
   return alpha
 end
 
-end # @muladd
+end # @muladd

src/solvers/dgsem_tree/indicators_2d.jl

@@ -229,6 +229,39 @@ 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, indicator_threaded, 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::Variable = 0.0
+    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]*0.25
+    end
+
+    if indicator_clamp.min <= mean <= indicator_clamp.max
+      alpha[element] =  1.0
+    else
+      alpha[element] = -1.0
+    end
+  end
+
+  return alpha
+end
+
 # this method is used when the indicator is constructed as for shock-capturing volume integrals
 # empty cache is default
 function create_cache(::Type{IndicatorNeuralNetwork},

src/solvers/dgsem_tree/indicators_3d.jl

@@ -246,5 +246,38 @@ 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, indicator_threaded, 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,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::Variable = 0.0
+    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]*0.125
+    end
+
+    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

@@ -408,6 +408,9 @@ isdir(outdir) && rm(outdir, recursive=true)
     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)
+
     equations = CompressibleEulerEquations2D(1.4)
     basis = LobattoLegendreBasis(3)
     indicator_neuralnetwork = IndicatorNeuralNetwork(