Skip to content

Commit ae35f82

Browse files
DanielDoehringranochabennibolm
authored
Single out u_mean computation into function (#2541)
* Single out `u_mean` computation into function * revisit * shorten * fmt * shorten * debug * Update src/callbacks_stage/positivity_zhang_shu_dg1d.jl Co-authored-by: Hendrik Ranocha <[email protected]> * re-order * node_vol * mesh type respect * real for treemesh * fix * Update src/callbacks_stage/positivity_zhang_shu_dg1d.jl Co-authored-by: Benjamin Bolm <[email protected]> * refactor * comment * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <[email protected]> --------- Co-authored-by: Hendrik Ranocha <[email protected]> Co-authored-by: Benjamin Bolm <[email protected]>
1 parent 28bd001 commit ae35f82

File tree

8 files changed

+70
-81
lines changed

8 files changed

+70
-81
lines changed

‎src/callbacks_stage/entropy_bounded_limiter_1d.jl‎

Lines changed: 1 addition & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -7,8 +7,6 @@
77

88
function limiter_entropy_bounded!(u, u_prev, exp_entropy_decrease_max,
99
mesh::AbstractMesh{1}, equations, dg::DGSEM, cache)
10-
@unpack weights = dg.basis
11-
1210
@threaded for element in eachelement(dg, cache)
1311
# Minimum exponentiated entropy within the current `element`
1412
# of the previous iterate `u_prev`
@@ -35,14 +33,7 @@ function limiter_entropy_bounded!(u, u_prev, exp_entropy_decrease_max,
3533
# Limiting only if entropy DECREASE below a user defined threshold is detected.
3634
d_exp_s_min < exp_entropy_decrease_max || continue
3735

38-
# Compute mean value
39-
u_mean = zero(get_node_vars(u, equations, dg, 1, element))
40-
for i in eachnode(dg)
41-
u_node = get_node_vars(u, equations, dg, i, element)
42-
u_mean += u_node * weights[i]
43-
end
44-
# Note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2
45-
u_mean = u_mean / 2^ndims(mesh)
36+
u_mean = compute_u_mean(u, element, mesh, equations, dg, cache)
4637

4738
entropy_change_mean = exp_entropy_change(pressure(u_mean, equations),
4839
density(u_mean, equations),

‎src/callbacks_stage/entropy_bounded_limiter_2d.jl‎

Lines changed: 2 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -7,9 +7,6 @@
77

88
function limiter_entropy_bounded!(u, u_prev, exp_entropy_decrease_max,
99
mesh::AbstractMesh{2}, equations, dg::DGSEM, cache)
10-
@unpack weights = dg.basis
11-
@unpack inverse_jacobian = cache.elements
12-
1310
@threaded for element in eachelement(dg, cache)
1411
# Minimum exponentiated entropy within the current `element`
1512
# of the previous iterate `u_prev`
@@ -35,18 +32,8 @@ function limiter_entropy_bounded!(u, u_prev, exp_entropy_decrease_max,
3532
# Detect if limiting is necessary.
3633
# Limiting only if entropy DECREASE below a user defined threshold is detected.
3734
d_exp_s_min < exp_entropy_decrease_max || continue
38-
# Compute mean value
39-
u_mean = zero(get_node_vars(u, equations, dg, 1, 1, element))
40-
total_volume = zero(eltype(u))
41-
for j in eachnode(dg), i in eachnode(dg)
42-
volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh,
43-
i, j, element)))
44-
u_node = get_node_vars(u, equations, dg, i, j, element)
45-
u_mean += u_node * weights[i] * weights[j] * volume_jacobian
46-
total_volume += weights[i] * weights[j] * volume_jacobian
47-
end
48-
# normalize with the total volume
49-
u_mean = u_mean / total_volume
35+
36+
u_mean = compute_u_mean(u, element, mesh, equations, dg, cache)
5037

5138
entropy_change_mean = exp_entropy_change(pressure(u_mean, equations),
5239
density(u_mean, equations),

‎src/callbacks_stage/entropy_bounded_limiter_3d.jl‎

Lines changed: 1 addition & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -7,9 +7,6 @@
77

88
function limiter_entropy_bounded!(u, u_prev, exp_entropy_decrease_max,
99
mesh::AbstractMesh{3}, equations, dg::DGSEM, cache)
10-
@unpack weights = dg.basis
11-
@unpack inverse_jacobian = cache.elements
12-
1310
@threaded for element in eachelement(dg, cache)
1411
# Minimum exponentiated entropy within the current `element`
1512
# of the previous iterate `u_prev`
@@ -35,18 +32,7 @@ function limiter_entropy_bounded!(u, u_prev, exp_entropy_decrease_max,
3532
# Detect if limiting is necessary. Avoid division by ("near") zero
3633
d_exp_s_min < exp_entropy_decrease_max || continue
3734

38-
# Compute mean value
39-
u_mean = zero(get_node_vars(u, equations, dg, 1, 1, 1, element))
40-
total_volume = zero(eltype(u))
41-
for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
42-
volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh,
43-
i, j, k, element)))
44-
u_node = get_node_vars(u, equations, dg, i, j, k, element)
45-
u_mean += u_node * weights[i] * weights[j] * weights[k] * volume_jacobian
46-
total_volume += weights[i] * weights[j] * weights[k] * volume_jacobian
47-
end
48-
# normalize with the total volume
49-
u_mean = u_mean / total_volume
35+
u_mean = compute_u_mean(u, element, mesh, equations, dg, cache)
5036

