Add positivity-limiter to the `AMRCallback` by patrickersing · Pull Request #2396 · trixi-framework/Trixi.jl

patrickersing · 2025-05-09T18:20:46Z

This PR introduces the possibility to pass a positivity-preserving limiter to the AMRCallback. The limiter is applied at the end of the coarsen!, refine! or adapt! step to ensure positivity after adaptation in the AMR callback. This generalizes the strategy from trixi-framework/TrixiShallowWater.jl#97 for arbitrary equations.

This tackles the first part of #2064. To close this issue we also need to implement positivity-limiting after the mortar projection.

github-actions · 2025-05-09T18:20:57Z

codecov · 2025-05-09T18:39:58Z

Codecov Report

❌ Patch coverage is 97.60000% with 6 lines in your changes missing coverage. Please review.
✅ Project coverage is 97.09%. Comparing base (e68e1ad) to head (f741aa3).

Files with missing lines	Patch %	Lines
src/callbacks_stage/positivity_zhang_shu_dg1d.jl	95.35%	2 Missing ⚠️
src/callbacks_stage/positivity_zhang_shu_dg2d.jl	95.35%	2 Missing ⚠️
src/callbacks_stage/positivity_zhang_shu_dg3d.jl	95.35%	2 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff            @@
##             main    #2396    +/-   ##
========================================
  Coverage   97.08%   97.09%            
========================================
  Files         612      612            
  Lines       47668    47906   +238     
========================================
+ Hits        46277    46510   +233     
- Misses       1391     1396     +5

Flag	Coverage Δ
unittests	`97.09% <97.60%> (+<0.01%)`	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

andrewwinters5000

Thanks for add this generic capability! It is much better than my hotfix from TrixiSW. My main question regards T8codeMesh

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>

DanielDoehring · 2025-05-11T18:36:39Z

Does it makes sense to call the limiter after refinement as well and not only after coarsening? Maybe make this also optional with a boolean variable? Normally, coarsening should be the troublemaker

andrewwinters5000 · 2025-05-12T04:56:41Z

Does it makes sense to call the limiter after refinement as well and not only after coarsening? Maybe make this also optional with a boolean variable? Normally, coarsening should be the troublemaker

There is no guarantee on that the sign of values after the refinement (interpolation) or coarsening (L2 projection in this case) when the solution is near vacuum.

patrickersing · 2025-05-12T07:56:37Z

Does it makes sense to call the limiter after refinement as well and not only after coarsening? Maybe make this also optional with a boolean variable? Normally, coarsening should be the troublemaker

For the elixir_euler_double_mach_amr.jl I checked that it is actually the refinement that causes positivity issues, if it runs without additional limiting.

patrickersing · 2025-05-12T08:31:45Z

Right now this only ensures positivity after the coarsening / refine step. I think in order to fix #2064 we would also need to add positivity-limiting after the mortar projection in prolong2mortars!. But for this we first need to find a good way to also access the positivity-limiter from the solver.
I would say we keep this as a draft for now, until we figure out a complete strategy how to implement this.

For the specialized SWE implementation in trixi-framework/TrixiShallowWater.jl#97 this additional limiting should not be necessary as we obtain positive values from the reconstruction strategy.

DanielDoehring · 2025-05-12T21:48:49Z

Right now this only ensures positivity after the coarsening / refine step. I think in order to fix #2064 we would also need to add positivity-limiting after the mortar projection in prolong2mortars!. But for this we first need to find a good way to also access the positivity-limiter from the solver. I would say we keep this as a draft for now, until we figure out a complete strategy how to implement this.

True - but this PR itself is already an improvement compared to the existing code. One option would be that we add the limiter! to the mortar and then supply such a mortar in the solver. But I agree that this would be something for a different PR.

I add this to the agenda of Tuesday's Trixi meeting.

patrickersing · 2025-05-13T07:35:48Z

True - but this PR itself is already an improvement compared to the existing code. One option would be that we add the limiter! to the mortar and then supply such a mortar in the solver. But I agree that this would be something for a different PR.

