# How to use Trixi Adaptive Mesh Refinement Indicators?

**URL:** <https://discourse.julialang.org/t/how-to-use-trixi-adaptive-mesh-refinement-indicators/132747>\
**Category:** Specific Domains\
**Tags:** trixi\
**Created:** [September 29, 2025, 1:59pm UTC](https://discourse.julialang.org/t/how-to-use-trixi-adaptive-mesh-refinement-indicators/132747 "2025-09-29T13:59:05Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![mra](https://avatars.discourse-cdn.com/v4/letter/m/d78d45/32.png) [@mra](https://discourse.julialang.org/u/mra)\
**Post date:** [September 29, 2025, 1:59pm UTC](https://discourse.julialang.org/t/how-to-use-trixi-adaptive-mesh-refinement-indicators/132747/1 "2025-09-29T13:59:05Z")

</div>

I have followed the tutorial at [15 Adaptive mesh refinement · Trixi.jl](https://trixi-framework.github.io/TrixiDocumentation/stable/tutorials/adaptive_mesh_refinement/) and have tried different indicators on this problem: `IndicatorMax`, `IndicatorLöhner` and `IndicatorHennemannGassner`.

`IndicatorLöhner` does the opposite of the refinement i.e. if refines anything but the interesting region.

`IndicatorHennemannGassner` does provide acceptable refinement on the tutorial problem, but fails in my own problem, where again it refines everything but the regions that should be refined.

What is the proper way of using those indicators so that they automatically detect the regions that need refinements (i.e. regions where refinements will improve the error the most)?

Modified tutorial code and plot:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/0/40b82bb1693224a108ad86881f69f23c96e2b325.png)

```julia
using Trixi
using OrdinaryDiffEq

advection_velocity = (0.2, -0.7)
equations = LinearScalarAdvectionEquation2D(advection_velocity)

# IDK why, but it is not provided in the tutorial
function initial_condition_gauss(x, t, equations::LinearScalarAdvectionEquation2D)
    u = exp(-1 * (x[1]^2 + x[2]^2))
    return SVector(u)
end

initial_condition = initial_condition_gauss
solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)

coordinates_min = (-5.0, -5.0)
coordinates_max = (5.0, 5.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level = 4,
                n_cells_max = 30_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

tspan = (0.0, 10.0)
ode = semidiscretize(semi, tspan);

#amr_indicator = IndicatorMax(semi, variable = first)
#amr_indicator = IndicatorHennemannGassner(semi, alpha_max=0.5, alpha_min=0.001, alpha_smooth=true, variable=first)
amr_indicator = IndicatorLoehner(semi, variable=first)

amr_controller = ControllerThreeLevel(semi, amr_indicator,
    base_level = 4,
    med_level = 5, med_threshold = 0.1,
    max_level = 6, max_threshold = 0.6
)

amr_callback = AMRCallback(semi, amr_controller,
                           interval = 5,
                           adapt_initial_condition = true,
                           adapt_initial_condition_only_refine = true)

stepsize_callback = StepsizeCallback(cfl = 0.9)

callbacks = CallbackSet(amr_callback);

ode_solver = Tsit5()
#ode_solver = CarpenterKennedy2N54(williamson_condition = false) # needs stepsize_callback

sol = solve(ode, ode_solver;
            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            ode_default_options()..., callback = callbacks);

using Plots
pd = PlotData2D(sol)
plot(pd)
plot!(getmesh(pd))

```

My problem - plot and code:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/d/5da65203d1cacc827df86b11993e93824cf81cd2.png)

```julia
using Trixi
using OrdinaryDiffEq
using Statistics
using Plots

function sigmoid(x)
    return 1.0 / (1.0 + exp(-x))
end

function ans_rect_pulse(x, t, a; s=16.0)
    if a*t < x <= a*t + 1
        # rising edge
        z = (x - a*t) - 1/2 # range from -1/2 to 1/2
        return sigmoid(z * s)
    elseif a*t + 1 < x <= a*t + 2
        # top
        return 1.0
    elseif a*t + 2 < x <= a*t + 3
        # falling edge
        z = (x - a*t) - 1/2 - 2 # range from -1/2 to 1/2
        return sigmoid(-z * s)
    else
        return 0.0
    end
end

function plot_ans(ans, t; x=0:0.01:10π)
    y = ans.(x, t)
    plot(x, y, label="analytical", xlabel="x", ylabel="u", title="Analytical Solution", size=(1200, 600), lw=1.5, color=:blue)
end

ode_solver = Tsit5()

tol = 1e-15
N = 4
K = 2^5

domain = (min = 0.0, max = 15)
tspan = (0.0, 1.0)

a = 5 # Advection velocity
ans_func = ans_rect_pulse
fig = plot_ans((x, t) -> ans_func(x, t, a), tspan[end]; x=domain.min:0.01:domain.max)

equations = LinearScalarAdvectionEquation1D(a)
coordinates_min = domain.min
coordinates_max = domain.max
mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level=Int(log2(K)), n_cells_max=10_000, periodicity=false)
solver = DGSEM(polydeg=N, surface_flux=flux_lax_friedrichs)

function initial_condition(x, t, equations::LinearScalarAdvectionEquation1D)
    u = ans_func.(x, t, a)
    return SVector(u)
end

function boundary_condition(x, t, equations::LinearScalarAdvectionEquation1D)
    u = ans_func(x[1], t, a)
    return SVector(u)
end

boundary_conditions = (
    x_neg = BoundaryConditionDirichlet(boundary_condition),
    #x_pos = BoundaryConditionDirichlet(boundary_condition)
    x_pos = Trixi.BoundaryConditionDoNothing(), # it seems that the for LinearScalarAdvectionEquation1D x_pos is ignored anyway
)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)
ode_trixi = semidiscretize(semi, tspan)

#amr_indicator = IndicatorMax(semi, variable = first)
#amr_indicator = IndicatorLoehner(semi, variable=first)
#amr_indicator = IndicatorHennemannGassner(semi, variable=first)
amr_indicator = IndicatorHennemannGassner(semi, alpha_max=0.5, alpha_min=0.001, alpha_smooth=true, variable=first)

Klog = Int(log2(K))
amr_controller = ControllerThreeLevel(semi, amr_indicator,
                                  base_level = Klog,
                                  med_level = Klog + 1, med_threshold = 0.1,
                                  max_level = Klog + 2, max_threshold = 0.3)

amr_callback = AMRCallback(semi, amr_controller,
                       interval = 1,
                       adapt_initial_condition = true,
                       adapt_initial_condition_only_refine = true)

summary_callback = SummaryCallback()
callbacks = CallbackSet(summary_callback, amr_callback)
callbacks_benchmark = CallbackSet()
    
@time sol_trixi = solve(ode_trixi, ode_solver; abstol=tol, reltol=tol, callback=callbacks)

pd = PlotData1D(sol_trixi)
plot!(fig, getmesh(pd), label=false, lw=0.5)
plot!(fig, pd; size=(1200, 600), label="Trixi", legend=:topright, color=:red, lw=1.5)
display(fig)

```

---

<div class="post-metadata">

**Author:** ![DanDoe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dandoe/32/52717_2.png) [@DanDoe](https://discourse.julialang.org/u/DanDoe)\
**Post date:** [October 2, 2025, 7:11am UTC](https://discourse.julialang.org/t/how-to-use-trixi-adaptive-mesh-refinement-indicators/132747/2 "2025-10-02T07:11:52Z")

</div>

`IndicatorMax` should give for this problem a decent result, once correctly set up.

See

```julia
using Trixi
using OrdinaryDiffEqLowStorageRK

advection_velocity = (0.2, -0.7)
equations = LinearScalarAdvectionEquation2D(advection_velocity)

initial_condition = initial_condition_gauss
solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)

coordinates_min = (-5.0, -5.0)
coordinates_max = (5.0, 5.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level = 4,
                n_cells_max = 30_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

tspan = (0.0, 10.0)
ode = semidiscretize(semi, tspan);

amr_indicator = IndicatorMax(semi, variable = first)

amr_controller = ControllerThreeLevel(semi, amr_indicator,
    base_level = 4,
    med_level = 5, med_threshold = 0.3,
    max_level = 6, max_threshold = 0.6
)

summary_callback = SummaryCallback()

# To see the grid distribution
analysis_callback = AnalysisCallback(semi, interval = 50,
                                      # Turn off error computation
                                     analysis_errors = Symbol[],
                                     analysis_integrals = ())

amr_callback = AMRCallback(semi, amr_controller,
                           interval = 5)

stepsize_callback = StepsizeCallback(cfl = 0.9)

callbacks = CallbackSet(summary_callback, analysis_callback,
                        amr_callback, stepsize_callback);

ode_solver = CarpenterKennedy2N54(williamson_condition = false) # needs stepsize_callback

sol = solve(ode, ode_solver;
            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            ode_default_options()..., callback = callbacks);

using Plots
pd = PlotData2D(sol)
plot(sol)
plot!(getmesh(pd))

```

which gives me

 ![Screenshot from 2025-10-02 09-11-25](https://global.discourse-cdn.com/julialang/original/3X/9/d/9d907d6b3fc4e70bedda086e8662eb24b81c6ae6.png)

---

<div class="post-metadata">

**Author:** ![mra](https://avatars.discourse-cdn.com/v4/letter/m/d78d45/32.png) [@mra](https://discourse.julialang.org/u/mra)\
**Post date:** [October 3, 2025, 7:54am UTC](https://discourse.julialang.org/t/how-to-use-trixi-adaptive-mesh-refinement-indicators/132747/3 "2025-10-03T07:54:10Z")

</div>

I am aware that `IndicatorMax` can be used to obtain the good result in this case (it is showed as an example in Trixi documentation), but that is not my issue.

I am specifically interested in

> [@mra](#):
>
> What is the proper way of using those indicators so that they automatically detect the regions that need refinements (i.e. regions where refinements will improve the error the most)?

`IndicatorMax` is basically manual indicator and it will only work on very simple problems. I am mostly interested in making those other indicators work since they promise they can automatically detect regions needing refinements, but given my minimal example they do not seem to work even on toy problems.

---

<div class="post-metadata">

**Author:** ![Benedict](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/benedict/32/208694_2.png) [@Benedict](https://discourse.julialang.org/u/Benedict)\
**Post date:** [October 9, 2025, 12:19pm UTC](https://discourse.julialang.org/t/how-to-use-trixi-adaptive-mesh-refinement-indicators/132747/4 "2025-10-09T12:19:28Z")

</div>

Identifying the most relevant regions is certainly the crucial part in adaptive schemes. For elliptic PDEs there are estimators for the error in the PDE solution, or even in derived cost functionals, with provable effectiveness and efficiency. However, for hyperbolic systems, such theoretical guarantees are, to my knowledge, still lacking and error indicators are mostly heuristic.

In Trixi.jl, `IndicatorMax` simply monitors the value of a chosen variable. In the example, the peak of the bell-shaped function is therefore nicely captured. `IndicatorLoehner` calculates a normalized, approximate second derivative of a chosen variable. The intended behavior can be seen in the original paper, where the front of some blast wave is tracked. In the example, it will therefore not capture the peak. (Regarding the refinement observed far away from the bell, my guess is that it is due to near-zero values of the observed quantity, which enter in the normalization.) `IndicatorHennemannGassner` is more specialized and dedicated to shock capturing.

All indicators are still generic in the sense that you have to pick or compute an appropriate quantity, which will naturally be problem-specific. A common approach would be to define your own variable

```julia-auto
function target(u, equations)
  # compute something from the current state
end

```

and pass this to the indicators.