Apply limiter after refinement/coarsening for T8codeMesh

examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl

@@ -36,7 +36,7 @@ initial_condition = initial_condition_weak_blast_wave
 # In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
 # Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
 # To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
-# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the 
+# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
 # `StepsizeCallback` (CFL-Condition) and less diffusion.
 surface_flux = FluxLaxFriedrichs(max_abs_speed_naive)
 volume_flux = flux_ranocha
@@ -97,7 +97,11 @@ 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! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6,
+                                                                                        5.0e-6),
+                                                                          variables = (Trixi.density,
+                                                                                       pressure)))
 
 stepsize_callback = StepsizeCallback(cfl = 0.5)
 

examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl

@@ -37,7 +37,7 @@ initial_condition = initial_condition_weak_blast_wave
 # In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
 # Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
 # To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
-# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the 
+# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
 # `StepsizeCallback` (CFL-Condition) and less diffusion.
 surface_flux = FluxLaxFriedrichs(max_abs_speed_naive)
 volume_flux = flux_ranocha
@@ -99,7 +99,11 @@ amr_controller = ControllerThreeLevel(semi, amr_indicator,
 amr_callback = AMRCallback(semi, amr_controller,
                            interval = 1,
                            adapt_initial_condition = false,
-                           adapt_initial_condition_only_refine = false)
+                           adapt_initial_condition_only_refine = false,
+                           limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6,
+                                                                                        5.0e-6),
+                                                                          variables = (Trixi.density,
+                                                                                       pressure)))
 
 stepsize_callback = StepsizeCallback(cfl = 0.5)
 

src/callbacks_step/amr.jl

@@ -869,7 +869,8 @@ end
 # Used when applying positivity limiter after refinement step
 # Saves every first id of the `2^ndims` child elements
 # All child elements can be addressed with `element_ids_new[i]:(element_ids_new[i] + 2^ndims - 1)`
-function compute_new_ids_refined_elements(elements_to_refine, mesh)
+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
@@ -883,7 +884,8 @@ 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)
+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)))
@@ -897,6 +899,56 @@ function compute_new_ids_coarsened_elements(elements_to_remove, mesh)
     return element_ids_new
 end
 
