Skip to content

Commit

Permalink
add limiter fix for curvilinear meshes (#46)
Browse files Browse the repository at this point in the history
  • Loading branch information
patrickersing authored May 22, 2024
1 parent da3af82 commit 117489e
Showing 1 changed file with 17 additions and 6 deletions.
23 changes: 17 additions & 6 deletions src/callbacks_stage/positivity_shallow_water_dg2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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).
Expand Down Expand Up @@ -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
Expand All @@ -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.
Expand Down

0 comments on commit 117489e

Please sign in to comment.