Skip to content

Commit

Permalink
fix weird formatting and add 'format: noindent' where missing. fix cr…
Browse files Browse the repository at this point in the history 
…ashing structured mesh run
  • Loading branch information
@andrewwinters5000
andrewwinters5000 committed Jul 5, 2023
1 parent a5b128f commit a2c12ad
Show file tree
Hide file tree
Showing 4 changed files with 152 additions and 142 deletions.
141 changes: 73 additions & 68 deletions src/callbacks_stage/positivity_shallow_water.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,79 +3,84 @@
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin
"""
PositivityPreservingLimiterShallowWater(; variables)
#! format: noindent

The limiter is specifically designed for the shallow water equations.
It is applied to all scalar `variables` in their given order
using the defined `threshold_limiter` from the [`ShallowWaterEquations1D`](@ref) struct
or the [`ShallowWaterEquations2D`](@ref) struct to determine the minimal acceptable values.
The order of the `variables` is important and might have a strong influence
on the robustness.
As opposed to the standard version of the [`PositivityPreservingLimiterZhangShu`](@ref),
nodes with a water height below the `threshold_limiter` are treated in a special way.
To avoid numerical problems caused by velocities close to zero,
the velocity is cut off, such that the node can be identified as "dry". The special feature of the
`ShallowWaterEquations` used here is that the bottom topography is stored as an additional
quantity in the solution vector `u`. However, the value of the bottom topography
should not be changed. That is why, it is not limited.
After the limiting process is applied to all degrees of freedom, for safety reasons,
the `threshold_limiter` is applied again on all the DG nodes in order to avoid water height below.
In the case where the cell mean value is below the threshold before applying the limiter,
there could still be dry nodes afterwards due to the logic of the limiter.
This fully-discrete positivity-preserving limiter is based on the work of
- Zhang, Shu (2011)
Maximum-principle-satisfying and positivity-preserving high-order schemes
for conservation laws: survey and new developments
[doi: 10.1098/rspa.2011.0153](https://doi.org/10.1098/rspa.2011.0153)
"""
struct PositivityPreservingLimiterShallowWater{N, Variables <: NTuple{N, Any}}
variables::Variables
end
"""
PositivityPreservingLimiterShallowWater(; variables)
function PositivityPreservingLimiterShallowWater(; variables)
PositivityPreservingLimiterShallowWater(variables)
end
The limiter is specifically designed for the shallow water equations.
It is applied to all scalar `variables` in their given order
using the defined `threshold_limiter` from the [`ShallowWaterEquations1D`](@ref) struct
or the [`ShallowWaterEquations2D`](@ref) struct to determine the minimal acceptable values.
The order of the `variables` is important and might have a strong influence
on the robustness.
function (limiter!::PositivityPreservingLimiterShallowWater)(u_ode, integrator,
semi::AbstractSemidiscretization,
t)
u = wrap_array(u_ode, semi)
@trixi_timeit timer() "positivity-preserving limiter" limiter_shallow_water!(u,
limiter!.variables,
mesh_equations_solver_cache(semi)...)
end
As opposed to the standard version of the [`PositivityPreservingLimiterZhangShu`](@ref),
nodes with a water height below the `threshold_limiter` are treated in a special way.
To avoid numerical problems caused by velocities close to zero,
the velocity is cut off, such that the node can be identified as "dry". The special feature of the
`ShallowWaterEquations` used here is that the bottom topography is stored as an additional
quantity in the solution vector `u`. However, the value of the bottom topography
should not be changed. That is why, it is not limited.
# Iterate over tuples in a type-stable way using "lispy tuple programming",
# similar to https://stackoverflow.com/a/55849398:
# Iterating over tuples of different functions isn't type-stable in general
# but accessing the first element of a tuple is type-stable. Hence, it's good
# to process one element at a time and replace iteration by recursion here.
# Note that you shouldn't use this with too many elements per tuple since the
# compile times can increase otherwise - but a handful of elements per tuple
# is definitely fine.
function limiter_shallow_water!(u, variables::NTuple{N, Any},
mesh,
equations::Union{ShallowWaterEquations1D,
ShallowWaterEquations2D},
solver, cache) where {N}
variable = first(variables)
remaining_variables = Base.tail(variables)
After the limiting process is applied to all degrees of freedom, for safety reasons,
the `threshold_limiter` is applied again on all the DG nodes in order to avoid water height below.
In the case where the cell mean value is below the threshold before applying the limiter,
there could still be dry nodes afterwards due to the logic of the limiter.
limiter_shallow_water!(u, equations.threshold_limiter, variable, mesh, equations,
solver, cache)
limiter_shallow_water!(u, remaining_variables, mesh, equations, solver, cache)
return nothing
end
This fully-discrete positivity-preserving limiter is based on the work of
- Zhang, Shu (2011)
Maximum-principle-satisfying and positivity-preserving high-order schemes
for conservation laws: survey and new developments
[doi: 10.1098/rspa.2011.0153](https://doi.org/10.1098/rspa.2011.0153)
"""
struct PositivityPreservingLimiterShallowWater{N, Variables <: NTuple{N, Any}}
variables::Variables
end

