Skip to content

Add positivity-limiter to the AMRCallback #2396

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 31 commits into
base: main
Choose a base branch
from
Open

Add positivity-limiter to the AMRCallback #2396

Show file tree
Hide file tree
Changes from 14 commits
Commits
Show all changes
31 commits
Select commit Hold shift + click to select a range
7b9a40d
add positivity limiter to the amr callback
patrickersing May 9, 2025
318bfad
remove duplicate limiter! argument
patrickersing May 9, 2025
f98cd62
Merge branch 'main' into amr_limiter
patrickersing May 9, 2025
e221e40
Update src/callbacks_step/amr.jl
patrickersing May 10, 2025
5faaf4d
Merge branch 'main' into amr_limiter
DanielDoehring Jun 28, 2025
d9b6c90
Merge branch 'main' into amr_limiter
bennibolm Jul 28, 2025
544ac69
Implement synchronized shifting for child elements
bennibolm Aug 29, 2025
24d8df0
Merge branch 'main' into amr_limiter
patrickersing Aug 29, 2025
2b7dc80
Merge branch 'main' into amr_limiter
patrickersing Aug 29, 2025
009496c
Merge branch 'main' into amr_limiter
bennibolm Sep 3, 2025
f214aea
Use `compute_u_mean`; Update docstring
bennibolm Sep 3, 2025
296d883
Use `compute_u_mean` in limiter after coarsening
bennibolm Sep 3, 2025
5605be8
Rename function
bennibolm Sep 3, 2025
15a8ff2
Implement suggestions
bennibolm Sep 3, 2025
0e82f4d
Apply suggestions from code review
DanielDoehring Sep 3, 2025
26abb91
Apply suggestions from code review
DanielDoehring Sep 3, 2025
ff35db1
Apply suggestions from code review
DanielDoehring Sep 3, 2025
c32a2e1
Apply suggestions from code review
DanielDoehring Sep 3, 2025
03d2936
Update src/callbacks_stage/positivity_zhang_shu_dg1d.jl
DanielDoehring Sep 3, 2025
07a1e0c
Apply suggestions from code review
DanielDoehring Sep 3, 2025
9861aee
Update src/callbacks_stage/positivity_zhang_shu_dg2d.jl
DanielDoehring Sep 3, 2025
29ab55a
Update src/callbacks_stage/positivity_zhang_shu_dg1d.jl
DanielDoehring Sep 3, 2025
00ca066
Update src/callbacks_stage/positivity_zhang_shu_dg1d.jl
DanielDoehring Sep 3, 2025
a5485f3
Apply suggestions from code review
DanielDoehring Sep 3, 2025
fdc1337
Rename for-variables; Remove comments
bennibolm Sep 4, 2025
b035070
Remove codes as suggested
bennibolm Sep 5, 2025
e042600
add comment for amr positivity-limiting
patrickersing Sep 8, 2025
6852972
Merge branch 'main' into amr_limiter
patrickersing Sep 8, 2025
7b8d2ab
Apply limiter after refinement/coarsening for T8codeMesh (#50)
bennibolm Sep 13, 2025
5c1a1a2
doc
DanielDoehring Sep 26, 2025
6d3cbfb
Merge branch 'main' into PatrickErsingLimitedMortars
DanielDoehring Sep 26, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 8 additions & 6 deletions examples/p4est_2d_dgsem/elixir_euler_double_mach_amr.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -131,6 +131,11 @@ save_solution = SaveSolutionCallback(interval = 100,
save_final_solution = true,
solution_variables = cons2prim)

# positivity limiter necessary for this example with strong shocks
positivity_limiter = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
variables = (Trixi.density,
pressure))

amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)

