indicators.jl

# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin
#! format: noindent

# Abstract supertype of indicators used for AMR, shock capturing, and
# adaptive volume-integral selection
abstract type AbstractIndicator end

function create_cache(typ::Type{IndicatorType},
                      semi) where {IndicatorType <: AbstractIndicator}
    return create_cache_indicator_for_amr(typ, mesh_equations_solver_cache(semi)...)
end

# this method is used when the indicator is constructed as for AMR
function create_cache_indicator_for_amr(typ::Type{IndicatorType},
                                        mesh, equations::AbstractEquations, dg::DGSEM,
                                        cache) where {IndicatorType <:
                                                      AbstractIndicator}
    return create_cache(typ, equations, dg.basis)
end

function get_element_variables!(element_variables, indicator::AbstractIndicator,
                                ::AbstractVolumeIntegralShockCapturing)
    element_variables[:indicator_shock_capturing] = indicator.cache.alpha
    return nothing
end

"""
    IndicatorHennemannGassner(equations::AbstractEquations, basis;
                              alpha_max=0.5,
                              alpha_min=0.001,
                              alpha_smooth=true,
                              variable)
    IndicatorHennemannGassner(semi::AbstractSemidiscretization;
                              alpha_max=0.5,
                              alpha_min=0.001,
                              alpha_smooth=true,
                              variable)

Indicator used for shock-capturing (when passing the `equations` and the `basis`)
or adaptive mesh refinement (AMR, when passing the `semi`).

See also [`VolumeIntegralShockCapturingHG`](@ref).

## References

- Hennemann, Gassner (2020)
  "A provably entropy stable subcell shock capturing approach for high order split form DG"
  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
"""
struct IndicatorHennemannGassner{RealT <: Real, Variable, Cache} <: AbstractIndicator
    alpha_max::RealT
    alpha_min::RealT
    alpha_smooth::Bool
    variable::Variable
    cache::Cache
end

# this method is used when the indicator is constructed as for shock-capturing volume integrals
function IndicatorHennemannGassner(equations::AbstractEquations, basis;
                                   alpha_max = 0.5,
                                   alpha_min = 0.001,
                                   alpha_smooth = true,
                                   variable)
    alpha_max, alpha_min = promote(alpha_max, alpha_min)
    cache = create_cache(IndicatorHennemannGassner, equations, basis)
    return IndicatorHennemannGassner{typeof(alpha_max), typeof(variable),
                                     typeof(cache)}(alpha_max,
                                                    alpha_min,
                                                    alpha_smooth,
                                                    variable,
                                                    cache)
end

# this method is used when the indicator is constructed as for AMR
function IndicatorHennemannGassner(semi::AbstractSemidiscretization;
                                   alpha_max = 0.5,
                                   alpha_min = 0.001,
                                   alpha_smooth = true,
                                   variable)
    alpha_max, alpha_min = promote(alpha_max, alpha_min)
    cache = create_cache(IndicatorHennemannGassner, semi)
    return IndicatorHennemannGassner{typeof(alpha_max), typeof(variable),
                                     typeof(cache)}(alpha_max,
                                                    alpha_min,
                                                    alpha_smooth,
                                                    variable,
                                                    cache)
end

function Base.show(io::IO, indicator::IndicatorHennemannGassner)
    @nospecialize indicator # reduce precompilation time

    print(io, "IndicatorHennemannGassner(")
    print(io, indicator.variable)
    print(io, ", alpha_max=", indicator.alpha_max)
    print(io, ", alpha_min=", indicator.alpha_min)
    print(io, ", alpha_smooth=", indicator.alpha_smooth)
    print(io, ")")
    return nothing
end

function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorHennemannGassner)
    @nospecialize indicator # reduce precompilation time
    setup = [
        "indicator variable" => indicator.variable,
        "max. α" => indicator.alpha_max,
        "min. α" => indicator.alpha_min,
        "smooth α" => (indicator.alpha_smooth ? "yes" : "no")
    ]
    summary_box(io, "IndicatorHennemannGassner", setup)
    return nothing
end