# terminate the type-stable iteration over tuples
function limiter_shallow_water!(u, variables::Tuple{},
mesh,
equations::Union{ShallowWaterEquations1D,
ShallowWaterEquations2D},
solver, cache)
nothing
end
function PositivityPreservingLimiterShallowWater(; variables)
PositivityPreservingLimiterShallowWater(variables)
end

include("positivity_shallow_water_dg2d.jl")
function (limiter!::PositivityPreservingLimiterShallowWater)(u_ode, integrator,
semi::AbstractSemidiscretization,
t)
u = wrap_array(u_ode, semi)
@trixi_timeit timer() "positivity-preserving limiter" limiter_shallow_water!(u,
limiter!.variables,
mesh_equations_solver_cache(semi)...)
end

# Iterate over tuples in a type-stable way using "lispy tuple programming",
# similar to https://stackoverflow.com/a/55849398:
# Iterating over tuples of different functions isn't type-stable in general
# but accessing the first element of a tuple is type-stable. Hence, it's good
# to process one element at a time and replace iteration by recursion here.
# Note that you shouldn't use this with too many elements per tuple since the
# compile times can increase otherwise - but a handful of elements per tuple
# is definitely fine.
function limiter_shallow_water!(u, variables::NTuple{N, Any},
mesh,
equations::Union{ShallowWaterEquations1D,
ShallowWaterEquations2D},
solver, cache) where {N}
variable = first(variables)
remaining_variables = Base.tail(variables)

limiter_shallow_water!(u, equations.threshold_limiter, variable, mesh, equations,
solver, cache)
limiter_shallow_water!(u, remaining_variables, mesh, equations, solver, cache)
return nothing
end

# terminate the type-stable iteration over tuples
function limiter_shallow_water!(u, variables::Tuple{},
mesh,
equations::Union{ShallowWaterEquations1D,
ShallowWaterEquations2D},
solver, cache)
nothing
end

include("positivity_shallow_water_dg2d.jl")
end # @muladd
138 changes: 71 additions & 67 deletions src/callbacks_stage/positivity_shallow_water_dg2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,82 +3,86 @@
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin
function limiter_shallow_water!(u, threshold::Real, variable,
mesh::AbstractMesh{2},
equations::ShallowWaterEquations2D, dg::DGSEM, cache)
@unpack weights = dg.basis

@threaded for element in eachelement(dg, cache)
# 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
#! format: noindent

# detect if limiting is necessary
value_min < threshold || continue
function limiter_shallow_water!(u, threshold::Real, variable,
mesh::AbstractMesh{2},
equations::ShallowWaterEquations2D, dg::DGSEM, cache)
@unpack weights = dg.basis

# compute mean value
u_mean = zero(get_node_vars(u, equations, dg, 1, 1, element))
for j in eachnode(dg), i in eachnode(dg)
u_node = get_node_vars(u, equations, dg, i, j, element)
u_mean += u_node * weights[i] * weights[j]
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 j in eachnode(dg), i in eachnode(dg)
u_node = get_node_vars(u, equations, dg, i, j, element)

# Cut off velocity in case that the water height is smaller than the threshold

h_node, h_v1_node, h_v2_node, b_node = u_node
h_mean, h_v1_mean, h_v2_mean, _ = u_mean # b_mean is not used as it must not be overwritten

if h_node <= threshold
h_v1_node = zero(eltype(u))
h_v2_node = zero(eltype(u))
h_v1_mean = zero(eltype(u))
h_v2_mean = zero(eltype(u))
end

u_node = SVector(h_node, h_v1_node, h_v2_node, b_node)
u_mean = SVector(h_mean, h_v1_mean, h_v2_mean, b_node)

# When velocities are cut off, the only averaged value is the water height,
# because the velocities are set to zero and this value is passed.
# Otherwise, the velocities are averaged, as well.
# Note that the auxiliary bottom topography variable `b` is never limited.
set_node_vars!(u, theta * u_node + (1 - theta) * u_mean,
equations, dg, i, j, element)
end
@threaded for element in eachelement(dg, cache)
# 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

# "Safety" application of the wet/dry thresholds over all the DG nodes
# on the current `element` after the limiting above in order to avoid dry nodes.
# If the value_mean < threshold before applying limiter, there
# could still be dry nodes afterwards due to logic of the limiting
for j in eachnode(dg), i in eachnode(dg)
u_node = get_node_vars(u, equations, dg, i, j, element)
# detect if limiting is necessary
value_min < threshold || continue

