From 2e1f82e75de7c61566680b846fa08da3fa81ac44 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 3 Sep 2025 15:19:56 +0200 Subject: [PATCH 1/2] Apply limiter after refinement/coarsening with T8codeMesh --- .../elixir_euler_weak_blast_wave_amr.jl | 8 +++- .../elixir_euler_weak_blast_wave_amr.jl | 8 +++- src/callbacks_step/amr.jl | 46 +++++++++++++++++++ src/callbacks_step/amr_dg2d.jl | 28 ++++++++++- src/callbacks_step/amr_dg3d.jl | 28 ++++++++++- 5 files changed, 112 insertions(+), 6 deletions(-) diff --git a/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl index f8f934395a1..fca36cf9c29 100644 --- a/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl +++ b/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) diff --git a/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl index 03b6dac3fb4..221f197d98c 100644 --- a/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl +++ b/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) diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index 1af67f54cd5..c7e45875e94 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -897,6 +897,52 @@ 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) + + @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, diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index a3931c40933..a2b7ba6224c 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -528,6 +528,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 @@ -606,7 +622,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 diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index de8fc2cd86a..01ce1a27669 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -448,6 +448,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 @@ -532,7 +548,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 From d87fdb43067bcf8f5acfb931b6cdfe31c537fdf2 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 9 Sep 2025 10:20:11 +0200 Subject: [PATCH 2/2] Implement suggestions --- src/callbacks_step/amr.jl | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index c7e45875e94..f19763b6678 100644 --- a/src/callbacks_step/amr.jl +++ b/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))) @@ -901,8 +903,12 @@ end # 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) + 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))