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)