From 4292a061d1fd4338585679cc6c79411f7fec74da Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 28 Jul 2025 11:50:11 +0200 Subject: [PATCH 01/17] Implement synchronized shifting for child elements --- src/callbacks_stage/positivity_zhang_shu.jl | 16 ++-- .../positivity_zhang_shu_dg2d.jl | 74 ++++++++++++++++++- src/callbacks_step/amr_dg2d.jl | 3 +- 3 files changed, 85 insertions(+), 8 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu.jl b/src/callbacks_stage/positivity_zhang_shu.jl index 79f9c2e9403..f65a3d2a291 100644 --- a/src/callbacks_stage/positivity_zhang_shu.jl +++ b/src/callbacks_stage/positivity_zhang_shu.jl @@ -40,9 +40,9 @@ end # Version used by the AMR callback function (limiter!::PositivityPreservingLimiterZhangShu)(u, mesh, equations, solver, - cache) + cache; refined_elements = []) limiter_zhang_shu!(u, limiter!.thresholds, limiter!.variables, mesh, equations, - solver, cache) + solver, cache; refined_elements = refined_elements) end # Iterate over tuples in a type-stable way using "lispy tuple programming", @@ -54,21 +54,25 @@ 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; + refined_elements = []) 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, + refined_elements = refined_elements) limiter_zhang_shu!(u, remaining_thresholds, remaining_variables, mesh, equations, - solver, cache) + solver, cache, + refined_elements = refined_elements) 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; + refined_elements = []) nothing end diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index 813dd65878b..01bb410e5c5 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -6,10 +6,82 @@ #! format: noindent function limiter_zhang_shu!(u, threshold::Real, variable, - mesh::AbstractMesh{2}, equations, dg::DGSEM, cache) + mesh::AbstractMesh{2}, equations, dg::DGSEM, cache; + refined_elements = []) @unpack weights = dg.basis @unpack inverse_jacobian = cache.elements + @trixi_timeit timer() "limiter! refined_elements" if !isempty(refined_elements) + element_id = 1 + theta_vec = ones(eltype(u), 2^ndims(mesh)) + u_mean_vec = zeros(eltype(u), size(u, 1), 2^ndims(mesh)) + + for element_id_old in 1:refined_elements[end] + if element_id_old in refined_elements + theta_vec .= one(eltype(u)) + u_mean_vec .= zero(eltype(u)) + # Iterate over the children of the current element + for new_element in 1:(2^ndims(mesh)) + new_element_id = element_id + new_element - 1 + + # determine minimum value + value_min = typemax(eltype(u)) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, new_element_id) + value_min = min(value_min, variable(u_node, equations)) + end + + # compute mean value + u_mean = zero(get_node_vars(u, equations, dg, 1, 1, new_element_id)) + total_volume = zero(eltype(u)) + for j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, + mesh, + i, j, + new_element_id))) + u_node = get_node_vars(u, equations, dg, i, j, new_element_id) + u_mean += u_node * weights[i] * weights[j] * volume_jacobian + total_volume += weights[i] * weights[j] * volume_jacobian + end + # normalize with the total volume + u_mean = u_mean / total_volume + u_mean_vec[:, new_element] .= u_mean + + # detect if limiting is necessary + value_min < threshold || continue + + # We compute the value directly with the mean values, as we assume that + # Jensen's inequality holds (e.g. pressure for compressible Euler equations). + value_mean = variable(u_mean, equations) + theta = (value_mean - threshold) / (value_mean - value_min) + theta_vec[new_element] = theta + end + theta_global = minimum(theta_vec) + theta_global < 1 || continue + + # Iterate again over the children to apply synchronized shifting + for new_element in 1:(2^ndims(mesh)) + new_element_id = element_id + new_element - 1 + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, new_element_id) + u_mean = get_node_vars(u_mean_vec, equations, dg, new_element) + set_node_vars!(u, + theta_global * u_node + + (1 - theta_global) * u_mean, + equations, dg, i, j, new_element_id) + end + end + + # Increment `element_id` on the refined mesh with the number + # of children, i.e., 4 in 2D + element_id += 2^ndims(mesh) + else + # Increment `element_id` on the unrefined mesh + element_id += 1 + end + end + end + @threaded for element in eachelement(dg, cache) # determine minimum value value_min = typemax(eltype(u)) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 188d8a0ef2c..48baeab401d 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -176,7 +176,8 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM # Apply the positivity limiter to the solution if limiter! !== nothing - limiter!(u, mesh, equations, dg, cache) + limiter!(u, mesh, equations, dg, cache; + refined_elements = elements_to_refine) end return nothing From abe6bc420a91279d380282384aa0b5cb7cea31d4 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 31 Jul 2025 16:43:28 +0200 Subject: [PATCH 02/17] Fix increase of `element_id` --- src/callbacks_stage/positivity_zhang_shu_dg2d.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index 01bb410e5c5..c4a4f3885de 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -18,11 +18,15 @@ function limiter_zhang_shu!(u, threshold::Real, variable, for element_id_old in 1:refined_elements[end] if element_id_old in refined_elements + # Increment `element_id` on the refined mesh with the number + # of children, i.e., 4 in 2D + element_id += 2^ndims(mesh) + theta_vec .= one(eltype(u)) u_mean_vec .= zero(eltype(u)) # Iterate over the children of the current element for new_element in 1:(2^ndims(mesh)) - new_element_id = element_id + new_element - 1 + new_element_id = element_id + new_element - 1 - 2^ndims(mesh) # determine minimum value value_min = typemax(eltype(u)) @@ -61,7 +65,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable, # Iterate again over the children to apply synchronized shifting for new_element in 1:(2^ndims(mesh)) - new_element_id = element_id + new_element - 1 + new_element_id = element_id + new_element - 1 - 2^ndims(mesh) for j in eachnode(dg), i in eachnode(dg) u_node = get_node_vars(u, equations, dg, i, j, new_element_id) u_mean = get_node_vars(u_mean_vec, equations, dg, new_element) @@ -71,10 +75,6 @@ function limiter_zhang_shu!(u, threshold::Real, variable, equations, dg, i, j, new_element_id) end end - - # Increment `element_id` on the refined mesh with the number - # of children, i.e., 4 in 2D - element_id += 2^ndims(mesh) else # Increment `element_id` on the unrefined mesh element_id += 1 From 31a4dbcbfa8ae1b49289496513ff0c77785e2051 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 6 Aug 2025 10:12:01 +0200 Subject: [PATCH 03/17] Revise implementation with mean value of parent element --- src/callbacks_stage/positivity_zhang_shu.jl | 19 ++- .../positivity_zhang_shu_dg2d.jl | 146 +++++++++--------- src/callbacks_step/amr_dg2d.jl | 25 ++- 3 files changed, 113 insertions(+), 77 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu.jl b/src/callbacks_stage/positivity_zhang_shu.jl index f65a3d2a291..0f56097383c 100644 --- a/src/callbacks_stage/positivity_zhang_shu.jl +++ b/src/callbacks_stage/positivity_zhang_shu.jl @@ -40,9 +40,12 @@ end # Version used by the AMR callback function (limiter!::PositivityPreservingLimiterZhangShu)(u, mesh, equations, solver, - cache; refined_elements = []) + cache; + refined_elements = nothing, + u_mean_refined_elements = nothing) limiter_zhang_shu!(u, limiter!.thresholds, limiter!.variables, mesh, equations, - solver, cache; refined_elements = refined_elements) + solver, cache; refined_elements = refined_elements, + u_mean_refined_elements = u_mean_refined_elements) end # Iterate over tuples in a type-stable way using "lispy tuple programming", @@ -55,24 +58,28 @@ end # is definitely fine. function limiter_zhang_shu!(u, thresholds::NTuple{N, <:Real}, variables::NTuple{N, Any}, mesh, equations, solver, cache; - refined_elements = []) where {N} + refined_elements = nothing, + u_mean_refined_elements = nothing) 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, - refined_elements = refined_elements) + refined_elements = refined_elements, + u_mean_refined_elements = u_mean_refined_elements) limiter_zhang_shu!(u, remaining_thresholds, remaining_variables, mesh, equations, solver, cache, - refined_elements = refined_elements) + refined_elements = refined_elements, + u_mean_refined_elements = u_mean_refined_elements) return nothing end # terminate the type-stable iteration over tuples function limiter_zhang_shu!(u, thresholds::Tuple{}, variables::Tuple{}, mesh, equations, solver, cache; - refined_elements = []) + refined_elements = nothing, + u_mean_refined_elements = nothing) nothing end diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index c4a4f3885de..645fcad6481 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -7,79 +7,22 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{2}, equations, dg::DGSEM, cache; - refined_elements = []) + refined_elements = nothing, + u_mean_refined_elements = nothing) @unpack weights = dg.basis @unpack inverse_jacobian = cache.elements - @trixi_timeit timer() "limiter! refined_elements" if !isempty(refined_elements) - element_id = 1 - theta_vec = ones(eltype(u), 2^ndims(mesh)) - u_mean_vec = zeros(eltype(u), size(u, 1), 2^ndims(mesh)) - - for element_id_old in 1:refined_elements[end] - if element_id_old in refined_elements - # Increment `element_id` on the refined mesh with the number - # of children, i.e., 4 in 2D - element_id += 2^ndims(mesh) - - theta_vec .= one(eltype(u)) - u_mean_vec .= zero(eltype(u)) - # Iterate over the children of the current element - for new_element in 1:(2^ndims(mesh)) - new_element_id = element_id + new_element - 1 - 2^ndims(mesh) - - # determine minimum value - value_min = typemax(eltype(u)) - for j in eachnode(dg), i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, j, new_element_id) - value_min = min(value_min, variable(u_node, equations)) - end - - # compute mean value - u_mean = zero(get_node_vars(u, equations, dg, 1, 1, new_element_id)) - total_volume = zero(eltype(u)) - for j in eachnode(dg), i in eachnode(dg) - volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, - mesh, - i, j, - new_element_id))) - u_node = get_node_vars(u, equations, dg, i, j, new_element_id) - u_mean += u_node * weights[i] * weights[j] * volume_jacobian - total_volume += weights[i] * weights[j] * volume_jacobian - end - # normalize with the total volume - u_mean = u_mean / total_volume - u_mean_vec[:, new_element] .= u_mean - - # detect if limiting is necessary - value_min < threshold || continue - - # We compute the value directly with the mean values, as we assume that - # Jensen's inequality holds (e.g. pressure for compressible Euler equations). - value_mean = variable(u_mean, equations) - theta = (value_mean - threshold) / (value_mean - value_min) - theta_vec[new_element] = theta - end - theta_global = minimum(theta_vec) - theta_global < 1 || continue - - # Iterate again over the children to apply synchronized shifting - for new_element in 1:(2^ndims(mesh)) - new_element_id = element_id + new_element - 1 - 2^ndims(mesh) - for j in eachnode(dg), i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, j, new_element_id) - u_mean = get_node_vars(u_mean_vec, equations, dg, new_element) - set_node_vars!(u, - theta_global * u_node + - (1 - theta_global) * u_mean, - equations, dg, i, j, new_element_id) - end - end - else - # Increment `element_id` on the unrefined mesh - element_id += 1 - end - end + if !isnothing(refined_elements) + @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." + + @trixi_timeit timer() "limiter_zhang_shu_refined_elements!" limiter_zhang_shu_refined_elements!(u, + threshold, + variable, + mesh, + equations, + dg, + refined_elements, + u_mean_refined_elements) end @threaded for element in eachelement(dg, cache) @@ -119,4 +62,67 @@ function limiter_zhang_shu!(u, threshold::Real, variable, return nothing end + +@inline function limiter_zhang_shu_refined_elements!(u, threshold::Real, variable, + mesh::AbstractMesh{2}, equations, + dg::DGSEM, + refined_elements, + u_mean_refined_elements) + @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." + + element_id_new = 1 + refined_element = 0 + for element_id_old in 1:refined_elements[end] + if element_id_old in refined_elements + refined_element += 1 + # Increment `element_id_new` on the refined mesh with the number + # of children, i.e., 4 in 2D + element_id_new += 2^ndims(mesh) + + u_mean = get_node_vars(u_mean_refined_elements, equations, dg, + refined_element) + + # We compute the value directly with the mean values, as we assume that + # Jensen's inequality holds (e.g. pressure for compressible Euler equations). + value_mean = variable(u_mean, equations) + + theta = one(eltype(u)) + + # Iterate over the children of the current element + for new_element in 1:(2^ndims(mesh)) + new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) + + # determine minimum value + value_min = typemax(eltype(u)) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, new_element_id) + value_min = min(value_min, variable(u_node, equations)) + end + + # detect if limiting is necessary + value_min < threshold || continue + + theta = min(theta, (value_mean - threshold) / (value_mean - value_min)) + end + theta < 1 || continue + + # Iterate again over the children to apply synchronized shifting + for new_element in 1:(2^ndims(mesh)) + new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) + + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, new_element_id) + set_node_vars!(u, + theta * u_node + (1 - theta) * u_mean, + equations, dg, i, j, new_element_id) + end + end + else + # Increment `element_id_new` on the unrefined mesh + element_id_new += 1 + end + end + + return nothing +end end # @muladd diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 48baeab401d..b549c63fa81 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -104,6 +104,28 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM GC.@preserve old_u_ode old_inverse_jacobian begin old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + # Compute mean values for the elements to be refined + (; weights) = dg.basis + 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 = zero(get_node_vars(old_u, equations, dg, 1, 1, old_element_id)) + total_volume = zero(eltype(old_u)) + for j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(get_inverse_jacobian(old_inverse_jacobian, + mesh, i, j, + old_element_id))) + u_node = get_node_vars(old_u, equations, dg, i, j, old_element_id) + u_mean += u_node * weights[i] * weights[j] * volume_jacobian + total_volume += weights[i] * weights[j] * volume_jacobian + end + # normalize with the total volume + u_mean = u_mean / total_volume + set_node_vars!(u_mean_refined_elements, u_mean, equations, dg, element) + end + if mesh isa P4estMesh # Loop over all elements in old container and scale the old solution by the Jacobian # prior to projection @@ -177,7 +199,8 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM # Apply the positivity limiter to the solution if limiter! !== nothing limiter!(u, mesh, equations, dg, cache; - refined_elements = elements_to_refine) + refined_elements = elements_to_refine, + u_mean_refined_elements = u_mean_refined_elements) end return nothing From 6be70f4221722a306125c58fabc2df7b7ea59e95 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 6 Aug 2025 10:12:41 +0200 Subject: [PATCH 04/17] Add 1D and 3D --- .../positivity_zhang_shu_dg1d.jl | 80 ++++++++++++++++++- .../positivity_zhang_shu_dg3d.jl | 80 ++++++++++++++++++- src/callbacks_step/amr_dg1d.jl | 17 ++++ src/callbacks_step/amr_dg3d.jl | 23 ++++++ 4 files changed, 198 insertions(+), 2 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl index 7797eb95b09..b45a60c6799 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl @@ -6,9 +6,24 @@ #! format: noindent function limiter_zhang_shu!(u, threshold::Real, variable, - mesh::AbstractMesh{1}, equations, dg::DGSEM, cache) + mesh::AbstractMesh{1}, equations, dg::DGSEM, cache; + refined_elements = nothing, + u_mean_refined_elements = nothing) @unpack weights = dg.basis + if !isnothing(refined_elements) + @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." + + @trixi_timeit timer() "limiter_zhang_shu_refined_elements!" limiter_zhang_shu_refined_elements!(u, + threshold, + variable, + mesh, + equations, + dg, + refined_elements, + u_mean_refined_elements) + end + @threaded for element in eachelement(dg, cache) # determine minimum value value_min = typemax(eltype(u)) @@ -42,4 +57,67 @@ function limiter_zhang_shu!(u, threshold::Real, variable, return nothing end + +@inline function limiter_zhang_shu_refined_elements!(u, threshold::Real, variable, + mesh::AbstractMesh{3}, equations, + dg::DGSEM, + refined_elements, + u_mean_refined_elements) + @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." + + element_id_new = 1 + refined_element = 0 + for element_id_old in 1:refined_elements[end] + if element_id_old in refined_elements + refined_element += 1 + # Increment `element_id_new` on the refined mesh with the number + # of children, i.e., 4 in 2D + element_id_new += 2^ndims(mesh) + + u_mean = get_node_vars(u_mean_refined_elements, equations, dg, + refined_element) + + # We compute the value directly with the mean values, as we assume that + # Jensen's inequality holds (e.g. pressure for compressible Euler equations). + value_mean = variable(u_mean, equations) + + theta = one(eltype(u)) + + # Iterate over the children of the current element + for new_element in 1:(2^ndims(mesh)) + new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) + + # determine minimum value + value_min = typemax(eltype(u)) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, new_element_id) + value_min = min(value_min, variable(u_node, equations)) + end + + # detect if limiting is necessary + value_min < threshold || continue + + theta = min(theta, (value_mean - threshold) / (value_mean - value_min)) + end + theta < 1 || continue + + # Iterate again over the children to apply synchronized shifting + for new_element in 1:(2^ndims(mesh)) + new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) + + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, new_element_id) + set_node_vars!(u, + theta * u_node + (1 - theta) * u_mean, + equations, dg, i, new_element_id) + end + end + else + # Increment `element_id_new` on the unrefined mesh + element_id_new += 1 + end + end + + return nothing +end end # @muladd diff --git a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl index 156abf35b4c..df43d8dd65a 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl @@ -6,10 +6,25 @@ #! format: noindent function limiter_zhang_shu!(u, threshold::Real, variable, - mesh::AbstractMesh{3}, equations, dg::DGSEM, cache) + mesh::AbstractMesh{3}, equations, dg::DGSEM, cache; + refined_elements = nothing, + u_mean_refined_elements = nothing) @unpack weights = dg.basis @unpack inverse_jacobian = cache.elements + if !isnothing(refined_elements) + @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." + + @trixi_timeit timer() "limiter_zhang_shu_refined_elements!" limiter_zhang_shu_refined_elements!(u, + threshold, + variable, + mesh, + equations, + dg, + refined_elements, + u_mean_refined_elements) + end + @threaded for element in eachelement(dg, cache) # determine minimum value value_min = typemax(eltype(u)) @@ -47,4 +62,67 @@ function limiter_zhang_shu!(u, threshold::Real, variable, return nothing end + +@inline function limiter_zhang_shu_refined_elements!(u, threshold::Real, variable, + mesh::AbstractMesh{3}, equations, + dg::DGSEM, + refined_elements, + u_mean_refined_elements) + @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." + + element_id_new = 1 + refined_element = 0 + for element_id_old in 1:refined_elements[end] + if element_id_old in refined_elements + refined_element += 1 + # Increment `element_id_new` on the refined mesh with the number + # of children, i.e., 4 in 2D + element_id_new += 2^ndims(mesh) + + u_mean = get_node_vars(u_mean_refined_elements, equations, dg, + refined_element) + + # We compute the value directly with the mean values, as we assume that + # Jensen's inequality holds (e.g. pressure for compressible Euler equations). + value_mean = variable(u_mean, equations) + + theta = one(eltype(u)) + + # Iterate over the children of the current element + for new_element in 1:(2^ndims(mesh)) + new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) + + # determine minimum value + value_min = typemax(eltype(u)) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, new_element_id) + value_min = min(value_min, variable(u_node, equations)) + end + + # detect if limiting is necessary + value_min < threshold || continue + + theta = min(theta, (value_mean - threshold) / (value_mean - value_min)) + end + theta < 1 || continue + + # Iterate again over the children to apply synchronized shifting + for new_element in 1:(2^ndims(mesh)) + new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, new_element_id) + set_node_vars!(u, + theta * u_node + (1 - theta) * u_mean, + equations, dg, i, j, k, new_element_id) + end + end + else + # Increment `element_id_new` on the unrefined mesh + element_id_new += 1 + end + end + + return nothing +end end # @muladd diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index 4f202523188..7eb420bd905 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -23,6 +23,23 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + # Compute mean values for the elements to be refined + (; weights) = dg.basis + 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 = zero(get_node_vars(old_u, equations, dg, 1, old_element_id)) + for i in eachnode(dg) + u_node = get_node_vars(old_u, equations, dg, i, element) + u_mean += u_node * weights[i] + end + # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2 + u_mean = u_mean / 2^ndims(mesh) + set_node_vars!(u_mean_refined_elements, u_mean, equations, dg, element) + end + # Get new list of leaf cells leaf_cell_ids = local_leaf_cells(mesh.tree) diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index f0eb5e65771..23d12b3e9f2 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -33,6 +33,29 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM GC.@preserve old_u_ode old_inverse_jacobian begin old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + # Compute mean values for the elements to be refined + (; weights) = dg.basis + 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 = zero(get_node_vars(old_u, equations, dg, 1, 1, 1, old_element_id)) + total_volume = zero(eltype(old_u)) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(get_inverse_jacobian(old_inverse_jacobian, + mesh, i, j, k, + old_element_id))) + u_node = get_node_vars(old_u, equations, dg, i, j, k, old_element_id) + u_mean += u_node * weights[i] * weights[j] * weights[k] * + volume_jacobian + total_volume += weights[i] * weights[j] * weights[k] * volume_jacobian + end + # normalize with the total volume + u_mean = u_mean / total_volume + set_node_vars!(u_mean_refined_elements, u_mean, equations, dg, element) + end + if mesh isa P4estMesh # Loop over all elements in old container and scale the old solution by the Jacobian # prior to projection From 8adec5d7d02891d3055b4f4ed9bb03a915d5ed0e Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 6 Aug 2025 10:46:07 +0200 Subject: [PATCH 05/17] Remove if and use dispatching --- src/callbacks_stage/positivity_zhang_shu.jl | 8 ++++++ .../positivity_zhang_shu_dg1d.jl | 25 ++++++++----------- .../positivity_zhang_shu_dg2d.jl | 23 ++++++++--------- .../positivity_zhang_shu_dg3d.jl | 23 ++++++++--------- 4 files changed, 39 insertions(+), 40 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu.jl b/src/callbacks_stage/positivity_zhang_shu.jl index 0f56097383c..dcb769b54bc 100644 --- a/src/callbacks_stage/positivity_zhang_shu.jl +++ b/src/callbacks_stage/positivity_zhang_shu.jl @@ -83,6 +83,14 @@ function limiter_zhang_shu!(u, thresholds::Tuple{}, variables::Tuple{}, nothing end +@inline function limiter_zhang_shu_refined_elements!(u, threshold::Real, variable, + mesh, equations, + dg::DGSEM, + refined_elements::Nothing, + u_mean_refined_elements) + nothing +end + include("positivity_zhang_shu_dg1d.jl") include("positivity_zhang_shu_dg2d.jl") include("positivity_zhang_shu_dg3d.jl") diff --git a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl index b45a60c6799..0a898df8b8d 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl @@ -11,18 +11,14 @@ function limiter_zhang_shu!(u, threshold::Real, variable, u_mean_refined_elements = nothing) @unpack weights = dg.basis - if !isnothing(refined_elements) - @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." - - @trixi_timeit timer() "limiter_zhang_shu_refined_elements!" limiter_zhang_shu_refined_elements!(u, - threshold, - variable, - mesh, - equations, - dg, - refined_elements, - u_mean_refined_elements) - end + @trixi_timeit timer() "limiter_zhang_shu_refined_elements!" limiter_zhang_shu_refined_elements!(u, + threshold, + variable, + mesh, + equations, + dg, + refined_elements, + u_mean_refined_elements) @threaded for element in eachelement(dg, cache) # determine minimum value @@ -59,10 +55,11 @@ function limiter_zhang_shu!(u, threshold::Real, variable, end @inline function limiter_zhang_shu_refined_elements!(u, threshold::Real, variable, - mesh::AbstractMesh{3}, equations, + mesh::AbstractMesh{1}, equations, dg::DGSEM, - refined_elements, + refined_elements::Vector{Int}, u_mean_refined_elements) + @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." element_id_new = 1 diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index 645fcad6481..d5bc54ffc46 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -12,18 +12,14 @@ function limiter_zhang_shu!(u, threshold::Real, variable, @unpack weights = dg.basis @unpack inverse_jacobian = cache.elements - if !isnothing(refined_elements) - @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." - - @trixi_timeit timer() "limiter_zhang_shu_refined_elements!" limiter_zhang_shu_refined_elements!(u, - threshold, - variable, - mesh, - equations, - dg, - refined_elements, - u_mean_refined_elements) - end + @trixi_timeit timer() "limiter_zhang_shu_refined_elements!" limiter_zhang_shu_refined_elements!(u, + threshold, + variable, + mesh, + equations, + dg, + refined_elements, + u_mean_refined_elements) @threaded for element in eachelement(dg, cache) # determine minimum value @@ -66,8 +62,9 @@ end @inline function limiter_zhang_shu_refined_elements!(u, threshold::Real, variable, mesh::AbstractMesh{2}, equations, dg::DGSEM, - refined_elements, + refined_elements::Vector{Int}, u_mean_refined_elements) + @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." element_id_new = 1 diff --git a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl index df43d8dd65a..ae4f5c72bc1 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl @@ -12,18 +12,14 @@ function limiter_zhang_shu!(u, threshold::Real, variable, @unpack weights = dg.basis @unpack inverse_jacobian = cache.elements - if !isnothing(refined_elements) - @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." - - @trixi_timeit timer() "limiter_zhang_shu_refined_elements!" limiter_zhang_shu_refined_elements!(u, - threshold, - variable, - mesh, - equations, - dg, - refined_elements, - u_mean_refined_elements) - end + @trixi_timeit timer() "limiter_zhang_shu_refined_elements!" limiter_zhang_shu_refined_elements!(u, + threshold, + variable, + mesh, + equations, + dg, + refined_elements, + u_mean_refined_elements) @threaded for element in eachelement(dg, cache) # determine minimum value @@ -66,8 +62,9 @@ end @inline function limiter_zhang_shu_refined_elements!(u, threshold::Real, variable, mesh::AbstractMesh{3}, equations, dg::DGSEM, - refined_elements, + refined_elements::Vector{Int}, u_mean_refined_elements) + @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." element_id_new = 1 From e50f54c26b7791a85cf1fd13872813da99a961cf Mon Sep 17 00:00:00 2001 From: bennibolm Date: Fri, 22 Aug 2025 10:25:39 +0200 Subject: [PATCH 06/17] Decrease theta by eps() --- src/callbacks_stage/positivity_zhang_shu_dg2d.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index d5bc54ffc46..46092456db1 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -103,6 +103,9 @@ end end theta < 1 || continue + # Make sure to really reach the threshold and not only by machine precision + theta -= eps(typeof(theta)) + # Iterate again over the children to apply synchronized shifting for new_element in 1:(2^ndims(mesh)) new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) From 27eebc2c4e9c9e30b2f9e0d0146d26f2fdbe0825 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 25 Aug 2025 13:42:33 +0200 Subject: [PATCH 07/17] Implement first suggestions --- src/callbacks_stage/positivity_zhang_shu_dg1d.jl | 10 +++++++--- src/callbacks_stage/positivity_zhang_shu_dg2d.jl | 7 ++++--- src/callbacks_stage/positivity_zhang_shu_dg3d.jl | 10 +++++++--- src/callbacks_step/amr_dg1d.jl | 4 +++- src/callbacks_step/amr_dg3d.jl | 4 +++- 5 files changed, 24 insertions(+), 11 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl index 0a898df8b8d..3b5d842f52e 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl @@ -62,8 +62,8 @@ end @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." - element_id_new = 1 - refined_element = 0 + element_id_new = 1 # Element id after refinement + refined_element = 0 # Index variable for the refined elements for element_id_old in 1:refined_elements[end] if element_id_old in refined_elements refined_element += 1 @@ -71,6 +71,7 @@ end # of children, i.e., 4 in 2D element_id_new += 2^ndims(mesh) + # Get the mean value from the parent element u_mean = get_node_vars(u_mean_refined_elements, equations, dg, refined_element) @@ -80,7 +81,7 @@ end theta = one(eltype(u)) - # Iterate over the children of the current element + # Iterate over the children of the current element to determine a joint limiting coefficient `theta` for new_element in 1:(2^ndims(mesh)) new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) @@ -98,6 +99,9 @@ end end theta < 1 || continue + # Make sure to really reach the threshold and not only by machine precision + theta -= eps(typeof(theta)) + # Iterate again over the children to apply synchronized shifting for new_element in 1:(2^ndims(mesh)) new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index 46092456db1..7607afa5ceb 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -67,8 +67,8 @@ end @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." - element_id_new = 1 - refined_element = 0 + element_id_new = 1 # Element id after refinement + refined_element = 0 # Index variable for the refined elements for element_id_old in 1:refined_elements[end] if element_id_old in refined_elements refined_element += 1 @@ -76,6 +76,7 @@ end # of children, i.e., 4 in 2D element_id_new += 2^ndims(mesh) + # Get the mean value from the parent element u_mean = get_node_vars(u_mean_refined_elements, equations, dg, refined_element) @@ -85,7 +86,7 @@ end theta = one(eltype(u)) - # Iterate over the children of the current element + # Iterate over the children of the current element to determine a joint limiting coefficient `theta` for new_element in 1:(2^ndims(mesh)) new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl index ae4f5c72bc1..11668e2b35f 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl @@ -67,8 +67,8 @@ end @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." - element_id_new = 1 - refined_element = 0 + element_id_new = 1 # Element id after refinement + refined_element = 0 # Index variable for the refined elements for element_id_old in 1:refined_elements[end] if element_id_old in refined_elements refined_element += 1 @@ -76,6 +76,7 @@ end # of children, i.e., 4 in 2D element_id_new += 2^ndims(mesh) + # Get the mean value from the parent element u_mean = get_node_vars(u_mean_refined_elements, equations, dg, refined_element) @@ -85,7 +86,7 @@ end theta = one(eltype(u)) - # Iterate over the children of the current element + # Iterate over the children of the current element to determine a joint limiting coefficient `theta` for new_element in 1:(2^ndims(mesh)) new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) @@ -103,6 +104,9 @@ end end theta < 1 || continue + # Make sure to really reach the threshold and not only by machine precision + theta -= eps(typeof(theta)) + # Iterate again over the children to apply synchronized shifting for new_element in 1:(2^ndims(mesh)) new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index 7eb420bd905..5d0f2b3fd54 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -92,7 +92,9 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, # Apply the positivity limiter to the solution if limiter! !== nothing - limiter!(u, mesh, equations, dg, cache) + limiter!(u, mesh, equations, dg, cache; + refined_elements = elements_to_refine, + u_mean_refined_elements = u_mean_refined_elements) end return nothing diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 23d12b3e9f2..277fec4d6f6 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -131,7 +131,9 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM # Apply the positivity limiter to the solution if limiter! !== nothing - limiter!(u, mesh, equations, dg, cache) + limiter!(u, mesh, equations, dg, cache; + refined_elements = elements_to_refine, + u_mean_refined_elements = u_mean_refined_elements) end return nothing From f927049a864df48afc401900a3030be10204ca76 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 25 Aug 2025 14:11:01 +0200 Subject: [PATCH 08/17] Compute mean within `limiter !== nothing` --- src/callbacks_step/amr_dg1d.jl | 30 +++++++++++++---------- src/callbacks_step/amr_dg2d.jl | 40 ++++++++++++++++-------------- src/callbacks_step/amr_dg3d.jl | 45 ++++++++++++++++++++-------------- 3 files changed, 65 insertions(+), 50 deletions(-) diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index 5d0f2b3fd54..54654f7260f 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -24,20 +24,24 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) # Compute mean values for the elements to be refined - (; weights) = dg.basis - 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 = zero(get_node_vars(old_u, equations, dg, 1, old_element_id)) - for i in eachnode(dg) - u_node = get_node_vars(old_u, equations, dg, i, element) - u_mean += u_node * weights[i] + # Only if limiter was passed + if limiter! !== nothing + (; weights) = dg.basis + 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 = zero(get_node_vars(old_u, equations, dg, 1, old_element_id)) + for i in eachnode(dg) + u_node = get_node_vars(old_u, equations, dg, i, element) + u_mean += u_node * weights[i] + end + # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2 + u_mean = u_mean / 2^ndims(mesh) + set_node_vars!(u_mean_refined_elements, u_mean, equations, dg, element) end - # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2 - u_mean = u_mean / 2^ndims(mesh) - set_node_vars!(u_mean_refined_elements, u_mean, equations, dg, element) end # Get new list of leaf cells diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index b549c63fa81..4430b99fada 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -105,25 +105,29 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) # Compute mean values for the elements to be refined - (; weights) = dg.basis - 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 = zero(get_node_vars(old_u, equations, dg, 1, 1, old_element_id)) - total_volume = zero(eltype(old_u)) - for j in eachnode(dg), i in eachnode(dg) - volume_jacobian = abs(inv(get_inverse_jacobian(old_inverse_jacobian, - mesh, i, j, - old_element_id))) - u_node = get_node_vars(old_u, equations, dg, i, j, old_element_id) - u_mean += u_node * weights[i] * weights[j] * volume_jacobian - total_volume += weights[i] * weights[j] * volume_jacobian + # Only if limiter was passed + if limiter! !== nothing + (; weights) = dg.basis + 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 = zero(get_node_vars(old_u, equations, dg, 1, 1, old_element_id)) + total_volume = zero(eltype(old_u)) + for j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(get_inverse_jacobian(old_inverse_jacobian, + mesh, i, j, + old_element_id))) + u_node = get_node_vars(old_u, equations, dg, i, j, old_element_id) + u_mean += u_node * weights[i] * weights[j] * volume_jacobian + total_volume += weights[i] * weights[j] * volume_jacobian + end + # normalize with the total volume + u_mean = u_mean / total_volume + set_node_vars!(u_mean_refined_elements, u_mean, equations, dg, element) end - # normalize with the total volume - u_mean = u_mean / total_volume - set_node_vars!(u_mean_refined_elements, u_mean, equations, dg, element) end if mesh isa P4estMesh diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 277fec4d6f6..7b70acc2020 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -34,26 +34,33 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) # Compute mean values for the elements to be refined - (; weights) = dg.basis - 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 = zero(get_node_vars(old_u, equations, dg, 1, 1, 1, old_element_id)) - total_volume = zero(eltype(old_u)) - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - volume_jacobian = abs(inv(get_inverse_jacobian(old_inverse_jacobian, - mesh, i, j, k, - old_element_id))) - u_node = get_node_vars(old_u, equations, dg, i, j, k, old_element_id) - u_mean += u_node * weights[i] * weights[j] * weights[k] * - volume_jacobian - total_volume += weights[i] * weights[j] * weights[k] * volume_jacobian + # Only if limiter was passed + if limiter! !== nothing + (; weights) = dg.basis + 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 = zero(get_node_vars(old_u, equations, dg, 1, 1, 1, + old_element_id)) + total_volume = zero(eltype(old_u)) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(get_inverse_jacobian(old_inverse_jacobian, + mesh, i, j, k, + old_element_id))) + u_node = get_node_vars(old_u, equations, dg, i, j, k, + old_element_id) + u_mean += u_node * weights[i] * weights[j] * weights[k] * + volume_jacobian + total_volume += weights[i] * weights[j] * weights[k] * + volume_jacobian + end + # normalize with the total volume + u_mean = u_mean / total_volume + set_node_vars!(u_mean_refined_elements, u_mean, equations, dg, element) end - # normalize with the total volume - u_mean = u_mean / total_volume - set_node_vars!(u_mean_refined_elements, u_mean, equations, dg, element) end if mesh isa P4estMesh From a30374f1f888d0a66aee2cf3677a36a3155d3ce4 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 25 Aug 2025 18:56:05 +0200 Subject: [PATCH 09/17] Add suggested comment --- src/callbacks_stage/positivity_zhang_shu_dg1d.jl | 9 +++++++++ src/callbacks_stage/positivity_zhang_shu_dg2d.jl | 9 +++++++++ src/callbacks_stage/positivity_zhang_shu_dg3d.jl | 9 +++++++++ 3 files changed, 27 insertions(+) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl index 3b5d842f52e..2179cc928e6 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl @@ -54,6 +54,15 @@ function limiter_zhang_shu!(u, threshold::Real, variable, return nothing end +# Modified version of the limiter used in the refinement step of the AMR callback. +# To ensure admissibility after the refinement step, we compute a joint +# limiting coefficient for all children elements and then limit against the +# admissible mean value of the parent element. +# This strategy is described in Remark 3 of the paper: +# - Arpit Babbar, Praveen Chandrashekar (2025) +# Lax-Wendroff flux reconstruction on adaptive curvilinear meshes with +# error based time stepping for hyperbolic conservation laws +# [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) @inline function limiter_zhang_shu_refined_elements!(u, threshold::Real, variable, mesh::AbstractMesh{1}, equations, dg::DGSEM, diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index 7607afa5ceb..f02a7afd596 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -59,6 +59,15 @@ function limiter_zhang_shu!(u, threshold::Real, variable, return nothing end +# Modified version of the limiter used in the refinement step of the AMR callback. +# To ensure admissibility after the refinement step, we compute a joint +# limiting coefficient for all children elements and then limit against the +# admissible mean value of the parent element. +# This strategy is described in Remark 3 of the paper: +# - Arpit Babbar, Praveen Chandrashekar (2025) +# Lax-Wendroff flux reconstruction on adaptive curvilinear meshes with +# error based time stepping for hyperbolic conservation laws +# [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) @inline function limiter_zhang_shu_refined_elements!(u, threshold::Real, variable, mesh::AbstractMesh{2}, equations, dg::DGSEM, diff --git a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl index 11668e2b35f..9e2c929cef4 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl @@ -59,6 +59,15 @@ function limiter_zhang_shu!(u, threshold::Real, variable, return nothing end +# Modified version of the limiter used in the refinement step of the AMR callback. +# To ensure admissibility after the refinement step, we compute a joint +# limiting coefficient for all children elements and then limit against the +# admissible mean value of the parent element. +# This strategy is described in Remark 3 of the paper: +# - Arpit Babbar, Praveen Chandrashekar (2025) +# Lax-Wendroff flux reconstruction on adaptive curvilinear meshes with +# error based time stepping for hyperbolic conservation laws +# [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) @inline function limiter_zhang_shu_refined_elements!(u, threshold::Real, variable, mesh::AbstractMesh{3}, equations, dg::DGSEM, From e8551a6f13d911787e70274c48af55d701bd29a9 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 25 Aug 2025 19:44:23 +0200 Subject: [PATCH 10/17] Restructure call of limiter for refined elements --- src/callbacks_stage/positivity_zhang_shu.jl | 30 ++++--------------- .../positivity_zhang_shu_dg1d.jl | 21 +++---------- .../positivity_zhang_shu_dg2d.jl | 21 +++---------- .../positivity_zhang_shu_dg3d.jl | 21 +++---------- src/callbacks_step/amr_dg1d.jl | 5 ++-- src/callbacks_step/amr_dg2d.jl | 5 ++-- src/callbacks_step/amr_dg3d.jl | 5 ++-- 7 files changed, 24 insertions(+), 84 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu.jl b/src/callbacks_stage/positivity_zhang_shu.jl index dcb769b54bc..7f9b4949cbc 100644 --- a/src/callbacks_stage/positivity_zhang_shu.jl +++ b/src/callbacks_stage/positivity_zhang_shu.jl @@ -40,12 +40,9 @@ end # Version used by the AMR callback function (limiter!::PositivityPreservingLimiterZhangShu)(u, mesh, equations, solver, - cache; - refined_elements = nothing, - u_mean_refined_elements = nothing) + cache, kwargs...) limiter_zhang_shu!(u, limiter!.thresholds, limiter!.variables, mesh, equations, - solver, cache; refined_elements = refined_elements, - u_mean_refined_elements = u_mean_refined_elements) + solver, cache, kwargs...) end # Iterate over tuples in a type-stable way using "lispy tuple programming", @@ -57,37 +54,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; - refined_elements = nothing, - u_mean_refined_elements = nothing) where {N} + mesh, equations, solver, cache, kwargs...) 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, - refined_elements = refined_elements, - u_mean_refined_elements = u_mean_refined_elements) + kwargs...) limiter_zhang_shu!(u, remaining_thresholds, remaining_variables, mesh, equations, - solver, cache, - refined_elements = refined_elements, - u_mean_refined_elements = u_mean_refined_elements) + solver, cache, kwargs...) return nothing end # terminate the type-stable iteration over tuples function limiter_zhang_shu!(u, thresholds::Tuple{}, variables::Tuple{}, - mesh, equations, solver, cache; - refined_elements = nothing, - u_mean_refined_elements = nothing) - nothing -end - -@inline function limiter_zhang_shu_refined_elements!(u, threshold::Real, variable, - mesh, equations, - dg::DGSEM, - refined_elements::Nothing, - u_mean_refined_elements) + mesh, equations, solver, cache, kwargs...) nothing end diff --git a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl index 2179cc928e6..c744f48a4a0 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl @@ -6,20 +6,9 @@ #! format: noindent function limiter_zhang_shu!(u, threshold::Real, variable, - mesh::AbstractMesh{1}, equations, dg::DGSEM, cache; - refined_elements = nothing, - u_mean_refined_elements = nothing) + mesh::AbstractMesh{1}, equations, dg::DGSEM, cache) @unpack weights = dg.basis - @trixi_timeit timer() "limiter_zhang_shu_refined_elements!" limiter_zhang_shu_refined_elements!(u, - threshold, - variable, - mesh, - equations, - dg, - refined_elements, - u_mean_refined_elements) - @threaded for element in eachelement(dg, cache) # determine minimum value value_min = typemax(eltype(u)) @@ -63,11 +52,9 @@ end # Lax-Wendroff flux reconstruction on adaptive curvilinear meshes with # error based time stepping for hyperbolic conservation laws # [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) -@inline function limiter_zhang_shu_refined_elements!(u, threshold::Real, variable, - mesh::AbstractMesh{1}, equations, - dg::DGSEM, - refined_elements::Vector{Int}, - u_mean_refined_elements) +function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{1}, + equations, dg::DGSEM, cache, + refined_elements::Vector{Int}, u_mean_refined_elements) @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index f02a7afd596..46329a31fad 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -6,21 +6,10 @@ #! format: noindent function limiter_zhang_shu!(u, threshold::Real, variable, - mesh::AbstractMesh{2}, equations, dg::DGSEM, cache; - refined_elements = nothing, - u_mean_refined_elements = nothing) + mesh::AbstractMesh{2}, equations, dg::DGSEM, cache) @unpack weights = dg.basis @unpack inverse_jacobian = cache.elements - @trixi_timeit timer() "limiter_zhang_shu_refined_elements!" limiter_zhang_shu_refined_elements!(u, - threshold, - variable, - mesh, - equations, - dg, - refined_elements, - u_mean_refined_elements) - @threaded for element in eachelement(dg, cache) # determine minimum value value_min = typemax(eltype(u)) @@ -68,11 +57,9 @@ end # Lax-Wendroff flux reconstruction on adaptive curvilinear meshes with # error based time stepping for hyperbolic conservation laws # [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) -@inline function limiter_zhang_shu_refined_elements!(u, threshold::Real, variable, - mesh::AbstractMesh{2}, equations, - dg::DGSEM, - refined_elements::Vector{Int}, - u_mean_refined_elements) +function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{2}, + equations, dg::DGSEM, cache, + refined_elements::Vector{Int}, u_mean_refined_elements) @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." diff --git a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl index 9e2c929cef4..15b3b2c907e 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl @@ -6,21 +6,10 @@ #! format: noindent function limiter_zhang_shu!(u, threshold::Real, variable, - mesh::AbstractMesh{3}, equations, dg::DGSEM, cache; - refined_elements = nothing, - u_mean_refined_elements = nothing) + mesh::AbstractMesh{3}, equations, dg::DGSEM, cache) @unpack weights = dg.basis @unpack inverse_jacobian = cache.elements - @trixi_timeit timer() "limiter_zhang_shu_refined_elements!" limiter_zhang_shu_refined_elements!(u, - threshold, - variable, - mesh, - equations, - dg, - refined_elements, - u_mean_refined_elements) - @threaded for element in eachelement(dg, cache) # determine minimum value value_min = typemax(eltype(u)) @@ -68,11 +57,9 @@ end # Lax-Wendroff flux reconstruction on adaptive curvilinear meshes with # error based time stepping for hyperbolic conservation laws # [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) -@inline function limiter_zhang_shu_refined_elements!(u, threshold::Real, variable, - mesh::AbstractMesh{3}, equations, - dg::DGSEM, - refined_elements::Vector{Int}, - u_mean_refined_elements) +function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{3}, + equations, dg::DGSEM, cache, + refined_elements::Vector{Int}, u_mean_refined_elements) @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index 54654f7260f..0a067753b1a 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -96,9 +96,8 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, # Apply the positivity limiter to the solution if limiter! !== nothing - limiter!(u, mesh, equations, dg, cache; - refined_elements = elements_to_refine, - u_mean_refined_elements = u_mean_refined_elements) + limiter!(u, mesh, equations, dg, cache, + elements_to_refine, u_mean_refined_elements) end return nothing diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 4430b99fada..112acc666e2 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -202,9 +202,8 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM # Apply the positivity limiter to the solution if limiter! !== nothing - limiter!(u, mesh, equations, dg, cache; - refined_elements = elements_to_refine, - u_mean_refined_elements = u_mean_refined_elements) + limiter!(u, mesh, equations, dg, cache, + elements_to_refine, u_mean_refined_elements) end return nothing diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 7b70acc2020..cd9565e08d6 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -138,9 +138,8 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM # Apply the positivity limiter to the solution if limiter! !== nothing - limiter!(u, mesh, equations, dg, cache; - refined_elements = elements_to_refine, - u_mean_refined_elements = u_mean_refined_elements) + limiter!(u, mesh, equations, dg, cache, + elements_to_refine, u_mean_refined_elements) end return nothing From a6bc4d2034a1978b7aacd0f888f23530ab7ee08b Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 28 Aug 2025 09:43:06 +0200 Subject: [PATCH 11/17] Implement suggestions --- src/callbacks_stage/positivity_zhang_shu.jl | 12 ++++++------ src/callbacks_stage/positivity_zhang_shu_dg1d.jl | 2 +- src/callbacks_stage/positivity_zhang_shu_dg2d.jl | 2 +- src/callbacks_stage/positivity_zhang_shu_dg3d.jl | 2 +- 4 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu.jl b/src/callbacks_stage/positivity_zhang_shu.jl index 7f9b4949cbc..e909c245b0b 100644 --- a/src/callbacks_stage/positivity_zhang_shu.jl +++ b/src/callbacks_stage/positivity_zhang_shu.jl @@ -40,9 +40,9 @@ end # Version used by the AMR callback function (limiter!::PositivityPreservingLimiterZhangShu)(u, mesh, equations, solver, - cache, kwargs...) + cache, args...) limiter_zhang_shu!(u, limiter!.thresholds, limiter!.variables, mesh, equations, - solver, cache, kwargs...) + solver, cache, args...) end # Iterate over tuples in a type-stable way using "lispy tuple programming", @@ -54,22 +54,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, kwargs...) 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, - kwargs...) + args...) limiter_zhang_shu!(u, remaining_thresholds, remaining_variables, mesh, equations, - solver, cache, kwargs...) + 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, kwargs...) + mesh, equations, solver, cache, args...) nothing end diff --git a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl index c744f48a4a0..f9727fc6c8f 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl @@ -51,7 +51,7 @@ end # - Arpit Babbar, Praveen Chandrashekar (2025) # Lax-Wendroff flux reconstruction on adaptive curvilinear meshes with # error based time stepping for hyperbolic conservation laws -# [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) +# [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{1}, equations, dg::DGSEM, cache, refined_elements::Vector{Int}, u_mean_refined_elements) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index 46329a31fad..47f1f965913 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -56,7 +56,7 @@ end # - Arpit Babbar, Praveen Chandrashekar (2025) # Lax-Wendroff flux reconstruction on adaptive curvilinear meshes with # error based time stepping for hyperbolic conservation laws -# [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) +# [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{2}, equations, dg::DGSEM, cache, refined_elements::Vector{Int}, u_mean_refined_elements) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl index 15b3b2c907e..5da37420af7 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl @@ -56,7 +56,7 @@ end # - Arpit Babbar, Praveen Chandrashekar (2025) # Lax-Wendroff flux reconstruction on adaptive curvilinear meshes with # error based time stepping for hyperbolic conservation laws -# [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) +# [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{3}, equations, dg::DGSEM, cache, refined_elements::Vector{Int}, u_mean_refined_elements) From 3f9e72c4ae493ee4929d461b3665b94d5be3603d Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 28 Aug 2025 12:05:48 +0200 Subject: [PATCH 12/17] Implement thread parallel version; Add timer --- .../positivity_zhang_shu_dg1d.jl | 92 +++++++++---------- .../positivity_zhang_shu_dg2d.jl | 92 +++++++++---------- .../positivity_zhang_shu_dg3d.jl | 92 +++++++++---------- src/callbacks_step/amr_dg2d.jl | 7 +- 4 files changed, 133 insertions(+), 150 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl index f9727fc6c8f..b3225959128 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl @@ -58,60 +58,54 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{1}, @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." - element_id_new = 1 # Element id after refinement - refined_element = 0 # Index variable for the refined elements - for element_id_old in 1:refined_elements[end] - if element_id_old in refined_elements - refined_element += 1 - # Increment `element_id_new` on the refined mesh with the number - # of children, i.e., 4 in 2D - element_id_new += 2^ndims(mesh) - - # Get the mean value from the parent element - u_mean = get_node_vars(u_mean_refined_elements, equations, dg, - refined_element) - - # We compute the value directly with the mean values, as we assume that - # Jensen's inequality holds (e.g. pressure for compressible Euler equations). - value_mean = variable(u_mean, equations) - - theta = one(eltype(u)) - - # Iterate over the children of the current element to determine a joint limiting coefficient `theta` - for new_element in 1:(2^ndims(mesh)) - new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) - - # determine minimum value - value_min = typemax(eltype(u)) - for i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, new_element_id) - value_min = min(value_min, variable(u_node, equations)) - end - - # detect if limiting is necessary - value_min < threshold || continue - - theta = min(theta, (value_mean - threshold) / (value_mean - value_min)) + # Precompute list with new element ids after refinement + # Only save the last element id of the child elements and address all with (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + element_ids_new = copy(refined_elements) + for i in eachindex(element_ids_new) + # Each refined element increases the ids of all following elements by 2^ndims(mesh) - 1 + for j in i:length(element_ids_new) + element_ids_new[j] += 2^ndims(mesh) - 1 + end + end + + @threaded for i in eachindex(element_ids_new) + # Get the mean value from the parent element + u_mean = get_node_vars(u_mean_refined_elements, equations, dg, i) + + # We compute the value directly with the mean values, as we assume that + # Jensen's inequality holds (e.g. pressure for compressible Euler equations). + value_mean = variable(u_mean, equations) + + theta = one(eltype(u)) + + # Iterate over the children of the current element to determine a joint limiting coefficient `theta` + for new_element_id in (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + # determine minimum value + value_min = typemax(eltype(u)) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, new_element_id) + value_min = min(value_min, variable(u_node, equations)) end - theta < 1 || continue - # Make sure to really reach the threshold and not only by machine precision - theta -= eps(typeof(theta)) + # detect if limiting is necessary + value_min < threshold || continue + + theta = min(theta, (value_mean - threshold) / (value_mean - value_min)) + end + + theta < 1 || continue - # Iterate again over the children to apply synchronized shifting - for new_element in 1:(2^ndims(mesh)) - new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) + # Make sure to really reach the threshold and not only by machine precision + theta -= eps(typeof(theta)) - for i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, new_element_id) - set_node_vars!(u, - theta * u_node + (1 - theta) * u_mean, - equations, dg, i, new_element_id) - end + # Iterate again over the children to apply joint shifting + for new_element_id in (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, new_element_id) + set_node_vars!(u, + theta * u_node + (1 - theta) * u_mean, + equations, dg, i, new_element_id) end - else - # Increment `element_id_new` on the unrefined mesh - element_id_new += 1 end end diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index 47f1f965913..44b8705342d 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -63,60 +63,54 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{2}, @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." - element_id_new = 1 # Element id after refinement - refined_element = 0 # Index variable for the refined elements - for element_id_old in 1:refined_elements[end] - if element_id_old in refined_elements - refined_element += 1 - # Increment `element_id_new` on the refined mesh with the number - # of children, i.e., 4 in 2D - element_id_new += 2^ndims(mesh) - - # Get the mean value from the parent element - u_mean = get_node_vars(u_mean_refined_elements, equations, dg, - refined_element) - - # We compute the value directly with the mean values, as we assume that - # Jensen's inequality holds (e.g. pressure for compressible Euler equations). - value_mean = variable(u_mean, equations) - - theta = one(eltype(u)) - - # Iterate over the children of the current element to determine a joint limiting coefficient `theta` - for new_element in 1:(2^ndims(mesh)) - new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) - - # determine minimum value - value_min = typemax(eltype(u)) - for j in eachnode(dg), i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, j, new_element_id) - value_min = min(value_min, variable(u_node, equations)) - end - - # detect if limiting is necessary - value_min < threshold || continue - - theta = min(theta, (value_mean - threshold) / (value_mean - value_min)) + # Precompute list with new element ids after refinement + # Only save the last element id of the child elements and address all with (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + element_ids_new = copy(refined_elements) + for i in eachindex(element_ids_new) + # Each refined element increases the ids of all following elements by 2^ndims(mesh) - 1 + for j in i:length(element_ids_new) + element_ids_new[j] += 2^ndims(mesh) - 1 + end + end + + @threaded for i in eachindex(element_ids_new) + # Get the mean value from the parent element + u_mean = get_node_vars(u_mean_refined_elements, equations, dg, i) + + # We compute the value directly with the mean values, as we assume that + # Jensen's inequality holds (e.g. pressure for compressible Euler equations). + value_mean = variable(u_mean, equations) + + theta = one(eltype(u)) + + # Iterate over the children of the current element to determine a joint limiting coefficient `theta` + for new_element_id in (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + # determine minimum value + value_min = typemax(eltype(u)) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, new_element_id) + value_min = min(value_min, variable(u_node, equations)) end - theta < 1 || continue - # Make sure to really reach the threshold and not only by machine precision - theta -= eps(typeof(theta)) + # detect if limiting is necessary + value_min < threshold || continue + + theta = min(theta, (value_mean - threshold) / (value_mean - value_min)) + end + + theta < 1 || continue - # Iterate again over the children to apply synchronized shifting - for new_element in 1:(2^ndims(mesh)) - new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) + # Make sure to really reach the threshold and not only by machine precision + theta -= eps(typeof(theta)) - for j in eachnode(dg), i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, j, new_element_id) - set_node_vars!(u, - theta * u_node + (1 - theta) * u_mean, - equations, dg, i, j, new_element_id) - end + # Iterate again over the children to apply joint shifting + for new_element_id in (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, new_element_id) + set_node_vars!(u, + theta * u_node + (1 - theta) * u_mean, + equations, dg, i, j, new_element_id) end - else - # Increment `element_id_new` on the unrefined mesh - element_id_new += 1 end end diff --git a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl index 5da37420af7..5f85f1ce3c3 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl @@ -63,60 +63,54 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{3}, @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." - element_id_new = 1 # Element id after refinement - refined_element = 0 # Index variable for the refined elements - for element_id_old in 1:refined_elements[end] - if element_id_old in refined_elements - refined_element += 1 - # Increment `element_id_new` on the refined mesh with the number - # of children, i.e., 4 in 2D - element_id_new += 2^ndims(mesh) - - # Get the mean value from the parent element - u_mean = get_node_vars(u_mean_refined_elements, equations, dg, - refined_element) - - # We compute the value directly with the mean values, as we assume that - # Jensen's inequality holds (e.g. pressure for compressible Euler equations). - value_mean = variable(u_mean, equations) - - theta = one(eltype(u)) - - # Iterate over the children of the current element to determine a joint limiting coefficient `theta` - for new_element in 1:(2^ndims(mesh)) - new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) - - # determine minimum value - value_min = typemax(eltype(u)) - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, j, k, new_element_id) - value_min = min(value_min, variable(u_node, equations)) - end - - # detect if limiting is necessary - value_min < threshold || continue - - theta = min(theta, (value_mean - threshold) / (value_mean - value_min)) + # Precompute list with new element ids after refinement + # Only save the last element id of the child elements and address all with (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + element_ids_new = copy(refined_elements) + for i in eachindex(element_ids_new) + # Each refined element increases the ids of all following elements by 2^ndims(mesh) - 1 + for j in i:length(element_ids_new) + element_ids_new[j] += 2^ndims(mesh) - 1 + end + end + + @threaded for i in eachindex(element_ids_new) + # Get the mean value from the parent element + u_mean = get_node_vars(u_mean_refined_elements, equations, dg, i) + + # We compute the value directly with the mean values, as we assume that + # Jensen's inequality holds (e.g. pressure for compressible Euler equations). + value_mean = variable(u_mean, equations) + + theta = one(eltype(u)) + + # Iterate over the children of the current element to determine a joint limiting coefficient `theta` + for new_element_id in (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + # determine minimum value + value_min = typemax(eltype(u)) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, new_element_id) + value_min = min(value_min, variable(u_node, equations)) end - theta < 1 || continue - # Make sure to really reach the threshold and not only by machine precision - theta -= eps(typeof(theta)) + # detect if limiting is necessary + value_min < threshold || continue + + theta = min(theta, (value_mean - threshold) / (value_mean - value_min)) + end + + theta < 1 || continue - # Iterate again over the children to apply synchronized shifting - for new_element in 1:(2^ndims(mesh)) - new_element_id = element_id_new + new_element - 1 - 2^ndims(mesh) + # Make sure to really reach the threshold and not only by machine precision + theta -= eps(typeof(theta)) - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, j, k, new_element_id) - set_node_vars!(u, - theta * u_node + (1 - theta) * u_mean, - equations, dg, i, j, k, new_element_id) - end + # Iterate again over the children to apply joint shifting + for new_element_id in (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, new_element_id) + set_node_vars!(u, + theta * u_node + (1 - theta) * u_mean, + equations, dg, i, j, k, new_element_id) end - else - # Increment `element_id_new` on the unrefined mesh - element_id_new += 1 end end diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 112acc666e2..8e62d1423fd 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -202,8 +202,9 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM # Apply the positivity limiter to the solution if limiter! !== nothing - limiter!(u, mesh, equations, dg, cache, - elements_to_refine, u_mean_refined_elements) + @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, + elements_to_refine, + u_mean_refined_elements) end return nothing @@ -409,7 +410,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, # Apply the positivity limiter to the solution if limiter! !== nothing - limiter!(u, mesh, equations, dg, cache) + @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache) end return nothing From 2132f68966b7ef3451e3124a327f42ff76647cda Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 28 Aug 2025 13:47:06 +0200 Subject: [PATCH 13/17] Remove `@assert` --- src/callbacks_stage/positivity_zhang_shu_dg1d.jl | 1 - src/callbacks_stage/positivity_zhang_shu_dg2d.jl | 1 - src/callbacks_stage/positivity_zhang_shu_dg3d.jl | 1 - 3 files changed, 3 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl index b3225959128..e17ead7cb19 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl @@ -56,7 +56,6 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{1}, equations, dg::DGSEM, cache, refined_elements::Vector{Int}, u_mean_refined_elements) @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." - @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." # Precompute list with new element ids after refinement # Only save the last element id of the child elements and address all with (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index 44b8705342d..bbaeb6c227d 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -61,7 +61,6 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{2}, equations, dg::DGSEM, cache, refined_elements::Vector{Int}, u_mean_refined_elements) @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." - @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." # Precompute list with new element ids after refinement # Only save the last element id of the child elements and address all with (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] diff --git a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl index 5f85f1ce3c3..2d40a0c6d66 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl @@ -61,7 +61,6 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{3}, equations, dg::DGSEM, cache, refined_elements::Vector{Int}, u_mean_refined_elements) @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." - @assert maximum(refined_elements)==refined_elements[end] "The maximum element id in `refined_elements` must be equal to the last element id in the mesh." # Precompute list with new element ids after refinement # Only save the last element id of the child elements and address all with (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] From 8f3d782dbecaae2da269852a4c2192e37fd7fb6f Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 28 Aug 2025 13:50:53 +0200 Subject: [PATCH 14/17] Add function for coarsened elements --- .../positivity_zhang_shu_dg1d.jl | 55 ++++++++++++++++++ .../positivity_zhang_shu_dg2d.jl | 58 +++++++++++++++++++ .../positivity_zhang_shu_dg3d.jl | 58 +++++++++++++++++++ src/callbacks_step/amr_dg1d.jl | 8 ++- src/callbacks_step/amr_dg2d.jl | 3 +- src/callbacks_step/amr_dg3d.jl | 8 ++- 6 files changed, 183 insertions(+), 7 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl index e17ead7cb19..002ec9c4580 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl @@ -110,4 +110,59 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{1}, return nothing end + +# Modified version of the limiter used in the coarsening step of the AMR callback. +# To ensure admissibility after the coarsening step, we apply the limiter to +# the coarsened elements. +function limiter_zhang_shu!(u, threshold::Real, variable, + mesh::AbstractMesh{1}, equations, dg::DGSEM, cache, + removed_elements::Vector{Int}) + @assert length(removed_elements) % (2^ndims(mesh))==0 "The length of `removed_elements` must be a multiple of 2^ndims(mesh)." + + @unpack weights = dg.basis + @unpack inverse_jacobian = cache.elements + + # Precompute list with new element ids after refinement + element_ids_new = zeros(Int, div(length(removed_elements), 2^ndims(mesh))) + for i in eachindex(element_ids_new) + # The new element id is the id of the first child minus (2^ndims(mesh) - 1) times the number of already coarsened elements. + element_ids_new[i] = removed_elements[4 * (i - 1) + 1] - + (2^ndims(mesh) - 1) * (i - 1) + end + + # Apply limiter to coarsened elements + @threaded for element in element_ids_new + # determine minimum value + value_min = typemax(eltype(u)) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + value_min = min(value_min, variable(u_node, equations)) + end + + # detect if limiting is necessary + value_min < threshold || continue + + # compute mean value + u_mean = zero(get_node_vars(u, equations, dg, 1, element)) + total_volume = zero(eltype(u)) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + u_mean += u_node * weights[i] + end + # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2 + u_mean = u_mean / 2^ndims(mesh) + + # We compute the value directly with the mean values, as we assume that + # Jensen's inequality holds (e.g. pressure for compressible Euler equations). + value_mean = variable(u_mean, equations) + theta = (value_mean - threshold) / (value_mean - value_min) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + set_node_vars!(u, theta * u_node + (1 - theta) * u_mean, + equations, dg, i, element) + end + end + + return nothing +end end # @muladd diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index bbaeb6c227d..b55b068a2d1 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -115,4 +115,62 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{2}, return nothing end + +# Modified version of the limiter used in the coarsening step of the AMR callback. +# To ensure admissibility after the coarsening step, we apply the limiter to +# the coarsened elements. +function limiter_zhang_shu!(u, threshold::Real, variable, + mesh::AbstractMesh{2}, equations, dg::DGSEM, cache, + removed_elements::Vector{Int}) + @assert length(removed_elements) % (2^ndims(mesh))==0 "The length of `removed_elements` must be a multiple of 2^ndims(mesh)." + + @unpack weights = dg.basis + @unpack inverse_jacobian = cache.elements + + # Precompute list with new element ids after refinement + element_ids_new = zeros(Int, div(length(removed_elements), 2^ndims(mesh))) + for i in eachindex(element_ids_new) + # The new element id is the id of the first child minus (2^ndims(mesh) - 1) times the number of already coarsened elements. + element_ids_new[i] = removed_elements[4 * (i - 1) + 1] - + (2^ndims(mesh) - 1) * (i - 1) + end + + # Apply limiter to coarsened elements + @threaded for element in element_ids_new + # determine minimum value + value_min = typemax(eltype(u)) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + value_min = min(value_min, variable(u_node, equations)) + end + + # detect if limiting is necessary + value_min < threshold || continue + + # compute mean value + u_mean = zero(get_node_vars(u, equations, dg, 1, 1, element)) + total_volume = zero(eltype(u)) + for j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh, + i, j, element))) + u_node = get_node_vars(u, equations, dg, i, j, element) + u_mean += u_node * weights[i] * weights[j] * volume_jacobian + total_volume += weights[i] * weights[j] * volume_jacobian + end + # normalize with the total volume + u_mean = u_mean / total_volume + + # We compute the value directly with the mean values, as we assume that + # Jensen's inequality holds (e.g. pressure for compressible Euler equations). + value_mean = variable(u_mean, equations) + theta = (value_mean - threshold) / (value_mean - value_min) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + set_node_vars!(u, theta * u_node + (1 - theta) * u_mean, + equations, dg, i, j, element) + end + end + + return nothing +end end # @muladd diff --git a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl index 2d40a0c6d66..984afc1de29 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl @@ -115,4 +115,62 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{3}, return nothing end + +# Modified version of the limiter used in the coarsening step of the AMR callback. +# To ensure admissibility after the coarsening step, we apply the limiter to +# the coarsened elements. +function limiter_zhang_shu!(u, threshold::Real, variable, + mesh::AbstractMesh{3}, equations, dg::DGSEM, cache, + removed_elements::Vector{Int}) + @assert length(removed_elements) % (2^ndims(mesh))==0 "The length of `removed_elements` must be a multiple of 2^ndims(mesh)." + + @unpack weights = dg.basis + @unpack inverse_jacobian = cache.elements + + # Precompute list with new element ids after refinement + element_ids_new = zeros(Int, div(length(removed_elements), 2^ndims(mesh))) + for i in eachindex(element_ids_new) + # The new element id is the id of the first child minus (2^ndims(mesh) - 1) times the number of already coarsened elements. + element_ids_new[i] = removed_elements[4 * (i - 1) + 1] - + (2^ndims(mesh) - 1) * (i - 1) + end + + # Apply limiter to coarsened elements + @threaded for element in element_ids_new + # determine minimum value + value_min = typemax(eltype(u)) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + value_min = min(value_min, variable(u_node, equations)) + end + + # detect if limiting is necessary + value_min < threshold || continue + + # compute mean value + u_mean = zero(get_node_vars(u, equations, dg, 1, 1, 1, element)) + total_volume = zero(eltype(u)) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh, + i, j, k, element))) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + u_mean += u_node * weights[i] * weights[j] * weights[k] * volume_jacobian + total_volume += weights[i] * weights[j] * weights[k] * volume_jacobian + end + # normalize with the total volume + u_mean = u_mean / total_volume + + # We compute the value directly with the mean values, as we assume that + # Jensen's inequality holds (e.g. pressure for compressible Euler equations). + value_mean = variable(u_mean, equations) + theta = (value_mean - threshold) / (value_mean - value_min) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + set_node_vars!(u, theta * u_node + (1 - theta) * u_mean, + equations, dg, i, j, k, element) + end + end + + return nothing +end end # @muladd diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index 0a067753b1a..6ea547e51de 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -96,8 +96,9 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, # Apply the positivity limiter to the solution if limiter! !== nothing - limiter!(u, mesh, equations, dg, cache, - elements_to_refine, u_mean_refined_elements) + @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, + elements_to_refine, + u_mean_refined_elements) end return nothing @@ -248,7 +249,8 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, # Apply the positivity limiter to the solution if limiter! !== nothing - limiter!(u, mesh, equations, dg, cache) + @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, + elements_to_remove) end return nothing diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 8e62d1423fd..91064d1625f 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -410,7 +410,8 @@ function coarsen!(u_ode::AbstractVector, adaptor, # Apply the positivity limiter to the solution if limiter! !== nothing - @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache) + @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, + elements_to_remove) end return nothing diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index cd9565e08d6..cb8a1e49d35 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -138,8 +138,9 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM # Apply the positivity limiter to the solution if limiter! !== nothing - limiter!(u, mesh, equations, dg, cache, - elements_to_refine, u_mean_refined_elements) + @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, + elements_to_refine, + u_mean_refined_elements) end return nothing @@ -332,7 +333,8 @@ function coarsen!(u_ode::AbstractVector, adaptor, # Apply the positivity limiter to the solution if limiter! !== nothing - limiter!(u, mesh, equations, dg, cache) + @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, + elements_to_remove) end return nothing From aca2bb76cbc3fe3cb91742aa0d0981d4c2d68854 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 28 Aug 2025 13:58:58 +0200 Subject: [PATCH 15/17] Fix dimension --- src/callbacks_stage/positivity_zhang_shu_dg1d.jl | 2 +- src/callbacks_stage/positivity_zhang_shu_dg2d.jl | 2 +- src/callbacks_stage/positivity_zhang_shu_dg3d.jl | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl index 002ec9c4580..a168afd753c 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl @@ -126,7 +126,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable, element_ids_new = zeros(Int, div(length(removed_elements), 2^ndims(mesh))) for i in eachindex(element_ids_new) # The new element id is the id of the first child minus (2^ndims(mesh) - 1) times the number of already coarsened elements. - element_ids_new[i] = removed_elements[4 * (i - 1) + 1] - + element_ids_new[i] = removed_elements[2^ndims(mesh) * (i - 1) + 1] - (2^ndims(mesh) - 1) * (i - 1) end diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index b55b068a2d1..011f18b336d 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -131,7 +131,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable, element_ids_new = zeros(Int, div(length(removed_elements), 2^ndims(mesh))) for i in eachindex(element_ids_new) # The new element id is the id of the first child minus (2^ndims(mesh) - 1) times the number of already coarsened elements. - element_ids_new[i] = removed_elements[4 * (i - 1) + 1] - + element_ids_new[i] = removed_elements[2^ndims(mesh) * (i - 1) + 1] - (2^ndims(mesh) - 1) * (i - 1) end diff --git a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl index 984afc1de29..20f71a9c2ae 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl @@ -131,7 +131,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable, element_ids_new = zeros(Int, div(length(removed_elements), 2^ndims(mesh))) for i in eachindex(element_ids_new) # The new element id is the id of the first child minus (2^ndims(mesh) - 1) times the number of already coarsened elements. - element_ids_new[i] = removed_elements[4 * (i - 1) + 1] - + element_ids_new[i] = removed_elements[2^ndims(mesh) * (i - 1) + 1] - (2^ndims(mesh) - 1) * (i - 1) end From cfb168a690a1a45d571edecfe93ac1ecf6d8d65d Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 28 Aug 2025 15:03:13 +0200 Subject: [PATCH 16/17] Add auxiliary functions to compute new elements ids --- .../positivity_zhang_shu_dg1d.jl | 30 +++-------------- .../positivity_zhang_shu_dg2d.jl | 30 +++-------------- .../positivity_zhang_shu_dg3d.jl | 30 +++-------------- src/callbacks_step/amr.jl | 32 +++++++++++++++++++ src/callbacks_step/amr_dg1d.jl | 10 ++++-- src/callbacks_step/amr_dg2d.jl | 10 ++++-- src/callbacks_step/amr_dg3d.jl | 10 ++++-- 7 files changed, 71 insertions(+), 81 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl index a168afd753c..21bdf1acded 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl @@ -54,18 +54,8 @@ end # [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{1}, equations, dg::DGSEM, cache, - refined_elements::Vector{Int}, u_mean_refined_elements) - @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." - - # Precompute list with new element ids after refinement - # Only save the last element id of the child elements and address all with (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] - element_ids_new = copy(refined_elements) - for i in eachindex(element_ids_new) - # Each refined element increases the ids of all following elements by 2^ndims(mesh) - 1 - for j in i:length(element_ids_new) - element_ids_new[j] += 2^ndims(mesh) - 1 - end - end + element_ids_new::Vector{Int}, u_mean_refined_elements) + @assert length(element_ids_new)==size(u_mean_refined_elements, 2) "The length of `element_ids_new` must match the second dimension of `u_mean_refined_elements`." @threaded for i in eachindex(element_ids_new) # Get the mean value from the parent element @@ -78,7 +68,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{1}, theta = one(eltype(u)) # Iterate over the children of the current element to determine a joint limiting coefficient `theta` - for new_element_id in (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + for new_element_id in element_ids_new[i]:(element_ids_new[i] + 2^ndims(mesh) - 1) # determine minimum value value_min = typemax(eltype(u)) for i in eachnode(dg) @@ -98,7 +88,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{1}, theta -= eps(typeof(theta)) # Iterate again over the children to apply joint shifting - for new_element_id in (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + for new_element_id in element_ids_new[i]:(element_ids_new[i] + 2^ndims(mesh) - 1) for i in eachnode(dg) u_node = get_node_vars(u, equations, dg, i, new_element_id) set_node_vars!(u, @@ -116,20 +106,10 @@ end # the coarsened elements. function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{1}, equations, dg::DGSEM, cache, - removed_elements::Vector{Int}) - @assert length(removed_elements) % (2^ndims(mesh))==0 "The length of `removed_elements` must be a multiple of 2^ndims(mesh)." - + element_ids_new::Vector{Int}) @unpack weights = dg.basis @unpack inverse_jacobian = cache.elements - # Precompute list with new element ids after refinement - element_ids_new = zeros(Int, div(length(removed_elements), 2^ndims(mesh))) - for i in eachindex(element_ids_new) - # The new element id is the id of the first child minus (2^ndims(mesh) - 1) times the number of already coarsened elements. - element_ids_new[i] = removed_elements[2^ndims(mesh) * (i - 1) + 1] - - (2^ndims(mesh) - 1) * (i - 1) - end - # Apply limiter to coarsened elements @threaded for element in element_ids_new # determine minimum value diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index 011f18b336d..dc0730afad1 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -59,18 +59,8 @@ end # [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{2}, equations, dg::DGSEM, cache, - refined_elements::Vector{Int}, u_mean_refined_elements) - @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." - - # Precompute list with new element ids after refinement - # Only save the last element id of the child elements and address all with (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] - element_ids_new = copy(refined_elements) - for i in eachindex(element_ids_new) - # Each refined element increases the ids of all following elements by 2^ndims(mesh) - 1 - for j in i:length(element_ids_new) - element_ids_new[j] += 2^ndims(mesh) - 1 - end - end + element_ids_new::Vector{Int}, u_mean_refined_elements) + @assert length(element_ids_new)==size(u_mean_refined_elements, 2) "The length of `element_ids_new` must match the second dimension of `u_mean_refined_elements`." @threaded for i in eachindex(element_ids_new) # Get the mean value from the parent element @@ -83,7 +73,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{2}, theta = one(eltype(u)) # Iterate over the children of the current element to determine a joint limiting coefficient `theta` - for new_element_id in (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + for new_element_id in element_ids_new[i]:(element_ids_new[i] + 2^ndims(mesh) - 1) # determine minimum value value_min = typemax(eltype(u)) for j in eachnode(dg), i in eachnode(dg) @@ -103,7 +93,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{2}, theta -= eps(typeof(theta)) # Iterate again over the children to apply joint shifting - for new_element_id in (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + for new_element_id in element_ids_new[i]:(element_ids_new[i] + 2^ndims(mesh) - 1) for j in eachnode(dg), i in eachnode(dg) u_node = get_node_vars(u, equations, dg, i, j, new_element_id) set_node_vars!(u, @@ -121,20 +111,10 @@ end # the coarsened elements. function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{2}, equations, dg::DGSEM, cache, - removed_elements::Vector{Int}) - @assert length(removed_elements) % (2^ndims(mesh))==0 "The length of `removed_elements` must be a multiple of 2^ndims(mesh)." - + element_ids_new::Vector{Int}) @unpack weights = dg.basis @unpack inverse_jacobian = cache.elements - # Precompute list with new element ids after refinement - element_ids_new = zeros(Int, div(length(removed_elements), 2^ndims(mesh))) - for i in eachindex(element_ids_new) - # The new element id is the id of the first child minus (2^ndims(mesh) - 1) times the number of already coarsened elements. - element_ids_new[i] = removed_elements[2^ndims(mesh) * (i - 1) + 1] - - (2^ndims(mesh) - 1) * (i - 1) - end - # Apply limiter to coarsened elements @threaded for element in element_ids_new # determine minimum value diff --git a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl index 20f71a9c2ae..f8ce7c5d173 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl @@ -59,18 +59,8 @@ end # [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622) function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{3}, equations, dg::DGSEM, cache, - refined_elements::Vector{Int}, u_mean_refined_elements) - @assert length(refined_elements)==size(u_mean_refined_elements, 2) "The length of `refined_elements` must match the second dimension of `u_mean_refined_elements`." - - # Precompute list with new element ids after refinement - # Only save the last element id of the child elements and address all with (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] - element_ids_new = copy(refined_elements) - for i in eachindex(element_ids_new) - # Each refined element increases the ids of all following elements by 2^ndims(mesh) - 1 - for j in i:length(element_ids_new) - element_ids_new[j] += 2^ndims(mesh) - 1 - end - end + element_ids_new::Vector{Int}, u_mean_refined_elements) + @assert length(element_ids_new)==size(u_mean_refined_elements, 2) "The length of `element_ids_new` must match the second dimension of `u_mean_refined_elements`." @threaded for i in eachindex(element_ids_new) # Get the mean value from the parent element @@ -83,7 +73,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{3}, theta = one(eltype(u)) # Iterate over the children of the current element to determine a joint limiting coefficient `theta` - for new_element_id in (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + for new_element_id in element_ids_new[i]:(element_ids_new[i] + 2^ndims(mesh) - 1) # determine minimum value value_min = typemax(eltype(u)) for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) @@ -103,7 +93,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{3}, theta -= eps(typeof(theta)) # Iterate again over the children to apply joint shifting - for new_element_id in (element_ids_new[i] - 2^ndims(mesh) + 1):element_ids_new[i] + for new_element_id in element_ids_new[i]:(element_ids_new[i] + 2^ndims(mesh) - 1) for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) u_node = get_node_vars(u, equations, dg, i, j, k, new_element_id) set_node_vars!(u, @@ -121,20 +111,10 @@ end # the coarsened elements. function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{3}, equations, dg::DGSEM, cache, - removed_elements::Vector{Int}) - @assert length(removed_elements) % (2^ndims(mesh))==0 "The length of `removed_elements` must be a multiple of 2^ndims(mesh)." - + element_ids_new::Vector{Int}) @unpack weights = dg.basis @unpack inverse_jacobian = cache.elements - # Precompute list with new element ids after refinement - element_ids_new = zeros(Int, div(length(removed_elements), 2^ndims(mesh))) - for i in eachindex(element_ids_new) - # The new element id is the id of the first child minus (2^ndims(mesh) - 1) times the number of already coarsened elements. - element_ids_new[i] = removed_elements[2^ndims(mesh) * (i - 1) + 1] - - (2^ndims(mesh) - 1) * (i - 1) - end - # Apply limiter to coarsened elements @threaded for element in element_ids_new # determine minimum value diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index 196a5814f98..4ad8fac53d6 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -852,6 +852,38 @@ function original2refined(original_cell_ids, refined_original_cells, mesh) return shifted_cell_ids[original_cell_ids] end +# Auxiliary function to compute the new element ids for refined elements +# 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) + 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 removed elements +# Used when applying positivity limiter after coarsening step +function compute_new_ids_removed_elements(elements_to_remove, mesh) + @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 + """ ControllerThreeLevel(semi, indicator; base_level=1, med_level=base_level, med_threshold=0.0, diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index 6ea547e51de..591f2e7b0f1 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -96,8 +96,11 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, # Apply the positivity limiter to the solution if limiter! !== nothing + # Precompute list with new element ids after refinement + element_ids_new = compute_new_id_refined_elements(elements_to_refine, mesh) + @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, - elements_to_refine, + element_ids_new, u_mean_refined_elements) end @@ -249,8 +252,11 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, # Apply the positivity limiter to the solution if limiter! !== nothing + # Precompute list with new element ids after coarsening + element_ids_new = compute_new_id_removed_elements(elements_to_remove, mesh) + @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, - elements_to_remove) + element_ids_new) end return nothing diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 91064d1625f..cd518c2a925 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -202,8 +202,11 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM # Apply the positivity limiter to the solution if limiter! !== nothing + # Precompute list with new element ids after refinement + element_ids_new = compute_new_id_refined_elements(elements_to_refine, mesh) + @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, - elements_to_refine, + element_ids_new, u_mean_refined_elements) end @@ -410,8 +413,11 @@ function coarsen!(u_ode::AbstractVector, adaptor, # Apply the positivity limiter to the solution if limiter! !== nothing + # Precompute list with new element ids after coarsening + element_ids_new = compute_new_id_removed_elements(elements_to_remove, mesh) + @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, - elements_to_remove) + element_ids_new) end return nothing diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index cb8a1e49d35..779f0938252 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -138,8 +138,11 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM # Apply the positivity limiter to the solution if limiter! !== nothing + # Precompute list with new element ids after refinement + element_ids_new = compute_new_id_refined_elements(elements_to_refine, mesh) + @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, - elements_to_refine, + element_ids_new, u_mean_refined_elements) end @@ -333,8 +336,11 @@ function coarsen!(u_ode::AbstractVector, adaptor, # Apply the positivity limiter to the solution if limiter! !== nothing + # Precompute list with new element ids after coarsening + element_ids_new = compute_new_id_removed_elements(elements_to_remove, mesh) + @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, - elements_to_remove) + element_ids_new) end return nothing From c3d750f6d7da0828c43e5ee13cab660f94b386b6 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 28 Aug 2025 15:05:28 +0200 Subject: [PATCH 17/17] Fix function name --- src/callbacks_step/amr_dg1d.jl | 4 ++-- src/callbacks_step/amr_dg2d.jl | 4 ++-- src/callbacks_step/amr_dg3d.jl | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index 591f2e7b0f1..492b1efb233 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -97,7 +97,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, # Apply the positivity limiter to the solution if limiter! !== nothing # Precompute list with new element ids after refinement - element_ids_new = compute_new_id_refined_elements(elements_to_refine, mesh) + element_ids_new = compute_new_ids_refined_elements(elements_to_refine, mesh) @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, element_ids_new, @@ -253,7 +253,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, # Apply the positivity limiter to the solution if limiter! !== nothing # Precompute list with new element ids after coarsening - element_ids_new = compute_new_id_removed_elements(elements_to_remove, mesh) + element_ids_new = compute_new_ids_removed_elements(elements_to_remove, mesh) @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, element_ids_new) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index cd518c2a925..dce3b99d840 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -203,7 +203,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM # Apply the positivity limiter to the solution if limiter! !== nothing # Precompute list with new element ids after refinement - element_ids_new = compute_new_id_refined_elements(elements_to_refine, mesh) + element_ids_new = compute_new_ids_refined_elements(elements_to_refine, mesh) @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, element_ids_new, @@ -414,7 +414,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, # Apply the positivity limiter to the solution if limiter! !== nothing # Precompute list with new element ids after coarsening - element_ids_new = compute_new_id_removed_elements(elements_to_remove, mesh) + element_ids_new = compute_new_ids_removed_elements(elements_to_remove, mesh) @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, element_ids_new) diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 779f0938252..b09845fad3d 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -139,7 +139,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM # Apply the positivity limiter to the solution if limiter! !== nothing # Precompute list with new element ids after refinement - element_ids_new = compute_new_id_refined_elements(elements_to_refine, mesh) + element_ids_new = compute_new_ids_refined_elements(elements_to_refine, mesh) @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, element_ids_new, @@ -337,7 +337,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, # Apply the positivity limiter to the solution if limiter! !== nothing # Precompute list with new element ids after coarsening - element_ids_new = compute_new_id_removed_elements(elements_to_remove, mesh) + element_ids_new = compute_new_ids_removed_elements(elements_to_remove, mesh) @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache, element_ids_new)