Passing the limiter to the mortar would also be my favored option. If we do this, the limiter would be available to both the AMRCallback and the solver. The main question is then how we want to pass the limiter in the elixirs. We could either pass the limiter both to the AMRCallback and DGSEM, or only pass it to DGSEM (and then extract it inside the AMRCallback from the solver).

Passing it only once reduces the risk that someone forgets to set the limiter at both places. On the other hand passing it explicitly to the AMRCallback makes it clear that we apply the limiter for refinement / coarsening.

patrickersing · 2025-05-15T15:34:45Z

An option could be to store the limiter inside the mortar similar to how it was drafted for the EC mortars.
Then we could have a special mortar type:

struct LobattoLegendrePositivityPreservingMortarL2{RealT <: Real, NNODES,
                               ForwardMatrix <: AbstractMatrix{RealT},
                               ReverseMatrix <: AbstractMatrix{RealT}, Limiter} <:
       AbstractMortarL2{RealT}
    forward_upper::ForwardMatrix
    forward_lower::ForwardMatrix
    reverse_upper::ReverseMatrix
    reverse_lower::ReverseMatrix
    limiter::Limiter
end

DanielDoehring · 2025-06-05T15:15:02Z

@patrickersing so what I would suggest is to add a limiter to a mortar first and then see if we can extract the limter from the mortar (via semi via solver via mortar) in some nice way that still allows for passing a potentially different limiter to the AMR callback/adaptor. I am happy to help here!

patrickersing · 2025-06-05T16:01:21Z