amr_controller = ControllerThreeLevel(semi, amr_indicator,
Expand All @@ -141,18 +146,15 @@ amr_controller = ControllerThreeLevel(semi, amr_indicator,
amr_callback = AMRCallback(semi, amr_controller,
interval = 1,
adapt_initial_condition = true,
adapt_initial_condition_only_refine = true)
adapt_initial_condition_only_refine = true,
limiter! = positivity_limiter)

callbacks = CallbackSet(summary_callback,
analysis_callback, alive_callback,
save_solution,
amr_callback)

# positivity limiter necessary for this example with strong shocks
stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
variables = (Trixi.density, pressure))

###############################################################################
# run the simulation
sol = solve(ode, SSPRK43(stage_limiter!);
sol = solve(ode, SSPRK43(stage_limiter! = positivity_limiter);
ode_default_options()..., callback = callbacks);
14 changes: 8 additions & 6 deletions examples/p4est_2d_dgsem/elixir_euler_forward_step_amr.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,6 +136,11 @@ save_solution = SaveSolutionCallback(interval = 2000,
save_final_solution = true,
solution_variables = cons2prim)

# positivity limiter necessary for this example with strong shocks
positivity_limiter = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
variables = (Trixi.density,
pressure))

amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)

amr_controller = ControllerThreeLevel(semi, amr_indicator,
Expand All @@ -146,19 +151,16 @@ amr_controller = ControllerThreeLevel(semi, amr_indicator,
amr_callback = AMRCallback(semi, amr_controller,
interval = 5,
adapt_initial_condition = true,
adapt_initial_condition_only_refine = true)
adapt_initial_condition_only_refine = true,
limiter! = positivity_limiter)

callbacks = CallbackSet(summary_callback,
analysis_callback, alive_callback,
save_solution,
amr_callback)

# positivity limiter necessary for this example with strong shocks
stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
variables = (Trixi.density, pressure))

###############################################################################
# run the simulation
sol = solve(ode, SSPRK43(stage_limiter!);
sol = solve(ode, SSPRK43(stage_limiter! = positivity_limiter);
maxiters = 999999, ode_default_options()...,
callback = callbacks);
16 changes: 9 additions & 7 deletions examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,12 @@ save_solution = SaveSolutionCallback(interval = 1000,
save_final_solution = true,
solution_variables = cons2prim)

# positivity limiter necessary for this example with strong shocks. Very sensitive
# to the order of the limiter variables, pressure must come first.
positivity_limiter = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-7, 1.0e-6),
variables = (pressure,
Trixi.density))

amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)

amr_controller = ControllerThreeLevel(semi, amr_indicator,
Expand All @@ -126,19 +132,15 @@ amr_controller = ControllerThreeLevel(semi, amr_indicator,
amr_callback = AMRCallback(semi, amr_controller,
interval = 1,
adapt_initial_condition = true,
adapt_initial_condition_only_refine = true)
adapt_initial_condition_only_refine = true,
limiter! = positivity_limiter)

callbacks = CallbackSet(summary_callback,
analysis_callback, alive_callback,
save_solution,
amr_callback)

# positivity limiter necessary for this example with strong shocks. Very sensitive
# to the order of the limiter variables, pressure must come first.
stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-7, 1.0e-6),
variables = (pressure, Trixi.density))

###############################################################################
# run the simulation
sol = solve(ode, SSPRK43(stage_limiter!);
sol = solve(ode, SSPRK43(stage_limiter! = positivity_limiter);
ode_default_options()..., callback = callbacks);
16 changes: 11 additions & 5 deletions examples/tree_1d_dgsem/elixir_euler_positivity.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,16 +82,23 @@ save_solution = SaveSolutionCallback(interval = 100,
save_final_solution = true,
solution_variables = cons2prim)

positivity_limiter = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
variables = (Trixi.density,
pressure))

amr_indicator = IndicatorLöhner(semi,
variable = density_pressure)

amr_controller = ControllerThreeLevel(semi, amr_indicator,
base_level = 4,
med_level = 0, med_threshold = 0.1, # med_level = current level
max_level = 6, max_threshold = 0.3)

