From fed45fc36927d9536ee14ff2ef4c964cebef1ed7 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Fri, 17 May 2024 09:08:31 +0200 Subject: [PATCH] add limiter fix for curvilinear meshes --- .../positivity_shallow_water_dg2d.jl | 23 ++++++++++++++----- 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/src/callbacks_stage/positivity_shallow_water_dg2d.jl b/src/callbacks_stage/positivity_shallow_water_dg2d.jl index 9acdc4d..bea301f 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg2d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg2d.jl @@ -12,6 +12,7 @@ function limiter_shallow_water!(u, threshold::Real, variable, equations::ShallowWaterEquationsWetDry2D, dg::DGSEM, cache) @unpack weights = dg.basis + @unpack inverse_jacobian = cache.elements Trixi.@threaded for element in eachelement(dg, cache) # determine minimum value @@ -26,12 +27,16 @@ function limiter_shallow_water!(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(Trixi.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). @@ -96,6 +101,7 @@ function limiter_shallow_water!(u, threshold::Real, variable, equations::ShallowWaterMultiLayerEquations2D, dg::DGSEM, cache) @unpack weights = dg.basis + @unpack inverse_jacobian = cache.elements Trixi.@threaded for element in eachelement(dg, cache) # Limit layerwise @@ -112,12 +118,17 @@ function limiter_shallow_water!(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(Trixi.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. # The waterheight `h` is limited independently in each layer.