function (indicator_hg::IndicatorHennemannGassner)(u, mesh, equations, dg::DGSEM, cache;
                                                   kwargs...)
    @unpack alpha_smooth = indicator_hg
    @unpack alpha, alpha_tmp = indicator_hg.cache
    # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR?
    #       Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)`
    #       or just `resize!` whenever we call the relevant methods as we do now?
    resize!(alpha, nelements(dg, cache))
    if alpha_smooth
        resize!(alpha_tmp, nelements(dg, cache))
    end

    # magic parameters
    # TODO: Are there better values for Float32?
    RealT = real(dg)
    threshold = 0.5f0 * 10^(convert(RealT, -1.8) * nnodes(dg)^convert(RealT, 0.25))
    o_0001 = convert(RealT, 0.0001)
    parameter_s = log((1 - o_0001) / o_0001)

    @threaded for element in eachelement(dg, cache)
        # This is dispatched by mesh dimension.
        # Use this function barrier and unpack inside to avoid passing closures to
        # Polyester.jl with `@batch` (`@threaded`).
        # Otherwise, `@threaded` does not work here with Julia ARM on macOS.
        # See https://github.com/JuliaSIMD/Polyester.jl/issues/88.
        calc_indicator_hennemann_gassner!(indicator_hg, threshold, parameter_s, u,
                                          element, mesh, equations, dg, cache)
    end

    if alpha_smooth
        apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache)
    end

    return alpha
end

"""
    IndicatorLöhner (equivalent to IndicatorLoehner)

    IndicatorLöhner(equations::AbstractEquations, basis;
                    f_wave=0.2, variable)
    IndicatorLöhner(semi::AbstractSemidiscretization;
                    f_wave=0.2, variable)

AMR indicator adapted from a FEM indicator by Löhner (1987), also used in the
FLASH code as standard AMR indicator.
The indicator estimates a weighted second derivative of a specified variable locally.

When constructed to be used for AMR, pass the `semi`. Pass the `equations`,
and `basis` if this indicator should be used for shock capturing.

## References

- Löhner (1987)
  "An adaptive finite element scheme for transient problems in CFD"
  [doi: 10.1016/0045-7825(87)90098-3](https://doi.org/10.1016/0045-7825(87)90098-3)
- [https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p62/node59.html#SECTION05163100000000000000](https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p62/node59.html#SECTION05163100000000000000)
"""
struct IndicatorLöhner{RealT <: Real, Variable, Cache} <: AbstractIndicator
    f_wave::RealT # TODO: Taal documentation
    variable::Variable
    cache::Cache
end

# this method is used when the indicator is constructed as for shock-capturing volume integrals
function IndicatorLöhner(equations::AbstractEquations, basis;
                         f_wave = 0.2, variable)
    cache = create_cache(IndicatorLöhner, equations, basis)
    return IndicatorLöhner{typeof(f_wave), typeof(variable), typeof(cache)}(f_wave,
                                                                            variable,
                                                                            cache)
end

# this method is used when the indicator is constructed as for AMR
function IndicatorLöhner(semi::AbstractSemidiscretization;
                         f_wave = 0.2, variable)
    cache = create_cache(IndicatorLöhner, semi)
    return IndicatorLöhner{typeof(f_wave), typeof(variable), typeof(cache)}(f_wave,
                                                                            variable,
                                                                            cache)
end

function Base.show(io::IO, indicator::IndicatorLöhner)
    @nospecialize indicator # reduce precompilation time

    print(io, "IndicatorLöhner(")
    print(io, "f_wave=", indicator.f_wave, ", variable=", indicator.variable, ")")
    return nothing
end

function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorLöhner)
    @nospecialize indicator # reduce precompilation time

    if get(io, :compact, false)
        show(io, indicator)
    else
        setup = [
            "indicator variable" => indicator.variable,
            "f_wave" => indicator.f_wave
        ]
        summary_box(io, "IndicatorLöhner", setup)
    end
end

"""
    IndicatorLoehner

Same as [`IndicatorLöhner`](@ref).
"""
const IndicatorLoehner = IndicatorLöhner

# dirty Löhner estimate, direction by direction, assuming constant nodes
@inline function local_löhner_estimate(um::Real, u0::Real, up::Real,
                                       löhner::IndicatorLöhner)
    num = abs(up - 2 * u0 + um)
    den = abs(up - u0) + abs(u0 - um) +
          löhner.f_wave * (abs(up) + 2 * abs(u0) + abs(um))
    return num / den
end