amr_callback = AMRCallback(semi, amr_controller,
interval = 2,
adapt_initial_condition = true,
adapt_initial_condition_only_refine = true)
adapt_initial_condition_only_refine = true,
limiter! = positivity_limiter)

stepsize_callback = StepsizeCallback(cfl = 0.5)

Expand All @@ -100,12 +107,11 @@ callbacks = CallbackSet(summary_callback,
save_solution,
amr_callback, stepsize_callback)

stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
variables = (Trixi.density, pressure))

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition = false);
sol = solve(ode,
CarpenterKennedy2N54(stage_limiter! = positivity_limiter,
williamson_condition = false);
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
ode_default_options()..., callback = callbacks);
14 changes: 8 additions & 6 deletions examples/tree_2d_dgsem/elixir_euler_astro_jet_amr.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,11 @@ save_solution = SaveSolutionCallback(interval = 5000,
save_final_solution = true,
solution_variables = cons2prim)

# positivity limiter necessary for this tough example
positivity_limiter = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
variables = (Trixi.density,
pressure))

amr_indicator = IndicatorHennemannGassner(semi,
alpha_max = 1.0,
alpha_min = 0.0001,
Expand All @@ -102,18 +107,15 @@ amr_controller = ControllerThreeLevelCombined(semi, amr_indicator, indicator_sc,
amr_callback = AMRCallback(semi, amr_controller,
interval = 1,
adapt_initial_condition = true,
adapt_initial_condition_only_refine = true)
adapt_initial_condition_only_refine = true,
limiter! = positivity_limiter)

callbacks = CallbackSet(summary_callback,
analysis_callback, alive_callback,
amr_callback, save_solution)

# positivity limiter necessary for this tough example
stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
variables = (Trixi.density, pressure))

###############################################################################
# run the simulation
# use adaptive time stepping based on error estimates, time step roughly dt = 1e-7
sol = solve(ode, SSPRK43(stage_limiter!);
sol = solve(ode, SSPRK43(stage_limiter! = positivity_limiter);
ode_default_options()..., callback = callbacks);
14 changes: 8 additions & 6 deletions examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,11 @@ save_solution = SaveSolutionCallback(interval = 1000,
save_final_solution = true,
solution_variables = cons2prim)

# positivity limiter necessary for this tough example
positivity_limiter = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
variables = (Trixi.density,
pressure))

# Simulation also feasible without AMR: AMR reduces CPU time by a factor of about 2
amr_indicator = IndicatorHennemannGassner(semi,
alpha_max = 1.0,
Expand All @@ -105,18 +110,15 @@ amr_controller = ControllerThreeLevelCombined(semi, amr_indicator, indicator_sc,
amr_callback = AMRCallback(semi, amr_controller,
interval = 1,
adapt_initial_condition = true,
adapt_initial_condition_only_refine = true)
adapt_initial_condition_only_refine = true,
limiter! = positivity_limiter)

callbacks = CallbackSet(summary_callback,
analysis_callback, alive_callback,
amr_callback, save_solution)

# positivity limiter necessary for this tough example
stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
variables = (Trixi.density, pressure))

###############################################################################
# run the simulation
# use adaptive time stepping based on error estimates, time step roughly dt = 5e-3
sol = solve(ode, SSPRK43(stage_limiter!);
sol = solve(ode, SSPRK43(stage_limiter! = positivity_limiter);
ode_default_options()..., callback = callbacks);
16 changes: 11 additions & 5 deletions examples/tree_2d_dgsem/elixir_euler_positivity.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,16 +84,23 @@ save_solution = SaveSolutionCallback(interval = 100,
save_final_solution = true,
solution_variables = cons2prim)

positivity_limiter = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
variables = (Trixi.density,
pressure))

amr_indicator = IndicatorLöhner(semi,
variable = density_pressure)

amr_controller = ControllerThreeLevel(semi, amr_indicator,
base_level = 4,
med_level = 0, med_threshold = 0.1, # med_level = current level
max_level = 6, max_threshold = 0.3)

