From b98b6f07c3b4242eaa4ffe5d17dd143a5a267e60 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 17 May 2024 08:21:24 +0200 Subject: [PATCH] Mean value in Zhang-Shu limiter on curved meshes (#1945) * solution mean value computation must take metric terms into account on mesh types other than TreeMesh * introduce helper function to compute mean value on 2d element * apply formatter * introduce helper function in 3d * simplify mean value computation. no need for a separate function * fix setting total_volume to be zero * put calc_mean_value functions back in * apply formatter * add generic get_inverse_jacobian function and remove other specializations * swap order of element and indices for consistency * revert ordering. see if tests actually run * fix index order * fix mistake in TreeMesh version * clarify the tree mesh version --- src/callbacks_stage/positivity_zhang_shu_dg2d.jl | 11 ++++++++--- src/callbacks_stage/positivity_zhang_shu_dg3d.jl | 11 ++++++++--- src/solvers/dgsem_structured/dg.jl | 8 ++++++++ src/solvers/dgsem_tree/dg.jl | 6 ++++++ 4 files changed, 30 insertions(+), 6 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index b37ed9c49d5..813dd65878b 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -8,6 +8,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{2}, equations, dg::DGSEM, cache) @unpack weights = dg.basis + @unpack inverse_jacobian = cache.elements @threaded for element in eachelement(dg, cache) # determine minimum value @@ -22,12 +23,16 @@ function limiter_zhang_shu!(u, threshold::Real, variable, # 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] + u_mean += u_node * weights[i] * weights[j] * volume_jacobian + total_volume += weights[i] * weights[j] * volume_jacobian end - # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2 - u_mean = u_mean / 2^ndims(mesh) + # 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). diff --git a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl index 773a236d831..156abf35b4c 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl @@ -8,6 +8,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{3}, equations, dg::DGSEM, cache) @unpack weights = dg.basis + @unpack inverse_jacobian = cache.elements @threaded for element in eachelement(dg, cache) # determine minimum value @@ -22,12 +23,16 @@ function limiter_zhang_shu!(u, threshold::Real, variable, # 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] + u_mean += u_node * weights[i] * weights[j] * weights[k] * volume_jacobian + total_volume += weights[i] * weights[j] * weights[k] * volume_jacobian end - # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2 - u_mean = u_mean / 2^ndims(mesh) + # 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). diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index 00e321fba65..424ed5e1e7e 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -76,6 +76,14 @@ end end end +@inline function get_inverse_jacobian(inverse_jacobian, + mesh::Union{StructuredMesh, StructuredMeshView, + UnstructuredMesh2D, P4estMesh, + T8codeMesh}, + indices...) + return inverse_jacobian[indices...] +end + include("containers.jl") include("dg_1d.jl") include("dg_2d.jl") diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index ef9a42b4c1a..0993b3c9b85 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -42,6 +42,12 @@ function volume_jacobian(element, mesh::TreeMesh, cache) return inv(cache.elements.inverse_jacobian[element])^ndims(mesh) end +@inline function get_inverse_jacobian(inverse_jacobian, mesh::TreeMesh, + indices...) + element = last(indices) + return inverse_jacobian[element] +end + # Indicators used for shock-capturing and AMR include("indicators.jl") include("indicators_1d.jl")