"""
    IndicatorMax(equations::AbstractEquations, basis; variable)
    IndicatorMax(semi::AbstractSemidiscretization; variable)

A simple indicator returning the maximum of `variable` in an element.
When constructed to be used for AMR, pass the `semi`. Pass the `equations`,
and `basis` if this indicator should be used for shock capturing.

If an AMR indicator depending not only on a solution variable but also the
space and time is desired, consider using [`IndicatorNodalFunction`](@ref)
instead.
"""
struct IndicatorMax{Variable, Cache <: NamedTuple} <: AbstractIndicator
    variable::Variable
    cache::Cache
end

# this method is used when the indicator is constructed as for shock-capturing volume integrals
function IndicatorMax(equations::AbstractEquations, basis;
                      variable)
    cache = create_cache(IndicatorMax, equations, basis)
    return IndicatorMax{typeof(variable), typeof(cache)}(variable, cache)
end

# this method is used when the indicator is constructed as for AMR
function IndicatorMax(semi::AbstractSemidiscretization;
                      variable)
    cache = create_cache(IndicatorMax, semi)
    return IndicatorMax{typeof(variable), typeof(cache)}(variable, cache)
end

function Base.show(io::IO, indicator::IndicatorMax)
    @nospecialize indicator # reduce precompilation time

    print(io, "IndicatorMax(")
    print(io, "variable=", indicator.variable, ")")
    return nothing
end

function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorMax)
    @nospecialize indicator # reduce precompilation time

    if get(io, :compact, false)
        show(io, indicator)
    else
        setup = [
            "indicator variable" => indicator.variable
        ]
        summary_box(io, "IndicatorMax", setup)
    end
end

@doc raw"""
    IndicatorEntropyChange(; maximum_entropy_increase::Real = 0.0)

This indicator checks the difference in mathematical [`entropy`](@ref) (``S``) due to the application
of a volume integral (VI) compared to the true/analytical entropy evolution
(without any dissipation inside the element).
In particular, the indicator computes
```math
\Delta S = \dot{S}_\mathrm{VI} - \dot{S}_\text{true} =
\int_{\Omega_m}
\frac{\partial S}{\partial \boldsymbol{u}} \cdot \dot{\boldsymbol u}_\mathrm{VI}
\mathrm{d} \Omega_m
-
\int_{\partial \Omega_m}
\boldsymbol{\psi} \cdot \hat{\boldsymbol{n}}
\mathrm{d} \partial \Omega_m
```
for the currently processed element/cell ``m``.
Here, ``\dot{\boldsymbol u}_\mathrm{VI}`` is the change in the DG right-hand-side due to the volume integral only.
``\dot{S}_\text{true}`` is the true entropy evolution, which can be computed from the
entropy potential ``\boldsymbol{\psi}`` (see also [`entropy_potential`](@ref)).

This is discussed in more detail in
- Chen, Shu (2017)
  "Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules for hyperbolic conservation laws"
  [DOI: 10.1016/j.jcp.2017.05.025](https://doi.org/10.1016/j.jcp.2017.05.025)
- Lin, Chan (2024)
  "High order entropy stable discontinuous Galerkin spectral element methods through subcell limiting"
  [DOI: 10.1016/j.jcp.2023.112677](https://doi.org/10.1016/j.jcp.2023.112677)

For ``\Delta S < \sigma \leq 0`` with ``\sigma`` being set to `maximum_entropy_increase`,
the e.g. [`VolumeIntegralWeakForm`](@ref) is more entropy-diffusive than the true entropy change
(which could be recovered with the [`VolumeIntegralFluxDifferencing`](@ref) and an
entropy-conservative flux such as [`flux_ranocha`](@ref), for instance).

If ``\sigma > 0`` is set, i.e., `maximum_entropy_increase > 0`, the indicator allows for
limited entropy increase, thereby allowing to use e.g. the cheaper weak-form volume integral
even in slightly entropy-producing situations to reduce computational cost.

Supposed to be used in conjunction with [`VolumeIntegralAdaptive`](@ref) which then selects
a volume integral for every cell/element ``m``.

The logic behind this indicator is similar to the "companion" scheme
approach proposed in Chapter 5 of

- Carpenter, Fisher, Nielsen, and Frankel (2014)
  "Entropy Stable Spectral Collocation Schemes for the Navier-Stokes Equations: Discontinuous Interfaces"
  [DOI: 10.1137/130932193](https://doi.org/10.1137/130932193)

Here, we thus equip e.g. the flux-differencing volume integral with a "companion" weak-form
volume integral.
However, usage of the entropy potential allows for comparison with the true entropy change.

!!! note
    This indicator is **not implemented as an AMR indicator**, i.e., it is **not
    possible** to employ this as the `indicator` in [`ControllerThreeLevel`](@ref),
    for instance.
"""
struct IndicatorEntropyChange{RealT <: Real} <:
       AbstractIndicator
    maximum_entropy_increase::RealT

    function IndicatorEntropyChange(; maximum_entropy_increase = 0.0)
        return new{typeof(maximum_entropy_increase)}(maximum_entropy_increase)
    end