amr_callback = AMRCallback(semi, amr_controller,
interval = 2,
adapt_initial_condition = true,
adapt_initial_condition_only_refine = true)
adapt_initial_condition_only_refine = true,
limiter! = positivity_limiter)

stepsize_callback = StepsizeCallback(cfl = 0.8)

Expand All @@ -102,12 +109,11 @@ callbacks = CallbackSet(summary_callback,
save_solution,
amr_callback, stepsize_callback)

stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
variables = (Trixi.density, pressure))

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition = false);
sol = solve(ode,
CarpenterKennedy2N54(stage_limiter! = positivity_limiter,
williamson_condition = false);
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
ode_default_options()..., callback = callbacks);
16 changes: 11 additions & 5 deletions examples/tree_3d_dgsem/elixir_euler_blob_amr.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,16 +92,23 @@ save_solution = SaveSolutionCallback(interval = 200,
save_final_solution = true,
solution_variables = cons2prim)

positivity_limiter = PositivityPreservingLimiterZhangShu(thresholds = (1.0e-4, 1.0e-4),
variables = (Trixi.density,
pressure))

amr_indicator = IndicatorLöhner(semi,
variable = Trixi.density)

amr_controller = ControllerThreeLevel(semi, amr_indicator,
base_level = 1,
med_level = 0, med_threshold = 0.1, # med_level = current level
max_level = 6, max_threshold = 0.3)

amr_callback = AMRCallback(semi, amr_controller,
interval = 3,
adapt_initial_condition = false,
adapt_initial_condition_only_refine = true)
adapt_initial_condition_only_refine = true,
limiter! = positivity_limiter)

stepsize_callback = StepsizeCallback(cfl = 1.7)

Expand All @@ -110,12 +117,11 @@ callbacks = CallbackSet(summary_callback,
save_solution,
amr_callback, stepsize_callback)

stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (1.0e-4, 1.0e-4),
variables = (Trixi.density, pressure))

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition = false);
sol = solve(ode,
CarpenterKennedy2N54(stage_limiter! = positivity_limiter,
williamson_condition = false);
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
ode_default_options()..., callback = callbacks);
16 changes: 12 additions & 4 deletions src/callbacks_stage/positivity_zhang_shu.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,13 @@ function (limiter!::PositivityPreservingLimiterZhangShu)(u_ode, integrator,
return nothing
end

# Version used by the AMR callback
function (limiter!::PositivityPreservingLimiterZhangShu)(u, mesh, equations, solver,
cache, args...)
limiter_zhang_shu!(u, limiter!.thresholds, limiter!.variables, mesh, equations,
solver, cache, args...)
end

# Iterate over tuples in a type-stable way using "lispy tuple programming",
# similar to https://stackoverflow.com/a/55849398:
# Iterating over tuples of different functions isn't type-stable in general
Expand All @@ -49,21 +56,22 @@ end
# compile times can increase otherwise - but a handful of elements per tuple
# is definitely fine.
function limiter_zhang_shu!(u, thresholds::NTuple{N, <:Real}, variables::NTuple{N, Any},
mesh, equations, solver, cache) where {N}
mesh, equations, solver, cache, args...) where {N}
threshold = first(thresholds)
remaining_thresholds = Base.tail(thresholds)
variable = first(variables)
remaining_variables = Base.tail(variables)

limiter_zhang_shu!(u, threshold, variable, mesh, equations, solver, cache)
limiter_zhang_shu!(u, threshold, variable, mesh, equations, solver, cache,
args...)
limiter_zhang_shu!(u, remaining_thresholds, remaining_variables, mesh, equations,
solver, cache)
solver, cache, args...)
return nothing
end

# terminate the type-stable iteration over tuples
function limiter_zhang_shu!(u, thresholds::Tuple{}, variables::Tuple{},
mesh, equations, solver, cache)
mesh, equations, solver, cache, args...)
nothing
end

Expand Down
Loading
Loading