How to use Trixi Adaptive Mesh Refinement Indicators?

Specific Domains
1

I have followed the tutorial at 15 Adaptive mesh refinement · Trixi.jl 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
image721×596 54.3 KB

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
image1786×892 37.2 KB

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)
2

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

See

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
Screenshot from 2025-10-02 09-11-251084×940 78.1 KB

3

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

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.

4

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

function target(u, equations)
  # compute something from the current state
end

and pass this to the indicators.

1 Like