h, h_v1, h_v2, b = u_node
# compute mean value
u_mean = zero(get_node_vars(u, equations, dg, 1, 1, element))
for j in eachnode(dg), i in eachnode(dg)
u_node = get_node_vars(u, equations, dg, i, j, element)
u_mean += u_node * weights[i] * weights[j]
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 j in eachnode(dg), i in eachnode(dg)
u_node = get_node_vars(u, equations, dg, i, j, element)

if h <= threshold
h = threshold
h_v1 = zero(eltype(u))
h_v2 = zero(eltype(u))
end
# Cut off velocity in case that the water height is smaller than the threshold

u_node = SVector(h, h_v1, h_v2, b)
h_node, h_v1_node, h_v2_node, b_node = u_node
h_mean, h_v1_mean, h_v2_mean, _ = u_mean # b_mean is not used as it must not be overwritten

set_node_vars!(u, u_node, equations, dg, i, j, element)
if h_node <= threshold
h_v1_node = zero(eltype(u))
h_v2_node = zero(eltype(u))
h_v1_mean = zero(eltype(u))
h_v2_mean = zero(eltype(u))
end

u_node = SVector(h_node, h_v1_node, h_v2_node, b_node)
u_mean = SVector(h_mean, h_v1_mean, h_v2_mean, b_node)

# When velocities are cut off, the only averaged value is the water height,
# because the velocities are set to zero and this value is passed.
# Otherwise, the velocities are averaged, as well.
# Note that the auxiliary bottom topography variable `b` is never limited.
set_node_vars!(u, theta * u_node + (1 - theta) * u_mean,
equations, dg, i, j, element)
end
end

# "Safety" application of the wet/dry thresholds over all the DG nodes
# on the current `element` after the limiting above in order to avoid dry nodes.
# If the value_mean < threshold before applying limiter, there
# could still be dry nodes afterwards due to logic of the limiting
@threaded for element in eachelement(dg, cache)
for j in eachnode(dg), i in eachnode(dg)
u_node = get_node_vars(u, equations, dg, i, j, element)

h, h_v1, h_v2, b = u_node

if h <= threshold
h = threshold
h_v1 = zero(eltype(u))
h_v2 = zero(eltype(u))
end

u_node = SVector(h, h_v1, h_v2, b)

return nothing
set_node_vars!(u, u_node, equations, dg, i, j, element)
end
end

return nothing
end
end # @muladd
9 changes: 5 additions & 4 deletions src/equations/shallow_water_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,13 +56,14 @@ References for the SWE are many but a good introduction is available in Chapter
struct ShallowWaterEquations2D{RealT <: Real} <: AbstractShallowWaterEquations{2, 4}
gravity::RealT # gravitational constant
H0::RealT # constant "lake-at-rest" total water height
threshold_limiter::RealT # Threshold to use in `PositivityPreservingLimiterShallowWater` on water height,
# as a (small) shift on the initial condition and cutoff before the
# next time step.
# `threshold_limiter` used in `PositivityPreservingLimiterShallowWater` on water height,
# as a (small) shift on the initial condition and cutoff before the next time step.
# Default is 500*eps() which in double precision is ≈1e-13.
threshold_wet::RealT # Threshold to be applied on water height to define when the flow is "wet"
threshold_limiter::RealT
# `threshold_wet` applied on water height to define when the flow is "wet"
# before calculating the numerical flux.
# Default is 5*eps() which in double precision is ≈1e-15.
threshold_wet::RealT
end

# Allow for flexibility to set the gravitational constant within an elixir depending on the
Expand Down
6 changes: 3 additions & 3 deletions test/test_structured_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ isdir(outdir) && rm(outdir, recursive=true)
end

@trixi_testset "elixir_advection_coupled.jl" begin
@test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_coupled.jl"),
@test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_coupled.jl"),
l2 = [7.816742843181738e-6, 7.816742843196112e-6],
linf = [6.314906965543265e-5, 6.314906965410039e-5],
coverage_override = (maxiters=10^5,))
Expand Down Expand Up @@ -272,8 +272,8 @@ isdir(outdir) && rm(outdir, recursive=true)

@trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl" begin
@test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_wet_dry.jl"),
l2 = [0.019731646454946954, 6.682390396110676e-14, 1.1925157549257398e-14, 0.07715172600379547],
linf = [0.5000000000001029, 3.959125826488326e-13, 4.556539260222467e-14, 2.0],
l2 = [0.019731646454942086, 1.0694532773278277e-14, 1.1969913383405568e-14, 0.0771517260037954],
linf = [0.4999999999998892, 6.067153702623552e-14, 4.4849667259339357e-14, 1.9999999999999993],
tspan = (0.0, 0.25))
end

Expand Down

0 comments on commit a2c12ad

Please sign in to comment.