+# Auxiliary function to compute the new element ids for refined and coarsened elements
+# Used when applying positivity limiter after refinement and coarsening step for T8codeMesh
+function compute_new_ids_refined_coarsened_elements(difference, mesh::T8codeMesh,
+                                                    old_nelems, new_nelems)
+    T8_CHILDREN = 2^ndims(mesh) # number of children elements
+
+    # The `difference` vector has the length `old_nelems`.
+    # It contains a `1` at the specific index (= element id) if the element was refined,
+    # and `T8_CHILDREN` consecutive `-1`s if those elements were coarsened.
+    # Non-refined/coarsened elements have a `0` at their specific index.
+    @assert sum(difference .< 0) % T8_CHILDREN==0 "$(difference .< 0)"
+
+    refined_element_ids_new = Vector{Int}(undef, sum(difference .> 0))
+    coarsened_element_ids_new = Vector{Int}(undef,
+                                            div(sum(difference .< 0), T8_CHILDREN))
+    refined_index = 0
+    coarsened_index = 0
+
+    # Local element indices.
+    old_index = 1
+    new_index = 1
+    while old_index <= old_nelems && new_index <= new_nelems
+        if difference[old_index] > 0 # Refine.
+            refined_index += 1
+            refined_element_ids_new[refined_index] = new_index
+
+            # Increment `old_index` on the original mesh and the `new_index`
+            # on the refined mesh with the number of children, i.e., T8_CHILDREN = 4
+            old_index += 1
+            new_index += T8_CHILDREN
+        elseif difference[old_index] < 0 # Coarsen.
+            coarsened_index += 1
+            coarsened_element_ids_new[coarsened_index] = new_index
+
+            # Increment `old_index` on the original mesh with the number of children
+            # (T8_CHILDREN = 4 in 2D) and the `new_index` by one for the single
+            # coarsened element
+            old_index += T8_CHILDREN
+            new_index += 1
+        else # No changes.
+            # No refinement / coarsening occurred, so increment element index
+            # on each mesh by one
+            old_index += 1
+            new_index += 1
+        end
+    end
+
+    return refined_element_ids_new, coarsened_element_ids_new
+end
+
 """
     ControllerThreeLevel(semi, indicator; base_level=1,
                                           med_level=base_level, med_threshold=0.0,

src/callbacks_step/amr_dg2d.jl

@@ -527,6 +527,22 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations,
     GC.@preserve old_u_ode begin
         old_u = wrap_array(old_u_ode, mesh, equations, dg, cache)
 
+        # Compute mean values for the elements to be refined
+        # Only if limiter was passed
+        if limiter! !== nothing
+            elements_to_refine = findall(difference .> 0)
+            u_mean_refined_elements = Matrix{eltype(u_ode)}(undef,
+                                                            nvariables(equations),
+                                                            length(elements_to_refine))
+            for element in eachindex(elements_to_refine)
+                old_element_id = elements_to_refine[element]
+                # compute mean value
+                u_mean = compute_u_mean(old_u, old_element_id,
+                                        mesh, equations, dg, cache)
+                set_node_vars!(u_mean_refined_elements, u_mean, equations, dg, element)
+            end
+        end
+
         # Loop over all elements in old container and scale the old solution by the Jacobian
         # prior to interpolation or projection
         for old_element_id in 1:old_nelems
@@ -605,7 +621,17 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations,
 
     # Apply the positivity limiter to the solution
     if limiter! !== nothing
-        limiter!(u, mesh, equations, dg, cache)
+        # Precompute list with new element ids after refinement and coarsening
+        refined_element_ids_new, coarsened_element_ids_new = compute_new_ids_refined_coarsened_elements(difference,
+                                                                                                        mesh,
+                                                                                                        old_nelems,
+                                                                                                        new_nelems)
+        # Refined elements
+        limiter!(u, mesh, equations, dg, cache, refined_element_ids_new,
+                 u_mean_refined_elements)
+
+        # Coarsened elements
+        limiter!(u, mesh, equations, dg, cache, coarsened_element_ids_new)
     end
 
     return nothing

src/callbacks_step/amr_dg3d.jl

@@ -447,6 +447,22 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations,
     GC.@preserve old_u_ode begin
         old_u = wrap_array(old_u_ode, mesh, equations, dg, cache)
 
+        # Compute mean values for the elements to be refined
+        # Only if limiter was passed
+        if limiter! !== nothing
+            elements_to_refine = findall(difference .> 0)
+            u_mean_refined_elements = Matrix{eltype(u_ode)}(undef,
+                                                            nvariables(equations),
+                                                            length(elements_to_refine))
+            for element in eachindex(elements_to_refine)
+                old_element_id = elements_to_refine[element]
+                # compute mean value
+                u_mean = compute_u_mean(old_u, old_element_id,
+                                        mesh, equations, dg, cache)
+                set_node_vars!(u_mean_refined_elements, u_mean, equations, dg, element)
+            end
+        end
+
         # Loop over all elements in old container and scale the old solution by the Jacobian
         # prior to interpolation or projection
         for old_element_id in 1:old_nelems
@@ -531,7 +547,17 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations,
 
     # Apply the positivity limiter to the solution
     if limiter! !== nothing
-        limiter!(u, mesh, equations, dg, cache)
+        # Precompute list with new element ids after refinement and coarsening
+        refined_element_ids_new, coarsened_element_ids_new = compute_new_ids_refined_coarsened_elements(difference,
+                                                                                                        mesh,
+                                                                                                        old_nelems,
+                                                                                                        new_nelems)
+        # Refined elements
+        limiter!(u, mesh, equations, dg, cache,
+                 refined_element_ids_new, u_mean_refined_elements)
+
+        # Coarsened elements
+        limiter!(u, mesh, equations, dg, cache, coarsened_element_ids_new)
     end
 
     return nothing