5137
entropy_change_mean = exp_entropy_change(pressure(u_mean, equations),
5238
density(u_mean, equations),

‎src/callbacks_stage/positivity_zhang_shu_dg1d.jl‎

Lines changed: 1 addition & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -7,8 +7,6 @@
77

88
function limiter_zhang_shu!(u, threshold::Real, variable,
99
mesh::AbstractMesh{1}, equations, dg::DGSEM, cache)
10-
@unpack weights = dg.basis
11-
1210
@threaded for element in eachelement(dg, cache)
1311
# determine minimum value
1412
value_min = typemax(eltype(u))
@@ -20,14 +18,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable,
2018
# detect if limiting is necessary
2119
value_min < threshold || continue
2220

23-
# compute mean value
24-
u_mean = zero(get_node_vars(u, equations, dg, 1, element))
25-
for i in eachnode(dg)
26-
u_node = get_node_vars(u, equations, dg, i, element)
27-
u_mean += u_node * weights[i]
28-
end
29-
# note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2
30-
u_mean = u_mean / 2^ndims(mesh)
21+
u_mean = compute_u_mean(u, element, mesh, equations, dg, cache)
3122

3223
# We compute the value directly with the mean values, as we assume that
3324
# Jensen's inequality holds (e.g. pressure for compressible Euler equations).

‎src/callbacks_stage/positivity_zhang_shu_dg2d.jl‎

Lines changed: 1 addition & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -7,9 +7,6 @@
77

88
function limiter_zhang_shu!(u, threshold::Real, variable,
99
mesh::AbstractMesh{2}, equations, dg::DGSEM, cache)
10-
@unpack weights = dg.basis
11-
@unpack inverse_jacobian = cache.elements
12-
1310
@threaded for element in eachelement(dg, cache)
1411
# determine minimum value
1512
value_min = typemax(eltype(u))
@@ -21,18 +18,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable,
2118
# detect if limiting is necessary
2219
value_min < threshold || continue
2320

24-
# compute mean value
25-
u_mean = zero(get_node_vars(u, equations, dg, 1, 1, element))
26-
total_volume = zero(eltype(u))
27-
for j in eachnode(dg), i in eachnode(dg)
28-
volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh,
29-
i, j, element)))
30-
u_node = get_node_vars(u, equations, dg, i, j, element)
31-
u_mean += u_node * weights[i] * weights[j] * volume_jacobian
32-
total_volume += weights[i] * weights[j] * volume_jacobian
33-
end
34-
# normalize with the total volume
35-
u_mean = u_mean / total_volume
21+
u_mean = compute_u_mean(u, element, mesh, equations, dg, cache)
3622

3723
# We compute the value directly with the mean values, as we assume that
3824
# Jensen's inequality holds (e.g. pressure for compressible Euler equations).

‎src/callbacks_stage/positivity_zhang_shu_dg3d.jl‎

Lines changed: 1 addition & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -7,9 +7,6 @@
77

88
function limiter_zhang_shu!(u, threshold::Real, variable,
99
mesh::AbstractMesh{3}, equations, dg::DGSEM, cache)
10-
@unpack weights = dg.basis
11-
@unpack inverse_jacobian = cache.elements
12-
1310
@threaded for element in eachelement(dg, cache)
1411
# determine minimum value
1512
value_min = typemax(eltype(u))
@@ -21,18 +18,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable,
2118
# detect if limiting is necessary
2219
value_min < threshold || continue
2320