@patrickersing so what I would suggest is to add a limiter to a mortar first and then see if we can extract the limiter from the mortar (via semi via solver via mortar) in some nice way that still allows for passing a potentially different limiter to the AMR callback/àdaptor`. I am happy to help here!

That sounds like a good plan! If we manage to set the mortar limiter by default, but keep the flexibility to pass different limiters to the AMR callback that would be great.

Going forward I would then keep this PR dedicated to the AMR callback, and create a separate PR to implement the new mortar type.

bennibolm · 2025-07-28T09:09:53Z

I'm currently working on IDP limiting at mortars, and I'm interested in using the Shu limiter to ensure positivity when transferring the solution while refinement and coarsening. So, I only need this positivity limiter for AMR.
I tried some things in the current state of this PR and also implemented a version where all 4 (or 8 in 3D) refined child elements are shifted towards the mean value by the same theta. This (hopefully) preserves conservation when the current version fails.
I will probably set up a PR for this in the next days and discuss it when @patrickersing is back from vacation.
Things I noticed:

In many simulations, the refinement and coarsening produced negative value_mean (and therefore negative thetas) directly while initialization of the amr callback. I had to adapt the initial_refinement_level such that this doesn't happen.
I haven't found an example where the implementation with simultaneous shifting (=same theta for child elements) preserves conservation and the one without doesn't. I'll test some more setups.

DanielDoehring · 2025-07-28T13:08:42Z

In many simulations, the refinement and coarsening produced negative value_mean (and therefore negative thetas) directly while initialization of the amr callback.

Interesting, this sounds like AMR is happening " too late" in some sense.

I haven't found an example where the implementation without simultaneous shifting preserves conservation

What do you mean by "simultaneous shifting"?

bennibolm · 2025-07-28T13:23:05Z

In many simulations, the refinement and coarsening produced negative value_mean (and therefore negative thetas) directly while initialization of the amr callback.

Interesting, this sounds like AMR is happening " too late" in some sense.

With "initialization of the amr callback" I really mean the refinement at the very start of the simulation to get to max_level. So, there is no "too late" here.
One setup here is for instance initial_refinement_level=4 and for the amr controller base_level=4, max_level=7. So, there are at least three round of refinement necessary here.
Additionally, we still have the sharp discontinuous from the initial condition.
So, I think this really is expected, I guess.

I haven't found an example where the implementation without simultaneous shifting preserves conservation

What do you mean by "simultaneous shifting"?

I added a remark above to clarify this. Simultaneous shifting is when I use the same (minimum) theta for all just refined child elements:
The workflow then is:

go through all 4 (or 8) child elements and compute the theta as always
get the minimum theta of all child elements
use this minimum theta to shift all solutions of the child elements toward the mean value as always with: theta * u_node + (1-theta) * u_mean

DanielDoehring · 2025-07-28T13:30:25Z

In many simulations, the refinement and coarsening produced negative value_mean (and therefore negative thetas) directly while initialization of the amr callback.

With "initialization of the amr callback" I really mean the refinement at the very start of the simulation to get to max_level. So, there is no "too late" here. One setup here is for instance initial_refinement_level=4 and for the amr controller base_level=4, max_level=7. So, there are at least three round of refinement necessary here. Additionally, we still have the sharp discontinuous from the initial condition. So, I think this really is expected, I guess.

Ah I somehow missed that you were really talking about the init, my bad. Thanks for the explanations, though 👍

patrickersing · 2025-08-05T07:32:35Z

Thanks a lot for taking the initiative on this @bennibolm!
Hopefully this approach solves the present positivity issues and preserves conservation.

In many simulations, the refinement and coarsening produced negative value_mean (and therefore negative thetas) directly while initialization of the amr callback. I had to adapt the initial_refinement_level such that this doesn't happen.

Does this appear with or without the simultaneous shifting technique? I also think without any fixes this is to be expected due to the reinterpolation.

I haven't found an example where the implementation with simultaneous shifting (=same theta for child elements) preserves conservation and the one without doesn't. I'll test some more setups.

So both approaches lose conservation? Do you already have an idea in which step we lose the conservation? How big are the conservation errors that you observe?

bennibolm · 2025-08-21T07:42:48Z

My first implementation in patrickersing#47 is ready. Now, I first compute the (nonnegative) mean solution of elements that are going to be refined before the refinement process and then pass that information to the limiter. Then, all the 4 child elements are shifted w.r.t. this mean and with the same theta.
I'm already using this implementation in a branch of mine where I use mortars with subcell limiting. There it works very good.
Unfortunately, I'm struggling with finding a good example without subcell limiting. What were elixirs and setups where you found positivity issues with the solution transfer during the refinement/coarsening process? Feel free to comment here or directly in the PR.

bennibolm · 2025-08-21T07:47:50Z

In many simulations, the refinement and coarsening produced negative value_mean (and therefore negative thetas) directly while initialization of the amr callback. I had to adapt the initial_refinement_level such that this doesn't happen.

This is fixed by using the mean value from the original parent element. So before reinterpolation.

I haven't found an example where the implementation with simultaneous shifting (=same theta for child elements) preserves conservation and the one without doesn't. I'll test some more setups.

This was because I chose a too small initial refinement level at the start. The solution integrals change slightly during the initial refinement process. When I update those numbers after this initialization prozess I get conservation.

* Implement synchronized shifting for child elements * Fix increase of `element_id` * Revise implementation with mean value of parent element * Add 1D and 3D * Remove if and use dispatching * Decrease theta by eps() * Implement first suggestions * Compute mean within `limiter !== nothing` * Add suggested comment * Restructure call of limiter for refined elements * Implement suggestions * Implement thread parallel version; Add timer * Remove `@assert` * Add function for coarsened elements * Fix dimension * Add auxiliary functions to compute new elements ids * Fix function name

* Apply limiter after refinement/coarsening with T8codeMesh * Implement suggestions

DanielDoehring

What is the current status here? Are we good apart from coverage?

patrickersing · 2025-09-26T07:27:14Z

Yes, apart from the missing test coverage the feature should be complete.
What is missing to ensure positivity with AMR (and therefore fix #2064) is the limiting step after the mortar projection. However, this will be done in a separate PR.

DanielDoehring · 2025-09-26T07:28:34Z

Yes, apart from the missing test coverage the feature should be complete. What is missing to ensure positivity with AMR (and therefore fix #2396) is the limiting step after the mortar projection. However, this will be done in a separate PR.

Nice, yeah that is what I also recall 👍

bennibolm · 2025-11-06T12:59:53Z

Coming back to this PR, do you think it would be better to first merge the implementation in 2D? There, we would "only" need a simulation where the actual shifting after coarsening is used to reach full coverage.
In 1D and 3D, we need more. See above #2396 (comment)
@DanielDoehring @patrickersing

DanielDoehring · 2025-11-06T13:06:44Z

Coming back to this PR, do you think it would be better to first merge the implementation in 2D? There, we would "only" need a simulation where the actual shifting after coarsening is used to reach full coverage.

Yeah let's carve this out!

patrickersing · 2025-11-11T10:22:16Z

Coming back to this PR, do you think it would be better to first merge the implementation in 2D? There, we would "only" need a simulation where the actual shifting after coarsening is used to reach full coverage. In 1D and 3D, we need more. See above #2396 (comment) @DanielDoehring @patrickersing

Sounds good! What do we do with the uncovered lines of the 1D / 3D implementation?
Should we comment them out / remove them or accept the uncovered lines for now and open an Issue to keep track?

ranocha

As discussed in the Trixi.jl meeting today, please add a simple test for this implementation in all dimensions, e.g., just initializing a discontinuity within one element, refine it, and fix the arsing negative values.

bennibolm · 2025-11-17T10:23:12Z

As discussed in the Trixi.jl meeting today, please add a simple test for this implementation in all dimensions, e.g., just initializing a discontinuity within one element, refine it, and fix the arsing negative values.

I created simple test cases. But (at least to my understanding) I need to set up equations, solver, semidiscretization to be able to refine and coarsen properly. So, it's not that short anymore (50 lines per dimension or about 85 for all three since many lines can be used for 1D, 2D and 3D).
What would be the best place to put/save/call these test cases? Should I simply add it to the unit test file? Or alternatively, put everything into one test file and then call it? Something else?

DanielDoehring · 2025-11-17T10:28:44Z

I created simple test cases. But (at least to my understanding) I need to set up equations, solver, semidiscretization to be able to refine and coarsen properly. So, it's not that short anymore (50 lines per dimension or about 85 for all three since many lines can be used for 1D, 2D and 3D). What would be the best place to put/save/call these test cases? Should I simply add it to the unit test file? Or alternatively, put everything into one test file and then call it? Something else?

Thanks for taking care of this!
Yeah I think there is no way around setting up equation, mesh etc. 👍
I think test_unit is a decent place for this.

DanielDoehring · 2025-11-17T19:20:31Z

+# Set `limiter! = positivity_limiter` to apply the positivity-preserving limiter after 
+# coarsening and refinement steps.
 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)


I would probably keep only the positivity limiter in the "positivity" elixirs - the other elixirs were running before as well, right?

So, you would for instance remove it from all elixirs using HennemannGassner shock capturing? In the end, I don't really care. @patrickersing did you add the limiter to all the elixirs with a reason?

Hm yeah I would maybe keep this at a minimum for coverage, but not sure

I added the limiter to all elixirs that were already using the PositivityPreservingLimiterZhangShu with AMR. I think it makes sense to keep the AMR limiter for these examples, since the authors anticipated positvitiy issues.
The only exception I noticed are the t8code examples. I think I just added the limiter here for testing purpose. Here, we could remove the limiter and instead introduce a unit test.

All elixirs did run before, though for some we observed positivity violations.

Sounds reasonable to me. But as I said, I don't have a strong opinion on this.

I think, the positivity issues occurred mostly directly in the beginning when refining the initial condition.

DanielDoehring · 2025-11-17T19:39:11Z

+                                                           refined_original_cells,
+                                                           limiter!)


Somewhat unrelated to this PR, but this will probably also not really be tested.

Actually, I am not sure if this gets tested at all, as for this the Euler Gravity stuff with AMR would need to be called:

Trixi.jl/src/semidiscretization/semidiscretization_euler_gravity.jl

Lines 592 to 609 in ff6ba9b

@inline function (amr_callback::AMRCallback)(u_ode,

semi::SemidiscretizationEulerGravity,

t, iter; kwargs...)

passive_args = ((semi.cache.u_ode,

mesh_equations_solver_cache(semi.semi_gravity)...),)

has_changed = amr_callback(u_ode, mesh_equations_solver_cache(semi.semi_euler)...,

semi, t, iter;

kwargs..., passive_args = passive_args)

if has_changed

new_length = length(semi.cache.u_ode)

resize!(semi.cache.du_ode, new_length)

resize!(semi.cache.u_tmp1_ode, new_length)

resize!(semi.cache.u_tmp2_ode, new_length)

end

return has_changed

end

DanielDoehring · 2025-11-17T19:39:39Z

+function compute_new_ids_refined_elements(elements_to_refine,
+                                          mesh::Union{TreeMesh, P4estMesh})
+    element_ids_new = copy(elements_to_refine)
+    for i in eachindex(element_ids_new)
+        # Each refined element increases all ids of the following elements by 2^ndims - 1
+        for j in (i + 1):length(element_ids_new)
+            element_ids_new[j] += 2^ndims(mesh) - 1
+        end
+    end
+
+    return element_ids_new
+end
+
+# Auxiliary function to compute the new element ids for coarsened elements
+# Used when applying positivity limiter after coarsening step
+function compute_new_ids_coarsened_elements(elements_to_remove,
+                                            mesh::Union{TreeMesh, P4estMesh})
+    @assert length(elements_to_remove) % (2^ndims(mesh))==0 "The length of `elements_to_remove` must be a multiple of 2^ndims(mesh)."
+
+    element_ids_new = zeros(Int, div(length(elements_to_remove), 2^ndims(mesh)))
+    for i in eachindex(element_ids_new)
+        # New element id is the id of the first child
+        # Additionally, all following ids decrease by (2^ndims - 1) per coarsened element
+        element_ids_new[i] = elements_to_remove[2^ndims(mesh) * (i - 1) + 1] -
+                             (2^ndims(mesh) - 1) * (i - 1)
+    end
+
+    return element_ids_new
+end


I wonder if these could somehow be re-used for the other AMR stuff

add positivity limiter to the amr callback

7b9a40d

patrickersing and others added 2 commits May 9, 2025 21:22

remove duplicate limiter! argument

318bfad

Merge branch 'main' into amr_limiter

f98cd62

andrewwinters5000 reviewed May 9, 2025

View reviewed changes

Comment thread src/callbacks_stage/positivity_zhang_shu.jl

Comment thread src/callbacks_step/amr.jl Outdated

Comment thread src/callbacks_step/amr_dg2d.jl

Update src/callbacks_step/amr.jl

e221e40

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>

DanielDoehring added the enhancement New feature or request label May 12, 2025

Merge branch 'main' into amr_limiter

5faaf4d

bennibolm mentioned this pull request Jul 28, 2025

Implement synchronized shifting for child elements patrickersing/Trixi.jl#47

Merged

Merge branch 'main' into amr_limiter

d9b6c90

bennibolm and others added 3 commits August 29, 2025 08:37

Merge branch 'main' into amr_limiter

24d8df0

Merge branch 'main' into amr_limiter

2b7dc80

Apply limiter after refinement/coarsening for T8codeMesh (#50)

7b8d2ab

* Apply limiter after refinement/coarsening with T8codeMesh * Implement suggestions

DanielDoehring reviewed Sep 26, 2025

View reviewed changes

DanielDoehring marked this pull request as ready for review September 26, 2025 07:28

DanielDoehring added 2 commits September 26, 2025 15:22

doc

5c1a1a2

Merge branch 'main' into PatrickErsingLimitedMortars

6d3cbfb

bennibolm reviewed Nov 6, 2025

View reviewed changes

Comment thread src/callbacks_step/amr_dg2d.jl Outdated

ranocha requested changes Nov 13, 2025

View reviewed changes

bennibolm marked this pull request as draft November 17, 2025 08:37

Add simple tests

b338193

bennibolm requested changes Nov 17, 2025

View reviewed changes

Comment thread test/test_unit.jl Outdated

Update test/test_unit.jl

9864c9d

bennibolm requested a review from DanielDoehring November 17, 2025 16:58

DanielDoehring reviewed Nov 17, 2025

View reviewed changes

DanielDoehring requested changes Nov 17, 2025

View reviewed changes

Update src/callbacks_step/amr.jl

f99c917

DanielDoehring reviewed Nov 17, 2025

View reviewed changes

bennibolm and others added 6 commits November 18, 2025 10:59

Remove dispatching

8dbb521

Use Array{T, 2} instead of Matrix{T}

aebbea7

Add timer to computation of parent's mean value

e29c9da

Merge branch 'main' into amr_limiter

89e1aeb

Merge branch 'main' into amr_limiter

5fb37aa

Fix tests (Set periodicity and boundary_conditions)

f741aa3

	@inline function (amr_callback::AMRCallback)(u_ode,
	semi::SemidiscretizationEulerGravity,
	t, iter; kwargs...)
	passive_args = ((semi.cache.u_ode,
	mesh_equations_solver_cache(semi.semi_gravity)...),)
	has_changed = amr_callback(u_ode, mesh_equations_solver_cache(semi.semi_euler)...,
	semi, t, iter;
	kwargs..., passive_args = passive_args)

	if has_changed
	new_length = length(semi.cache.u_ode)
	resize!(semi.cache.du_ode, new_length)
	resize!(semi.cache.u_tmp1_ode, new_length)
	resize!(semi.cache.u_tmp2_ode, new_length)
	end

	return has_changed
	end

Conversation

patrickersing commented May 9, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

github-actions Bot commented May 9, 2025

Review checklist

Purpose and scope

Code quality

Documentation

Testing

Performance

Verification

Uh oh!

codecov Bot commented May 9, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

andrewwinters5000 left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

DanielDoehring commented May 11, 2025

Uh oh!

andrewwinters5000 commented May 12, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

patrickersing commented May 12, 2025

Uh oh!

patrickersing commented May 12, 2025

Uh oh!

DanielDoehring commented May 12, 2025

Uh oh!

patrickersing commented May 13, 2025

Uh oh!

patrickersing commented May 15, 2025

Uh oh!

DanielDoehring commented Jun 5, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

patrickersing commented Jun 5, 2025

Uh oh!

bennibolm commented Jul 28, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

DanielDoehring commented Jul 28, 2025

Uh oh!

bennibolm commented Jul 28, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

DanielDoehring commented Jul 28, 2025

Uh oh!

patrickersing commented Aug 5, 2025

Uh oh!

bennibolm commented Aug 21, 2025

Uh oh!

bennibolm commented Aug 21, 2025

Uh oh!

DanielDoehring left a comment

Choose a reason for hiding this comment

Uh oh!

patrickersing commented Sep 26, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

DanielDoehring commented Sep 26, 2025

Uh oh!

bennibolm commented Nov 6, 2025

Uh oh!

DanielDoehring commented Nov 6, 2025

Uh oh!

Uh oh!

patrickersing commented Nov 11, 2025

Uh oh!

ranocha left a comment

Choose a reason for hiding this comment

Uh oh!

bennibolm commented Nov 17, 2025

Uh oh!

DanielDoehring commented Nov 17, 2025

patrickersing commented May 9, 2025 •

edited

Loading

codecov Bot commented May 9, 2025 •

edited

Loading

andrewwinters5000 commented May 12, 2025 •

edited

Loading

DanielDoehring commented Jun 5, 2025 •

edited

Loading

bennibolm commented Jul 28, 2025 •

edited

Loading

bennibolm commented Jul 28, 2025 •

edited

Loading

patrickersing commented Sep 26, 2025 •

edited

Loading