end

function Base.show(io::IO, indicator::IndicatorEntropyChange)
    @nospecialize indicator # reduce precompilation time

    print(io, "IndicatorEntropyChange(")
    print(io, "maximum_entropy_increase=", indicator.maximum_entropy_increase, ")")

    return nothing
end

function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorEntropyChange)
    @nospecialize indicator # reduce precompilation time

    if get(io, :compact, false)
        show(io, indicator)
    else
        setup = [
            "maximum_entropy_increase" => indicator.maximum_entropy_increase
        ]
        summary_box(io, "IndicatorEntropyChange", setup)
    end
end

"""
    IndicatorEntropyCorrection(equations::AbstractEquations, basis;
                               scaling=true)

Indicator used for entropy correction using subcell FV schemes, where the
blending is determined so that the volume integral entropy production is the
same or more than that of an entropy-conservative (EC) scheme.

This is intended to guide the convex blending of a `volume_integral_default`
(for example, [`VolumeIntegralWeakForm`](@ref)) and `volume_integral_stabilized`
(for example, [`VolumeIntegralPureLGLFiniteVolume`](@ref) with an entropy stable
finite volume flux).

The parameter `scaling ≥ 1` in [`IndicatorEntropyCorrection`](@ref) scales the DG-FV blending
parameter ``\\alpha``(see the [tutorial on shock-capturing](https://trixi-framework.github.io/TrixiDocumentation/stable/tutorials/shock_capturing/#Shock-capturing-with-flux-differencing))
by a constant, increasing the amount of the subcell FV added in (up to 1, i.e., pure subcell FV).
This can be used to add shock capturing-like behavior. Note though that ``\\alpha`` is computed
here from the entropy defect, **not** using [`IndicatorHennemannGassner`](@ref).

The use of `IndicatorEntropyCorrection` requires either
`entropy_potential(u, orientation, equations)` for TreeMesh, or
`entropy_potential(u, normal_direction, equations)` for other mesh types
to be defined.

"""
struct IndicatorEntropyCorrection{Cache, ScalingT} <: AbstractIndicator
    cache::Cache
    scaling::ScalingT # either Bool or Real
end

# this method is used when the indicator is constructed as for shock-capturing volume integrals
function IndicatorEntropyCorrection(equations::AbstractEquations,
                                    basis::LobattoLegendreBasis;
                                    scaling = true) # true = 1 in floating point multiplication
    cache = create_cache(IndicatorEntropyCorrection, equations, basis)
    return IndicatorEntropyCorrection{typeof(cache), typeof(scaling)}(cache, scaling)
end

# this method is used when the indicator is constructed as for
# shock-capturing volume integrals.
function create_cache(::Type{IndicatorEntropyCorrection},
                      equations::AbstractEquations{NDIMS, NVARS},
                      basis::LobattoLegendreBasis) where {NDIMS, NVARS}
    uEltype = real(basis)
    AT = Array{uEltype, NDIMS + 1}

    # container for elementwise volume integrals
    volume_integral_values_threaded = AT[AT(undef, NVARS,
                                            ntuple(_ -> nnodes(basis), NDIMS)...)
                                         for _ in 1:Threads.maxthreadid()]

    # stores the blending coefficients
    alpha = Vector{uEltype}()

    return (; alpha, volume_integral_values_threaded)
end

function Base.show(io::IO, indicator::IndicatorEntropyCorrection)
    @nospecialize indicator # reduce precompilation time
    print(io, "IndicatorEntropyCorrection")
    return nothing
end

function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorEntropyCorrection)
    @nospecialize indicator # reduce precompilation time
    summary_box(io, "IndicatorEntropyCorrection")
    return nothing