24-
# compute mean value
25-
u_mean = zero(get_node_vars(u, equations, dg, 1, 1, 1, element))
26-
total_volume = zero(eltype(u))
27-
for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
28-
volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh,
29-
i, j, k, element)))
30-
u_node = get_node_vars(u, equations, dg, i, j, k, element)
31-
u_mean += u_node * weights[i] * weights[j] * weights[k] * volume_jacobian
32-
total_volume += weights[i] * weights[j] * weights[k] * volume_jacobian
33-
end
34-
# normalize with the total volume
35-
u_mean = u_mean / total_volume
21+
u_mean = compute_u_mean(u, element, mesh, equations, dg, cache)
3622

3723
# We compute the value directly with the mean values, as we assume that
3824
# Jensen's inequality holds (e.g. pressure for compressible Euler equations).

‎src/meshes/tree_mesh.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -252,5 +252,7 @@ end
252252
isperiodic(mesh::TreeMesh) = isperiodic(mesh.tree)
253253
isperiodic(mesh::TreeMesh, dimension) = isperiodic(mesh.tree, dimension)
254254

255+
Base.real(::TreeMesh{NDIMS, TreeType, RealT}) where {NDIMS, TreeType, RealT} = RealT
256+
255257
include("parallel_tree_mesh.jl")
256258
end # @muladd

‎src/solvers/dgsem/dgsem.jl‎

Lines changed: 61 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,7 @@ function DGSEM(RealT, polydeg::Integer,
5151
end
5252

5353
# This API is no longer documented, and we recommend avoiding its public use.
54-
function DGSEM(polydeg, surface_flux = flux_central,
54+
function DGSEM(polydeg::Integer, surface_flux = flux_central,
5555
volume_integral = VolumeIntegralWeakForm())
5656
DGSEM(Float64, polydeg, surface_flux, volume_integral)
5757
end
@@ -71,4 +71,64 @@ end
7171
@inline polydeg(dg::DGSEM) = polydeg(dg.basis)
7272

7373
Base.summary(io::IO, dg::DGSEM) = print(io, "DGSEM(polydeg=$(polydeg(dg)))")
74+
75+
# `compute_u_mean` used in:
76+
# (Stage-) Callbacks `EntropyBoundedLimiter` and `PositivityPreservingLimiterZhangShu`
77+
78+
# positional arguments `mesh` and `cache` passed in to match signature of 2D/3D functions
79+
@inline function compute_u_mean(u::AbstractArray{<:Any, 3}, element,
80+
mesh::AbstractMesh{1}, equations, dg::DGSEM, cache)
81+
@unpack weights = dg.basis
82+
83+
u_mean = zero(get_node_vars(u, equations, dg, 1, element))
84+
for i in eachnode(dg)
85+
u_node = get_node_vars(u, equations, dg, i, element)
86+
u_mean += u_node * weights[i]
87+
end
88+
# normalize with the total volume
89+
# note that the reference element is [-1,1], thus the weights sum to 2
90+
return 0.5f0 * u_mean
91+
end
92+
93+
@inline function compute_u_mean(u::AbstractArray{<:Any, 4}, element,
94+
mesh::AbstractMesh{2}, equations, dg::DGSEM, cache)
95+
@unpack weights = dg.basis
96+
@unpack inverse_jacobian = cache.elements
97+
98+
node_volume = zero(real(mesh))
99+
total_volume = zero(node_volume)
100+
101+
u_mean = zero(get_node_vars(u, equations, dg, 1, 1, element))
102+
for j in eachnode(dg), i in eachnode(dg)
103+
volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh,
104+
i, j, element)))
105+
node_volume = weights[i] * weights[j] * volume_jacobian
106+
total_volume += node_volume
107+
108+
u_node = get_node_vars(u, equations, dg, i, j, element)
109+
u_mean += u_node * node_volume
110+
end
111+
return u_mean / total_volume # normalize with the total volume
112+
end
113+
114+
@inline function compute_u_mean(u::AbstractArray{<:Any, 5}, element,
115+
mesh::AbstractMesh{3}, equations, dg::DGSEM, cache)
116+
@unpack weights = dg.basis
117+
@unpack inverse_jacobian = cache.elements
118+
119+
node_volume = zero(real(mesh))
120+
total_volume = zero(node_volume)
121+
122+
u_mean = zero(get_node_vars(u, equations, dg, 1, 1, 1, element))
123+
for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
124+
volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh,
125+
i, j, k, element)))
126+
node_volume = weights[i] * weights[j] * weights[k] * volume_jacobian
127+
total_volume += node_volume
128+
129+
u_node = get_node_vars(u, equations, dg, i, j, k, element)
130+
u_mean += u_node * node_volume
131+
end
132+
return u_mean / total_volume # normalize with the total volume
133+
end
74134
end # @muladd

0 commit comments

Comments
 (0)