end

"""
    IndicatorEntropyCorrectionShockCapturingCombined(; indicator_shock_capturing,
                                                       indicator_entropy_correction)

Indicator used for entropy correction using subcell FV schemes, where the blending
is taken to be the maximum between a blending determined by shock capturing
(`indicator_shock_capturing`) and a blending determined so that the volume integral
entropy production is the same or more than that of an EC scheme (`indicator_entropy_correction`).

This is intended to guide the convex blending of a `volume_integral_default` (for
example, [`VolumeIntegralWeakForm`](@ref)) and `volume_integral_stabilized` (for
example, [`VolumeIntegralPureLGLFiniteVolume`](@ref) with an entropy stable finite
volume flux).

The use of `IndicatorEntropyCorrectionShockCapturingCombined` requires either
`entropy_potential(u, orientation, equations)` for TreeMesh, or
`entropy_potential(u, normal_direction, equations)` for other mesh types
to be defined.
"""
struct IndicatorEntropyCorrectionShockCapturingCombined{IndicatorEC, IndicatorSC} <:
       AbstractIndicator
    indicator_entropy_correction::IndicatorEC
    indicator_shock_capturing::IndicatorSC
end

function IndicatorEntropyCorrectionShockCapturingCombined(; indicator_shock_capturing,
                                                          indicator_entropy_correction)
    return IndicatorEntropyCorrectionShockCapturingCombined(indicator_entropy_correction,
                                                            indicator_shock_capturing)
end

function Base.show(io::IO, indicator::IndicatorEntropyCorrectionShockCapturingCombined)
    @nospecialize indicator # reduce precompilation time
    print(io, "IndicatorEntropyCorrectionShockCapturingCombined(")
    print(io, indicator.indicator_entropy_correction)
    print(io, ", ")
    print(io, indicator.indicator_shock_capturing |> typeof |> nameof)
    print(io, ")")
    return nothing
end

function Base.show(io::IO, ::MIME"text/plain",
                   indicator::IndicatorEntropyCorrectionShockCapturingCombined)
    @nospecialize indicator # reduce precompilation time

    if get(io, :compact, false)
        show(io, indicator)
    else
        setup = [
            "indicator EC" => indicator.indicator_entropy_correction,
            "indicator SC" => indicator.indicator_shock_capturing |> typeof |> nameof
        ]
        summary_box(io, "IndicatorEntropyCorrectionShockCapturingCombined", setup)
    end
    return nothing
end

"""
    IndicatorNodalFunction(f)

Create an AMR indicator from a solution, space and time dependent indicator function `f(u, x, t)`.
The function `f` is evaluated at the nodal points. The maximum of `f` over all nodes in each element is used as indicator for the element.
The function can be solution independent allowing for user-defined mesh refinement/coarsening varying in space and time, similar to the `refinement_patches` keyword for the [`TreeMesh`](@ref).

If the function `f` depends only on the solution, you can also use the [`IndicatorMax`](@ref) instead.
"""
struct IndicatorNodalFunction{F, Cache} <: AbstractIndicator
    indicator_function::F
    cache::Cache
end

# this method is used when the indicator is constructed as for AMR
function IndicatorNodalFunction(indicator_function, semi::AbstractSemidiscretization)
    cache = create_cache(IndicatorNodalFunction, semi)
    return IndicatorNodalFunction{typeof(indicator_function), typeof(cache)}(indicator_function,
                                                                             cache)
end

# this method is directly used when the indicator is constructed as for shock-capturing volume integrals
# and by the dimension-independent method called for AMR
function create_cache(::Type{IndicatorNodalFunction},
                      equations::AbstractEquations, basis::LobattoLegendreBasis)
    uEltype = real(basis)
    alpha = Vector{uEltype}()

    return (; alpha)
end

function Base.show(io::IO, indicator::IndicatorNodalFunction)
    @nospecialize indicator # reduce precompilation time
    print(io, "IndicatorNodalFunction")
    print(io, indicator.indicator_function)
    return nothing
end

function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorNodalFunction)
    @nospecialize indicator # reduce precompilation time
    if get(io, :compact, false)
        show(io, indicator)
    else
        setup = [
            "Indicator Function" => indicator.indicator_function
        ]
        summary_box(io, "IndicatorNodalFunction", setup)
    end
end
end # @muladd